#!/usr/bin/perl -W

# UNTESTED!!!

# calculate linear regression for a set of (x,Y) data
# http://www.statsoftinc.com/textbook/stmulreg.html

# The idea is to find the (a,b) coefficients for
# formula y = f(x) = a + b * x 
# which provides the smallest possible sum( (f(x) - y_measured_at(x))**2 )
# In other words: the straight line that best fits the measured values.

# possible uses: verify calculated edgetime_corr data
# against measured values (x=edgetime, y=measured jumps in ADCcount units)

# sequences for appr. b==2 :
# X-file:
# a=01 04 06 07 08 0a 10 12 13 15 17
# Y-file:
# a=02 08 0c 0f 0f 14 20 24 26 29 2e

# diff_table.pl is cloned for the administrative stuff
# we should really go object oriented and apply code-reuse...

if (@ARGV != 2) {
	print "reading 'rpm   kpa   ignition_advance' log from STDIN\n";
	read_advfile();
    print "other usage:\n";
    print "\t linreg.pl X-table Y-table\n\n";
    exit;
}

$tblX = read_table($ARGV[0]);
$tblY = read_table($ARGV[1]);

foreach $n (sort keys %$tblX) {
    if (defined($tblY->{$n})) {
	print "\n---- set of data: $n \n";
	print linreg_line($tblX->{$n}, $tblY->{$n}) . "\n";
    } else {
	print "'$n' not found in file $ARGV[1]\n";
    }
}
#------------------
# call before any calculations
sub stats_reset
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	$st->{sumX} = 0;
	$st->{sumXY} = 0;
	$st->{sumXX} = 0;
	$st->{sumY} = 0;
	$st->{N} = 0;
	$st->{sumErrsq} = 0;

	# calculated data:
	$st->{denom} = undef;
	$st->{Anom} = undef;
	$st->{Bnom} = undef;
	$st->{A} = undef;
	$st->{B} = undef;
}
#-------------
# call for each (x,y) pairs
sub stats_xy
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	my $X = shift;
	my $Y = shift;

	$st->{sumX} += $X;
	$st->{sumXY} += $X * $Y;
	$st->{sumXX} += $X * $X;
	$st->{sumY} += $Y;
	$st->{N} ++;
}

#================
sub hate0
{
	my $de = shift;
	my $epsilon = 0.00001;
	return 1 if (($de < $epsilon) && ($de > 0-$epsilon));	# we hate divby0
	return $de;
}
#-----------------
# calculates the common denominator of "a" and "b" (in y = a + b*x)
# aware that the formulas are from my mother, I did not verify :-)
sub linreg_denom
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	return $st->{denom} if(defined $st->{denom});
	my $de = $st->{sumX} * $st->{sumX};
	$de -= $st->{N} * $st->{sumXX};
	my $zz = hate0($de);
	$st->{denom} = $zz;
	return $zz;
}
#-----------------
# calculates the nominator of "a" (in y = a + b*x)
sub linreg_Anom
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	$st->{Anom}= $st->{sumX} * $st->{sumXY};
	$st->{Anom} -= $st->{sumXX} * $st->{sumY};
	$st->{A} = $st->{Anom} / linreg_denom($st);
	return $st->{Anom};
}

#-----------------
# calculates the nominator of "b" (in y = a + b*x)
sub linreg_Bnom
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	$st->{Bnom}= $st->{sumX} * $st->{sumY} - $st->{N} * $st->{sumXY};
	$st->{B} = $st->{Bnom} / linreg_denom($st);
	return $st->{Bnom};
}

#-----------------
# calculate the sum(error ** 2)
sub linreg_err
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	my $X = shift;
	my $Y = shift;

	my $err = linreg_val($st, $X) - $Y;
	$st->{sumErrsq} += $err * $err;
}

# read a Y value (from the fitted line) belonging to an X
sub linreg_val
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	my $X = shift;
	print "where is st->{A}\n" if(! defined $st->{A});
	print "where is st->{B}\n" if(! defined $st->{B});
	print "where is X\n" if(! defined $X);
	return $st->{A} + $st->{B} * $X;
}

sub linreg_statistics
{
	my $st = shift;	# we get a hash reference (use it as a struct to store statistical data)
	my $str = "samplecount=$st->{N}\ny=$st->{A} + $st->{B} * x\nsumerrsquare=$st->{sumErrsq}\n";
	return $str;
}
#----------------
# fits a straight line on data input as Xarray and Yarray
sub linreg_line
{
    my $line1 = shift @_;
    my $line2 = shift @_;

    @values1 = split / /, $line1;
    @values2 = split / /, $line2;
    
	$stat->{yep} = 1;	# hash reference 
	stats_reset($stat);
	#print $stat->{sumXX};

    my $diff = "";
    for ($i=0; $i<=$#values1; $i++) {
#	$diff .= hex($values1[$i])-hex($values2[$i]) . " ";
#	printf("%3d ", hex($values2[$i])-hex($values1[$i]));
	stats_xy($stat, hex($values1[$i]), hex($values2[$i]));
#	printf("%02X ", hex($values1[$i])-hex($values2[$i]));
    }

	# calculate statistics
	# linreg_denom($stat); # will be called automagically
	linreg_Anom($stat);
	linreg_Bnom($stat);

	# calculate error
    	for ($i=0; $i<=$#values1; $i++) {
		linreg_err($stat, hex($values1[$i]), hex($values2[$i]));
	}

    	return linreg_statistics($stat);
}


#-----------------
# reads the lines of a file and stores key = values
# in the hash given as parameter
sub read_table
{
    my $fname = shift @_;
    my $tbl;	# hash reference

    open(FILE, "<$fname") or die "can't open $fname";

    while($l = <FILE>) {
	chomp $l;

	if ($l =~ /^\s*($|\/\/)/) {  # skip lines consisting of comments only
	    next;
	}
    
	if ($l =~ /(\w+)\s*=\s*(.+)/) {  # search for key=values
	    $array_name = $1;
	    $values = $2;
	    $tbl->{$array_name} = $values;
	}
    }
    
    close FILE;
    return $tbl;
}


# =================================
# reading logfile and fitting ignition advance
sub read_advfile
{
	my $r;
	my $k;
	# some wired variables from mik's tables:
	@kpabins=(0x14,0x1E,0x28,0x32,0x3C,0x4B,0x5A,0x64);
	@rpmorigbins=(0x05,0x0A,0x0F,0x14,0x1C,0x24,0x2C,0x34);

	# orig resolution above was 100 RPM
	@rpmbins = midvals(\@rpmorigbins, 100);
	@kpamids = midvals(\@kpabins, 1);

	# now read the input
	my $li;
	while($li = <STDIN>) {
		my ($rpm, $kpa, $adv);
		# chomp $l;
		if(($rpm, $kpa, $adv) = ($li =~ /(\d+)\t(\d+)\t(\d+)/)){
			# next if(!defined $adv);
			rpm_kpa_adv($rpm, $kpa, $adv / 2);
		}
	}	# while

	print_adv_stats();

}	# read_advfile()

# ==================
sub midvals
{
	my $aref = shift;
	my $mult = shift;
	my @a;
	my $i;
	my $prev;
	foreach $i (@$aref){	# we want the threshold halfway between the columns 
		if(defined $prev){
			push @a, (0.5 * $mult * ($prev + $i));
		}
		$prev = $i;
	}
	return @a;
}
# =================================
sub rpm_kpa_adv
{
	my ($rpm, $kpa, $adv) = @_;

    return if ($adv < 1 || $adv > 90);
	return if ($rpm < 200);

	$r = search_table($rpm, \@rpmbins);
	$k = search_table($kpa, \@kpabins);

		# the values that don't fit to a the line well
		# (actually rpmcolumn, currently!)
		# must be eliminated (alternatively downweighted?)
	my $rorig = 100*$rpmorigbins[$r];
	my $rpm_err = $rorig - $rpm;
	simple_ave($rpm, $kpa, $adv);

	return if(($rpm_err * $rpm_err) > (200*200));	# allow +-200 RPM

	my $key = $r."_".$k;
	setdebug($key);
	
	print "data: $rpm($rorig), $kpa, $adv\n" if $debug;
	
		# we'll have a lotsof regression-arrays, $ad{$key} is one
	if(!defined $ad{$key}){
		$ad{$key}->{yep} = 1;	# hash reference 
		stats_reset($ad{$key});
		print "wow, new loadsite: $key\n" if $debug;
	}
	# yes, use the values for statistics:
	print "data: $key, $kpa, $adv\n" if $debug;
	stats_xy($ad{$key}, $kpa, $adv);	# how is adv in function of kpa?
}	# rpm_kpa_adv()

# only debug a few loadsites
sub setdebug
{
	my $key = shift;
	if($key eq "1_2xxx"){
		$debug = 1;
	} else {
		$debug = undef;
	}
}
#=================
# simple averaging
sub simple_ave
{
	my ($rpm, $kpa, $adv) = @_;
	my $r;
	my $k;
	# we use the midpoints for searching !!
	$r = search_table($rpm, \@rpmbins);
	$k = search_table($kpa, \@kpamids);
	
	# drop data if too far away from grid
	my $rpm_err = 100*$rpmorigbins[$r] - $rpm;
	return if(($rpm_err * $rpm_err) > (150*150));	# allow +-150 RPM

	my $kpa_err = $kpabins[$k] - $kpa;
	return if(($kpa_err * $kpa_err) > (7*7));	# allow +-7 kPa

	my $key = $r."_".$k;
	$ave{$key}->{N}++;
	$ave{$key}->{sum} += $adv;
}
# ================
# results for above
sub calc_ave
{
	my $st = shift;	# statistics for loadsite
	my $adv = shift; # default return value, if no data found
	return (0,$adv) if(! defined $st);
	if(defined $st->{N}){
		return ($st->{N}, myround($st->{sum} / hate0($st->{N})));
	} else {
		return (0, $adv);
	}
}
# =================================
sub print_adv_stats
{
	my $k;
	my $r;
	my $limp=10;
	# print the results:
	for($k=0; $k<=8; $k++){	# yes, an extra line 8 was also calculated! 
		my $ad1 = $limp;
		my $ada = $limp; 
		my $ad2 = $limp;	# limphome for no datapoints
		my $k2 = ($k==0 ? $k : ($k-1)); # previous line (except at 0)
		my $kpa = $kpabins[$k2]; # what was the kpa?
		my $l = "n[$k2]($kpa):";
		my $a = $l;
		my $l2 = $l;

		for($r=0; $r<=7; $r++){
			my $key = $r."_".$k;
			my $key2 = $r."_".$k2;
			my $N; my $Na; my $N2;	# matching loadsites k and k2
			
			# the same kpa value calced from
			# the loadsite above, and loadsite below:
			($N2, $ad2) = calc_loadsite_adv($ad{$key2}, $kpa, $ad2);
			($Na, $ada) = calc_ave($ave{$key2}, $ada);
			($N, $ad1) = calc_loadsite_adv($ad{$key}, $kpa, $ad1);
			
			$l2 .= "\t$ad2($N2)";
			$a  .= "\t$ada($Na)";
			$l	.= "\t$ad1($N)";
			
		}
		print "$l2\n$a\n$l\n\n";
	}
}	# print_adv_stats()

# ==================
# returns the number of data found, and the calculated 'y' value
sub calc_loadsite_adv
{
	my $st = shift;	# statistics for loadsite
	my $x = shift;	# the x input
	my $adv = shift; # default return value, if no data found
	return (0,$adv) if(! defined $st);
	
	if(!defined $st->{denom}){
		# calculate statistics once
		# linreg_denom($ad{$key}); # automagically called
		linreg_Anom($st);
		linreg_Bnom($st);
	}
	return ($st->{N}, myround(linreg_val($st, $x)));
}
#===========
sub myround 
{
	my $ret = shift;
	$ret *= 10;
	return 0.1 * int($ret+0.5);
}
# =================================
# returns the index (of the array) which element is
# greater or equal than the given threshold
# for 8 element array, returns 0..8 ?
sub search_table
{
	my $val = shift;
	my $aref = shift;	# array reference that we traverse
	my $i = 0;
	foreach $e (@$aref){
		return $i if($e >= $val);
		$i++;
	}
	return $i;
}