# Modul (soubor mypoly.py) s naší implementací polyval a polyfit
import numpy
import scipy
from numba import njit

# Vyčíslení polynomu p[0]*x**n + ... + p[n]
@njit
def polyval(p,x):
  result=numpy.empty(x.size)
  for nn in range(result.size):
    result[nn]=p[0]
    for n in range(1,p.size): result[nn]=result[nn]*x[nn]+p[n]
  return result

# Interpolační polynom p[0]*x**n + ... + p[n]
@njit
def polyfit2(x,y):
  n=x.size-1
  V=numpy.empty((n+1,n+1))         # Vandermondova matice
  # for i in range(n+1):           # pomocí cyklů
  #   V[i,0]=1
  #   for j in range(1,n+1): V[i,j]=V[i,j-1]*x[i]
  V[:,0]=1                         # totéž vektorově
  for j in range(1,n+1): V[:,j]=V[:,j-1]*x[:]
  c=numpy.linalg.solve(V,y)        # řešení algebraické soustavy
  return c[-1::-1]                 # zrcadlení vektoru

# ... druhá část souboru následuje v dalším odstavci
# ... druhá část souboru mypoly.py s naší implementací funkce polyfit

# Regresní polynom p[0]*x**m + ... + p[m], m=deg
@njit
def polyfit(x,y,deg):
  m=deg; n=x.size-1
  # if m==n: return polyfit2(x,y)  # interpolační polynom
  xx=numpy.ones(n+1); sumx=numpy.empty(m*2+1); sumy=numpy.empty(m+1)
  for i in range(m*2+1):
    sumx[i]=sum(xx)
    if i<=m: sumy[i]=sum(y*xx)
    xx*=x
  G=numpy.empty((m+1,m+1))         # Gramova matice
  for i in range(m+1):
    for j in range(m+1): G[i,j]=sumx[i+j]
  c=numpy.linalg.solve(G,sumy)     # řešení algebraické soustavy
  return c[-1::-1]                 # zrcadlení vektoru

# Regresní polynom p[0]*x**m + ... + p[m], m=deg
# varianta s funkcí numpy.linalg.lstsq() a Vandermondovou maticí
@njit
def polyfita(x,y,deg):
  m=deg; n=x.size-1
  V=numpy.empty((n+1,m+1))         # Vandermondova matice
  V[:,0]=1                         # totéž vektorově
  for j in range(1,m+1): V[:,j]=V[:,j-1]*x[:]
  c, res, rank, s= numpy.linalg.lstsq(V,y)
  return c[-1::-1]                 # zrcadlení vektoru

# Regresní polynom p[0]*x**m + ... + p[m], m=deg
# varianta s chybami (sig je pole standardních
# odchylek pro jednotlivé hodnoty y
@njit
def polyfite(x,y,sig,deg):
  m=deg; n=x.size-1
  sig2= 1.0/sig**2
  # if m==n: return polyfit2(x,y)  # interpolační polynom
  xx=numpy.ones(n+1); sumx=numpy.empty(m*2+1); sumy=numpy.empty(m+1)
  for i in range(m*2+1):
    sumx[i]= numpy.dot(xx,sig2)
    if i<=m: sumy[i]=numpy.dot(y*xx,sig2)
    xx*=x
  G=numpy.empty((m+1,m+1))         # Gramova matice
  for i in range(m+1):
    for j in range(m+1): G[i,j]=sumx[i+j]
  c=numpy.linalg.solve(G,sumy)     # řešení algebraické soustavy
  return c[-1::-1]                 # zrcadlení vektoru