// Copyright (C) 2010 - 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

//
// This is a demonstration of Scilab optimization features.
//

///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////
//
// Non linear least squares
//

function dy = myModel(t, y, a, b)
    // The right-hand side of the
    // Ordinary Differential Equation.
    dy(1) = -a*y(2) + y(1) + t^2 + 6*t + b
    dy(2) = b*y(1) - a*y(2) + 4*t + (a+b)*(1-t^2)
 endfunction
 function f = myDifferences(x, t, y_exp)
    // Returns the difference between the
    // simulated differential
    // equation and the experimental data.
    a = x(1)
    b = x(2)
    y0 = y_exp(1,:)
    t0 = 0
    y_calc=ode(y0',t0,t,list(myModel,a,b))
    diffmat = y_calc' - y_exp
    // Maxe a column vector
    f = diffmat(:)
 endfunction
 //
 // 1. Experimental data
 t = [0 1 2 3 4 5 6]';
 y_exp(:,1) = [-1 2 11 26 47 74 107]';
 y_exp(:,2) = [ 1 3 09 19 33 51 73]';
 //
 // 2. Optimize
 a = 0.1;
 b = 0.4;
 x0 = [a;b];
 [fopt,xopt,gopt]=leastsq(list(myDifferences,t,y_exp), x0)

 // Comparing before/after optimization.

function plotmyODE(a,b,t0,y0,t,y_exp,mytitle)
    y_calc=ode(y0',t0,t,list(myModel,a,b));
    subplot(2,1,1);
    plot(t,y_calc(1,:),"r*");
    plot(t,y_exp(:,1),"bo-");
    legend(["Data","ODE"]);
    xtitle(mytitle,"t","Y1");
    subplot(2,1,2);
    plot(t,y_calc(2,:),"r*");
    plot(t,y_exp(:,2),"bo-");
    legend(["Data","ODE"]);
    xtitle(mytitle,"t","Y2");
    // Put transparent marks
    h = gcf();
    for i = 1 : size(h.children,"*")
        for j = 2 : 3
            h.children(i).children(j).children.mark_background=0;
        end
    end
endfunction


// Before optimization
a=x0(1);
b=x0(2);
y0 = y_exp(1,:);
t0 = 0;
mytitle = "Before Optimization";
scf();
plotmyODE(a,b,t0,y0,t,y_exp,mytitle);

// After optimization
a=xopt(1);
b=xopt(2);
y0 = y_exp(1,:);
t0 = 0;
mytitle = "After Optimization";
scf();
plotmyODE(a,b,t0,y0,t,y_exp,mytitle);

///////////////////////////////////////////////////////////
//
// Genetic algorithms
//
function y = rastriginV ( x1 , x2 )
    // Vectorized function for contouring.
    y = x1.^2 + x2.^2-cos(12*x1)-cos(18*x2)
endfunction
function y = rastrigin ( x )
    // Non-vectorized function for optimization.
    y = rastriginV ( x(1) , x(2) )
endfunction
function y= rastriginBinary ( x )
    BinLen = 8
    lb = [-1 -1]';
    ub = [1 1]';
    tmp = convert_to_float(x, BinLen, ub, lb)
    y = rastrigin (tmp)
endfunction
PopSize = 100;
Proba_cross = 0.7;
Proba_mut = 0.1;
NbGen = 10;
Log = %T;
ga_params = init_param();
ga_params = add_param(ga_params,"minbound",[-1 -1]');
ga_params = add_param(ga_params,"maxbound",[1 1]');
ga_params = add_param(ga_params,"dimension",2);
ga_params = add_param(ga_params,"binary_length",8);
ga_params = add_param(ga_params,"crossover_func",crossover_ga_binary);
ga_params = add_param(ga_params,"mutation_func",mutation_ga_binary);
ga_params = add_param(ga_params,"codage_func",coding_ga_binary);
ga_params = add_param(ga_params,"multi_cross",%T);
[pop_opt,fobj_pop_opt, pop_init , fobj_init] = optim_ga(rastriginBinary, PopSize, ..
  NbGen, Proba_mut, Proba_cross, Log,ga_params);
// 3. Compute contouring values.
x1 = linspace(-1,1,100); 
x2 = linspace(-1,1,100); 
[XX1,XX2]=meshgrid(x1,x2);
Z = rastriginV ( XX1 , XX2  );

// 4. Plot the results.

scf();
plot3d(x1,x2,Z');
xtitle('Ratrigin Function - 3D Plot','x1','x2');

scf();
drawlater;
// Plot initial population
subplot(2,1,1);
xset("fpf"," ");
contour ( x1 , x2 , Z' , 5 );
xtitle("Initial Population","X1","X2");
for i=1:length(pop_init)
  plot(pop_init(i)(1),pop_init(i)(2),"r.");
end
// Plot optimum population
subplot(2,1,2);
xset("fpf"," ");
contour ( x1 , x2 , Z' , 5 );
xtitle("Optimum Population","X1","X2");
for i=1:length(pop_opt)
  plot(pop_opt(i)(1),pop_opt(i)(2),"g.");
end
drawnow;

///////////////////////////////////////////////////////////