(*tenhle textik je o tzv. linearni regesi, o graficich a zobrazovani dat*)
(*"linearni regrese" tomu rikame mimojine proto, ze hledame koeficienty linearni kombinace zvolenych funkci , koeficienty takove, aby vysledna linearni kombinace co nejlepe (ve zvolene norme) pasovala na zadana data*)
(*jak to udelat v pripade, ze pozadovana nahrada neni linearni kombinace zadanych /zvolenych/ funkci, tak to si ukazeme pozdeji:-)*)

(*hlavne je to nalezni vzorcu... ale je zde i jak udelat interpolacni funkci a zobrazit ji*)

(*uvedomme si: interpolaci pouzivame, kdyz nemame vzorec a verime namerenym datum*)

(*proklad pouzivame kdyz verime hodne vzorci a ne az tak datum*)

(*jak poznat co pouzit v inzenyrske praxi? Hlavne znalosti proslusne teorie, napr. nabijeni kondenzatoru pred rezistor s konstantnim odporem ze zdroje konstantniho napeti bude exponenciela, tedy pouzijeme vzorecek a regresi, pro zavislost hustoty vody na teplote pri konstantnim tlaku pouzijeme tabulky a interpolaci a podobne.*)

(*vycisteni promennych*)
Quiet@Remove["Global`*"];
(*nastaveni adresare, kde je tento soubor, jako defaultniho *)
SetDirectory[NotebookDirectory[]]

(*tady si jen vytvorime data, se kterymi budeme pracovat.... vlastne jen vam spis chceme ukazat, jak je tvorit a exportovat... ty listy si muzete i napsat rucne, napriklad
data= {{-8.`,2.76},{-6.5`,1.197},{-5.`,0.3579`},{-3.5`,-0.28616},{-2.`,-0.32},{-0.5`,0.3646},{1.`,2.1415},{2.5`,3.51},{4.`,5.2787}} atd...*)
(*nastavime sadu pseudonahodnych cisel, delame to proto, ze obvykle program ladime, tedy spoustime mnohokrat s upravami a kdybychom teohle neudelali, dostavali bychom ruzna cisla a nevedeli bychom na prvni pohled, jestli se zmena stala tim, ze jsme zmenili kod, nebo proto, ze se zmenily vstupni hodnoty*)
(*temi pseudonahodnymi cisly simulujeme chyby mereni, sumy a rusive signaly a podobne*)
SeedRandom[1];
(*rn definujeme pomoci SetFelayed, tedy ":=", pri kazdem pouziti se vygeneruje nove nahodne cislo*)
rn:=0.5*RandomReal[{-1,1}];
(*vygenerujeme si dvojice hodnot (v grafu "{x,y}"), vyraz jsme zvolili nahodne, nechteli jsme opakovat notoricky se vyskytujici polynomy*)
data={#,1+0.5*#+0.1((#)^2)+0.7*Tanh[#]+rn}&/@Range[-8,5,1.5];
(*nez s daty neco zacneme delat tak pokud to je mozno si je prohledneme, klidne obycejnym ListPlotem bez popisu os, je to pro nas, ne vystup ven, kde by se to netolerovalo, odstranite-li strednik, uvidite je*)
ListPlot[data];
(*abychom si pozdeji ukazali, jak se data importuji, tak si je vyexportujeme*)
(*Wolfram pozna pozadovany format podle pripony*)
Export["datacsv.csv",data]
Export["dataxls.xls",data]

(*nez psat nazev souboru, je lepsi a spolehlivejsi ho z kopirovat z vystupu FileNames[]... jen nesmime byt hlupaci a  mit vsechno na Plose, pak vystup FileNames[] obsahuje mnoho jmen souboru a moc to nepomuze...*)
(*FileNames[] , coz muzeme cist "FileNames bez parametru" zobrazi nazvy souboru v defaultnim adresari, kdybychom napsali jako parametr cestu k adresari, dostaneme soubory podle zadane cesty*)
FileNames[]

dat1=Import["datacsv.csv"];
(*format Excelu predpoklada vice sesitu... a kazdy odpovida jednomu Listu co vraci Import, proto bereme Import["dataxls.xls"][[1]]  , protoze bereme prvni (a jediny) prvek importovaneho listu... no proste si to zkuste:-)*)
dat2=Import["dataxls.xls"][[1]];
dat1-dat2//Union
(*cili je to stejne, namisto abychom si prohlizeli data a koukali, jestli mame shodne hodnoty, tak ty listy po prvcich odecteme a Union provede sjednoceni listu a poskrtani duplicit*)
(*jelikoz vysly same dvojice {0.,0.} tak nam to vrati {{0.,0.}} a vidime, ze data jsou shodna*)

(*mame tedy data ve tvaru listu dvojic {(Subscript[x, i]),f[(Subscript[x, i])]}, kde predpokladame (tady vime, protoze jsme si to vymysleli), ze f[x_]:=a+b*x+c*(x^2)+d*Tanh[x]*)
(*tedy vime, ze jde o linerni kombinaci:
- jednicky (je "a krat")
-ikska (je "be krat")
-(x^2) (je "ce krat")
-Tanh[x] (je "d krat")
*)
(*proto bude na konci syntaxe 
prokladFit[x_]=Fit[dat2, {1,x,(x^2),Tanh[x]},x]
*)
(*{1,x,(x^2),Tanh[x]} je baze linearniho prostoru funkci (jasne, zde spis vyrazu... ale dovolme si byt nepresni, jako matematici:-))*)

f[x_]:=a+b*x+c*(x^2)+d*Tanh[x];

(*Δ je odchylka, bod:{xMerene_,fMerene_} je jedna dvojice {x,y} z merenych ci jinak ziskanych dat*)
(*vsimnete si, ze pro ciselne {x,y}={1,1} je linearni formou symbolickych parametru {a,b,c,d}*)
(*pojmenovali jsem si parametr bod v  Δ[bod:{xMerene_,fMerene_}]:=fMerene-f[xMerene] castecne z pleziru, castecne pro pochopeni... dulezite to neni, ale verte mi, ze kdybychom pozdeji delali nejake NestList a podobne, tenhle format promenne nam usnadni premysleni a programovani*)
Δ[bod:{xMerene_,fMerene_}]:=fMerene-f[xMerene];
Δ[{1,1}]

(*odchylka Δ[{xi,yi}] muze byt kladna i zaporna, takze zadna jeji hodnota nemuze byt kriteriem optimality, tedy minima celkove odchylky... proste odchylky "nahoru" by se vyrusily s odchylkami "dolu" a vysledna funkce by se nam vubec nelibila*)
(*zavedeme tedy funkci norm, ktera bude vracet...takovou hodnotu, aby, chapeme-li ciselnou hodnotu Δ[{xi,yi}] jako vektor, splnovala pozadavky na tzv. normu vektoru... viz. napr. skvele texty kolegy doc. Habaly*)
(*nam staci, ze norm[x]≥0 pro kazde x, norm[x]==0 prave jen pro x==0*)
(*zvolime nejjednodussi norm, tedy eukleidovskou druhou mocninu*)
(*ale Abs by fungovalo taky a mnoho dalsich norem... ale jak uvidime, vysledne rovnice bro minimum by nebyly linearni a take jde o to, co chceme, jaky pomerny vliv maji meale odchylky vyci velkym *)
(*z algebry vite, ze eukleidovska norma popisuje svet, jak ho bezne vidime a pokud nejsou vzdalenosti tak velike, ze by prostor prestal byt nezakriveny, staci to a odpovida vizualni zkusenosti*)
(*a ostatne takhle to puvodne genius Gauss odvodil*)

norm=((#)^2)&;(*to je ta druha mocnina*)

(*takze normovana odchylka se spocte takhle*)
normovanaOdchylka[bod:{xMerene_,fMerene_}]:=norm[Δ[bod]];
(*s touto normou je to pro {xMerene,fMerene} konkretni kvadraticka forma koeficientu {a,b,c,d}, kazme si to na {xMerene,fMerene}={1,1}*)
normovanaOdchylka[{1,1}]

(*namapuji-li funkci normovanaOdchylka na vstupni data, treba dat2, dostanu list normovanych odchylek... kdyz date strednik pryc, muzete si ho prohlednout *)
normovaneOdchylky=normovanaOdchylka/@dat2;

(*tyto normovane odchylky secteme*)
(*a zjednodusime to..*)
soucetNormovanychOdchylek=Simplify@Total[normovaneOdchylky]

(*ted kousek programovani... jasne, vime, ze jde o nalezeni koeficientu {a,b,c,d}... jenze co kdyz bychom menili modelovou funkci a i jmena koeficientu? Radeji si je najdeme z promenne soucetNormovanychOdchylek*)
(*zkuste si odtranit Union@ ...*)
(*posledni parametr Cases, tedy {0,∞} znamena "prohledej tu stromovou strukturu v pozadi celou"*)
proms=Union@Cases[soucetNormovanychOdchylek,_Symbol,{0,∞}]

(*minimum hledame (jde-li to, ono to tak skoro nikdy krome skolskych pripadu a prave takovychhle prokladu nejde, ale to az v magistru v predmetu IAP, kdyz spolu zustaneme na svete a v EEM:-)) tak, ze derivujeme podle promennych a derivaci davame rovnou nule*)
(*udelame si tedy funkci, ktera vrati rovnici, ze derivace podle parametru tehle funkce je rovna nule*)
rceMin:=0==D[soucetNormovanychOdchylek,#]&;

(*fumkci rceMin namapujeme na promenne proms a vysledek nazveme rceMinParams*)
rceMinParams=rceMin/@proms

(*tady jsme nechali vypis... vidime totiz, ze volba normy jako druhe mocniny vede na soustavu linearnich rovnic... a pritom popisuje nase eukleidovske videni sveta... docela zajimave, ale spis tema do hospody... ale ty jsou zavrene:-)*)
(*tak je vyresime*)
resMin=Solve[rceMinParams][[1]]
(*zase beru jen vnitrni list, Solve je pripraveno na vice reseni...*)

(*ted pouzijeme nalezene koeficienty*)
(*takhle bychom to udelali funkce, nicmene neprakticky, protoze musime kopirovat cely vyraz, potencialni zdroj omylu*)
fVyslednaHloupe[x_]:=a+b*x+c*(x^2)+d*Tanh[x]/.resMin;
(*kdyz napiseme funkci, ktera ma vracet pro ciselny parametr cislo, je dobre si ji zavolat a vyzkouset si to*)
fVyslednaHloupe[1]

(*programatorsky lepsi, vyraz jsmem zadavali jen jednou, ted se na nej odkazeme...*)
fVysledna[x_]:=f[x]/.resMin;
fVysledna[1]

(*kdykoli delame nejake proklady, vzdycky se puvodni body a nalezeny proklad zobrazime v jednom grafu*)
(*neni spatne si zobrazit nejen proklad, ale i interpolaci dat*)
(*navic z te interpolace snadno ziskame i meze pro Plot*)
interpolace=Interpolation[dat2,InterpolationOrder->3];

(*vidite, ze objekt interpolace obsahuje informaci, kde plati... zkuste TreeForm a uvidite... my to vyuzijeme pro zminene meze grafu*)
{xmin,xmax}=interpolace[[1,1]]

(*zapis {xmin,xmax}=interpolace[[1,1]] je jen pro usporu jednoho radku a nejakych zavorek... ve Wolframu funguje
{a,b}={aSet,bSet}
*)

(*zobrazime si puvodni (merene) body*) 
(*opet jsme to nejprve zkusili bez stredniku a pak strednik pridali*)
(*grafiku jsme si ulozili do promenne plotBody, bude se to hodit az pomoci Show[gr1,gr2,..] budeme zobrazovat vice grafu v jednom obrazku*)
(*a opet, jako vzdy, nejprve jsme si ji zobrazili bez stredniku a ten az kdyz se nam graf libil jsme odstranili*)
plotBody=ListPlot[dat2,AxesLabel->{"x","f[x]"},PlotStyle->{Red,PointSize[0.01]},GridLines->Automatic,PlotRange->All];

(*interpolace prochazi body a mezi body je vycislena jako hodnota polynomu radu InterpolationOrder*)
(*tedy je-li nastaveno InterpolationOrder->order, tak se kazdych order+1 bodu proklada polynomem radu order, vime, ze polynom 1. radu je jednoznacne definovan dvema ruznymi body atd.*)
(*o stredniku plati samosebou stale totez*)
(*vsimnete si, ze nepopisujeme osy a podobne... kdyz bude Show mit jako prvni parametr grafiku plotBody, pouziji se popisy aos a podobne tak, jak je nastaveno v plotBody*)
plotInterpolace=Plot[interpolace[x],{x,xmin,xmax},PlotStyle->Black];

(*plot prokladu*)
plotProkladu=Plot[fVysledna[x],{x,xmin,xmax},PlotStyle->Blue];

(*ted se uz jen podivame na vsechno dohromady*)
grafik=Show[plotBody,plotInterpolace,plotProkladu]
(*vidime, ze interpolace prochazi body, proklad nikoli: v pripade interpolace verime bodum, v pripade prokladu vzorci*)

(*no a abychom to mohli prezentovat na pocitacich, kde neni Wolfram, exportujeme si grafik*)
(*pro carovou grafiku se hodi format png a obecne dobre funguje pdf*)
Export["grafik.pdf",grafik,ImageResolution->300]

(*poznamenejme ale, ze ve Wolframu to pokud volime norm=((#)^2)&, nemuzime to delat tak slozite, uz je to pro nas pripraveno*)
(*Chop odstrani "skoronuly" a vidime, ze vysledek je stejny*)
(*SetDelayed, tedy ":="  nefunguje, zkuste si  prokladFit[x_]:=Fit[dat2, {1,x,(x^2),Tanh[x]},x]... nemuze fungovat, protoze se neprovede proklad v okamziku definice... kdy je "x" symbolicka promenna a jde proklad delat. V okamziku volani ma "x" konkretni hodnotu a nejde ho chapat jako promennou... atd... neni to uplne jednoduche, doporucujeme ten zakladni vyukovy text na powerwiki.cz:-)*)
prokladFit[x_]=Fit[dat2, {1,x,(x^2),Tanh[x]},x]
fVysledna[x]-%//Chop

(*V minulem jsme ukazali, jak najit nejlepsi linearni kombinaci zadanych funkci, tak aby ve smyslu dane normy prokladala zadana data*)

(*No... je to super, pokud nam staci eukleidovska norma, tedy norm=((#)^2)&  a spokojime se s tim, ze hledame koeficienty linearni kombinace zadanych bazovych funkci, ziskali jsme resni jednoznacne a rychle: reseni soustavy linearnich rovnic je dost dobre zvladnuta problematika a jde prevest na male maticove soustavy atd. *)

(*Jenze obcas to takhle nejde... napriklad nabijeni/vybijeni kondenzatoru s konstantni kapacitou pres rezistor s konstantnim odporem a nabiji/vybiji se z/do zdroje napeti exponencialne podle vzorce:
u[t_]:=u0+(Uz-u0)*(1-Exp[((-t)/(R*c))])*)

Quiet@Remove["Global`*"];
gridlines=Automatic;
R=100;
μF=(10^(-6));
c=4.7μF;
Uz=10;
τ=R*c;
u0=3;
uMer[t_]:=u0+(Uz-u0)*(1-Exp[((-t)/τ)]);
Plot[uMer[t],{t,0,5τ},AxesOrigin->{0,0},AxesLabel->{"t [s]","u [V]"},PlotStyle->{Blue, Thickness[0.004]},GridLines->gridlines]

(*osciloskopem ale sejmeme jinaci data, svetem leta ruseni elektricke a magneticke, to vas (snad) nauci kolegove... no, bude to vypadat nejak takhle*)
tmax=6.τ;
SeedRandom[1];
rn:=0.2*RandomReal[{-1,1}];
dat2={#,uMer[#]+rn}&/@Range[0,tmax,(tmax/200.)];
plOsc=ListPlot[dat2,AxesLabel->{"t [s]","u [V]"},PlotStyle->{Yellow, PointSize[0.007]},GridLines->gridlines,Background->GrayLevel[0.2],Joined->True]

(*verime vzorecku u[t_]:=u0+(uinf-u0)*(1-Exp[((-t)/τ)]), ale nezname hodnoty {u0, uinf, τ}*)
(*{u0m,Uzm,τm} jsou "merene" hodnoty, proto je oznacuji _<>"m"*)
u[t_]:=u0m+(Uzm-u0m)*(1-Exp[((-t)/(τm))]);
Δ[bod:{tMerene_,uMerene_}]:=uMerene-u[tMerene];
Δ[{1,1}]

(*odchylka Δ[{xi,yi}] muze byt kladna i zaporna, takze zadna jeji hodnota nemuze byt kriteriem optimality, tedy minima celkove odchylky... proste odchylky "nahoru" by se vyrusily s odchylkami "dolu" a vysledna funkce by se nam vubec nelibila*)
(*zavedeme tedy funkci norm, ktera bude vracet...takovou hodnotu, aby, chapeme-li ciselnou hodnotu Δ[{xi,yi}] jako vektor, splnovala pozadavky na tzv. normu vektoru... viz. napr. skvele texty kolegy doc. Habaly*)
(*nam staci, ze norm[x]≥0 pro kazde x, norm[x]==0 prave jen pro x==0*)
(*zvolime nejjednodussi norm, tedy eukleidovskou druhou mocninu*)
(*ale Abs by fungovalo taky a mnoho dalsich norem... jde o to, co chceme:-)*)
(*jelikoz nebudeme odkazani na hledani minima resenim rovnic pro nulovost derivaci podle koeficientu, mame ve volbe normy velkou volnost*)

norm=((#)^2)&;(*to je ta druha mocnina*)
(*takze normovana odchylka se spocte takhle*)
normovanaOdchylka[bod:{tMerene_,uMerene_}]:=norm[Δ[bod]];
(*s touto normou je to kdadraticka forma koeficientu {u0m,Uzm,τm}*)
normovanaOdchylka[{1,1}]

(*namapujeme-li funkci normovanaOdchylka na vstupni data, treba dat2, dostaneme list normovanych odchylek... kdyz date strednik pryc, muzete si ho prohlednout *)
normovaneOdchylky=normovanaOdchylka/@dat2;
(*zcela analogicky linearnimu pripadu*)
soucetNormovanychOdchylek=Total[normovaneOdchylky];
(*ted kousek programovani... jasne, vime, ze jde o koeficientu {u0m,Uzm,τm}... jenze co kdyz bychom menili modelovou funkci a i jmena koeficientu? Radeji si je najdeme z promenne soucetNormovanychOdchylek*)
(*zkuste si odtranit Union@ ...*)
(*posledni parametr Cases, tedy {0,∞} znamena "prohledej tu stromovou strukturu v pozadi celou"*)
(*klidne se na tu syntaxi vykaslete, pokud se nenudite... takhle jsem se zbavil "cisel" typi Pi a E a podobne, ktera jsou typu Symbol... *)
(*((Head[#]==Symbol)&&!NumericQ[#]&) je boolska funkce, zkuste si jizavolat na ruzne parametry   a _?boolskaFce dela omezeni patternu na to, pro co boolskaFce vraci True*)
(*!NumericQ je klasicky negace boolske funkce NumericQ*)
proms=Union@Cases[soucetNormovanychOdchylek,_?((Head[#]==Symbol)&&!NumericQ[#]&),{0,∞}]

(*najdeme minimum, vraci strukturu {hodnotaMinima,{prom1->val1...}}*)
min=NMinimize[soucetNormovanychOdchylek,proms]

(*potrebujeme je list dosazeni z hledani minima, tedy druhy prvek vystupu NMinimize, min[[2]]*)
prokladNapeti[t_]:=u[t]/.min[[2]];
plotProklad=Plot[prokladNapeti[t],{t,0,6τ},AxesOrigin->{0,0},AxesLabel->{"t [s]","u [V]"},PlotStyle->{Red, Thickness[0.004]},GridLines->gridlines];
Show[plOsc,plotProklad]

(*takhle to delat nemusime, mame pripravenou funkci FindFit*)
model=u0m+(Uzm-u0m)*(1-Exp[((-t)/(τm))]);
fit=FindFit[dat2,model,proms,t]

(*skonci to stejne... ale v predchozim postupu to mame mnohem vice pod kontrolou, v ruce...*)
plotFindFit=Plot[model/.fit,{t,0,6τ},AxesOrigin->{0,0},AxesLabel->{"t [s]","u [V]"},PlotStyle->{Red, Thickness[0.004]},GridLines->gridlines];
Show[plOsc,plotFindFit]