remove_polynomials.pdl

=head2 remove_polynomials - hi2 pipeline component

=for ref

This fits a polynomial to the time series of each pixel the given cube in the pipeline data structure, and
subtracts it.  There are some options:

=over 3

=item POLY_PRESMOOTH - if set, this does a rowwise running-median followed by a columnwise running-median operation to remove stars and such, before generating the fit.

=item MIN_SMOOTH_SIZE - the data are shifted to their minimum value at the end of the operation.  The MIN_SMOOTH_SIZE performs a rowwise/columnwise median on the minimum value plane before the subtraction.

=item POLY_ORDER - the order of the polynomial (default is 3 -- cubic)

=item POLY_CUBE - if present, this is the name of the cube to modify.  Default is CEL, which modifies CEL_CUBE.

=item POLY_TRIM_PERCENTILE - if present, this trims some percentile of the max and min values from each time series before the fit.

=item POLY_DOUBLESTUFF - if present, this trims the <n> percentile of outlier points that might be driving the fit away from the true background, and refits. 

=back


=cut
use PDL::Slatec;
use PDL::NiceSlice;
use strict;
sub remove_polynomials {
    my $hash = shift;
    my $opt = shift;
    my $inst_cube = shift || 0;
    my $VERSION = "2.3 1-Mar-2013";

    print "into remove_polynomials...\n";

    $opt = {} unless defined($opt);
    
    $opt->{POLY_PRESMOOTH} = 0 unless(defined($opt->{POLY_PRESMOOTH}));
    $opt->{MIN_SMOOTH_SIZE} = 11 unless(defined($opt->{MIN_SMOOTH_SIZE}));
    $opt->{POLY_ORDER} = 3 unless(exists($opt->{POLY_ORDER}));
    $opt->{POLY_CUBE} = "CEL" unless(exists($opt->{POLY_CUBE}));
    $opt->{POLY_TRIM_PERCENTILE} = 0 unless(exists($opt->{POLY_TRIM_PERCENTILE}));
    $opt->{POLY_DOUBLESTUFF} = 0 unless(exists($opt->{POLY_DOUBLESTUFF}));

    my $pc = $opt->{POLY_CUBE}."_CUBE";
    my $cube = $hash->{$pc};
    my $n = $cube->dim(2);

    my $c2;

    # If we're presmoothing, then presmooth.
    if($opt->{POLY_PRESMOOTH}) {
	print "presmoothing for polynomial fit";
	$c2 = PDL->new_from_specification($cube->mv(2,0)->dims);
	for my $i(0..$c2->dim(0)-1) {
	    $c2->(($i)) .= $cube->(:,:,($i))->med2d(ones(1,$opt->{POLY_PRESMOOTH}))->med2d(ones($opt->{POLY_PRESMOOTH},1));
	    print ".";
	}
	print "\n";
    } else {
	$c2 = $cube->mv(2,0);  # c2 is (time,x,y)
    }

    print "fitting polynomials (order=$opt->{POLY_ORDER})...\n";
    
    # Break the polyfit up into pieces to avoid hogging too much temp memory at once
    my $pieces = 10;
    print "10 easy pieces...\n";
    my $ii;
    for $ii(0..9){
	print "$ii...";
	my $min = int($ii * $c2->dim(2)/10);
	my $max = int(($ii+1) * $c2->dim(2)/10)-1;
	if($max >= $c2->dim(2)-1) {
	    $max = $c2->dim(2)-1;
	}
	my $c3 = $c2->(:,:,$min:$max)->copy; # c3 is (time, x, y)
	$c3->badflag(1);

	if($opt->{POLY_TRIM_PERCENTILE}) {
	    my $ptt = $opt->{POLY_TRIM_PERCENTILE};
	    my $maxtrim = (ref($ptt) eq 'ARRAY') ? $ptt->[1] : $ptt;
	    my $mintrim = (ref($ptt) eq 'ARRAY') ? $ptt->[0] : $ptt;
	    my $rejrow = (ref($ptt) eq 'ARRAY') ? ($ptt->[2] // 0.1) : 0.1;
	    
	    my $c3i = $c3->qsorti;  # BAD gets sorted at the end
	    my $idex = ndcoords($c3i);
	    $idex->((0)) .= $c3i;
	    my $c3_xf = $c3->indexND($idex);

	    my $n_keep = (1 - $maxtrim/100) * $c3i->isgood->sumover->(*1);
	    my $n_rej =  ($mintrim/100) * $c3i->isgood->sumover->(*1);

	    my $badrow = ($n_keep - $n_rej) < $rejrow;

	    $c3_xf->where((xvals($c3i) < $n_rej) | (xvals($c3i) >= $n_keep) | $badrow) .= $c3->badvalue;

	} 


	my $eps = 0;
	print "c3 is ".join("x",$c3->dims)."\n";
	my($ndeg, $r, $ierr, $a) = polyfit( xvals($c3->dim(0)), $c3, ones($c3->dim(0)), $opt->{POLY_ORDER}, $eps );
	undef $ierr;
	undef $a;
	undef $ndeg;

	if($opt->{POLY_DOUBLESTUFF}) {
	    my $c4_num = (  ($opt->{POLY_DOUBLESTUFF} / 100) * $c3->dim(0) );
	    print "c3 is ".join("x",$c3->dims)."\n";
	    print "remove_polynomials: POLY_DOUBLESTUFF=$opt->{POLY_DOUBLESTUFF}: rejecting worst $c4_num points from each timeseries...\n";
	    our $c4 = (-($c3 - $r)->abs)->qsorti;
	    my $c3foo = $c3->index1d($c4);

	    $c3foo->(0:$c4_num) .= $c3->badvalue;
	    our @z;
	    our @r;
	    push(@r,$r);
	    push(@z,$c3);
	    ($ndeg, $r, $ierr, $a) = polyfit( xvals($c3->dim(0)), $c3, ones($c3->dim(0)), $opt->{POLY_ORDER}, $eps );
	    undef $ierr;
	    undef $a;
	    undef $ndeg;
	}

	$r -= $r->minimum->(*1);
	$r->where($r->isbad) .= 0;
	$cube->mv(2,0)->(:,:,$min:$max) -= $r;
    }
    print "\n";

    unless(exists($hash->{log})) {
	$hash->{log} = [];
    } 
    push(@{$hash->{log}}, "  remove_polynomials v$VERSION");
    
    return $hash;
}