#!/usr/bin/env python
# coding: utf-8

# <div style="font-size:35px"> Metode najmanjih kvadrata </div><br>

# Često se pri raznim mjerenjima pojavljuje veliki broj podataka. Da bi lakše protumačili te podatke pribjegava se grafičkom predstavljanju tih podataka. Na primjer pri izvodjenju nekog eksperimenta dobijeni su podaci 
# $$(x_0,y_0),\,(x_1,y_1),\dots,(x_7,y_7).$$
# 
# Ovakve podatke moguće je predstaviti grafikom, na primjer kao na sljedećoj slici.
# 
# 
# <center style="font-size:20px"><img src="Slike/najmanji1.png" width=800px>  
#     <caption > Podaci iz prethodne tabele predstavljeni grafički</caption> 
# </center>
# 
# Na prvi pogled grafik linearne funkcije najbolje bi odgovara ovim eksperimentalnim podacima. Sljedeći korak bi bio da odredimo linearnu funkciju 
# 
# $$y=ax+b,$$  
# 
# odnosno koeficijente $a$ i $b.$   Postavlja se pitanje: Koja je to prava koja je najbliže prolazi predstavljenim tačkama $(x_0,y_0)-(x_7,y_7)?$ ili u opštem slučaju tačkama $(x_0,y_0)-(x_n,y_n)?$ Pretpostavimo da su $a$ i $b$ tražene vrijednosti. Neka je $(x_k,y_k)$ jedna od datih tačaka i neka kroz ovu tačku ne prolazi prava $y=ax+b,$ tada je $y_k\neq ax_k+b,$ tj. $y_k-a x_k-b\neq 0.$  
# 
# 
# Da bi izračunali $a$ i $b$ posmatrajmo funkciju 
# $$\varphi(a,b)=\sum_{k=0}^{n}(ax_k+b-y_k)^2. \:(\star)$$ 
# 
# 
# Ako greška ima normalnu statističku raspodjelu, onda minimizacijom funkcije $\varphi$ dobijamo najbolje procjena za $a$ i $b.$  Računanje parametara $a$ i $b$ na ovakav način, odnosno linearne funkcije kojom ćemo aproksimirati naziva se $l_2$ aproksimacija. <br><br>
# 
# 
# <div style="border-left: 2px solid black;border-right: 1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom: 1px solid #E6E6E6;background-color:#FAFAFA; padding: 5px 5px 5px 15px" >
# <b> Napomena.</b>  Norma $l_p$ definisana je na sljedeći način 
# $$ || x ||_p=\left\{ \sum_{i=1}^{n}|x_i|^p \right\}^{1/p},\:1\leqslant p<\infty,$$
# za neki vektor $x=(x_1,\ldots,x_n)^T.$
# </div><br>
# 
# Odredimo sada minimum funkcije $\varphi.$ Potrebni uslovi za minimum funkcije $\varphi$ su
# 
# $$ \frac{\partial \varphi}{\partial a}=0,\quad  \frac{\partial \varphi}{\partial b}=0,$$ 
# 
# računajući prethodne uslove dobijamo
# 
# \begin{align*}
#  \sum_{k=0}^{n}{2(ax_k+b-y_k)x_k}=&0\\
#  \sum_{k=0}^{n}{2(ax_k+b-y_k)}=&0 ,  
# \end{align*}
# 
# Prethodne jednačine su po nepoznatim $a$ i $b$ i  možemo ih napisati u obliku
# 
# 
# \begin{align}
# \left(\sum_{k=0}^{n}{x^2_k}\right)a+   \left(\sum_{k=0}^{n}{x_k}\right)b=&\sum_{k=0}^{n}{x_k y_k}\\ 
# \left(\sum_{k=0}^{n}{x_k}\right)a+   (n+1)b=&\sum_{k=0}^{n}{y_k}, 
# \end{align}
# 
# 
# gdje je $\sum_{k=0}^{n}{1}=n+1.$
# 
# 
# Rješenje prethodnog sistema je 
# 
# 
# 
# \begin{align*}
#  a=& \frac{1}{d}\left[(n+1)\sum_{k=0}^{n}{x_k y_k}-\left(\sum_{k=0}^{n}{x_k}\right)\left(\sum_{k=0}^{n} {y_k}\right)  \right] \\
#  b=&\frac{1}{d}\left[\left(\sum_{k=0}^{n} {x^2_k}\right)\left(\sum_{k=0}^{n} {y_k}\right) 
#            -\left(\sum_{k=0}^{n}{x_k}\right)\left(\sum_{k=0}^{n} {x_ky_k}\right)  \right], 
# \end{align*}
# 
# 
# gdje je $$d=(n+1)\sum_{k=0}^{n}{x^2_k}-\left(\sum_{k=0}^{n}{x_k^2 }\right)^2.$$<br>
# 
# 
# 
# 
# 
# <div style="border-left:2px solid black;border-right:1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom:1px solid #E6E6E6; background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b> Zadatak.</b>  Izračunati za podatke 
#    $(0.5,2.25),(1,3.7),\,(2,4.1),\,(2.5,4.3),\,(3,5);$  linearnu funkciju $y=ax+b $ metodom najmanjih kvadrata.
# </div><br>

# In[1]:


from IPython import get_ipython;  # briše varijable, funkcije ...
get_ipython().magic('reset -sf')  #

import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import pyplot
from matplotlib import rc

pyplot.style.use('tableau-colorblind10')
newparams={'figure.figsize':(15,8),'font.size': 13,'axes.grid': True,'lines.markersize':10, 'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

xdata=[.5,1,2,2.5,3]
ydata=[2.25,3.7,4.1,4.3,5]

#xdata=np.linspace(1,15,200)
#def f1(x):
#    return 4*np.sin(3*x)*np.cos(14*x)+x
#ydata=f1(xdata)
broj=len(xdata)

def lin(x):
    xzbir=0; yzbir=0; xyzbir=0; xkvadrat=0; d=0; a=0; b=0;   
    for i in range(broj):
        xzbir=xzbir+xdata[i]
        yzbir=yzbir+ydata[i]
        xyzbir=xyzbir+xdata[i]*ydata[i]
        xkvadrat=xkvadrat+(xdata[i])**2
    
    d=(broj)*xkvadrat-xzbir**2    
    a=((broj)*xyzbir-xzbir*yzbir)/d
    b=(xkvadrat*yzbir-xzbir*xyzbir)/d
    print('Koeficijenti su a=%f i'%a, 'b=%f.'%b)
    return a*x+b

iks=np.linspace(np.min(xdata)-.25,np.max(xdata)+.25,1000)
fig=pyplot.figure()
pyplot.plot(xdata,ydata,'o',markersize=10)
pyplot.plot(iks,lin(iks),linewidth=3)
pyplot.legend(('Podaci','Linearna funkcija'))
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/najmanji2a.pdf')
pyplot.close(fig)
pyplot.show()


# <center style="font-size:20px"><img src="Slike/najmanji2a.png" width=1400px>  
#     <caption> Konstrukcija linearne funkcije za uzorak od 5 tačaka </caption> 
# </center>
# 
# <center style="font-size:20px"><img src="Slike/najmanji2.png" width=1400px>  
#     <caption> Konstrukcija linearne funkcije za uzorak od 200 tačaka </caption> 
# </center><br>

# <div style="font-size:30px"> Generalni slučaj </div>

# Umjesto računanja linearne funkcije $y=ax+b,$ metoda najmanjih kvadrata može se poopštiti na računanje funkcije $$y=\sum_{j=0}^{m}{c_jg_j(x)}$$ kojom ćemo predstaviti podatke koji su nam na raspolaganju. Funkcije $g_0,\ldots,g_m$ su poznate, fiksirane i linearno nezavisne.  Koeficijente  $c_0,\ldots,c_m$ ponovo treba izračunati koristeći prethodno izloženi princip najmanjih kvadrata. Drugim riječima, ponovo su poznate kordinate tačaka  $(x_k,y_k),\,i=0,\ldots,n;$ koje smo dobili mjerenjem, eksperimentalno ili na neki drugi način, i cilj nam je izračunati koeficijente $c_0,\ldots,c_m;$ tako da funkcija  
# \begin{equation}
# \varphi(c_0,\ldots,c_m)=\sum_{k=0}^{n}{\left[ \sum_{j=0}^{m}{c_jg_j(x_k)-y_k}\right]^2 }  
# \label{funkcija1}
# \end{equation}
# ima što je moguće manju  vrijednost. Potrebni uslovi za dobijanje minimuma funkcije $\varphi$ su 
# $$ \frac{\partial \varphi}{\partial c_i}=0,\:i=0,\ldots,m$$
# gdje su 
# $$\frac{\partial \varphi}{\partial c_i}=\sum_{k=0}^{n}{2\left[\sum_{j=0}^{m}{c_jg_j(x_k)-y_k} \right]}g_i(x_k),\:i=0,\ldots,m.$$
# Na kraju dobijamo
# $$ \sum_{j=0}^{m}{\left[ \sum_{k=0}^{n}{g_i(x_k)g_j(x_k)}\right]c_j}= \sum_{k=0}^{n}{y_kg_i(x_k)},\:i=0,\ldots,m.$$

# <div style="border-left:2px solid black;border-right:1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom:1px solid #E6E6E6; background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b> Zadatak.</b>  Izračunati za podatke $(-1.5,2.3),\, (-1,1.1),\,(-0.5,0.24),\,(0,0.01),\,(0.5,0.26),\,(1,1.2),\,(1.5,2.25)$ kvadratnu funkciju $y=c_2x^2+c_1x+c_0$ metodom najmanjih kvadrata.
# </div><br>

# In[8]:


from IPython import get_ipython;  # briše varijable, funkcije ...
get_ipython().magic('reset -sf')  #

import numpy as np
from matplotlib import pyplot
from numpy import linalg
from numpy import polynomial
from matplotlib import rc

pyplot.style.use('tableau-colorblind10')
newparams={'figure.figsize':(15,8),'font.size': 13,'axes.grid': True,'lines.markersize':10, 'lines.linewidth':2,'figure.dpi':300}
pyplot.rcParams.update(newparams)

xdata=np.array([-1.5,-1,-.5,0,.5 ,1,1.5])
ydata=np.array([2.3,1.1,.24,.01,.26,1.2,2.25])

brojTacaka=len(xdata)
brojBaznih=3

matricaBaznih=np.zeros((brojBaznih,brojTacaka))
matricaSistema=np.zeros((brojBaznih,brojBaznih))
vektorSlob=np.zeros(brojBaznih)

def bazne(x,stepen):
    return x**stepen

for i in range(brojBaznih):
    for j in range(brojTacaka):
        matricaBaznih[i,j]=bazne(xdata[j],i)
        
for i in range(brojBaznih):
    for j in range(brojBaznih):
        for k in range(brojTacaka):
            matricaSistema[i,j]=matricaSistema[i,j]+matricaBaznih[i,k]*matricaBaznih[j,k]

for i in range(brojBaznih):
    for k in range(brojTacaka):
        vektorSlob[i]=vektorSlob[i]+ydata[k]*matricaBaznih[i,k]      
        
koeficijenti=linalg.solve(matricaSistema,vektorSlob)        

pol=np.polynomial.Polynomial(koeficijenti)

iks=np.linspace(min(xdata),max(xdata),250)
y=pol(iks)
fig=pyplot.subplot()
pyplot.plot(xdata,ydata,'o',markersize=10)
pyplot.plot(iks,y)  
pyplot.legend(('Podaci','polinom $p_{10}$'),loc='upper center');
pyplot.savefig('/home/samir/Dropbox/Python/Zadaci/Slike/najmanji4d.png')
#pyplot.xlabel(' Broj tačaka je %i'%len(xdata));
#pyplot.show()
pyplot.close()


#  <center style="font-size:20px"><img src="Slike/najmanji3a.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_2(x)=c_2x²+c_1x+c_0$ za uzorak od 7 tačaka </caption> 
# </center>
# <center style="font-size:20px"><img src="Slike/najmanji4a.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_2(x)=c_2x²+c_1x+c_0$ za uzorak od 200 tačaka </caption> 
# </center>
# <center style="font-size:20px"><img src="Slike/najmanji4b.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_3(x)=c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka </caption> 
# </center>
# <center style="font-size:20px"><img src="Slike/najmanji4c.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_4(x)=c_4x^4+c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka</caption> 
# </center>
# <center style="font-size:20px"><img src="Slike/najmanji4e.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_5(x)=c_5x^5+c_4x^4+c_3x^3+c_2x^2+c_1x+c_0$ za uzorak od 200 tačaka </caption> 
# </center>
# 
# 
# 
# 
# <center style="font-size:20px"><img src="Slike/najmanji4d.png" width=1400px>  
#     <caption> Konstrukcija polinoma $p_{10}$ za uzorak od 500 tačaka</caption> 
# </center>

# <div style="border-left:2px solid black;border-right:1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom:1px solid #E6E6E6; background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b> Napomena 1.</b> Kao bazne funkcije umjesto $1,x,\ldots,x^n$ mogu se koristiti i druge funkcije.</div><br><br>
# 
# <div style="border-left:2px solid black;border-right:1px solid #E6E6E6;border-top: 1px solid #E6E6E6;border-bottom:1px solid #E6E6E6; background-color:#FAFAFA;padding: 5px 5px 5px 15px" >
# <b> Napomena 2.</b> Ovdje je predstavljena najjednostavija metoda. Druge metode i detaljnija analiza može se vidjeti na primjer u  [J97] ili [TB97].    </div><br>

# <div style="font-size:35px"> Literatura </div><br>
# 
# [CD04]  W. Cheney, D. Kincaid, <i>Numerical Mathematics and Computing</i>, Brooks/Cole-Thompson Learning, USA, 2004. ISBN: 053489937.<br><br>
# [J97] J.W. Demmel, <i>Applied Numerical Linear Algebra</i>, SIAM, USA, 1997. ISBN: 0-89871-389-7.
# 
# [TB97] L.N.Trefethen, D.Bau III, <i>Numerical Linear Algebra</i>, SIAM, USA, 1997. ISBN: 0-89871-361-7.