#!/usr/bin/env perl

use warnings;
use strict;

## report the program and machine
use Sys::Hostname;
my $host = hostname();
print "\n\n";
print "Starting script $0 on $host\n\n";

use vars qw( $VERSION );
$VERSION = '0.01';

use IPC::Cmd qw( can_run run );
use PS::IPP::Metadata::Config;
use PS::IPP::Metadata::List qw( parse_md_list );
use Statistics::Descriptive;

use PS::IPP::Config qw(
    $PS_EXIT_SUCCESS
    $PS_EXIT_UNKNOWN_ERROR
    $PS_EXIT_SYS_ERROR
    $PS_EXIT_CONFIG_ERROR
    $PS_EXIT_PROG_ERROR
    $PS_EXIT_DATA_ERROR
    $PS_EXIT_TIMEOUT_ERROR
    );
my $ipprc = PS::IPP::Config->new(); # IPP configuration

use Getopt::Long qw( GetOptions :config auto_help auto_version gnu_getopt );
use Pod::Usage qw( pod2usage );

my ($det_id, $iter, $det_type, $camera, $filter, $workdir, $dbname, $no_update, $no_op);
GetOptions(
	   'det_id|d=s'        => \$det_id,
	   'iteration=s'       => \$iter,
	   'det_type|t=s'      => \$det_type,
	   'camera=s'          => \$camera,
	   'filter=s'          => \$filter,
	   'workdir|w=s'       => \$workdir, # Working directory for output files
	   'dbname|d=s'        => \$dbname, # Database name
	   'no-update'         => \$no_update,
	   'no-op'             => \$no_op,
) or pod2usage( 2 );

pod2usage( -msg => "Unknown option: @ARGV", -exitval => 2 ) if @ARGV;
pod2usage(
    -msg => "Required options: --det_id --iteration --det_type --camera",
    -exitval => 3,
) unless defined $det_id
    and defined $iter
    and defined $det_type
    and defined $camera;


# Look for programs we need
my $missing_tools;
my $dettool = can_run('dettool') or (warn "Can't find dettool" and $missing_tools = 1);
if ($missing_tools) { 
    warn("Can't find required tools.");
    exit($PS_EXIT_CONFIG_ERROR); 
}

my $mdcParser = PS::IPP::Metadata::Config->new;	# Parser for metadata config files

# Get list of component files
my $exposures;			# Array of exposures
{
    my $command = "$dettool -residexp -det_id $det_id -iteration $iter"; # Command to run
    $command .= " -dbname $dbname" if defined $dbname;
    my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	run(command => $command, verbose => 1);
    unless ($success) {
	$error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	&my_die("Unable to perform dettool -residexp: $error_code", $det_id, $iter, $error_code);
    }

    my $metadata = $mdcParser->parse(join "", @$stdout_buf) or
	&my_die("Unable to parse metadata config doc", $det_id, $iter, $PS_EXIT_PROG_ERROR);
    # XXX need to test that the result list is 0 length and return with a different error
    $exposures = parse_md_list($metadata) or
	&my_die("Unable to parse metadata list", $det_id, $iter, $PS_EXIT_PROG_ERROR);
}

my @expTags;			# Array of exposure IDs
my @means;			# Array of means
my @stdevs;			# Array of standard deviations
my @meanStdevs;			# Array of standard deviations of the mean (normalised by the mean)
my @variances;			# Array of variances
my @accept;			# Array of accept flags
my @include;			# Array of include flags
foreach my $exposure (@$exposures) {
    &my_die("Unable to find exposure id.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{exp_tag};
    &my_die("Unable to find mean.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg};
    &my_die("Unable to find stdev.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg_stdev};
    &my_die("Unable to find mean stdev.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{bg_mean_stdev};
    &my_die("Unable to find accept.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{accept};
    &my_die("Unable to find include.\n", $det_id, $iter, $PS_EXIT_SYS_ERROR) unless defined $exposure->{include};
    push @expTags, $exposure->{exp_tag};
    push @means, $exposure->{bg};
    push @stdevs, $exposure->{bg_stdev};
    # XXX why are we using mean_stdev / bd?
    # push @meanStdevs, $exposure->{bg_mean_stdev} / $exposure->{bg};
    push @meanStdevs, $exposure->{bg_mean_stdev};
    push @variances, $exposure->{bg_stdev}**2;
    push @accept, $exposure->{accept};
    push @include, $exposure->{include};
}
my $meanStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
$meanStats->add_data(@means);
my $stdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
$stdevStats->add_data(@stdevs);
my $meanStdevStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
$meanStdevStats->add_data(@meanStdevs);

$ipprc->define_camera($camera);
# Rejection thresholds
my $reject_mean = rejection_limit( 'ENSEMBLE.MEAN', $det_type, $filter );
my $reject_stdev = rejection_limit( 'ENSEMBLE.STDEV', $det_type, $filter );
my $reject_meanstdev = rejection_limit( 'ENSEMBLE.MEANSTDEV', $det_type, $filter );

print "Ensemble mean " . $meanStats->mean() . " +/- " . $meanStats->standard_deviation .
    ", stdev " . $stdevStats->mean() . " +/- " . $stdevStats->standard_deviation() . "\n\n";

# Go through again to do rejection, and update the database for each exposure
my $numChanges = 0;		# Number of exposures with changed status
my $numReject = 0;		# Number of exposures rejected
my $command;
my $reject;
my $expTag;
for (my $i = 0; $i < scalar @means; $i++) {
    $expTag = $expTags[$i];	# Exposure ID
    $command = "$dettool -updateresidexp -det_id $det_id -iteration $iter -exp_tag $expTag"; # Command to run
    $command .= " -dbname $dbname" if defined $dbname;
    
    $reject = 0;		# Reject this exposure?

    if (not $accept[$i]) {
	# Rejected this at an earlier stage
	print "Rejecting $expTag based on earlier determination.\n";
	$reject = 1;
	goto UPDATE;
    }

    # Cop-out if we're not doing any operations
    if ($no_op) {
	# Make sure something gets rejected (just once!), just so that
	# we can trace the full range of the workflow
	if ($i == 0 and $iter == 0) {
	    $reject = 1;
	}
	goto UPDATE;
    }
    
    if ($reject_mean and defined $meanStats->standard_deviation() ) {
	my $dMean = abs($means[$i] - $meanStats->mean()) ;
	if ($dMean > ($reject_mean * $meanStats->standard_deviation())) {
	    print "Rejecting $expTag based on outlier mean value: " .
		"$means[$i] --> $dMean vs " . $reject_mean * $meanStats->standard_deviation() . "\n";
	    $reject = 1;
	    goto UPDATE;
	}
    } else {
	print "no rejection for exposure mean\n";
    }

    if ($reject_stdev and defined $stdevStats->standard_deviation() > 0) {
	my $dMean = abs($stdevs[$i] - $stdevStats->mean());
	if ($dMean > ($reject_stdev * $stdevStats->standard_deviation())) {
	    print "Rejecting $expTag based on outlier stdev: " .
		$stdevs[$i] . " --> $dMean sigma vs " . $reject_stdev . "\n";
	    $reject = 1;
	    goto UPDATE;
	}
    } else {
	print "no rejection for exposure stdev\n";
    }
    
    if ($reject_meanstdev and defined $meanStdevStats->standard_deviation() > 0) {
	my $dMean = abs($meanStdevs[$i] - $meanStdevStats->mean());
	if ($dMean > ($reject_meanstdev * $meanStdevStats->standard_deviation())) {
	    print "Rejecting $expTag based on outlier mean stdev: " .
		$meanStdevs[$i] . " --> $dMean sigma vs " . $reject_meanstdev . "\n";
	    $reject = 1;
	    goto UPDATE;
	}
    } else {
	print "no rejection for exposure mean stdev\n";
    }

  UPDATE:
    if ($reject) {
	$command .= ' -reject';
	$numReject++;
    }
    
    # Check for status changes
    if ((not $include[$i] and not $reject) or ($include[$i] and $reject)) {
	print "Status of $expTag has changed.\n";
	$numChanges++;
    }

    unless ($no_update) {
	# Update 
	my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	    run(command => $command, verbose => 1);
	unless ($success) {
	    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	    &my_die("Unable to perform dettool -updateresidexp: $error_code", $det_id, $iter, $error_code);
	}
    }
}

# Decide if the current is sufficient to use as a master, and if we can stop iterating
my $master = 1;			# This is good enough for a master
my $stop = 1;			# Stop iterating

if ($numChanges > 0) {
    $master = 0;
    $stop = 0;
}

# Rejecting everything --- stop before something bad happens!
if ($numReject == scalar @means) {
    $master = 0;
    $stop = 1;
}

print "Master: $master\n";
print "Stop: $stop\n";

# Put results into the database
unless ($no_update) {

    {
	# Add the summary
	my $varianceStats = Statistics::Descriptive::Sparse->new(); # Statistics calculator
	$varianceStats->add_data(@variances);
	
	my $bg = ($meanStats->mean() or 'NAN');
	my $bg_stdev = (sqrt( $varianceStats->mean() ) or 'NAN');
	my $bg_mean_stdev = ($meanStats->standard_deviation() or 'NAN');
	
	my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter";
	$command .= " -bg $bg -bg_stdev $bg_stdev -bg_mean_stdev $bg_mean_stdev";
	$command .= " -accept" if $master;
	$command .= " -dbname $dbname" if defined $dbname;
	
	my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	    run(command => $command, verbose => 1);
	unless ($success) {
	    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	    warn("Unable to perform dettool -adddetrunsummary: $error_code");
	    exit($error_code);
	}
    }

    # Re-run processing if required
    {
	my $command = "$dettool -updatedetrun -det_id $det_id";
	if ($stop) {
	    $command .= ' -state stop';
	} else {
	    $command .= ' -again';
	}
	
	$command .= " -dbname $dbname" if defined $dbname;
	
	my ( $success, $error_code, $full_buf, $stdout_buf, $stderr_buf ) =
	    run(command => $command, verbose => 1);
	unless ($success) {
	    $error_code = (($error_code >> 8) or $PS_EXIT_PROG_ERROR);
	    warn("Unable to perform dettool -updatedetrun: $error_code");
	    exit($error_code);
	}
    }
}

END {
    my $status = $?;
    system("sync") == 0
        or die "failed to execute sync: $!" ;
    $? = $status;
}


sub my_die
{
    my $msg = shift; # Warning message on die
    my $det_id = shift;		# Detrend identifier
    my $iter = shift;		# Iteration
    my $exit_code = shift; # Exit code to add

    warn($msg);
    if ($det_id and $iter and not $no_update) {
	my $command = "$dettool -adddetrunsummary -det_id $det_id -iteration $iter -code $exit_code";
	$command .= " -dbname $dbname" if defined $dbname;
        system ($command);
    }
    exit $exit_code;
}



# Retrieve the requested rejection limit, dying unless extant
sub rejection_limit
{
    my $name = shift;		# Rejection limit to 
    my $type = shift;		# Type of exposure
    my $filter = shift;		# Filter

    my $value = $ipprc->rejection( $name, $det_type, $filter );
    if (not defined $value) {
	$filter = "(no filter)" unless defined $filter;
	die "Unable to determine $name rejection limit for $det_type with $filter.\n";
    }

    return $value;
}


__END__
