import csv
import sys
import numpy as np
from numpy import sin, cos, tan, pi
from pprint import pprint


from tklib.tkutils import pint, pfloat, safe_getelement, check_attributes, merge_attributes, joinf, save_csv



def convolve_func(x, x0, whalf):
    if func_type == 'gauss':
        return Gaussian(x, x0, Wa)
    else:
        return Lorentzian(x, x0, Wa)
    
def Ridge(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, is_print = False):
    n = len(x)
    Si  = np.empty([m, 1])
    Sij = np.empty([m, m])

    def cal_w(j):
        if j <= nGmax:
            wj = w[j]
        elif n - 1 - j < nGmax:
            wj = w[n-1-j]
        else:
            wj = w[nGmax - 1]
        return wj
        
    for l in range(0, m):
        v = 0.0
        for i in range(0, n):
            v += y[i] * lsqfunc(x[l], x[i])
#        Si[l, 0] = v
        Si[l, 0] = v / cal_w(l)

    for j in range(0, m):
        for l in range(j, m):
            v = 0.0
            for i in range(0, n):
                v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
#            Sij[j, l] = Sij[l, j] = v
            Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)
        Sij[j, j] += alpha

    if is_print:
        print("Vector and Matrix:")
        print("Si=")
        print(Si)
        print("Sij=")
        print(Sij)
        print("")

    ci = np.linalg.inv(Sij) @ Si
    ci = ci.transpose().tolist()
    return ci[0]

def Smoothing_Penalty_Regression_by_base(x, y, m, lsqfunc, alpha = 0.0, nsmooth = 1, is_print = False):
    n = len(x)
    Si  = np.empty([m, 1])
    Sij = np.empty([m, m])

    for l in range(m):
        Si[l, 0] = sum([y[i] * lsqfunc(l, x[i]) for i in range(n)])

    for j in range(m):
        for l in range(j, m):
            v = sum([lsqfunc(j, x[i]) * lsqfunc(l, x[i]) for i in range(n)])
            Sij[j, l] = Sij[l, j] = v

    nband = int(nsmooth / 2 + 1.0001)
    use_ridge = 1
    if use_ridge > 1:
        for j in range(0, m):
            Sij[j, j] += alpha
    elif use_ridge > 0:
        for j in range(0, m):
            j0 = j - nband
            j1 = j + nband
            if j0 < 0:
                j0 = 0
            if j1 >= m:
                j1 = m
            for k in range(j0, j):
                Sij[j, k] -= alpha
            for k in range(j+1, j1):
                Sij[j, k] -= alpha

    if is_print:
        print("Vector and Matrix:")
        print("Si=")
        print(Si)
        print("Sij=")
        print(Sij)
        print("")

    ci = np.linalg.inv(Sij) @ Si
    ci = ci.transpose().tolist()

    return ci[0], Si, Sij


def Smoothing_Penalty_Regression(x, y, m, lsqfunc, w, nGmax, alpha = 0.0, nsmooth = 1, is_print = False):
    n = len(x)
    Si  = np.empty([m, 1])
    Sij = np.empty([m, m])

    def cal_w(j):
        if j <= nGmax:
            wj = w[j]
        elif n - 1 - j < nGmax:
            wj = w[n-1-j]
        else:
            wj = w[nGmax - 1]
        return wj

    for l in range(0, m):
        v = 0.0
        for i in range(0, n):
            v += y[i] * lsqfunc(x[l], x[i])
#        Si[l, 0] = v
        Si[l, 0] = v / cal_w(l)

    for j in range(0, m):
        for l in range(j, m):
            v = 0.0
            for i in range(0, n):
                v += lsqfunc(x[j], x[i]) * lsqfunc(x[l], x[i])
#            Sij[j, l] = Sij[l, j] = v
            Sij[j, l] = Sij[l, j] = v / cal_w(j) / cal_w(l)

    nband = int(nsmooth / 2 + 1.0001)
    use_ridge = 1
    if use_ridge > 1:
        for j in range(0, m):
            Sij[j, j] += alpha
    elif use_ridge > 0:
        for j in range(0, m):
            j0 = j - nband
            j1 = j + nband
            if j0 < 0:
                j0 = 0
            if j1 >= m:
                j1 = m
            for k in range(j0, j):
                Sij[j, k] -= alpha
            for k in range(j+1, j1):
                Sij[j, k] -= alpha

    if is_print:
        print("Vector and Matrix:")
        print("Si=")
        print(Si)
        print("Sij=")
        print(Sij)
        print("")

    ci = np.linalg.inv(Sij) @ Si
    ci = ci.transpose().tolist()
    return ci[0]


def mlsq(x, y, m, iPrint = 0):
        """
            LSQ for m-th order polynomial
            x, y: data to be fitted
            m: order of polynomial
            iPrint: output control [0|1]
        """

        n = len(x)
        Si  = np.empty([m+1, 1])
        Sij = np.empty([m+1, m+1])

        for l in range(0, m+1):
             Si[l, 0] = sum([y[i] * pow(x[i], l) for i in range(n)])

        for j in range(0, m+1):
            for l in range(j, m+1):
                v = sum([pow(x[i], j+l) for i in range(n)])
                Sij[j, l] = Sij[l, j] = v

        if iPrint == 1:
            print("tkoptimize.mslq:: Vector and Matrix:") 
            print("Si=")
            pprint(Si)
            print("Sij=")
            pprint(Sij)
            print("")

        ci = np.linalg.inv(Sij) @ Si
        ci = ci.transpose().tolist()

        return ci[0], Si, Sij

def mlsq_general(x, y, m, lsqfunc, iPrint = 0):
        """
            LSQ for linear combinationof m functions
            x, y: data to be fitted
            m: number of functions
            lsqfunc(i, x): functions
            iPrint: output control [0|1]
        """

        n = len(x)
        Si  = np.empty([m, 1])
        Sij = np.empty([m, m])
#        print("n=", n, "  m=", m)

        for l in range(m):
            Si[l, 0] = sum([y[i] * lsqfunc(l, x[i]) for i in range(n)])

        for j in range(m):
            for l in range(j, m):
                v = sum([lsqfunc(j, x[i]) * lsqfunc(l, x[i]) for i in range(n)])
                Sij[j, l] = Sij[l, j] = v

        if iPrint == 1:
            print("tkoptimize.mslq_general:: Vector and Matrix:") 
            print("Si=")
            pprint(Si)
            print("Sij=")
            pprint(Sij)
            print("")

        ci = np.linalg.inv(Sij) @ Si
        ci = ci.transpose().tolist()

        return ci[0], Si, Sij

