%
% MATLAB for DUMMIES 11-12
% Exercice 7
%
% Solution detaillee
%  Vincent Legat, Damien Bouvy, Denis-Gabriel Caprace, Gilles De Neyer, 
%  Kevin Degraux, Mohammad Iranmanesh, Christophe Lorant, Remy Moulin, 
%  Matthew Philippe, Karim Slaoui, Olivier Thiry, Laurent Van Eesbeeck,
%  Thomas Vanderstraeten, Antoine Weynants
%
% -------------------------------------------------------------------------
%
function matlab7()
    global count 

    close all;
    functions = {@easy,@difficult}; 
    for i=1:length(functions) 
        u = functions{i}; 

        fprintf('\n -----------------------------------------------------------------------------');    
        count = 0; 
        [x,y,z] = minimumSearch(u);
        fprintf('\n ===> (x,y,z) = %14.7e %14.7e %14.7e  ',x,y,z)
        fprintf('\n Count %6d : Obtained value                                = %14.7e ',count,u(x,y,z))

        count = 0; 
        xi = -2:.1:2;
        [X Y Z] = meshgrid(xi,xi,xi); 
        U = u(X,Y,Z);
        figure; slice(X,Y,Z,U,[-1.2,.8,2],[2],[-2])
        fprintf('\n Count %6d : Global minimum with a stupid greedy method    = %14.7e \n',count,min(min(min(U))));
        fprintf('\n -----------------------------------------------------------------------------\n');


    end
end


function [x,y,z] = minimumSearch(u)
%MINIMUMSEARCH Minimizes the function u with maximum 200 evaluations
%   [X Y Z] = MINIMSEARCH(U) finds the coordinates of the point (X,Y,Z)
%   in the box [-2,2]x[-2,2]x[-2,2] of the space.
%
%   Based on an random search and and Newton-Raphson algorithm
%
    global count
%
% On accede a la variable "count" afin d'effectuer un suivi des
% iterations...  Mais on ne triche pas hein :-)
%
% ======== Principe global de la solution ===============
%
% On effectue une recherche aleatoire de 135 coups
% Ensuite, on effectue 4 iterations de NR avec 16 evaluations
%    ==> 135 + 4 * 16 = 199
% La derniere evaluation sert a estimer la valeur finale du minimum
% Oui : il faut bien la calculer
%
% Il y a evidemment plein d'autres algorithmes possibles.
% Cette solution montre qu'une option SIMPLE existait :-)
% Globalement, c'est une excellente option qui concurrence des trucs
% vachement plus tordus.
%
% Pour etre pedagogique, on optimise aussi la methode de NR 
% en faisant passer le cout d'une iteration de 21 evaluations 
% a 16 evaluations par rapport a l'exemple fourni.
%
% Pour les fans d'algorithmes plus compliques (style Nelder-Mead :-), je
% laisse le soin de les introduire a mon ami Francois Glineur.
% C'est plus complique et pas toujours plus efficace au passage...
%
% =======================================================
%
% -1- Recherche aleatoire avec 135 points
%
    xstart = -2 + 4*rand(135,1);
    ystart = -2 + 4*rand(135,1);
    zstart = -2 + 4*rand(135,1);
    [minf i] = min(u(xstart,ystart,zstart));
    fprintf('\n Count %6d : Obtained value after random search            = %14.7e ',count,minf)
%
% -2- Utilisation du minimum trouve pour lancer NR
%
%     On peut mimiser le nombre d'evaluations en conservant u(x,y,z) et
%     f(x,y,z) et en evitant de les recalculer successivement !!!
%     Le code est un fifrelin moins joli, mais plus efficace :-)
%
    x = xstart(i);
    y = ystart(i);
    z = zstart(i);
    del = 10e-6; tol = 10e-3; delta = tol+1;
    nmax = 4;n = 0;
    while abs(norm(delta)) >= tol && n <= nmax
        n = n+1;
        floc = f(x,y,z);
        uest = uloc;
        dfdx(1,:) = (f(x+del,y,z) - floc)/del; 
        dfdx(2,:) = (f(x,y+del,z) - floc)/del; 
        dfdx(3,:) = (f(x,y,z+del) - floc)/del; 
        delta = -floc / dfdx;
        x = x + delta(1);
        y = y + delta(2);
        z = z + delta(3);
        fprintf('\n Count %6d : NR-iteration %1d : error = %9.2e, objective = %14.7e ',count,n,max(abs(delta)),uest)
    end


    function f = f(x,y,z)
        uloc = u(x,y,z);
        f(1) = (u(x+del,y,z) - uloc)/del; 
        f(2) = (u(x,y+del,z) - uloc)/del; 
        f(3) = (u(x,y,z+del) - uloc)/del; 
    end      


end

%
% Fonctions a minimiser : le test a ete modifie pour ceux qui pensaient 
% qu'on pouvait evaluer davantage de points en incluant un scalaire en x 
% et des vecteurs en y et en z : bande de petits taquins va !
%

function u = difficult(x,y,z)
    global count 
    count = count + max([prod(size(x)) prod(size(y)) prod(size(y))]);
    u = sin(3*y-x.^2+1)+3 *cos(2*y.^4-2*x)+ 4 *sin(3*y-z.^3+1) + 2* cos(2*y.^2-2*z);
end


function u = easy(x,y,z)
    global count 
    count = count + max([prod(size(x)) prod(size(y)) prod(size(y))]);
    u = (x-0.1).^2 + (y - 0.4).^2 + (z - 0.3).^2;
end