###############################################################################
## TIEDOSTO: visualisointeja.m
##
## Numeerisen optimoinnin havainnollistamisia kurssille MATH2040 kevaalla 2010
## Vaasan yliopistossa.
##
## Tarvittavat itsemaaritellyt funktiot:
##    newton_stat
##      num_grad
##      num_hesse
##    jyrkka_max
##      num_grad
##     
## (copyleft) 2010 Tommi Sottinen <tommi.sottinen(miuku)uwasa.fi>
###############################################################################

## Seuraava "tyhja" komento estaa Octavea tulkitsemasta tata tiedostoa 
## funktiotiedoksi.
1;

###############################################################################
## Esimerkkifunktio
###############################################################################

function y = kumpuja(x) 
  y = cos(x(1))*cos(x(2))/(x(1)^2+x(2)^2+5); 
##y = exp(-x(1)^4-x(2)^2);  ## Talla Neton nayttaa menevan pieleen
endfunction

###############################################################################
## Visualisointi
###############################################################################

## Esimerkkifunktio samplataan matriisiksi

n=32;                                            # Samplauspisteiden lkm.  
x1_min = -3.25;                                  # Sampauspisteiden rajat
x1_max =  3.25; 
x2_min = -3.25;
x2_max =  3.25;

x1_hila = linspace(x1_min,x1_max,n);             # Hilat lineaarisesti
x2_hila = linspace(x2_min,x2_max,n);
z=zeros(length(x1_hila),length(x2_hila));        # Samplausmatriisin alustus

for i=1:n
  for j=1:n
    z(i,j) = kumpuja([x1_hila(i),x2_hila(j)]');  # Itse samplaus
  endfor
endfor


###############################################################################
## Havainnollistetaan esimerkkifunktiota
###############################################################################
## hold off;                                     # Vapautetaan kuva (Paranoiaa)
## surf(x1_hila,x2_hila,z);                      # Pinnan z kuva hilassa
###############################################################################

## Piirretaan esimerkkifunktion tasa-arvokayrasto, mutta sita
## ennen varmuuden vuoksi vapautetaan/pyyhitaan kuvapohja 
hold off;                     
contour(x1_hila,x2_hila,z);

## Piirretaan jyrkimman nousun optimointi
hold on;                                         # Kaikki yhteen kuvaan
axis([x1_min, x1_max,  x2_min, x2_max]);         # Kiinnitetaan nakyma 

m = 7;                                           # Iteraatioiden maara
s = 1.0;                                         # Hypyn (suhteellinen) koko

##Iteraation alkupiste (Newton parempi)
x0 = [0.7 0.9]';                                 
##Iteraation alkupiste (Newton lahtee "vaaraan" suuntaan) 
##x0 = [1 1.5]';                                   
##Iteraation alkupiste (Newton "ampuu yli")
##x0 = [0.3 0.7]';                                 

[x,p] = newton_stat(@kumpuja,x0,1,m,s);          # Newtonin menetelma
[y,q] = jyrkka_max(@kumpuja,x0,1,m,s);           # Jyrkan nousun menetelma

## Satunnaistuksia visualisointia varten
p = p + 0.02*rand(2,m);
q = q + 0.02*rand(2,m);

## Piirretaan Newtonin polku punaisilla o-viivoilla
plot(p(1,:),p(2,:),'r-o');                       
## Piirretaan jyrkimman nousun polku sinisilla *-viivoilla
plot(q(1,:),q(2,:),'b-*');                       
hold off;                                        # Vapautetaan kuvapohja