package Stat;
# assume a, b, and h have same length
# h = [ Häufigkeit, Gewicht ]
# N = Anzahl der Beobachtungen
# a = [ Einzelwerte der Beobachtungen ]

sub nsort {
    { $a <=> $b }
}
sub prod($$) {
    my $a = shift;
    my $b = shift;

    if (!defined($b)) {
	return $a;
    }

    my @arr = @$a;
    for (my $ix = 0; $ix < @arr; $ix++) {
	$arr[$ix] *= $$b[$ix];
    }

    [ @arr ];
}

sub sum($) {
    my $r = shift;
    my $sum = 0;

    for (@$r) {
	$sum += $_;
    }
    $sum;
}

sub am($$) { # Arithmetisches Mittel
    my $a = shift;
    my $h = shift;

    my $am = 0;
    if (defined($h)) {
	$am = sum(prod($a, $h)) / sum($h);
    } else {
	$am = sum($a) / @$a;
    }

    $am;
}

sub gm($$) { # Geometrisches Mittel
    # assuming all @$a are > 0
    my $a = shift;
    my $h = shift;

    my $gm = 1;
    my $N;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $gm *= $arr[$ix]**$$h[$ix];
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $gm *= $arr[$ix];
	}
    }

    $gm = $gm**(1/$hs);
}

sub qm($$) { # Quadratisches Mittel
    my $a = shift;
    my $h = shift;

    my $qm = 0;
    my $N;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $qm += $arr[$ix]*$arr[$ix]*$$h[$ix];
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $qm += $arr[$ix]*$arr[$ix];
	}
    }
    $qm /= $N;
}

sub hm($$) { # Harmonisches Mittel
    my $a = shift;
    my $h = shift;

    my $hm = 0;
    my $N;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $hm += $$h[$ix]/$arr[$ix];
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    $hm += 1/$arr[$ix];
	}
    }
    $hm = $N/$hm;
}

sub zw($$) { # Zentralwert, Median
    my $a = shift;
    my $h = shift;
    
    my $zw = 0;
    my $N;
    if (defined($h)) {
	my @arr;
	$N = sum($h);
	for (my $ix = 0; $ix < @$a; $ix++) {
	    push @arr, [ $$a[$ix], $$h[$ix] ];
	}
	my @b = sort { $$a[0] <=> $$b[0] } @arr;
	my $ho = 0;
	for ($ix = 0; $ix < @b; $ix++) {
	    my $hi = $b[$ix][1];
	    my $hn = $ho + $hi*$b[$i][0];
	    if ($hn < $N/2) {
		next;
	    } elsif ($hn == $N/2) {
		# TODO
	    } else {
		# TODO
		$zw = $b[$ix][0];
		last;
	    }
	}
	#....
    } else { # it's like all h is 1's
	my @arr = sort nsort @$a;
	$N = @arr;
	my $n = int($N / 2);
	if ( ($N % 2) != 0 ) { # $N is odd
	    $zw = $arr[$n];
	} else { # $N is even
	    $zw = ( $arr[$n] + $arr[$n] ) / 2;
	}
    }
    $zw;
}

sub ls($$) { # Durchschnittliche Abweichung, Standardabw., lineare Streuung
    my $a = shift;
    my $h = shift;

    my $am = am($a,$h); # can also be median
    my $ls;
    my $N;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $ls += $$h[$ix]*$x*$x;
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $ls += $x*$x;
	}
    }
    $ls = sqrt($ls/($N-1));
}

sub qs($$) { # Mittlere quadratische Abweichung, Standardavw., quadratische Streuung
    my $a = shift;
    my $h = shift;

    my $am = am($a,$h);
    my $qs;
    my $N;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $qs += $$h[$ix]*$x*$x;
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $qs += $x*$x;
	}
    }
    $qs = sqrt($qs/($N-1));

    [ $N, $am, $qs ];
}

sub stat($$) {
    my $a = shift;
    my $h = shift;

    my $am = am($a,$h);
    my $qs;
    my $N;
    my $min;
    my $max;
    if (defined($h)) {
	my @arr = @$a;
	$N = sum($h);
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $qs += $$h[$ix]*$x*$x;
	}
    } else { # it's like all h is 1's
	my @arr = @$a;
	$N = @arr;
	for (my $ix = 0; $ix < @arr; $ix++) {
	    my $x = $arr[$ix] - $am;
	    $qs += $x*$x;
	}
    }
    if ($N > 1) {
	$qs = sqrt($qs/($N-1));
    } else {
	$qs = -1;
    }
    my @arr = sort nsort @$a;
    $min = $arr[0];
    $max = $arr[$#arr];
    [ $N, $min, $am, $max, $qs ];
}

sub stat_linreg($$$) {
    my $xr = shift;
    my $yr = shift;
    my $h = shift;

    my @x = @$xr;
    my @y = @$yr;
    if (@x != @y || @x == 0) {
	printf "Err (stat_linreg): arrays of lengths %d %d\n", @x+0, @y+0;
	return [0, 0];
    }
    my $a0 = 0;
    my $a1 = 0;

    my $N  = @x+0;
    my $x  = 0;
    my $y  = 0;
    my $xx = 0;
    my $xy = 0;
    if (defined($h)) {
	$N = sum($h);
	for (my $ix = 0; $ix < @x; $ix++) {
	    $x  += $$h[$ix] * $x[$ix];
	    $y  += $$h[$ix] * $y[$ix];
	    $xx += $$h[$ix] * $x[$ix] * $x[$ix];
	    $xy += $$h[$ix] * $x[$ix] * $y[$ix];
	}
    } else {
	for (my $ix = 0; $ix < @x; $ix++) {
	    $x  += $x[$ix];
	    $y  += $y[$ix];
	    $xx += $x[$ix] * $x[$ix];
	    $xy += $x[$ix] * $y[$ix];
	}
    }
    #printf "%.1f %.1f %.1f %.1f\n", $x, $y, $xx, $xy;

    #
    my $b = $N * $xx - $x * $x;
    if ($b == 0) {
	printf "Err (stat_linreg): divide by zero\n";
	return [0, 0];
    }
    $a0 = ($y  * $xx - $xy * $x) / $b;
    $a1 = ($xy * $N  - $y  * $x) / $b;

    ($a0, $a1);
}


1;
