import sys
import numpy as np
from numpy import sin, cos, tan, sinh, cosh, tanh, exp, log, sqrt, pi
from pprint import pprint

from tklib.tkutils import replaceif


def inverse_quadratic(x, y, i):
    xa = x[i-1]
    fa = y[i-1]
    xb = x[i]
    fb = y[i]
    xc = x[i+1]
    fc = y[i+1]
    return xa * fb * fc / (fa - fb) / (fa - fc) \
         + xb * fa * fc / (fb - fa) / (fb - fc) \
         + xc * fa * fb / (fc - fa) / (fc - fb)   

def lagrange(x, y, x0, *, IsPrint = 0):
    ndata = len(x)
    if ndata < 3:
        if IsPrint == 1:
            print("Error: ndata [{}] must be larger than 2".format(ndata))
        exit()
    w = 0.0;
    for i in range(0, ndata):
        w1 = 1.0
        for j in range(0, ndata):
            if i != j:
                w1 *= (x0 - x[j]) / (x[i] - x[j])
        w += w1 * y[i];
    return w

def linear(x, y, x0):
    if x0 < x[0]:
        return y[0]
    n = len(x)
    if x[n-2] <= x0:
        return y[n-1]
    xstep = x[1] - x[0]
    i = int((x0 - x[0]) / xstep)
    yi = y[i] + (y[i+1] - y[i]) / (x[i+1] - x[i]) * (x0 - x[i])
    return yi

def fft(y):
    return np.fft.fft(y)


"""
sub Interpolate
{
    my ($pX, $pY, $x) = @_;
    return $pY->[0] if($x <= $pX->[0]);
    my $nData = @$pX;
    return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

    my $idx  = FindIndexByBinarySearch($pX, $x);
    my $idx1 = $idx+1;
    my $k = ($pY->[$idx1] - $pY->[$idx]) / ($pX->[$idx1] - $pX->[$idx]);
    return $pY->[$idx] + $k * ($x - $pX->[$idx]);
}


sub FindIndexByBinarySearch
{
    my ($pX, $x) = @_;
    return 0 if($x <= $pX->[0]);
    my $nData = @$pX;
    return $nData - 1 if($pX->[$nData-1] <= $x);

    my $idx0 = 0;
    my $idx1 = $nData - 1;
    my $x0 = $pX->[$idx0];
    my $x1 = $pX->[$idx1];
    while(1) {
        my $idxh = int( ($idx0 + $idx1) /2 );
        my $xh   = $pX->[$idxh];
#print "x: $x ($x0,$xh,$x1) => ";
#        return $idxh if(($x1 - $x) * ($xh - $x) == 0.0);
        if(($x1 - $x) * ($xh - $x) <= 0.0) {
            $idx0 = $idxh;
            $x0 = $pX->[$idx0];
        }
        else {
#        elsif(($xh - $x) * ($x0 - $x) < 0.0) {
            $idx1 = $idxh;
            $x1 = $pX->[$idx1];
        }
#        else {
#            return undef;
#        }
#print "($x0 [$idx0],$x1 [$idx1])\n";
        last if($idx1 - $idx0 <= 1);
    }
#exit if($x > 1140);
    return $idx0;
}


def interpolate_simpson(x, y, x0):
    return $pY->[0] if($x <= $pX->[0]);
    my $nData = @$pX;
    return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

    my $idx  = FindIndexByBinarySearch($pX, $x);
    $idx = 1        if($idx <= 1);
    $idx = $nData-2 if($idx >= $nData-3);
#my $idxm1=$idx-1;
#my $idxp1=$idx+1;
#print "idx=$idx\tx=$x ($pX->[$idxm1],$pX->[$idx],$pX->[$idxp1])\n";
    my $y = InterpolateByLagrange( [@$pX[$idx-1..$idx+1]], [@$pY[$idx-1..$idx+1]], $x);
#print "\t\ty=$y ($pY->[$idxm1],$pY->[$idx],$pY->[$idxp1])\n";
    return $y;
}

sub InterpolateByCubicSpline
{
    my ($pX, $pY, $x) = @_;

    my $nData = @$pX;
    return undef           if($nData < 3);
    return $pY->[0]        if($x <= $pX->[0]);
    return $pY->[$nData-1] if($x >= $pX->[$nData-1]);

    for(my $i = 0 ; $i < $nData-1 ; $i++) {
        if($pX->[$i+1] <= $pX->[$i]) {
my $i1 = $i+1;
print "Warning in Algorism::InterpolateByCubicSpline: $i: d[i+1] must be > d[i] ($pX->[$i1] < $pX->[$i1])\n";
            return undef;
        }
    }

    my ($ph, $psp) = subspl($pX, $pY);

    for(my $i = 1 ; $i < $nData ; $i++) {
        if($pX->[$i-1] <= $x and $x < $pX->[$i]) {
            my $sm  = $psp->[$i-1];
            my $si  = $psp->[$i];
            my $hi  = $ph->[$i];
            my $hi2 = $hi * $hi;
            my $dxp = $pX->[$i] - $x;
            my $dxm = $x - $pX->[$i-1];
            my $v = $sm * $dxp * $dxp * $dxp
                  + $si * $dxm * $dxm * $dxm
                  + (6.0 * $pY->[$i-1] - $hi2 * $sm) * $dxp
                  + (6.0 * $pY->[$i]   - $hi2 * $si) * $dxm;
            return $v / $hi / 6.0;
        }
    }
    return undef;
}

sub subspl
{
    my ($pX, $pY) = @_;

    my $nData = @$pX;
    my (@h, @sp, @du, @dc);
    $du[$nData-1] = 0.0;
    $dc[$nData-1] = 0.0;
    $sp[$nData-1] = 0.0;
    for(my $i = 1 ; $i < $nData ; $i++) {
        $h[$i] = $pX->[$i] - $pX->[$i-1];
    }
    for(my $i = 1 ; $i < $nData-1 ; $i++) {
        my $hip1 = $pX->[$i+1] - $pX->[$i];
        $sp[$i] = 6.0 * (($pY->[$i+1] - $pY->[$i]) / ($h[$i] * $hip1)
                - ($pY->[$i] - $pY->[$i-1]) / ($h[$i] * $h[$i]));
        $du[$i] = $h[$i+1] / $h[$i];
        $dc[$i] = 2.0 * (1.0 + $du[$i]);
    }
    $dc[$nData-1] += 1.0;
    $dc[$nData-2] += $h[$nData-1] / $h[$nData-2];
    $du[1] /= $dc[1];
    $sp[1] /= $dc[1];
    for(my $i = 2 ; $i < $nData ; $i++) {
#print "$i: $dc[$i], $du[$i-1]\n";
        my $g = 1.0 / ($dc[$i] - $du[$i-1]);
        $du[$i] *= $g;
        $sp[$i] = ($sp[$i] - $sp[$i-1]) * $g;
    }
    for(my $i = $nData-2 ; $i >= 1 ; $i--) {
#print "$i: $sp[$i+1], $du[$i]\n";
        $sp[$i] -= $sp[$i+1] * $du[$i];
    }
    $sp[0] = $sp[1];
    $sp[$nData-1] = $sp[$nData-2];
    return (\@h, \@sp);
}

"""