#!/usr/bin/perl
#
# testParams.pl
# Program to run through the parameters for STAC.


@REJECT = (3.0); # List of "reject" parameters
@FRAC = (0.5); # List of "frac" parameters
@GRAD = (0.6); # List of "grad" parameters
@SEEING = (2.5); # List of seeing
$ITERATIONS = 1; # Number of iterations for each parameter combination

$RUN = "rm -f testout.fits* chi2_*.fits test_[0-3].fits.\{err,grad,mask,rejmap,shift*}; ./stac -v testout.fits test_0.fits test_1.fits test_2.fits test_3.fits";		# Command with which to run the STAC program
$RUNREJECT = "-k";		# Rejection parameter flag
$RUNFRAC = "-f";		# Frac parameter flag
$RUNGRAD = "-G";		# Grad parameter flag
$CREATEDIR = "/home/mithrandir/price/images/"; # Directory in which the image creation routines reside
$CREATE = "./mkfakeimages.pl -diag"; # Command for the image creation
$CREATESEEING = "-seeing";	# Switch for seeing in the image creation
$CREATENOCRS = "-numcrs 0";	# Switch for no CRs in the image creation

$WORKINGDIR = `pwd`;		# Working directory (current directory)
chomp($WORKINGDIR);

# Iterate over the parameter space
open OUTPUT, "> testParams.out";
foreach $seeing (@SEEING) {
    foreach $reject (@REJECT) {
	foreach $frac (@FRAC) {
	    foreach $grad (@GRAD) {
		@trueRejected = (); # Number of true CRs rejected in each iteration
		@falseRejected = (); # Number of false CRs rejected in each iteration
		@dBonB = (); # Relative difference between input and output slope
		
		# Do each iteration
		for ($i=1;$i<=$ITERATIONS;$i++) {
		    
		    ### Run the program on the data with CRs
		    
		    # Make the data
		    chdir $CREATEDIR;
		    print "$CREATE $CREATESEEING $seeing\n";
		    system "$CREATE $CREATESEEING $seeing";
		    system "mv -f test_[0-3]\{.fits,.fits.map,.cr\} $WORKINGDIR";
		    chdir $WORKINGDIR;
		    
		    # Run the program and count the number of CRs
		    print "$RUN $RUNREJECT $reject $RUNFRAC $frac $RUNGRAD $grad |\n";
		    open STAC, "$RUN $RUNREJECT $reject $RUNFRAC $frac $RUNGRAD $grad |";
		    $numCRs = 0;	# Number of CRs masked
		    while (<STAC>) {
			print;
			if (/(\d+) pixels masked in image \d+/) {
			    $numCRs += $1;
			}
		    }
		    close STAC;
		    printf "---> %d pixels masked in total.\n", $numCRs;

		    my @trueFlux = (); # Flux of pixels that are real CRs
		    my @falseFlux = ();	# Flux of pixels that were masked, but not CRs
		    for (my $i = 0; $i < 4; $i++) {
			# Read list of CRs injected into image
			open CR, "test_$i.cr";
			my %originalCRs = ();
			while (<CR>) {
			    my ($x,$y,$flux) = split;
			    $originalCRs{"$x,$y"} = $flux;
			}
			close CR;

			# Read list of CRs masked by program
			open CR, "test_$i.fits.cr";
			while (<CR>) {
			    my ($x,$y,$junk,$flux,$frac,$grad) = split;
			    # See if it was injected
			    if (exists($originalCRs{"$x,$y"})) {
				push @trueFlux, $flux;
			    } else {
				push @falseFlux, $flux;
			    }
			}
		    }

		    # Number of true and false CRs
		    push @trueRejected, scalar @trueFlux;
		    push @falseRejected, scalar @falseFlux;

		    print "True: " . scalar @trueFlux . "   False: " . scalar @falseFlux . "\n";

		    # Get input slope
		    system "cat test_[0-3].cr | sort -k 3 -n | awk \'{num++; print num,\$3;}\' > test_in.cr";
		    open FIT, "guessline.pl test_in.cr |";
		    while (<FIT>) {
			if (/b = (\S+)/) {
			    $inSlope = $1;
			    last;
			}
		    }
		    close FIT;

		    # Get output slope
		    open CR, "| sort -n | awk \'{num++; print num,\$1;}\' > test_out.cr";
		    foreach $flux (@trueFlux) {
			print CR $flux . "\n";
		    }
		    close CR;
		    
		    open FIT, "guessline.pl test_out.cr |";
		    while (<FIT>) {
			if (/b = (\S+)/) {
			    $outSlope = $1;
			    last;
			}
		    }
		    close FIT;
		    print "inSlope: $inSlope\noutSlope: $outSlope\n";

		    push @dBonB, abs($outSlope - $inSlope) / $inSlope;
		    
		}
		
		# Calculate the mean
		$trueMean = mean(\@trueRejected);
		$falseMean = mean(\@falseRejected);
		$trueMedian = median(\@trueRejected);
		$falseMedian = median(\@falseRejected);
		$dBonBMean = mean(\@dBonB);
		$dBonBMedian = median(\@dBonB);
		
		printf "$seeing\t$reject\t$frac\t$grad\t%.1f\t%.1f\t%.1f\t%.1f\t%.3f\t%.3f\n",$trueMean,$trueMedian,$falseMean,$falseMedian,$dBonBMean,$dBonBMedian;
		printf OUTPUT "$seeing\t$reject\t$frac\t$grad\t%.1f\t%.1f\t%.1f\t%.1f\t%.3f\t%.3f\n",$trueMean,$trueMedian,$falseMean,$falseMedian,$dBonBMean,$dBonBMedian;
	    }
	}
    }
}
close OUTPUT;




# Return the mean
sub mean
{
    my ($arrayRef) = @_;

    my @array = @$arrayRef;
    my $mean = 0;
    foreach my $value (@array) {
	$mean += $value;
    }
    return $mean / (scalar @array);
}

# Return the median
sub median
{
    my ($arrayRef) = @_;
    my @array = @$arrayRef;
    # Insertion sort
    my @sort = ();
    foreach my $value (@array) {
	for ($i=0;$i<=$#sort;$i++) {
	    if ($sort[$i] < $value) {
		last;
	    }
	}
	splice @sort, $i, 0, $value;
    }
    # Median is the middle one.
    if (scalar @sort % 2 == 1) {
	return $sort[(scalar @sort  - 1) / 2];
    } else {
	return ($sort[(scalar @sort)/2] + $sort[(scalar @sort)/2 - 1])/2;
    }
}
