#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 25 17:37:54 2024

@author: laurentdumas
"""
import numpy as np
import matplotlib.pyplot as plt
phi1 = lambda h: h
phi2 = lambda h: h**3
phi3 = lambda h: np.exp(-h**2)
#C
def C(x, phi):
    n = len(x)
    C = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            C[i][j] = phi(abs(x[i]-x[j]))
    return C

#P
def P(x, k): 
    n = len(x)
    P = np.zeros((n, k))
    for i in range(n):
        for j in range(k):
            P[i,j] = x[i]**j
    return P

#M
def M(x, phi,C, k):
    n = len(x)
    M = np.zeros((n+k, n+k))
    M[0 : n , 0 : n]  = C(x,phi)
    M[0 : n , n : n+k]= P(x,k)
    M[n:n+k , 0 : n]  = P(x,k).T
    return M

#alpha & a
def coeff(M,U):
    u = np.zeros(len(M))
    u[:len(U)] = U
    X = np.linalg.solve(M, u)
    alpha = X[:len(U)]
    a = X[len(U):]
    return alpha, a

# u
def u(X, x, U, phi,k):
    def d(x,a):
        s=0.
        for i in range (0,len(a)):
            s+= a[i]*x**i
        return s
    def w(X,x,alpha, phi):
        s=0.
        for i in range (len(alpha)):
            s+= alpha[i]*phi(abs(x-X[i]))
        return s
    Mat = M(x, phi, C,k)
    alpha, a = coeff(Mat, U)
    return d(X, a) + w(x, X, alpha, phi)
I = np.array([1,2,5,7,8,11,13,16])
x=I
U = np.array([4,6,9,10,8,7,9,5])
E = np.linspace(0, max(x)+1, 500)
#I = np.linspace(0, 4*np.pi, 10)
#U = np.sin(x)
#U=1./(1.+x**2)
#E = np.linspace(0, 4*np.pi, 300)


u1 = u(E, I, U, phi1, 2)
u2 = u(E, I, U, phi2, 2)
u3 = u(E, I, U, phi3, 2)
plt.plot(E, u1, label="h")
plt.scatter(x,U)
plt.show()
plt.plot(E, u2, label="h3")
plt.scatter(x,U)
plt.show()
plt.plot(E, u3, label="exp(-h2)")
plt.scatter(x,U)
plt.show()