#!/usr/bin/perl -w

# Mostly copyright Greg Baker (gregb@ifost.org.au) 2008. Portions
# used with permission by Peter Acklam.
#
# Keep the copyright attribution(s); other than that feel free to
# use Greg's work here in any way you want. If you want to tweak
# Peter Acklam's code you'll need to look at the (open) license for
# his stuff on CPAN.

=head1 NAME

observe

=head1 SYNOPSIS

C<observe [--datafile=...] --record var=value [ var=value ... ]>

C<observe [--datafile=...] --stats var [var ...]>

C<observe [--datafile=...] --probability var=value>

C<observe [--datafile=...] --trend var>

C<observe [--datafile=...] --predict var=value>

C<observe [--datafile=...] --list>

C<observe [--datafile=...] --version>

C<observe [--datafile=...] --graph var [var...]>

=head1 OPTIONS

=over

=item C<--datafile> ...

Choose a specific datafile to read from, overriding the environment
variable OBSERVE_DATAFILE.

If neither C<--datafile> nor C<OBSERVE_DATAFILE> are set, then:

=over

=item On MS-Windows, the default is F<.observed> in whatever directory B<observe> is running from.

=item On Unix, when running as 
if running as I<root> or as a user able to write to F</var/.observed>
then the default is F</var/.observed>

=item On Unix for other users, F<.observed> in the users home directory.

=back

To do. Other options need to be documented.

=cut


use strict;
use English;
use Fcntl qw{:flock :seek};
use Time::Local;
use File::Basename;
use POSIX qw(strftime);

my ($datafile);
if ($#ARGV > -1 and $ARGV[0] =~ m:^[-/ ]*datafile:) {
  if ($ARGV[0] =~ m:^[-/ ]*datafile *= *([^ ].*)$:) {
    $datafile = $1;
    shift @ARGV;
  } elsif ($ARGV[0] =~ m:^[-/ ]*datafile *$:) {
    $datafile = $ARGV[1];
    shift @ARGV;
    shift @ARGV;
  } else {
    die "Datafile argument #1 ($ARGV[0]) incomprehensible. Stopped";
  }
} elsif (exists($ENV{OBSERVE_DATAFILE})) {
 $datafile = $ENV{OBSERVE_DATAFILE};
} elsif ($OSNAME =~ /^MSWin/) {
  # I can't figure out how to do anything better than this.
  $datafile = ".observed";
} elsif ($> == 0 || -w '/var/.observed') {
 $datafile = '/var/.observed';
} elsif (exists ($ENV{HOME})) {
 # how else to figure out what my home directory is? getpwent?
 $datafile = $ENV{HOME} . "/.observed";
} else  {
 $datafile = ".observed";
}
my ($arg);
my ($line);
my ($timestamp);

sub stats {
 my ($arg) = shift;
 my ($timestamp);
 my ($timestampname);
 my (@fields);
 my (%varvals);
 my ($count) = 0;
 my ($sum) = 0.0;
 my ($sum_of_squares) = 0.0;
 my ($most_recent);
 my $minimum = undef;
 my $maximum = undef;
 my $earliest_observation_time = undef;
 my $latest_observation_time = undef;
 open(DATAFILE,$datafile) || die "Can't read data";
 flock(DATAFILE,LOCK_SH);
 while (<DATAFILE>) {
   chomp;
   ($timestampname,$timestamp,@fields) = split(/\t/);
   %varvals = @fields;
   if (exists $varvals{$arg}) {
     $count++;
     $sum += $varvals{$arg};
     $sum_of_squares += ($varvals{$arg} ** 2);
     $most_recent = $varvals{$arg};
     $minimum = $varvals{$arg} if !defined $minimum || $minimum > $varvals{$arg};
     $maximum = $varvals{$arg} if !defined $maximum || $maximum < $varvals{$arg};
     $earliest_observation_time = $timestamp if !defined $earliest_observation_time || $earliest_observation_time > $timestamp;
     $latest_observation_time = $timestamp if !defined $latest_observation_time || $latest_observation_time < $timestamp;
   }
 }
 flock(DATAFILE,LOCK_UN);
 close(DATAFILE);
 return (0,undef,undef,undef) if ($count == 0);
 return ($count,($sum / $count),sqrt(($sum_of_squares - ($sum ** 2 /$count))  /$count),$most_recent,$minimum,$maximum,$earliest_observation_time,$latest_observation_time);
}

sub varnames {
  my %known_varnames;
  my $varname;
  open(DATAFILE,$datafile) || die "Can't read datafile $datafile";
  flock(DATAFILE,LOCK_SH);
  while (<DATAFILE>) {
    chomp;
    my ($timestampname,$timestamp,@fields) = split(/\t/);
    my %varvals = @fields;
    foreach $varname (keys %varvals) { $known_varnames{$varname} = 1; }
  }
  flock(DATAFILE,LOCK_UN);
  close(DATAFILE);
  return keys %known_varnames;
}

sub phi {
 my ($arg) = shift;
 my ($erfarg) = $arg / sqrt(2.0);
 my ($erfresult) = Math::SpecFun::Erf::erf($erfarg ); # See below
 my ($answer) = 0.5 * (1.0 + $erfresult);
 return $answer;
}
 
sub linear_regression {
  my ($independent) = shift;
  my ($dependent) = shift;
  my (@fields);
  my (%varvals);
  my ($independent_sum) = 0.0;
  my ($dependent_sum) = 0.0;
  my ($count) = 0;
  open (DATAFILE,$datafile) || die "Can't read data";
  flock(DATAFILE,LOCK_SH);
  while (<DATAFILE>) {
    chomp;
    @fields = split(/\t/);
    %varvals = @fields;
    if (exists($varvals{$independent}) && exists($varvals{$dependent})) {
      $count++;
      $independent_sum += $varvals{$independent};
      $dependent_sum += $varvals{$dependent};
    }
  }
  seek(DATAFILE,0,SEEK_SET);
  my ($dependent_mean) = $dependent_sum / $count;
  my ($independent_mean) = $independent_sum / $count;
  my ($slope_sum) = 0.0;
  my ($slope_divisor) = 0.0;
  while (<DATAFILE>) {
    chomp;
    @fields = split(/\t/);
    %varvals = @fields;
    if (exists($varvals{$independent}) && exists($varvals{$dependent})) {
      $slope_sum += ($varvals{$independent} - $independent_mean) *
                    ($varvals{$dependent} - $dependent_mean);
      $slope_divisor += ($varvals{$independent} - $independent_mean) ** 2;
    }  
  }
  flock(DATAFILE,LOCK_UN);
  close(DATAFILE);
  my ($slope) = $slope_sum / $slope_divisor;
  my ($intercept) = $dependent_mean - $slope * $independent_mean;
  return ($slope,$intercept);
}

sub graph {
  my ($data_point_count) = shift;
  my ($units_of_time) = shift;
  my ($time_units) = shift;
  my ($start_time) = shift;
  my ($end_time) = shift;
  my (@varnames_to_plot) = @_;
  my %durations = ('second' => 1, 'minute' => 60, 'hour' => 3600,
		   'day' => 86400, 'week' => '604800', 'year' => 31536000);
  my %format_strings = ('second' => '%H:%M:%S', 'minute' => '%H:%M:%S',
			'hour' => '%a %e %H:%M',
			'day' => '%a %b %e',
			'week' => '%e %b %Y',
			'month' => '%e %b %Y',
			'year' => '%b %Y');
  # I should add another layer below the x legend for second hour and day
  # to give the month
  my $v;
  if ($time_units eq 'epoch') { 
      $start_time = undef; 
  } elsif ($time_units eq 'month') { 
      $start_time = MiniDateManip::subtract_months($end_time,$units_of_time);
  } elsif (exists $durations{$time_units}) {
    $start_time = $end_time - $units_of_time * $durations{$time_units};
  } else {
    die "Confused by time unit $time_units";
  }

  my %summaries;
  foreach $v (@varnames_to_plot) {
      $summaries{$v} = new DataSummary($start_time,$end_time,$data_point_count);
  }

  my $have_data = 0;
  open(DATAFILE,$datafile) || die "Can't read datafile $datafile";
  flock(DATAFILE,LOCK_SH);
  my @legend_texts;
 DATAFILE_LINE:
  while (<DATAFILE>) {
    chomp;
    my ($timestampname,$timestamp,@fields) = split(/\t/);
    my %varvals = @fields;
    my $relevant = 0;
    next DATAFILE_LINE unless grep(defined $varvals{$_} ,@varnames_to_plot);
    if (!defined $start_time) { 
	$start_time = $timestamp; 
	foreach $v (values %summaries) { $v->set_start_time($start_time); }
    }
    next DATAFILE_LINE if $timestamp < $start_time;

    foreach $v (@varnames_to_plot) {
      $summaries{$v}->store($timestamp,$varvals{$v});
    }
    $have_data = 1;
  }
  flock(DATAFILE,LOCK_UN);
  close(DATAFILE);
  return "chd=s:_" unless $have_data;

  my $subremum = undef;
  my $supremum = undef;
  foreach $v (values %summaries) {
      my $m = $v->maximum_anywhere();
      $supremum = $m if !defined $supremum or $m > $supremum;
      $m = $v->minimum_anywhere();
      $subremum = $m if !defined $subremum or $m < $subremum;
  }

  foreach $v (@varnames_to_plot) {
      push(@legend_texts,html_encode($v));
  }

  ($subremum,$supremum) = SimplifyRange::min_max_pair($subremum,$supremum);
  my $encoder = new GoogleChart::SimpleEncoder($subremum,$supremum);
  my $step_remum = SimplifyRange::increment_around($supremum-$subremum);
  my @y_axis_labels;
  my $remum;
  for($remum = $subremum;$remum <= $supremum; $remum += $step_remum) {
      push(@y_axis_labels,$remum);
  }
  my $y_axis_string = join("|",@y_axis_labels);

  my @datastring;
  foreach $v (values %summaries) {
      @datastring = (@datastring,$v->google_chart($encoder));
  }
  
  if (!defined $time_units or $time_units eq 'epoch') {
      my $tu;
      foreach $tu (qw{second minute hour day week year}) {
	  if ($end_time - $start_time > $durations{$tu} * 6) {
	      $time_units = $tu;
	  }
      }
  }
  my @times;
  my $t;
  my $time_step = $durations{$time_units};
  while ($end_time - $start_time > 8 * $time_step) {
      $time_step += $durations{$time_units};
      if ($time_units eq 'day' && $time_step >= 4 * $durations{'day'}) {
	  $time_step = $durations{'week'};
      }
  }
  for($t=$start_time;$t<=$end_time;$t+=$time_step) {
  #  for($t=$start_time;$t<=$end_time;$t+=$end_time-$start_time) {
      my @l = localtime $t;
      my $s = strftime($format_strings{$time_units},@l);
      push(@times,html_encode($s));
  }
  my $time_string = join("|",@times);
  my $legend_string = join("|",@legend_texts);
  my $chart_title = File::Basename::basename($datafile);
  $chart_title =~ s/^\.//;
  if ($chart_title =~ /observed/ && $#legend_texts == 0) {
      $chart_title = $legend_string;
      $legend_string = ""; 
  } else {
      $legend_string = 'chdl=' . $legend_string. "&";
  }
  $chart_title = html_encode($chart_title);
  return 'chd=s:'.join(",",@datastring)."&chxt=x,y&chxl=0:|$time_string|1:|$y_axis_string&${legend_string}chtt=$chart_title";
}


sub html_encode {
    my $x = shift;
    $x =~ s/%/%25/g;
    $x =~ s/\+/%2b/g;
    $x =~ s/ /+/g;
    return $x;
}


sub usage {
 die "Usage: $0 [datafile=...] record|stats|probability|trend|predict|list|version|graph";
}
&usage() unless ($#ARGV > -1);

if ($ARGV[0] =~ m:^[- /]*record:i) {
 $timestamp = scalar(time());
 $line = "TIMESTAMP\t$timestamp";
 foreach $arg (@ARGV[1..$#ARGV]) {
   die "Confused by $arg" unless $arg =~ /([^=\t]*) *= *([\d.]+)/;
   $line .= "\t" . $1 . "\t" . $2;
 }
 # probably should lock here, or something
 open(DATAFILE,">>$datafile") || die "Can't record data.";
 flock(DATAFILE,LOCK_EX);
 print DATAFILE "$line\n";
 flock(DATAFILE,LOCK_UN);
 close(DATAFILE);
 exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*stats:i) {
 foreach $arg (@ARGV[1..$#ARGV]) {
   my ($count,$mean,$stddev,$most_recent,$minimum,$maximum,$earliest_time,$latest_time) = stats($arg);
   print "count($arg) = $count\n";
   print "mean($arg) = $mean\n" if defined $mean;
   print "stddev($arg) = $stddev\n" if defined $stddev;
   print "latest($arg) = $most_recent\n" if defined $most_recent;
   print "minimum($arg) = $minimum\n" if defined $minimum;
   print "maximum($arg) = $maximum\n" if defined $maximum;
 }
 exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*probability:i) {
 # Print how many standard deviations away from normal this is.
 die "Don't understand result" unless $ARGV[1] =~ /([^=\t]*) *= *([\d.]+)/;
 my ($variable) = $1;
 my ($val) = $2;
 my ($count,$mean,$stddev,$most_recent,$minimum,$maximum,$earliest_time,$latest_time) = stats($variable);
 my ($stddevs_away) = (($val - $mean) / $stddev);
 my ($probability) = abs(Math::SpecFun::Erf::erf($stddevs_away));
 print "Probability of a measurement closer to the mean is $probability.\n";
 my ($nice_prob) = (1-$probability) * 100.0;
 my ($displayed_prob);
 my ($i);
 for ($i=2;$i<8;$i++)  {
  $displayed_prob = sprintf ("%.${i}f",$nice_prob);
  if ($displayed_prob != 0.0) { last; }
 }
 print "$displayed_prob% chance event.\n";
 exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*trend:i) {
 my ($arg) = $ARGV[1];
 my ($slope,$intercept) = linear_regression('TIMESTAMP',$arg);
 print "Formula: $arg = $slope * (time in seconds since epoch) + $intercept\n";
 print "$arg is increasing at ".($slope * 86400). " units per day\n";
 print "$arg is increasing at ".($slope * 7 * 86400). " units per week\n";
 print "$arg is increasing at ".($slope * 365 * 86400). " units per year\n";
 exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*predict:i) {
 die "Didn't understand $ARGV[1]" unless $ARGV[1] =~ /([^=\t]*) *= *([\d.]+)/;
 my ($arg) = $1;
 my ($target) = $2;
 my ($slope,$intercept) = linear_regression($arg,'TIMESTAMP');
 my ($predicted_time) = $slope * $target + $intercept;
 my ($ltime) = scalar(localtime($predicted_time));
 my ($now) = scalar(time);
 my ($soon) = $predicted_time - $now;
 print "$arg is expected to be $target at $ltime\n";
 print "That is $soon seconds from now\n";
 print "That is ".($soon / 86400)." days from now\n";
 print "That is ".($soon / (86400*7))." weeks from now\n";
 print "That is ".($soon / (86400*365))." years from now\n";
 exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*list:i) {
  my @varnames = &varnames();
  my $varname;
  foreach $varname (@varnames) {
    print "$varname\n";
  }
  exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*graph:i) {
  shift(@ARGV);
  &usage() unless $#ARGV > -1;
  my $units_of_time = 1;
  if ($ARGV[0] =~ /^ *([\d.]+) *$/) { $units_of_time = $1; shift(@ARGV); }
  my $time_units = "epoch";
 TIME_PERIOD:
  foreach my $period (qw{second minute hour day week month year}) {
    if ($ARGV[0] =~ /^ *${period}s? *$/i) {
      $time_units = $period;
      shift(@ARGV);
      last TIME_PERIOD;
    }
  }
  &usage() unless $#ARGV > -1;
  my $variable;
  my $number_of_points_to_plot = 0;
  my $very_earliest_time = undef;
  my $very_latest_time = undef;

  foreach $variable(@ARGV) {
   my ($count,$mean,$stddev,$most_recent,$minimum,$maximum,$earliest_time,$latest_time) = stats($variable);
   if ($count > $number_of_points_to_plot) {
      $number_of_points_to_plot = $count;
   }
   if (!defined $very_earliest_time  || $very_earliest_time > $earliest_time) {
       $very_earliest_time = $earliest_time;
   }

   if (!defined $very_latest_time  || $very_latest_time < $latest_time) {
       $very_latest_time = $latest_time;
   }
  }
  if ($number_of_points_to_plot > 400) {
   $number_of_points_to_plot = 400;
   # no point in plotting more points than we have pixels on screen
  } 
  my $google_chart_data = &graph($number_of_points_to_plot,$units_of_time,$time_units,$very_earliest_time,$very_latest_time,@ARGV);
  print "http://chart.apis.google.com/chart?chs=600x400&cht=lc&${google_chart_data}\n";
  exit(0);
} elsif ($ARGV[0] =~ m:^[- /]*version:i) {
 print '$Id: observe 171 2008-10-03 14:17:44Z gregb $';
 print "\n";
} else {
  &usage();
}


#####################################################################
#####################################################################
# The rest of this file is just statistical math code so we can calculate
# probabilities. Thank you to  Peter J. Acklam.

# Author:      Peter J. Acklam
# Time-stamp:  2004-01-21 11:42:33 +0100
# E-mail:      pjacklam@online.no
# URL:         http://home.online.no/~pjacklam

package Math::SpecFun::Erf;
require 5.000;
require Exporter;

use strict;
use vars qw($VERSION @ISA @EXPORT_OK %EXPORT_TAGS);

$VERSION = '0.02';
@ISA = qw(Exporter);
@EXPORT_OK = qw(erf erfc erfcx erfinv erfcinv erfcxinv);
%EXPORT_TAGS = ( all => [ @EXPORT_OK ] );

########################################################################
## Internal functions.
########################################################################

sub calerf {
    my ($arg, $result, $jint) = @_;
    local $[ = 1;
#------------------------------------------------------------------
#
# This packet evaluates  erf(x),  erfc(x),  and  exp(x*x)*erfc(x)
#   for a real argument  x.  It contains three FUNCTION type
#   subprograms: ERF, ERFC, and ERFCX (or DERF, DERFC, and DERFCX),
#   and one SUBROUTINE type subprogram, CALERF.  The calling
#   statements for the primary entries are:
#
#                   Y=ERF(X)     (or   Y=DERF(X)),
#
#                   Y=ERFC(X)    (or   Y=DERFC(X)),
#   and
#                   Y=ERFCX(X)   (or   Y=DERFCX(X)).
#
#   The routine  CALERF  is intended for internal packet use only,
#   all computations within the packet being concentrated in this
#   routine.  The function subprograms invoke  CALERF  with the
#   statement
#
#          CALL CALERF(ARG,RESULT,JINT)
#
#   where the parameter usage is as follows
#
#      Function                     Parameters for CALERF
#       call              ARG                  Result          JINT
#
#     ERF(ARG)      ANY REAL ARGUMENT         ERF(ARG)          0
#     ERFC(ARG)     ABS(ARG) < XBIG           ERFC(ARG)         1
#     ERFCX(ARG)    XNEG < ARG < XMAX         ERFCX(ARG)        2
#
#   The main computation evaluates near-minimax approximations
#   from "Rational Chebyshev approximations for the error function"
#   by W. J. Cody, Math. Comp., 1969, PP. 631-638.  This
#   transportable program uses rational functions that theoretically
#   approximate  erf(x)  and  erfc(x)  to at least 18 significant
#   decimal digits.  The accuracy achieved depends on the arithmetic
#   system, the compiler, the intrinsic functions, and proper
#   selection of the machine-dependent constants.
#
#*******************************************************************
#*******************************************************************
#
# Explanation of machine-dependent constants
#
#   XMIN   = the smallest positive floating-point number.
#   XINF   = the largest positive finite floating-point number.
#   XNEG   = the largest negative argument acceptable to ERFCX;
#            the negative of the solution to the equation
#            2*exp(x*x) = XINF.
#   XSMALL = argument below which erf(x) may be represented by
#            2*x/sqrt(pi)  and above which  x*x  will not underflow.
#            A conservative value is the largest machine number X
#            such that   1.0 + X = 1.0   to machine precision.
#   XBIG   = largest argument acceptable to ERFC;  solution to
#            the equation:  W(x) * (1-0.5/x**2) = XMIN,  where
#            W(x) = exp(-x*x)/[x*sqrt(pi)].
#   XHUGE  = argument above which  1.0 - 1/(2*x*x) = 1.0  to
#            machine precision.  A conservative value is
#            1/[2*sqrt(XSMALL)]
#   XMAX   = largest acceptable argument to ERFCX; the minimum
#            of XINF and 1/[sqrt(pi)*XMIN].
#
#   Approximate values for some important machines are:
#
#                          XMIN       XINF        XNEG     XSMALL
#
#  CDC 7600      (S.P.)  3.13E-294   1.26E+322   -27.220  7.11E-15
#  CRAY-1        (S.P.)  4.58E-2467  5.45E+2465  -75.345  7.11E-15
#  IEEE (IBM/XT,
#    SUN, etc.)  (S.P.)  1.18E-38    3.40E+38     -9.382  5.96E-8
#  IEEE (IBM/XT,
#    SUN, etc.)  (D.P.)  2.23D-308   1.79D+308   -26.628  1.11D-16
#  IBM 195       (D.P.)  5.40D-79    7.23E+75    -13.190  1.39D-17
#  UNIVAC 1108   (D.P.)  2.78D-309   8.98D+307   -26.615  1.73D-18
#  VAX D-Format  (D.P.)  2.94D-39    1.70D+38     -9.345  1.39D-17
#  VAX G-Format  (D.P.)  5.56D-309   8.98D+307   -26.615  1.11D-16
#
#
#                          XBIG       XHUGE       XMAX
#
#  CDC 7600      (S.P.)  25.922      8.39E+6     1.80X+293
#  CRAY-1        (S.P.)  75.326      8.39E+6     5.45E+2465
#  IEEE (IBM/XT,
#    SUN, etc.)  (S.P.)   9.194      2.90E+3     4.79E+37
#  IEEE (IBM/XT,
#    SUN, etc.)  (D.P.)  26.543      6.71D+7     2.53D+307
#  IBM 195       (D.P.)  13.306      1.90D+8     7.23E+75
#  UNIVAC 1108   (D.P.)  26.582      5.37D+8     8.98D+307
#  VAX D-Format  (D.P.)   9.269      1.90D+8     1.70D+38
#  VAX G-Format  (D.P.)  26.569      6.71D+7     8.98D+307
#
#*******************************************************************
#*******************************************************************
#
# Error returns
#
#  The program returns  ERFC = 0      for  ARG >= XBIG;
#
#                       ERFCX = XINF  for  ARG < XNEG;
#      and
#                       ERFCX = 0     for  ARG >= XMAX.
#
#
# Intrinsic functions required are:
#
#     ABS, AINT, EXP
#
#
#  Author: W. J. Cody
#          Mathematics and Computer Science Division
#          Argonne National Laboratory
#          Argonne, IL 60439
#
#  Latest modification: March 19, 1990
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#------------------------------------------------------------------
    my ($i);
    my ($x, $del, $xden, $xnum, $y, $ysq);
#------------------------------------------------------------------
#  Mathematical constants
#------------------------------------------------------------------
    my ($four, $one, $half, $two, $zero) = (4, 1, 0.5, 2, 0);
    my $sqrpi = 5.6418958354775628695e-1;
    my $thresh = 0.46875;
    my $sixten = 16;
#------------------------------------------------------------------
#  Machine-dependent constants
#------------------------------------------------------------------
    my ($xinf, $xneg, $xsmall) = (1.79e308, -26.628, 1.11e-16);
    my ($xbig, $xhuge, $xmax) = (26.543, 6.71e7, 2.53e307);
#------------------------------------------------------------------
#  Coefficients for approximation to  erf  in first interval
#------------------------------------------------------------------
    my @a = (3.16112374387056560e00, 1.13864154151050156e02,
             3.77485237685302021e02, 3.20937758913846947e03,
             1.85777706184603153e-1);
    my @b = (2.36012909523441209e01, 2.44024637934444173e02,
             1.28261652607737228e03, 2.84423683343917062e03);
#------------------------------------------------------------------
#  Coefficients for approximation to  erfc  in second interval
#------------------------------------------------------------------
    my @c = (5.64188496988670089e-1, 8.88314979438837594e0,
             6.61191906371416295e01, 2.98635138197400131e02,
             8.81952221241769090e02, 1.71204761263407058e03,
             2.05107837782607147e03, 1.23033935479799725e03,
             2.15311535474403846e-8);
    my @d = (1.57449261107098347e01, 1.17693950891312499e02,
             5.37181101862009858e02, 1.62138957456669019e03,
             3.29079923573345963e03, 4.36261909014324716e03,
             3.43936767414372164e03, 1.23033935480374942e03);
#------------------------------------------------------------------
#  Coefficients for approximation to  erfc  in third interval
#------------------------------------------------------------------
    my @p = (3.05326634961232344e-1, 3.60344899949804439e-1,
             1.25781726111229246e-1, 1.60837851487422766e-2,
             6.58749161529837803e-4, 1.63153871373020978e-2);
    my @q = (2.56852019228982242e00, 1.87295284992346047e00,
             5.27905102951428412e-1, 6.05183413124413191e-2,
             2.33520497626869185e-3);
#------------------------------------------------------------------
    $x = $arg;
    $y = abs($x);
    if ($y <= $thresh) {
#------------------------------------------------------------------
#  Evaluate  erf  for  |X| <= 0.46875
#------------------------------------------------------------------
        $ysq = $zero;
        if ($y > $xsmall) { $ysq = $y * $y }
        $xnum = $a[5]*$ysq;
        $xden = $ysq;
        for (my $i = 1 ; $i <= 3 ; ++$i) {
            $xnum = ($xnum + $a[$i]) * $ysq;
            $xden = ($xden + $b[$i]) * $ysq;
        }
        $$result = $x * ($xnum + $a[4]) / ($xden + $b[4]);
        if ($jint != 0) { $$result = $one - $$result }
        if ($jint == 2) { $$result = exp($ysq) * $$result }
        goto x800;
#------------------------------------------------------------------
#  Evaluate  erfc  for 0.46875 <= |X| <= 4.0
#------------------------------------------------------------------
    } elsif ($y <= $four) {
        $xnum = $c[9]*$y;
        $xden = $y;
        for (my $i = 1 ; $i <= 7 ; ++$i) {
            $xnum = ($xnum + $c[$i]) * $y;
            $xden = ($xden + $d[$i]) * $y;
        }
        $$result = ($xnum + $c[8]) / ($xden + $d[8]);
        if ($jint != 2) {
            $ysq = int($y*$sixten)/$sixten;
            $del = ($y-$ysq)*($y+$ysq);
            $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result;
        }
#------------------------------------------------------------------
#  Evaluate  erfc  for |X| > 4.0
#------------------------------------------------------------------
    } else {
        $$result = $zero;
        if ($y >= $xbig) {
            if (($jint != 2) || ($y >= $xmax)) { goto x300 }
            if ($y >= $xhuge) {
                $$result = $sqrpi / $y;
                goto x300;
            }
        }
        $ysq = $one / ($y * $y);
        $xnum = $p[6]*$ysq;
        $xden = $ysq;
        for (my $i = 1 ; $i <= 4 ; ++$i) {
            $xnum = ($xnum + $p[$i]) * $ysq;
            $xden = ($xden + $q[$i]) * $ysq;
        }
        $$result = $ysq *($xnum + $p[5]) / ($xden + $q[5]);
        $$result = ($sqrpi -  $$result) / $y;
        if ($jint != 2) {
            $ysq = int($y*$sixten)/$sixten;
            $del = ($y-$ysq)*($y+$ysq);
            $$result = exp(-$ysq*$ysq) * exp(-$del) * $$result;
        }
    }
#------------------------------------------------------------------
#  Fix up for negative argument, erf, etc.
#------------------------------------------------------------------
  x300:
    if ($jint == 0) {
        $$result = ($half - $$result) + $half;
        if ($x < $zero) { $$result = -$$result }
    } elsif ($jint == 1) {
        if ($x < $zero) { $$result = $two - $$result }
    } else {
        if ($x < $zero) {
            if ($x < $xneg) {
                $$result = $xinf;
            } else {
                $ysq = int($x*$sixten)/$sixten;
                $del = ($x-$ysq)*($x+$ysq);
                $y = exp($ysq*$ysq) * exp($del);
                $$result = ($y+$y) - $$result;
            }
        }
    }
  x800:
    return 1;
#---------- Last card of CALERF ----------
}

sub erf {
    my $x = @_ ? $_[0] : $_;
#--------------------------------------------------------------------
#
# This subprogram computes approximate values for erf(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, January 8, 1985
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#--------------------------------------------------------------------
    my $result;
    my $jint = 0;
    calerf($x, \$result, $jint);
    my $erf = $result;
    return $erf;
#---------- Last card of ERF ----------
}

########################################################################
## User functions.
########################################################################

sub erfc {
    my $x = @_ ? $_[0] : $_;
#--------------------------------------------------------------------
#
# This subprogram computes approximate values for erfc(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, January 8, 1985
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#--------------------------------------------------------------------
    my ($result);
    my $jint = 1;
    calerf($x, \$result, $jint);
    my $erfc = $result;
    return $erfc;
#---------- Last card of ERFC ----------
}

sub erfcx {
    my $x = @_ ? $_[0] : $_;
#------------------------------------------------------------------
#
# This subprogram computes approximate values for exp(x*x) * erfc(x).
#   (see comments heading CALERF).
#
#   Author/date: W. J. Cody, March 30, 1987
#
# Translation to Perl by Peter J. Acklam, December 3, 1999
#
#------------------------------------------------------------------
    my ($result);
    my $jint = 2;
    calerf($x, \$result, $jint);
    my $erfcx = $result;
    return $erfcx;
#---------- Last card of ERFCX ----------
}

sub erfinv {
    my $y = @_ ? $_[0] : $_;

    return  0             if $y == 0;
    return  erfcinv(1-$y) if $y >  0.5;
    return -erfcinv(1+$y) if $y < -0.5;

    #
    # Halley's rational 3rd order method:
    #   u <- f(x)/f'(x)
    #   v <- f''(x)/f'(x)
    #   x <- x - u/(1-u*v/2)
    #
    # Here:
    #   f(x) = erf(x) - y
    #   f'(x) = 2/sqrt(pi)*exp(-x*x)
    #   f''(x) = -4/sqrt(pi)*x*exp(-x*x)
    #
    my $x = 0;
    my $dx;
    my $c = .88622692545275801364908374167055;  # sqrt(pi)/2
    my $eps = 5e-15;
    do {
        my $f = erf($x) - $y;
        my $u = $c*$f*exp($x*$x);
        $dx = -$u/(1+$u*$x);
        $x += $dx;
    } until abs($dx/$x) <= $eps;
    return $x;
}

sub erfcinv {
    my $y = @_ ? $_[0] : $_;

    return 0 if $y == 1;

    my $flag = $y > 1;
    $y = 2 - $y if $flag;

    #
    # Halley's rational 3rd order method:
    #   u <- f(x)/f'(x)
    #   v <- f''(x)/f'(x)
    #   x <- x - u/(1-u*v/2)
    #
    # Here:
    #   f(x) = erfc(x) - y
    #   f'(x) = -2/sqrt(pi)*exp(-x*x)
    #   f''(x) = 4/sqrt(pi)*x*exp(-x*x)
    #
    my $x = 0;
    my $dx;
    my $c = -.88622692545275801364908374167055; # sqrt(pi)/2
    my $eps = 5e-15;
    do {
        my $u = $c*(erfcx($x) - $y*exp($x*$x));
        $dx = -$u/(1+$u*$x);
        $x += $dx;
    } until abs($dx/$x) <= $eps;

    return $flag ? -$x : $x;
}

sub erfcxinv {
    my $y = @_ ? $_[0] : $_;

    return 0 if $y == 1;

    #
    # Halley's 3rd order method:
    #   u <- f(x)/f'(x)
    #   v <- f''(x)/f'(x)
    #   x <- x - u/(1-u*v/2)
    #
    # Here:
    #   f(x) = erfcx(x) - y
    #   f'(x) = 2*(x*erfcx(x)-1/sqrt(pi));
    #   f''(x) = (2+4*x*x)*erfcx(x) - 4*x/sqrt(pi);
    #
    my $x = 0;
    my $dx;
    my $c = .56418958354775628694807945156079;  # 1/sqrt(pi)
    my $d = 2.2567583341910251477923178062432;  # 4/sqrt(pi)
    my $eps = 5e-15;
    do {
        my $f = erfcx($x) - $y;
        my $df = 2*($x*erfcx($x)-$c);
        my $ddf = (2+4*$x*$x)*erfcx($x) - $x*$d;
        my $u = $f/$df;
        my $v = $ddf/$df;
        $dx = -$u/(1-$u*$v/2);
        $x += $dx;
    } until abs($dx/$x) <= $eps;
    return $x;
}

package GoogleChart::SimpleEncoder;

sub new {
  my $class = shift;
  my $self = {};
  $self->{"lower_bound"} = shift;
  $self->{"upper_bound"} = shift;
  my %translations = 
     (0 => 'A', 1 => 'B', 2 => 'C', 3 => 'D', 4 => 'E', 5 => 'F', 6 => 'G', 
      7 => 'H', 8 => 'I', 9 => 'J', 10 => 'K',11 => 'L', 12 => 'M', 
      13 => 'N', 14 => 'O', 15 => 'P', 16 => 'Q', 17 => 'R', 18 => 'S', 
      19 => 'T', 20 => 'U', 21 => 'V', 22 => 'W', 23 => 'X', 24 => 'Y', 
      25 => 'Z', 26 => 'a', 27 => 'b', 28 => 'c', 29 => 'd', 30 => 'e', 
      31 => 'f', 32 => 'g', 33 => 'h', 34 => 'i', 35 => 'j', 36 => 'k', 
      37 => 'l', 38 => 'm', 39 => 'n', 40 => 'o', 41 => 'p', 42 => 'q', 
      43 => 'r', 44 => 's', 45 => 't', 46 => 'u', 47 => 'v', 48 => 'w', 
      49 => 'x', 50 => 'y', 51 => 'z', 52 => '0', 53 => '1', 54 => '2', 
      55 => '3', 56 => '4', 57 => '5', 58 => '6', 59 => '7', 60 => '8', 
      61 => '9');
  $self->{"translations"} = \%translations;
  $self->{"granule"} = (1.0 * $self->{"upper_bound"} - $self->{"lower_bound"}) / 61;
  bless $self,$class;
  return $self;
}

sub encode {
  my $self = shift;
  my $lower_bound = $self->{"lower_bound"};
  my $upper_bound = $self->{"upper_bound"};
  my $value = shift;
  return "_" unless defined $value;
  return "_" if $value > $upper_bound;
  return "_" if $value < $lower_bound;
  my $granule = $self->{"granule"};
  my $which_granule = int(0.5 + (($value - $lower_bound) / $granule));
  return $self->{"translations"}{$which_granule};
}


package SimplifyRange;

sub increment_around {
    my $x = shift;
    if ($x < 0.0) { $x = -$x; }
    if ($x <= 10.0) { return 1; }
    if ($x > 10.0) { return 10.0 * &increment_around($x * 0.1); }
    if ($x < 1.0) { return 0.1 * &increment_around($x * 10.0); }
}

sub down {
  my $number = shift;
  die unless defined $number;
  if ($number < 0) { return - &up(-$number); }
  if ($number >= 10) { return 10 * &down(int($number) / 10); }
  return 0.0 if $number == 0.0;
  if ($number <= 0.1) { return 0.1 * &down($number * 10);  }
  return int($number);
}
sub up {
  my $number = shift;
  die unless defined $number;
  if ($number < 0) { return - &down(-$number); }
  if ($number >= 10) { return 10 * &up(int($number) / 10); }
  if ($number <= 0.1) { return 0.1 * &up($number * 10); }
  if ($number == int($number)) { return $number; }
  return int($number)+1;
}

sub min_max_pair {
  my $lower = shift;
  my $upper = shift;
  $upper = &up($upper);
  $lower = &down($lower);
  if ($lower < $upper / 2 and $lower > 0) {
    $lower = 0.0;
  }
  return ($lower,$upper);
}

package DataSummary;

sub new {
  my $class = shift;
  my $self = {};
  $self->{"start_time"} = shift;
  $self->{"end_time"} = shift;
  my $buckets = shift;
  $self->{"buckets"} = $buckets;
  $self->{"min"} = [];
  $self->{"max"} = [];
  $self->{"count"} = [];
  $self->{"sum"} = [];
  my $i;
  for ($i=0;$i<$buckets;$i++) {
    $self->{"min"}->[$i] = undef;
    $self->{"max"}->[$i] = undef;
    $self->{"count"}->[$i] = 0;
    $self->{"sum"}->[$i] = 0;
  }
  bless $self,$class;
  return $self;
}

sub set_start_time {
  my $self = shift;
  $self->{"start_time"} = shift;
}

sub calculate_bucket_width {
  my $self = shift;
  return if exists $self->{"bucket_width"} and defined $self->{"bucket_width"};
  return unless exists $self->{"start_time"};
  return unless defined $self->{"start_time"};
  my $start_time = $self->{"start_time"};
  my $end_time = $self->{"end_time"};
  my $buckets = $self->{"buckets"};
  $self->{"bucket_width"} = (0.001 + $end_time - $start_time) / $buckets;
}

sub store {
  my $self = shift;
  my $timestamp = shift;
  $self->calculate_bucket_width();
  my $bucket_width = $self->{"bucket_width"};
  my $start_time = $self->{"start_time"};
  my $time_index = int (($timestamp - $start_time) / $bucket_width);
  my $buckets = $self->{"buckets"};
  die "Calculation error" if $time_index > $buckets;
  my $value = shift;
  my $min = $self->{"min"};
  my $max = $self->{"max"};
  my $count = $self->{"count"};
  my $sum = $self->{"sum"};
  $max->[$time_index] = $value unless defined $max->[$time_index];
  $min->[$time_index] = $value unless defined $min->[$time_index];
  if ($value > $max->[$time_index]) { $max->[$time_index] = $value;}
  if ($value < $min->[$time_index]) { $min->[$time_index] = $value;}
  $sum->[$time_index] += $value;
  $count->[$time_index] += 1;
}

sub maximum_anywhere {
  my $self = shift;
  my $max = $self->{"max"};
  my $v;
  my $best;
 BUCKET_MAXIMUM:
  foreach $v (@$max) { 
      next BUCKET_MAXIMUM unless defined $v;
      $best = $v if (!defined $best) || ($best < $v);
  }
  return $best;
}
sub minimum_anywhere {
  my $self = shift;
  my $min = $self->{"min"};
  my $v;
  my $best;
 BUCKET_MINIMUM:
  foreach $v (@$min) {
      next BUCKET_MINIMUM unless defined $v;
      $best = $v if !defined $best or $best > $v; 
  }
  return $best;
}

sub google_chart {
    my $self = shift;
    my $chart_encoder = shift;
    my $min = $self->{"min"};
    my $max = $self->{"max"};
    my $count = $self->{"count"};
    my $sum = $self->{"sum"};
    my $min_string = "";
    my $max_string = "";
    my $val_string = "";
    my $i;
    my $buckets = $self->{"buckets"};
  BUCKET:
    for($i=0;$i<$buckets;$i++) {
	if ($count->[$i] == 0) { 
	    $min_string .= "_"; $max_string .= "_"; $val_string .= "_";
	    next BUCKET;
	}
	$min_string .= $chart_encoder->encode($min->[$i]);
	$max_string .= $chart_encoder->encode($max->[$i]);
	$val_string .= $chart_encoder->encode(1.0 * $sum->[$i] / $count->[$i]);
    }
    return ($min_string,$val_string,$max_string);
}

package MiniDateManip;
# It would be much easier if we could rely on DateManip, but since it's
# not part of Perl core, and this is often running on bizarre things like
# the Perl that comes with HP OpenView, we do what we have to here.

sub subtract_months {
    my $now = shift;
    my $month_count = shift;
    my ($sec,$min,$hour,$mday,$month,$year,$wday,$yday,$dst) = localtime($now);
    $month -= $month_count;
    while ($month < 0) { $year--; $month += 12; }
    return Time::Local::timelocal($sec,$min,$hour,$mday,$month,$year,$wday,$yday,$dst);
}


package Main;