Update the ref now we use grand instead of rand.
[scilab.git] / scilab / modules / simulated_annealing / tests / unit_tests / optim_sa.dia.ref
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
3 // Copyright (C) 2010-2011 - DIGITEO - Michael Baudin
4 //
5 // This file must be used under the terms of the CeCILL.
6 // This source file is licensed as described in the file COPYING, which
7 // you should have received as part of this distribution.  The terms
8 // are also available at
9 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 // <-- CLI SHELL MODE -->
11 // <-- ENGLISH IMPOSED -->
12 // Test that we can run with default values for parameters
13 function y=test_func(x)
14   y=x^2
15 endfunction
16 rand("seed",0);
17 x0 = 10;
18 ItExt = 30;
19 ItInt = 50;
20 Log = %f;
21 ItMX = 100;
22 T_init = compute_initial_temp(x0, test_func, 0.8, ItMX);
23 [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history,iter] = optim_sa(x0, test_func, ItExt, ItInt, T_init, Log);
24 assert_checkalmostequal ( x_best , 0.0 , [] , 1.e-1 );
25 assert_checkalmostequal ( f_best , 0.0 , [] , 1.e-1 );
26 assert_checkequal ( size(mean_list) , [1 ItExt] );
27 assert_checkequal ( size(var_list) , [1 ItExt] );
28 assert_checkequal ( size(temp_list) , [1 ItExt] );
29 assert_checkequal ( size(f_history) , ItExt );
30 assert_checkequal ( size(x_history) , ItExt );
31 assert_checkequal ( iter>0 , %t );
32 ///////////////////////////////////////////
33 // Test that we can configure our own neighbour function
34 function f = quad ( x )
35   p = [4 3];
36   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
37 endfunction
38 // We produce a neighbor by adding some noise to each component of a given vector
39 function x_neigh = myneigh_func ( x_current, T , param)
40   nxrow = size(x_current,"r")
41   nxcol = size(x_current,"c")
42   sa_min_delta = -0.1*ones(nxrow,nxcol);
43   sa_max_delta = 0.1*ones(nxrow,nxcol);
44   x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;
45 endfunction
46 rand("seed",0);
47 x0          = [2 2];
48 Proba_start = 0.7;
49 It_Pre      = 100;
50 It_extern   = 50;
51 It_intern   = 50;
52 saparams = init_param();
53 saparams = add_param(saparams,"neigh_func", myneigh_func);
54 T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);
55 Log = %f;
56 [x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams);
57 assert_checkalmostequal ( x_opt , [4 3] ,  1.e-1 );
58 assert_checkalmostequal ( f_opt , 0 ,  [] , 1.e-1 );
59 ///////////////////////////////////////////
60 // Test that an additionnal parameter can be passed to the cost function
61 function f = quadp ( x , p )
62   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
63 endfunction
64 rand("seed",0);
65 x0 = [-1 -1];
66 p = [4 3];
67 T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);
68 [x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 30, 30, T0, %f);
69 assert_checkalmostequal ( x_opt , [4 3] ,  1.e-1 );
70 assert_checkalmostequal ( f_opt , 0 ,  [] , 1.e-1 );
71 ///////////////////////////////////////////
72 // Test with a plot function, which serves also as a stop function.
73 function f = quad ( x )
74   p = [4 3];
75   f = (x(1) - p(1))^2 + (x(2) - p(2))^2
76 endfunction
77 // See that the stop variable becomes true when the function value is near zero.
78 // The threshold is rather loose.
79 function stop = outfunc ( itExt , x_best , f_best , T , saparams )
80   [mythreshold,err] = get_param(saparams,"mythreshold",0);
81   v = format()
82   format("e",10)
83   sxbest = string(x_best)
84   mprintf ( "Iter = #%-4d, \t x_best=[%12s %12s], f_best = %12s, T = %12s\n", itExt , sxbest(1), sxbest(2) , string(f_best) , string(T) )
85   if ( v(1) == 0 ) then
86     format("e",v(2))
87   else
88     format("v",v(2))
89   end
90   stop = ( abs(f_best) < mythreshold )
91 endfunction
92 rand("seed",0);
93 x0 = [-1 -1];
94 saparams = init_param();
95 saparams = add_param(saparams,"output_func", outfunc );
96 saparams = add_param(saparams,"mythreshold", 1.e-1 );
97 T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);
98 // Notice that the number of external iterations is %inf, so 
99 // that the external loop never stops.
100 // This allows to check that the output function really allows to 
101 // stop the algorithm.
102 [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , iter ] = optim_sa(x0, quad , 1e6, 100, T0, %f, saparams);
103 Iter = #1   ,    x_best=[  -1.000D+00   -1.000D+00], f_best =    4.100D+01, T =    1.608D+00
104 Iter = #2   ,    x_best=[   1.078D-01    1.127D-01], f_best =    2.349D+01, T =    1.447D+00
105 Iter = #3   ,    x_best=[   6.984D-01    4.624D-01], f_best =    1.734D+01, T =    1.302D+00
106 Iter = #4   ,    x_best=[   8.860D-01    8.590D-01], f_best =    1.428D+01, T =    1.172D+00
107 Iter = #5   ,    x_best=[   1.012D+00    1.526D+00], f_best =    1.110D+01, T =    1.055D+00
108 Iter = #6   ,    x_best=[   2.233D+00    1.673D+00], f_best =    4.883D+00, T =    9.493D-01
109 Iter = #7   ,    x_best=[   2.884D+00    2.775D+00], f_best =    1.295D+00, T =    8.543D-01
110 Iter = #8   ,    x_best=[   3.356D+00    2.657D+00], f_best =    5.322D-01, T =    7.689D-01
111 Iter = #9   ,    x_best=[   3.356D+00    2.657D+00], f_best =    5.322D-01, T =    6.920D-01
112 Iter = #10  ,    x_best=[   3.356D+00    2.657D+00], f_best =    5.322D-01, T =    6.228D-01
113 Iter = #11  ,    x_best=[   3.985D+00    2.319D+00], f_best =    4.644D-01, T =    5.605D-01
114 Iter = #12  ,    x_best=[   4.048D+00    2.653D+00], f_best =    1.230D-01, T =    5.045D-01
115 Iter = #13  ,    x_best=[   4.048D+00    2.653D+00], f_best =    1.230D-01, T =    4.540D-01
116 Iter = #14  ,    x_best=[   3.942D+00    2.813D+00], f_best =    3.844D-02, T =    4.086D-01
117 assert_checkalmostequal ( x_best , [4 3] ,  1.e-1 );
118 assert_checkalmostequal ( f_best , 0 ,  [] , 1.e-1 );
119 assert_checkequal ( iter > 0 , %t );