#!/usr/bin/perl
#

use List::Util;

sub print_help {
    print STDERR <<EOF;
  usage:  average_variance_of_slope filename [xstart] [xend] [nsteps] [-q]
      This script returns the average and variance of the slope of f based
      on a two-column x-sorted datafile containing  x f(x)
      filename can be - for stdin (for usage in pipes)
      If xstart / xend are not one of the x values: linear interpolation
      The nstep for computing the variance defaults to 30 steps
                                      (or less if few data are available)
      -q : quiet output, no statistical information
  output: 2 numbers: the average and the variance
      created by Nick Manini & Andrea Stenco on 13-04-2016
                                    this version by Nick Manini, 11-09-2017
EOF
return;
}

$tmpfile="/tmp/ave_vas".$$;

$tin="unset";
$tfin="unset";
$step="unset";
$noise=1;

while($#ARGV+1){
    $p=shift(@ARGV);
    if($p eq "-h"){# print help
        print_help;
        exit;
    }
    if($p eq "-q"){ # quiet
        $noise=0;
        next;
    }
    if($#filename < 0 ){
	if ($p eq "-"){
	    push(@filename,$p);
	    next;
	}else{
	    if(-e $p){
		push(@filename,$p);
		next;
	    }else{
		die "non existing file: $p\n";
	    }
	}
    }
    if($tin eq "unset"){
	$tin=$p;
        next;
    }
    if($tfin eq "unset"){
	$tfin=$p;
        next;
    }
    if($step eq "unset"){
	$step=$p;
	$step>1 or die "better have at lest 2 steps!!\n";
        next;
    }
    die "error: \"$p\" is an unknown option\n";
}

if($filename[0] eq "-" || $#filename < 0 ){
    $file=$tmpfile;
    `cat > $file`;
}else{
    $file=$filename[0];
}

$_=`head -1 $file`;
s/^\s+//;
my @temp =  split(/\s+/);
$tinfile=shift(@temp)+0.;
$tbef=$tinfile;
$fbef=shift(@temp)+0.;
$_=`tail -1 $file`;
s/^\s+//;
@temp =  split(/\s+/);
$tfinfile=shift(@temp)+0.;
$lines=`awk '/^#/{c++}END{print NR-c}' $file`+0;
$lines>2 || die "average_variance_of_slope: No way to work with ",$lines," lines: terminating.\n";

# check input and set defaults:
if ($tin eq "unset" || $tin<$tinfile){
    $tin=$tinfile;
}
if ($tfin eq "unset" || $tfin>$tfinfile){
    $tfin=$tfinfile;
}
if($step eq "unset" || $step<2){
    $step=List::Util::min 30,$lines-1;
}

if($noise>0){
    print STDERR "average_variance_of_slope: file=",$file," tin=",$tin," tfin=",$tfin," nstep=",$step," ndata=",$lines,"\n";
}

$deltat=1.*($tfin-$tin)/$step;

#print "QUAA ",$deltat," ",$tfin," ",$tin," ",$step,"\n";

open(INPUT, $file) || die "Cannot open $file\n";

$s1=0;
$s2=0;
$count=0;
for($i=0;$i<$step;$i++){ 
    $ti=$tin+$i*$deltat;	
    $tf=$tin+($i+1)*$deltat;	
    LINE: while(<INPUT>){
	s/^\s*//;       # remove leading spaces
        if(!/^#/){	# ignore commented lines
	    @li = split('\s+');
	    $t=$li[0]+0.;
	    $f=$li[1]+0.;
	    if($t<=$ti){
		$tbef=$t;
		$fbef=$f;
	    }
	    last LINE if $t>$tf;
	    $taft=$t;
	    $faft=$f;
#	    print "QUUU ",$t," ",$tbef," ",$taft,"\n";
	}
    }
    if($taft>$tbef){
	$val=($faft-$fbef)/($taft-$tbef);
	#	print "QUII ",$deltat," ",$tbef," ",$taft," ",$fbef," ",$faft," ",$val,"\n";
	$s1+=$val;
	$s2+=$val*$val;
	$tbef=$taft;
	$fbef=$faft;
	$taft=$t;
	$faft=$f;
	$count++;
    }
#    print "QUXX ",$count," ",$i," ",$ti," ",$tf," ",$t,"\n";
}

$mean=$s1/List::Util::max 1,$count;
$diff=List::Util::max 0,$count*$s2 - $s1*$s1;
$std=sqrt($diff)/$count;
$sem=$std/sqrt($count); # standard error on the mean
# sem is an estimator of the standard deviation on the final mean,
# not on a measurament of an individual signal segment!!!

print $mean,"\t",$sem,"\n";

if (-e $tmpfile){
    `rm $tmpfile`;
}