#!/usr/bin/perl -w
use strict;
use warnings;
use File::Basename;
use Tie::Array::Sorted;

my $prog = basename($0);
die "usage: $prog alpha_up alpha_down alpha_var quantile fix_buf min_buf max_buf\n" if @ARGV != 7;
my ($alpha_down, $alpha_up, $alpha_var, $quantile, $fix_buf, $min_buf, $max_buf) = @ARGV;

my $decimals = 2;

my $upper_sds = ltqnorm($quantile);
warn "quantile: $quantile\nupper_sds: $upper_sds\n";

my $weight = 0;
my $weight_var = 0;
my $weight_fb = 0;
my $total = 0;
my $total_fb = 0;
my $total_variance = 0;
my $total_sd = 0;
my $n = 0;
my $prev_upper;
my $prev_upper_df;
my $prev_upper_fb;
my $prev_average = 0;
my $prev_average_variance = 0;

my $tot_ok = 0;
my $tot_ok_df = 0;
my $tot_ok_fb = 0;
my $tot_extra = 0;
my $tot_extra_df = 0;
my $tot_extra_fb = 0;
my $tot_v = 0;

my $min_n = 2;

my $compromise = 1;
my $prev_v = 0;
my $var_both_ways = 1;

my @difference_factors;
tie @difference_factors, "Tie::Array::Sorted", sub { $_[0] <=> $_[1] };

# start with high variance and high average for safety...?
#$total = 10;
#$total_variance = 10;
#$weight = 1;

my $alpha1_down = 1 - $alpha_down;
my $alpha1_up = 1 - $alpha_up;
my $alpha1_fb = $alpha1_down;
my $alpha1_var = 1 - $alpha_var;

my @head = qw(n v av diff sd upper buff upbuff df_i dfbuff upfb fixbuf ok ok_df ok_fb extra ext_df ext_fb);
print join("\t", @head), "\n";

# TODO make "df" somehow follow using exponential smoothing also?

while (defined (my $v = <STDIN>)) {
    chomp $v;
    ++$n;
    my $alpha1 = $v < $prev_average ? $alpha1_down : $alpha1_up;
    $weight = 1 + $alpha1 * $weight;
    $weight_fb = 1 + $alpha1_fb * $weight_fb;
    $total = $v + $alpha1 * $total;
    $total_fb = $v + $alpha1_fb * $total_fb;
    my $average = $total / $weight;
    my $average_fb = $total_fb / $weight_fb;
    my ($difference, $difference_factor);
    if ($prev_average) {
        $difference = $v - $prev_average;
        $difference_factor = $difference / $prev_average;
    } else {
        $difference = 0;
        $difference_factor = 0;
    }
    push @difference_factors, $difference_factor;
    my $diff2 = $difference**2;

    # TODO separate alpha1, weight for variance?
#    my $alpha1_var = 0.99; #$diff2 < $prev_average_variance ? $alpha1_down : $alpha1_up;
#    my $alpha1_var = $alpha1;
    if ($var_both_ways || $difference > 0) {
        $total_variance = $diff2 + $alpha1_var * $total_variance;
        $weight_var = 1 + $alpha1_var * $weight_var;
    }
    my $average_variance = $weight_var ? $total_variance / $weight_var : 0;
    my $sd = sqrt($average_variance);
    my $upper = $average + $sd * $upper_sds;
    my $buffer = $upper / $average - 1;
    $buffer = $min_buf if $buffer < $min_buf;
    $buffer = $max_buf if $buffer > $max_buf;
    $upper = $average * (1 + $buffer);
    my $okay = ($prev_upper||0) >= $v ? 1 : 0;
    my $df_i = int($n * $quantile + 1 - 1e-9); 
    if ($df_i >= $n) {
        $df_i = $n - 1;
    }
    my $df = $difference_factors[$df_i];
    my $df_buffer = $df;
    if ($compromise) {
        if ($buffer > $df_buffer) {
            $df_buffer = $buffer;  #($df_buffer + $buffer) / 2;
        }
    }
    $df_buffer = $min_buf if $df_buffer < $min_buf;
    $df_buffer = $max_buf if $df_buffer > $max_buf;
    my $upper_df = $average * (1+$df_buffer);
    my $okay_df = ($prev_upper_df||0) >= $v ? 1 : 0;
    my $extra_df = ($prev_upper_df||0) - $v;
    my $extra = ($prev_upper||0) - $v;

    my $upper_fb = $average_fb * (1+$fix_buf);
    my $okay_fb = ($prev_upper_fb||0) >= $v ? 1 : 0;
    my $extra_fb = ($prev_upper_fb||0) - $v;

    my @row = ($n, $v, $average, $difference, $sd, $upper, $buffer, $upper_df, $df_i, $df_buffer, $upper_fb, $fix_buf, $okay, $okay_df, $okay_fb, $extra, $extra_df, $extra_fb);
    print join("\t", map {fmt_num($_)} @row), "\n";
    $prev_upper = $upper;
    $prev_upper_df = $upper_df;
    $prev_upper_fb = $upper_fb;
    $prev_average = $average;
    $prev_average_variance = $average_variance;
    if ($n > $min_n) {
        $tot_v += $v;
        $tot_ok += $okay;
        $tot_ok_df += $okay_df;
        $tot_ok_fb += $okay_fb;
        $tot_extra += $extra if $extra > 0;
        $tot_extra_df += $extra_df if $extra_df > 0;
        $tot_extra_fb += $extra_fb if $extra_fb > 0;
    }
}

$n -= $min_n;

warn "n: ", $n, "\n";
warn "tot_ok: ", $tot_ok / $n, "\n";
warn "tot_ok_df: ", $tot_ok_df / $n, "\n";
warn "tot_ok_fb: ", $tot_ok_fb / $n, "\n";
warn "tot_extra frac: ", $tot_extra / $tot_v, "\n";
warn "tot_extra_df frac: ", $tot_extra_df / $tot_v, "\n";
warn "tot_extra_fb frac: ", $tot_extra_fb / $tot_v, "\n";

sub fmt_num {
    local ($_) = @_;
    $_ = sprintf("%.*f", $decimals, $_);
    if (/\./) {
        s/0+$//; s/\.$//;
    }
    return $_;
}

sub ltqnorm {
    #
    # Lower tail quantile for standard normal distribution function.
    #
    # This function returns an approximation of the inverse cumulative
    # standard normal distribution function.  I.e., given P, it returns
    # an approximation to the X satisfying P = Pr{Z <= X} where Z is a
    # random variable from the standard normal distribution.
    #
    # The algorithm uses a minimax approximation by rational functions
    # and the result has a relative error whose absolute value is less
    # than 1.15e-9.
    #
    # Author:      Peter John Acklam
    # Time-stamp:  2000-07-19 18:26:14
    # E-mail:      pjacklam@online.no
    # WWW URL:     http://home.online.no/~pjacklam

    my ($p) = @_;
    die "input argument must be in (0,1)\n" unless 0 < $p && $p < 1;

    # Coefficients in rational approximations.
    my @a = (-3.969683028665376e+01,  2.209460984245205e+02,
             -2.759285104469687e+02,  1.383577518672690e+02,
             -3.066479806614716e+01,  2.506628277459239e+00);
    my @b = (-5.447609879822406e+01,  1.615858368580409e+02,
             -1.556989798598866e+02,  6.680131188771972e+01,
             -1.328068155288572e+01 );
    my @c = (-7.784894002430293e-03, -3.223964580411365e-01,
             -2.400758277161838e+00, -2.549732539343734e+00,
              4.374664141464968e+00,  2.938163982698783e+00);
    my @d = ( 7.784695709041462e-03,  3.224671290700398e-01,
              2.445134137142996e+00,  3.754408661907416e+00);

    # Define break-points.
    my $plow  = 0.02425;
    my $phigh = 1 - $plow;

    # Rational approximation for lower region:
    if ( $p < $plow ) {
       my $q  = sqrt(-2*log($p));
       return ((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) /
               (((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1);
    }

    # Rational approximation for upper region:
    if ( $phigh < $p ) {
       my $q  = sqrt(-2*log(1-$p));
       return -((((($c[0]*$q+$c[1])*$q+$c[2])*$q+$c[3])*$q+$c[4])*$q+$c[5]) /
                (((($d[0]*$q+$d[1])*$q+$d[2])*$q+$d[3])*$q+1);
    }

    # Rational approximation for central region:
    my $q = $p - 0.5;
    my $r = $q*$q;
    return ((((($a[0]*$r+$a[1])*$r+$a[2])*$r+$a[3])*$r+$a[4])*$r+$a[5])*$q /
           ((((($b[0]*$r+$b[1])*$r+$b[2])*$r+$b[3])*$r+$b[4])*$r+1);
}