#!/usr/bin/perl -w

# This script reads Aquarius house-keeping telemetry files and
# computes statistics for each parameter (mean, median, minimum,
# maximum, count) over a fixed time interval in order to reduce
# the volume of the data so that plots may be generated more quickly.

# Norman Kuring		20-May-2011	Initial development
# Norman Kuring		 1-Jun-2011	Continued development
# Norman Kuring		 2-Jun-2011	Write output files to temporary
#					location, and only move them to their
#					final destination when they are
#					complete.
# Norman Kuring		15-Jun-2011	Also condense stored and real-time
#					Service Platform telemetry files.
# Norman Kuring		15-Jul-2011	Filter out data with bogus time stamps.
# Norman Kuring		15-Jul-2011	Apply leap-second correction to GPS
#					time stamps.
# Norman Kuring		20-Jul-2011	Fix logic to allow data from 3-day
#					stored telemetry files from CONAE and
#					to condense right up through the most
#					recent data available.
# Norman Kuring		 5-Aug-2011	Don't compare to the earliest day
#					condensed so far if nothing has been
#					condensed yet.
# Norman Kuring		29-Aug-2011	Skip L1AQL files just like L1A files.
# Norman Kuring		16-Nov-2011	Gather daily medians, extrema, counts.
# Norman Kuring		 7-Dec-2011	Change the way I determine earliest and
#					latest days condensed so far.
# Norman Kuring		13-Jun-2012	Add mean to the statistics.

use strict;
use Time::Local;
use File::Path;
use lib '/web/aquarius/scripts';
use GpsToUtc qw(&read_leap_secs &gps2utc);

$| = 1;		# ... so I can see log output as it happens.

# Initialize the GPS to UTC conversion function.
read_leap_secs;

my $mission_start = 991_754_785;	# 10 June 2011 15:26:10 UTC
my $now = $^T - timegm(0,0,0,6,0,80);	# I'm ignoring leap seconds when I set
					# the current time to the GPS epoch.

my %dir = (
AQ => '/web/data/Aquarius/DOWNLINKS/L1',
SP => '/web/data/Aquarius/HKT',
);

my %pattern = (
AQ => '^Q\d{8}_\d{6}_T\d{8}_\d{6}_\d{8}_\d{6}\.AQ_\w+\.txt$',
SP => '\d{8}_\d{6}_\d{8}_\d{6}\.SP_\w+\.txt$',
);

my @types = keys %dir;

my $condensed_dir = '/web/aquarius/metadata/condensed_hkt';

my %ctgts;	# category/time-stamp hash
my %count;

# Make a temporary directory for writing intermediate files and go there.
# Also make a subroutine to remove the directory when it's no longer needed.
my $cleanup;
{
  my $tmpdir = (split /\//,$0)[-1];
  $tmpdir =~ s/\.pl$//;
  $tmpdir = "/web/aquarius/metadata/$tmpdir.$^T.$$";
  mkpath($tmpdir,0,0755) or die "Could not create $tmpdir ($!)";
  chdir $tmpdir or die "Could not cd to $tmpdir ($!)";
  $cleanup = sub {
    chdir();
    system("/bin/rm -r $tmpdir")==0 or warn "Could not remove $tmpdir ($!)"
  };
}

foreach my $type (@types){
  my $earliest_day = earliest_day($condensed_dir,$type);
  my $latest_day   =   latest_day($condensed_dir,$type);

  my @years = read_dir($dir{$type},'^\d{4}$');
  foreach my $year (@years){
    my @mmdd = read_dir("$dir{$type}/$year",'^\d{4}$');
    foreach my $mmdd (@mmdd){

      my $nominal_day = yearmmdd2day($year,$mmdd);

      # Assume that telemetry data are coming in roughly chronologically
      # at the day level and ignore uncondensed data that are more than
      # a few days earlier than the most recently condensed data.
      # Presumably, those data will have already been condensed.
      next if $nominal_day < $latest_day - 4;

      # Keep the category/time-stamp hash purged of older entries
      # to keep memory usage down.
      foreach my $ctg (keys %ctgts){
        foreach my $ts (keys %{ $ctgts{$ctg} }){
          delete $ctgts{$ctg}{$ts} if int($ts/86400) < $nominal_day - 5;
        }
      }

      # Generate output files for days that are complete.
      # Do not, however, condense data from the first day encountered
      # during this processing run unless it happens to be the same
      # day or earlier than the earliest day already condensed.  This
      # is to prevent the script from replacing earlier-written, complete
      # data with partial results.
      my %done = earlier_than($nominal_day);
      my @done = sort {$a<=>$b} keys %done;
      delete $done{$done[0]} if $earliest_day > 0
                             && defined $done[0]
                             && $done[0] > $earliest_day;
      condense_data($type,%done) if %done;

      # Process the input (uncondensed) files.
      print "[",scalar localtime,
            "] Reading data from $dir{$type}/$year/$mmdd ...";
      my @files = grep !/AQ_L1A(?:QL)?\.txt$/,
                  read_dir("$dir{$type}/$year/$mmdd",$pattern{$type});
      foreach my $f (@files){
        my ($ctg) = $f =~ /${type}_(\w+)\.txt$/;
        my $path = "$dir{$type}/$year/$mmdd/$f";
        open IN,'<',$path or quit("Could not open $path ($!)",$cleanup);
        my $rec = 0;
        while(<IN>){
          $rec++;
          my ($ts,@v) = split ' ';

          # Filter out data with bogus time stamps.
          next if $ts < $mission_start || $ts > $now;

          # Eliminate records with duplicate time stamps.
          next if $ctgts{$ctg}{$ts}++;

          # Perform a consistency check on the input data.
          if(defined $count{$ctg}){
            if(@v != $count{$ctg}){
              quit(
              "Inconsistent number of parameter values in files of category, " .
              "$ctg. The initial count was $count{$ctg}, but line, $rec, " .
              "of file, $path, contains ",scalar @v," parameter values.\n",
              $cleanup
              );
            }
          }
          else{
            $count{$ctg} = @v;
          }

          # Store record groupings in temporary data files.
          my $day = int( gps2utc($ts)/86400 );
          store($day,$ctg,@v);
        }
        close IN;
      }
      print " done.\n";
    }
  }
  # Generate output files for the remainder of the data.
  my %done = earlier_than(int($^T/86400) + 1);	# Earlier than tomorrow
  my @done = sort {$a<=>$b} keys %done;
  delete $done{$done[0]} if $done[0] != $earliest_day;
  condense_data($type,%done) if %done;
}

&$cleanup;

print "[",scalar localtime,"] $0 finished.\n";


my $current_file;

sub store{
  my ($day,$ctg,@v) = @_;
  my $f = "$day.$ctg";
  if( not defined $current_file ){
    $current_file = $f;
    open TMP,'>>',$f or quit("Error opening $f for appending ($!)",$cleanup);
  }
  elsif( $f ne $current_file ){
    close TMP;
    $current_file = $f;
    open TMP,'>>',$f or quit("Error opening $f for appending ($!)",$cleanup);
  }

  print TMP "@v\n";
}

sub earlier_than{
  my $d = shift;
  my %done;
  my @f = read_dir('.','^\d+\.\w+$');
  foreach my $f (@f){
    my ($day,$ctg) = $f =~ /^(\d+)\.(\w+)$/;
    next unless $day < $d;
    $done{$day}{$ctg}++;
  }
  %done;
}

sub condense_data{
  my ($type,%done) = @_;
  foreach my $day (keys %done){
    my ($yyyy,$ddd) = day2yyyyddd($day);
    my $outdir = "$condensed_dir/$yyyy/$ddd";
    print "[",scalar localtime,"] Condensing data to $outdir ...";
    unless(-d $outdir){
      mkpath($outdir,0,0755) or quit("Could not create $outdir ($!)",$cleanup);
    }

    # For each telemetry category...
    foreach my $ctg (keys %{ $done{$day} }){
      my $out = "Q$yyyy$ddd.day.${type}_$ctg.txt";

      open OUT,'>',$out or quit("Could not create $out ($!)",$cleanup);

      my (@mean,@median,@min,@max,@n);

      if( defined $done{$day}{$ctg} ){

        # Read the data from a temporary file created earlier,
        # and then delete that file.
        my $f = "$day.$ctg";
        open TMP,'<',$f or quit("Could not open $f ($!)",$cleanup);
        my $rec;
        my @data;
        push @data, [split ' ',$rec] while( $rec = <TMP> );
        close TMP;
        unlink $f;

        my $nr = @data;		# number of recs in this time period
        my $np = @{ $data[0] };	# number of params in this category

        # For each parameter...
        for(my $p=0; $p<$np; $p++){

          my @population;

          # For each record in the time period...
          for(my $r=0; $r<$nr; $r++){
            my $val = $data[$r][$p];
            if( $val =~ /^[+\-]?(?:\d+\.?\d*|\.\d+)$/ ){
              push @population, $val;
            }
          }

          # The population size for each parameter is stored in @n.
          $n[$p] = @population;

          if($n[$p] > 0){
            @population = sort {$a <=> $b} @population;
            my $sum = 0;
            grep {$sum += $_; 0} @population;
            $mean[$p] = $sum/$n[$p];
            my $mid = int($n[$p]/2);
            $median[$p] = $n[$p] % 2        ?
                          $population[$mid] :
                         ($population[$mid - 1] + $population[$mid])/2;
            $min[$p] = $population[0];
            $max[$p] = $population[-1];
          }
          else{
            $mean[$p] = $median[$p] = $min[$p] = $max[$p] = 'NaN';
          }
        }
        # Write out the five statistics for each
        # parameter for this time period.
        print OUT "@n\n@mean\n@median\n@min\n@max\n";
      }
      else{
        # Output five records indicating that
        # no data are available for this period.
        print OUT "N\nN\nN\nN\nN\n";
      }

      close OUT;

      rename $out,"$outdir/$out"
      or quit("Could not rename $out to $outdir/$out ($!)",$cleanup);
    }
    print " done.\n";
  }
}

sub latest_day{
  my ($dir,$type) = @_;
  my @years = reverse read_dir($dir,'^\d{4}$');
  foreach my $y (@years){
    my @days = reverse read_dir("$dir/$y",'^\d{3}$');
    foreach my $d (@days){
      if( read_dir("$dir/$y/$d",'^Q\d{7}\.day\.' . $type . '_\w+\.txt$') ){
        return timegm(0,0,0,1,0,$y-1900)/86400 + $d - 1;
      }
    }
  }
  return 0;
}

sub earliest_day{
  my ($dir,$type) = @_;
  my @years = read_dir($dir,'^\d{4}$');
  foreach my $y (@years){
    my @days = read_dir("$dir/$y",'^\d{3}$');
    foreach my $d (@days){
      if( read_dir("$dir/$y/$d",'^Q\d{7}\.day\.' . $type . '_\w+\.txt$') ){
        return timegm(0,0,0,1,0,$y-1900)/86400 + $d - 1;
      }
    }
  }
  return 0;
}

sub read_dir{
  my ($dir,$pattern) = @_;
  opendir D,$dir or quit("Could not open $dir ($!)",$cleanup);
  my @results = sort grep /$pattern/ && -s "$dir/$_", readdir D;
  closedir D;
  @results;
}

sub yearmmdd2day{
  my ($year,$mmdd) = @_;
  my ($mm,$dd) = unpack "a2a2",$mmdd;
  timegm(0,0,0,$dd,$mm-1,$year-1900)/86400;
}

sub day2yyyyddd{
  my $day = shift;
  my ($yyyy,$ddd) = (gmtime 86400*$day)[5,7];
  ( sprintf("%4d",$yyyy+1900), sprintf("%03d",$ddd+1) );
}

sub quit{
  my ($err,$cu) = @_;
  &$cu;
  die $err;
}
