optim_sa: added output function 64/364/3
Michael Baudin [Mon, 3 May 2010 14:37:41 +0000 (16:37 +0200)]
Change-Id: Ia9f42b36f4d004b890e487ea6dc220e63deb1feb

scilab/modules/simulated_annealing/demos/SAIsing2ddemo.sce
scilab/modules/simulated_annealing/demos/SAdemo.sce
scilab/modules/simulated_annealing/help/en_US/neigh_func_default.xml
scilab/modules/simulated_annealing/help/en_US/optim_sa.xml
scilab/modules/simulated_annealing/macros/optim_sa.sci
scilab/modules/simulated_annealing/tests/unit_tests/optim_sa.dia.ref
scilab/modules/simulated_annealing/tests/unit_tests/optim_sa.tst

index 9b87dbe..366eef8 100755 (executable)
@@ -3,8 +3,8 @@
 //////////////////////////////////////////////////////
 
 my_handle = scf(100001);
-clf(my_handle,'reset');
-demo_viewCode('SAIsing2ddemo.sce');
+clf(my_handle,"reset");
+demo_viewCode("SAIsing2ddemo.sce");
 
 lines(0);
 old_funcprot = funcprot();
@@ -12,13 +12,13 @@ funcprot(0);
 
 // Loading the test problem
 
-path = get_absolute_file_path('SAIsing2ddemo.sce');
+path = get_absolute_file_path("SAIsing2ddemo.sce");
 
-getd(path + '/Ising');
+getd(path + "/Ising");
 
 // Loading the neighborhood function for the ising problem
 
-getd(path + '/.');
+getd(path + "/.");
 
 Proba_start = 0.8;
 It_intern = 1000;
@@ -31,13 +31,13 @@ Ising_Dim   = 10;
 Ising_Proba = 0.3;
 J       = 1.1;
 H       = 0.7;
-Magnet  = '%T';
-Connect = '%T';
+Magnet  = "%T";
+Connect = "%T";
 
 // For the minimization case, everything must be at -1 or 1 in the optimal solution
-deff('y=f(x)','y = ising2d(x,'+string(J)+','+string(H)+','+Magnet+','+Connect+')');
+deff("y=f(x)","y = ising2d(x,"+string(J)+","+string(H)+","+Magnet+","+Connect+")");
 // For the maximization case, we must have a checker board solution (+1 -1 +1 -1 ....)
-//deff('y=f(x)','y = - ising2d(x,'+string(J)+','+string(H)+','+Magnet+','+Connect+')');
+//deff("y=f(x)","y = - ising2d(x,"+string(J)+","+string(H)+","+Magnet+","+Connect+")");
 
 x0 = init_ising2d(Ising_Dim, Ising_Proba);
 
@@ -45,24 +45,24 @@ x0 = init_ising2d(Ising_Dim, Ising_Proba);
 // Simulated Annealing //
 /////////////////////////
 
-printf('SA: geometrical decrease temperature law\n');
+printf("SA: geometrical decrease temperature law\n");
 
 sa_params = init_param();
-sa_params = add_param(sa_params,'dimension',10);
-sa_params = add_param(sa_params,'proba',0.05);
-sa_params = add_param(sa_params,'neigh_func', neigh_func_ising2d); // Required because this operator is specific to the ising2d problem
-sa_params = add_param(sa_params,'accept_func', accept_func_default); // Optional
-sa_params = add_param(sa_params,'temp_law', temp_law_default); // Optional
-sa_params = add_param(sa_params,'alpha',alpha); // For the temperature decreasing law
+sa_params = add_param(sa_params,"dimension",10);
+sa_params = add_param(sa_params,"proba",0.05);
+sa_params = add_param(sa_params,"neigh_func", neigh_func_ising2d); // Required because this operator is specific to the ising2d problem
+sa_params = add_param(sa_params,"accept_func", accept_func_default); // Optional
+sa_params = add_param(sa_params,"temp_law", temp_law_default); // Optional
+sa_params = add_param(sa_params,"alpha",alpha); // For the temperature decreasing law
 
 T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
 
-printf('Initial temperature T0 = %f\n', T0);
+printf("Initial temperature T0 = %f\n", T0);
 
 [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log, sa_params);
 
-printf('optimal solution:\n'); disp(x_opt);
-printf('value of the objective function = %f\n', f_opt);
+printf("optimal solution:\n"); disp(x_opt);
+printf("value of the objective function = %f\n", f_opt);
 
 plot_ising2d(x_opt);
 
index 42fa9bb..367e4c0 100755 (executable)
@@ -3,8 +3,8 @@
 //////////////////////////////////////////////////////
 
 my_handle = scf(100001);
-clf(my_handle,'reset');
-demo_viewCode('SAdemo.sce');
+clf(my_handle,"reset");
+demo_viewCode("SAdemo.sce");
 
 lines(0);
 old_funcprot = funcprot();
@@ -27,7 +27,7 @@ function y = rastrigin(x)
   y = x(1)^2+x(2)^2-cos(12*x(1))-cos(18*x(2));
 endfunction
 
-func = 'rastrigin';
+func = "rastrigin";
 
 Proba_start = 0.8;
 It_intern = 1000;
@@ -43,47 +43,47 @@ DoHuang = %F;
 
 //////////////////////////////////////////
 
-Min = eval('min_bd_'+func+'()');
-Max = eval('max_bd_'+func+'()');
+Min = eval("min_bd_"+func+"()");
+Max = eval("max_bd_"+func+"()");
 x0  = (Max - Min).*rand(size(Min,1),size(Min,2)) + Min;
 
-deff('y=f(x)','y='+func+'(x)');
+deff("y=f(x)","y="+func+"(x)");
 
 /////////////////////////
 // Simulated Annealing //
 /////////////////////////
 
 if DoSA then
-  printf('SA: geometrical decrease temperature law\n');
+  printf("SA: geometrical decrease temperature law\n");
 
   sa_params = init_param();
-  sa_params = add_param(sa_params,'min_delta',-0.1*(Max-Min));
-  sa_params = add_param(sa_params,'max_delta', 0.1*(Max-Min));
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_default); // Optional
-  sa_params = add_param(sa_params,'accept_func', accept_func_default); // Optional
-  sa_params = add_param(sa_params,'temp_law', temp_law_default); // Optional
-  sa_params = add_param(sa_params,'min_bound',Min);
-  sa_params = add_param(sa_params,'max_bound',Max);
-  sa_params = add_param(sa_params,'alpha',alpha); // for the decreasing temperature law
+  sa_params = add_param(sa_params,"min_delta",-0.1*(Max-Min));
+  sa_params = add_param(sa_params,"max_delta", 0.1*(Max-Min));
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_default); // Optional
+  sa_params = add_param(sa_params,"accept_func", accept_func_default); // Optional
+  sa_params = add_param(sa_params,"temp_law", temp_law_default); // Optional
+  sa_params = add_param(sa_params,"min_bound",Min);
+  sa_params = add_param(sa_params,"max_bound",Max);
+  sa_params = add_param(sa_params,"alpha",alpha); // for the decreasing temperature law
 
   T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
-  printf('Initial temperature T0 = %f\n', T0);
+  printf("Initial temperature T0 = %f\n", T0);
 
   [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log = %T, sa_params);
 
-  printf('optimal solution:\n'); disp(x_opt);
-  printf('value of the objective function = %f\n', f_opt);
+  printf("optimal solution:\n"); disp(x_opt);
+  printf("value of the objective function = %f\n", f_opt);
 
   drawlater;
   subplot(2,1,1);
-  xtitle('Geometrical annealing','Iteration','Mean / Variance');
+  xtitle("Geometrical annealing","Iteration","Mean / Variance");
   t = 1:length(sa_mean_list);
-  plot(t,sa_mean_list,'r',t,sa_var_list,'g');
-  legend(['Mean','Variance']);
+  plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+  legend(["Mean","Variance"]);
   subplot(2,1,2);
-  xtitle('Temperature evolution','Iteration','Temperature');
+  xtitle("Temperature evolution","Iteration","Temperature");
   for i=1:length(t)-1
-    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],'k-');
+    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],"k-");
   end
   drawnow;
 end
@@ -93,40 +93,40 @@ end
 /////////
 
 if DoFSA then
-  printf('SA: the FSA algorithm\n');
+  printf("SA: the FSA algorithm\n");
 
   sa_params = init_param();
-  sa_params = add_param(sa_params,'min_delta',-0.1*(Max-Min));
-  sa_params = add_param(sa_params,'max_delta', 0.1*(Max-Min));
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_default); // Required for compute_initial_temp
-  sa_params = add_param(sa_params,'accept_func', accept_func_default); // Optional
-  sa_params = add_param(sa_params,'temp_law', temp_law_fsa); // Required to transform SA into FSA
-  sa_params = add_param(sa_params,'min_bound',Min);
-  sa_params = add_param(sa_params,'max_bound',Max);
+  sa_params = add_param(sa_params,"min_delta",-0.1*(Max-Min));
+  sa_params = add_param(sa_params,"max_delta", 0.1*(Max-Min));
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_default); // Required for compute_initial_temp
+  sa_params = add_param(sa_params,"accept_func", accept_func_default); // Optional
+  sa_params = add_param(sa_params,"temp_law", temp_law_fsa); // Required to transform SA into FSA
+  sa_params = add_param(sa_params,"min_bound",Min);
+  sa_params = add_param(sa_params,"max_bound",Max);
 
   T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
   
-  sa_params = remove_param(sa_params,'neigh_func');
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_fsa); // Required to transform SA into FSA
+  sa_params = remove_param(sa_params,"neigh_func");
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_fsa); // Required to transform SA into FSA
   
-  printf('Initial temperature T0 = %f\n', T0);
+  printf("Initial temperature T0 = %f\n", T0);
 
   [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log = %T, sa_params);
 
-  printf('optimal solution:\n'); disp(x_opt);
-  printf('value of the objective function = %f\n', f_opt);
+  printf("optimal solution:\n"); disp(x_opt);
+  printf("value of the objective function = %f\n", f_opt);
 
   scf();
   drawlater;
   subplot(2,1,1);
-  xtitle('FSA','Iteration','Mean / Variance');
+  xtitle("FSA","Iteration","Mean / Variance");
   t = 1:length(sa_mean_list);
-  plot(t,sa_mean_list,'r',t,sa_var_list,'g');
-  legend(['Mean','Variance']);
+  plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+  legend(["Mean","Variance"]);
   subplot(2,1,2);
-  xtitle('Temperature evolution','Iteration','Temperature');
+  xtitle("Temperature evolution","Iteration","Temperature");
   for i=1:length(t)-1
-    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],'k-');
+    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],"k-");
   end
   drawnow;
 end
@@ -136,44 +136,44 @@ end
 //////////
 
 if DoVFSA then
-  printf('SA: the VFSA algorithm\n');
+  printf("SA: the VFSA algorithm\n");
 
   sa_params = init_param();
-  sa_params = add_param(sa_params,'min_delta',-0.1*(Max-Min));
-  sa_params = add_param(sa_params,'max_delta', 0.1*(Max-Min));
-  sa_params = add_param(sa_params,'accept_func', accept_func_vfsa); // Required to transform SA into FSA
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_default); // Required for compute_initial_temp
-  sa_params = add_param(sa_params,'temp_law', temp_law_vfsa); // Required to transform SA into FSA
-  sa_params = add_param(sa_params,'type_accept', 'vfsa'); // Required to compute correctly the starting temperature for VFSA
-  sa_params = add_param(sa_params,'dimension', length(x0)); // Required to compute correctly the starting temperature for VFSA
-  sa_params = add_param(sa_params,'min_bound',Min);
-  sa_params = add_param(sa_params,'max_bound',Max);
+  sa_params = add_param(sa_params,"min_delta",-0.1*(Max-Min));
+  sa_params = add_param(sa_params,"max_delta", 0.1*(Max-Min));
+  sa_params = add_param(sa_params,"accept_func", accept_func_vfsa); // Required to transform SA into FSA
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_default); // Required for compute_initial_temp
+  sa_params = add_param(sa_params,"temp_law", temp_law_vfsa); // Required to transform SA into FSA
+  sa_params = add_param(sa_params,"type_accept", "vfsa"); // Required to compute correctly the starting temperature for VFSA
+  sa_params = add_param(sa_params,"dimension", length(x0)); // Required to compute correctly the starting temperature for VFSA
+  sa_params = add_param(sa_params,"min_bound",Min);
+  sa_params = add_param(sa_params,"max_bound",Max);
   
   T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
 
-  sa_params = remove_param(sa_params,'neigh_func');
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_vfsa); // Required to transform SA into VFSA
-  sa_params = remove_param(sa_params,'type_accept');
-  sa_params = add_param(sa_params,'type_accept', 'sa'); // We go back to the classical method for computing the starting temperature
+  sa_params = remove_param(sa_params,"neigh_func");
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_vfsa); // Required to transform SA into VFSA
+  sa_params = remove_param(sa_params,"type_accept");
+  sa_params = add_param(sa_params,"type_accept", "sa"); // We go back to the classical method for computing the starting temperature
 
-  printf('Initial temperature T0 = %f\n', T0);
+  printf("Initial temperature T0 = %f\n", T0);
 
   [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log = %T, sa_params);
 
-  printf('optimal solution:\n'); disp(x_opt);
-  printf('value of the objective function = %f\n', f_opt);
+  printf("optimal solution:\n"); disp(x_opt);
+  printf("value of the objective function = %f\n", f_opt);
 
   scf();
   drawlater;
   subplot(2,1,1);
-  xtitle('VFSA','Iteration','Mean / Variance');
+  xtitle("VFSA","Iteration","Mean / Variance");
   t = 1:length(sa_mean_list);
-  plot(t,sa_mean_list,'r',t,sa_var_list,'g');
-  legend(['Mean','Variance']);
+  plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+  legend(["Mean","Variance"]);
   subplot(2,1,2);
-  xtitle('Temperature evolution','Iteration','Temperature');
+  xtitle("Temperature evolution","Iteration","Temperature");
   for i=1:length(t)-1
-    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],'k-');
+    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],"k-");
   end
   drawnow;
 end
@@ -183,40 +183,40 @@ end
 /////////
 
 if DoCSA then
-  printf('SA: the CSA algorithm\n');
+  printf("SA: the CSA algorithm\n");
 
   sa_params = init_param();
-  sa_params = add_param(sa_params,'min_delta',-0.1*(Max-Min));
-  sa_params = add_param(sa_params,'max_delta', 0.1*(Max-Min));
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_default); // Required for compute_initial_temp
-  sa_params = add_param(sa_params,'accept_func', accept_func_default); // Optional
-  sa_params = add_param(sa_params,'temp_law', temp_law_csa); // Required to transform SA into CSA
-  sa_params = add_param(sa_params,'min_bound',Min);
-  sa_params = add_param(sa_params,'max_bound',Max);
+  sa_params = add_param(sa_params,"min_delta",-0.1*(Max-Min));
+  sa_params = add_param(sa_params,"max_delta", 0.1*(Max-Min));
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_default); // Required for compute_initial_temp
+  sa_params = add_param(sa_params,"accept_func", accept_func_default); // Optional
+  sa_params = add_param(sa_params,"temp_law", temp_law_csa); // Required to transform SA into CSA
+  sa_params = add_param(sa_params,"min_bound",Min);
+  sa_params = add_param(sa_params,"max_bound",Max);
 
   T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
   
-  sa_params = remove_param(sa_params,'neigh_func');
-  sa_params = add_param(sa_params,'neigh_func', neigh_func_csa); // Required to transform SA into CSA
+  sa_params = remove_param(sa_params,"neigh_func");
+  sa_params = add_param(sa_params,"neigh_func", neigh_func_csa); // Required to transform SA into CSA
 
-  printf('Initial temperature T0 = %f\n', T0);
+  printf("Initial temperature T0 = %f\n", T0);
 
   [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log = %T, sa_params);
 
-  printf('optimal solution:\n'); disp(x_opt);
-  printf('value of the objective function = %f\n', f_opt);
+  printf("optimal solution:\n"); disp(x_opt);
+  printf("value of the objective function = %f\n", f_opt);
 
   scf();
   drawlater;
   subplot(2,1,1);
-  xtitle('Classical simulated annealing','Iteration','Mean / Variance');
+  xtitle("Classical simulated annealing","Iteration","Mean / Variance");
   t = 1:length(sa_mean_list);
-  plot(t,sa_mean_list,'r',t,sa_var_list,'g');
-  legend(['Mean','Variance']);
+  plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+  legend(["Mean","Variance"]);
   subplot(2,1,2);
-  xtitle('Temperature evolution','Iteration','Temperature');
+  xtitle("Temperature evolution","Iteration","Temperature");
   for i=1:length(t)-1
-    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],'k-');
+    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],"k-");
   end
   drawnow;
 end
@@ -226,34 +226,34 @@ end
 ///////////
 
 if DoHuang then
-  printf('SA: the Huang annealing\n');
+  printf("SA: the Huang annealing\n");
 
   sa_params = init_param();
-  sa_params = add_param(sa_params,'min_delta',-0.1*(Max-Min));
-  sa_params = add_param(sa_params,'max_delta', 0.1*(Max-Min));
-  sa_params = add_param(sa_params,'temp_law', temp_law_huang);
-  sa_params = add_param(sa_params,'min_bound',Min);
-  sa_params = add_param(sa_params,'max_bound',Max);
+  sa_params = add_param(sa_params,"min_delta",-0.1*(Max-Min));
+  sa_params = add_param(sa_params,"max_delta", 0.1*(Max-Min));
+  sa_params = add_param(sa_params,"temp_law", temp_law_huang);
+  sa_params = add_param(sa_params,"min_bound",Min);
+  sa_params = add_param(sa_params,"max_bound",Max);
 
   T0 = compute_initial_temp(x0, f, Proba_start, It_Pre, sa_params);
-  printf('Initial temperature T0 = %f\n', T0);
+  printf("Initial temperature T0 = %f\n", T0);
 
   [x_opt, f_opt, sa_mean_list, sa_var_list, temp_list] = optim_sa(x0, f, It_extern, It_intern, T0, Log = %T, sa_params);
 
-  printf('optimal solution:\n'); disp(x_opt);
-  printf('value of the objective function = %f\n', f_opt);
+  printf("optimal solution:\n"); disp(x_opt);
+  printf("value of the objective function = %f\n", f_opt);
 
   scf();
   drawlater;
   subplot(2,1,1);
-  xtitle('Huang annealing','Iteration','Mean / Variance');
+  xtitle("Huang annealing","Iteration","Mean / Variance");
   t = 1:length(sa_mean_list);
-  plot(t,sa_mean_list,'r',t,sa_var_list,'g');
-  legend(['Mean','Variance']);
+  plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+  legend(["Mean","Variance"]);
   subplot(2,1,2);
-  xtitle('Temperature evolution','Iteration','Temperature');
+  xtitle("Temperature evolution","Iteration","Temperature");
   for i=1:length(t)-1
-    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],'k-');
+    plot([t(i) t(i+1)], [temp_list(i) temp_list(i)],"k-");
   end
   drawnow;
 end
index e753463..5473c07 100644 (file)
   <refsynopsisdiv>
     <title>Calling Sequence</title>
 
-    <synopsis>x_neigh = neigh_func_default(x_current,T,param)</synopsis>
+    <synopsis>
+      x_neigh = neigh_func_default(x_current,T)
+      x_neigh = neigh_func_default(x_current,T,param)
+    </synopsis>
   </refsynopsisdiv>
 
   <refsection>
@@ -48,7 +51,8 @@
         <term>T</term>
 
         <listitem>
-          <para>the current temperature</para>
+          <para>the current temperature. This parameter is ignored but is 
+          there to make all the neighbour function consistent.</para>
         </listitem>
       </varlistentry>
 
@@ -98,14 +102,14 @@ Proba_start = 0.7;
 It_Pre      = 100;
 It_extern   = 100;
 It_intern   = 1000;
-x_test = neigh_func_default(x0);
+x_test = neigh_func_default(x0,%nan);
 
 saparams = init_param();
 saparams = add_param(comp_t_params,'neigh_func', neigh_func_default);
 
 T0 = compute_initial_temp(x0, rastrigin, Proba_start, It_Pre, saparams);
-
-[x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log = %T);
+Log = %T;
+[x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log,saparams);
 
 printf('optimal solution:\n'); disp(x_opt);
 printf('value of the objective function = %f\n', f_opt);
index 239d837..cf3b7ad 100644 (file)
@@ -2,6 +2,7 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
+ * Copyright (C) 2010 - 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
@@ -17,7 +18,7 @@
           xmlns:ns4="http://www.w3.org/1999/xhtml"
           xmlns:mml="http://www.w3.org/1998/Math/MathML"
           xmlns:db="http://docbook.org/ns/docbook">
-  
+
   <refnamediv>
     <refname>optim_sa</refname>
 
   <refsynopsisdiv>
     <title>Calling Sequence</title>
 
-    <synopsis>[x_best,f_best,mean_list,var_list,f_history,temp_list,x_history] = optim_sa(x0,f,ItExt,ItInt,T0,Log,temp_law,param_temp_law,neigh_func,param_neigh_func)</synopsis>
+    <synopsis>
+      x_best = optim_sa(x0,f,ItExt,ItInt,T0,Log,temp_law,param_temp_law,neigh_func,param_neigh_func)
+      [x_best,f_best] = optim_sa(..)
+      [x_best,f_best,mean_list] = optim_sa(..)
+      [x_best,f_best,mean_list,var_list] = optim_sa(..)
+      [x_best,f_best,mean_list,var_list,f_history] = optim_sa(..)
+      [x_best,f_best,mean_list,var_list,f_history,temp_list] = optim_sa(..)
+      [x_best,f_best,mean_list,var_list,f_history,temp_list,x_history] = optim_sa(..)
+      [x_best,f_best,mean_list,var_list,f_history,temp_list,x_history,iter] = optim_sa(..)
+    </synopsis>
   </refsynopsisdiv>
 
   <refsection>
         <term>f</term>
 
         <listitem>
-          <para>the objective function to be optimized (the prototype if
-          f(x))</para>
+          <para>
+            the objective function to be optimized (the prototype if
+            f(x))
+          </para>
         </listitem>
       </varlistentry>
 
         <term>T0</term>
 
         <listitem>
-          <para>the initial temperature (see compute_initial_temp to compute
-          easily this temperature)</para>
+          <para>
+            the initial temperature (see compute_initial_temp to compute
+            easily this temperature)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>Log</term>
 
         <listitem>
-          <para>if %T, some information will be displayed during the run of
-          the simulated annealing</para>
+          <para>
+            if %T, some information will be displayed during the run of
+            the simulated annealing
+          </para>
         </listitem>
       </varlistentry>
 
         <term>temp_law</term>
 
         <listitem>
-          <para>the temperature decrease law (see temp_law_default for an
-          example of such a function)</para>
+          <para>
+            the temperature decrease law (see temp_law_default for an
+            example of such a function)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>param_temp_law</term>
 
         <listitem>
-          <para>a structure (of any kind - it depends on the temperature law
-          used) which is transmitted as a parameter to temp_law</para>
+          <para>
+            a structure (of any kind - it depends on the temperature law
+            used) which is transmitted as a parameter to temp_law
+          </para>
         </listitem>
       </varlistentry>
 
         <term>neigh_func</term>
 
         <listitem>
-          <para>a function which computes a neighbor of a given point (see
-          neigh_func_default for an example of such a function)</para>
+          <para>
+            a function which computes a neighbor of a given point (see
+            neigh_func_default for an example of such a function)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>param_neigh_func</term>
 
         <listitem>
-          <para>a structure (of any kind like vector, list, it depends on the
-          neighborhood function used) which is transmitted as a parameter to
-          neigh_func</para>
+          <para>
+            a structure (of any kind like vector, list, it depends on the
+            neighborhood function used) which is transmitted as a parameter to
+            neigh_func
+          </para>
         </listitem>
       </varlistentry>
 
         <term>mean_list</term>
 
         <listitem>
-          <para>the mean of the objective function value for each temperature
-          stage. A vector of float (optional)</para>
+          <para>
+            the mean of the objective function value for each temperature
+            stage. A vector of float (optional)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>var_list</term>
 
         <listitem>
-          <para>the variance of the objective function values for each
-          temperature stage. A vector of float (optional)</para>
+          <para>
+            the variance of the objective function values for each
+            temperature stage. A vector of float (optional)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>f_history</term>
 
         <listitem>
-          <para>the computed objective function values for each iteration.
-          Each input of the list corresponds to a temperature stage. Each
-          input of the list is a vector of float which gathers all the
-          objective function values computed during the corresponding
-          temperature stage - (optional)</para>
+          <para>
+            the computed objective function values for each iteration.
+            Each input of the list corresponds to a temperature stage. Each
+            input of the list is a vector of float which gathers all the
+            objective function values computed during the corresponding
+            temperature stage - (optional)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>temp_list</term>
 
         <listitem>
-          <para>the list of temperature computed for each temperature stage. A
-          vector of float (optional)</para>
+          <para>
+            the list of temperature computed for each temperature stage. A
+            vector of float (optional)
+          </para>
         </listitem>
       </varlistentry>
 
         <term>x_history</term>
 
         <listitem>
-          <para>the parameter values computed for each iteration. Each input
-          of the list corresponds to a temperature stage. Each input of the
-          list is a vector of input variables which corresponds to all the
-          variables computed during the corresponding temperature stage -
-          (optional - can slow down a lot the execution of optim_sa)</para>
+          <para>
+            the parameter values computed for each iteration. Each input
+            of the list corresponds to a temperature stage. Each input of the
+            list is a vector of input variables which corresponds to all the
+            variables computed during the corresponding temperature stage -
+            (optional - can slow down a lot the execution of optim_sa)
+          </para>
         </listitem>
       </varlistentry>
+
+      <varlistentry>
+        <term>iter</term>
+
+        <listitem>
+          <para>
+            a double, the actual number of external iterations in the
+            algorithm (optional).
+          </para>
+        </listitem>
+      </varlistentry>
+
     </variablelist>
   </refsection>
 
   <refsection>
     <title>Description</title>
 
-        <para>A Simulated Annealing optimization method.</para>
+    <para>A Simulated Annealing optimization method.</para>
+    <para>
+      Simulated annealing (SA) is a generic probabilistic meta-algorithm for the global optimization
+      problem, namely locating a good approximation to the global optimum of a given function in a
+      large search space. It is often used when the search space is discrete (e.g., all tours that
+      visit a given set of cities).
+    </para>
+
+    <para>
+      The current solver can find the solution of an optimization problem without constraints or 
+      with bound constraints. The bound constraints can be customized with the neighbour
+      function. This algorithm does not use the derivatives of the objective function.
+    </para>
+
+    <para>
+      The solver is made of Scilab macros, which enables a high-level programming model for
+      this optimization solver. The SA macros are based on the <literal>parameters</literal>
+      Scilab module for the management of the (many) optional parameters.
+    </para>
+
+    <para>
+      To use the SA algorithm, one should perform the following steps :
+      <itemizedlist>
+        <listitem>
+          configure the parameters with calls to <literal>init_param</literal> 
+          and <literal>add_param</literal> especially the neighbor function, the 
+          acceptance function, the temperature law,
+        </listitem>
+        <listitem>
+          compute an initial temperature with a call to <literal>compute_initial_temp</literal>,
+        </listitem>
+        <listitem>
+          find an optimum by using the <literal>optim_sa</literal> solver.
+        </listitem>
+      </itemizedlist>
+    </para>
+
+    <para>
+      The algorithm is based on an iterative update of two points :
+      <itemizedlist>
+        <listitem>
+          the current point is updated by taking into account the neighbour and the acceptance
+          functions,
+        </listitem>
+        <listitem>
+          the best point is the point which achieved the minimum of the objective function over the
+          iterations.
+        </listitem>
+      </itemizedlist>
+      While the current point is used internally to explore the domain, only the best point is returned
+      by the function.
+      The algorithm is based on an external loop and an internal loop. In the external loop,
+      the temperature is updated according to the temperature function. In the internal loop, the 
+      point is updated according to the neighbour function. A new point is accepted depending 
+      on its associated function value or the value of the acceptance function, which value 
+      depends on the current temperature and a uniform random number.
+    </para>
+
+    <para>
+      The acceptance of the new point depends on the output values produced
+      by the <literal>rand</literal> function. This implies that two consecutive 
+      calls to the <literal>optim_sa</literal> will not produce the same result.
+      In order to always get exactly the same results, please initialize the random number 
+      generator with a valid seed.
+    </para>
+
+    <para>
+      See the Demonstrations, in the "Optimization" section and "Simulated Annealing" subsection
+      for more examples.
+    </para>
+
   </refsection>
 
   <refsection>
-    <title>Examples</title>
+    <title>The objective function</title>
+    <para>
+      The objective function is expected to have the following header.
+    </para>
 
     <programlisting role="example">
-  <![CDATA[ 
+      <![CDATA[ 
+function y = f ( x )
+ ]]>
+</programlisting>
+
+    <para>
+      In the case where the objective function needs additionnal parameters,
+      the objective function can be defined as a list, where the first 
+      argument is the cost function, and the second argument is the 
+      additionnal parameter. See below for an example.
+    </para>
+
+  </refsection>
+
+<refsection>
+    <title>A simple example</title>
+    <para>
+      In the following example, we search the minimum of the
+      Rastriging function. This function has many local minimas, but only
+      one single global minimum located at x = (0,0), where the function value is
+      f(x) = -2. We use the simulated annealing algorithm with default settings
+      and the default neighbour function neigh_func_default.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
 function y = rastrigin(x)
   y = x(1)^2+x(2)^2-cos(12*x(1))-cos(18*x(2));
 endfunction
@@ -213,47 +357,211 @@ It_extern   = 100;
 It_intern   = 1000;
 x_test = neigh_func_default(x0);
 
-saparams = init_param();
-saparams = add_param(comp_t_params,'neigh_func', neigh_func_default);
-
-T0 = compute_initial_temp(x0, rastrigin, Proba_start, It_Pre, saparams);
+T0 = compute_initial_temp(x0, rastrigin, Proba_start, It_Pre);
 
-[x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log = %T);
+Log = %T;
+[x_opt, f_opt, sa_mean_list, sa_var_list] = optim_sa(x0, rastrigin, It_extern, It_intern, T0, Log);
 
-printf('optimal solution:\n'); disp(x_opt);
-printf('value of the objective function = %f\n', f_opt);
+printf("optimal solution:\n"); disp(x_opt);
+printf("value of the objective function = %f\n", f_opt);
 
 t = 1:length(sa_mean_list);
-plot(t,sa_mean_list,'r',t,sa_var_list,'g');
- ]]></programlisting>
+plot(t,sa_mean_list,"r",t,sa_var_list,"g");
+ ]]>
+    </programlisting>
   </refsection>
 
   <refsection>
-    <title>See Also</title>
+    <title>Configuring a neighbour function</title>
+    <para>
+      In the following example, we customize the
+      neighbourhood function. In order to pass this function to the
+      <literal>optim_sa</literal> function, we setup a parameter where the
+      <literal>"neigh_func"</literal> key is associated with our particular neighbour function.
+      The neighbour function can be customized at will, provided that the
+      header of the function is the same. The particular implementation shown
+      below is the same, in spirit, as the <literal>neigh_func_default</literal>
+      function.
+    </para>
 
-    <simplelist type="inline">
-      <member><link linkend="compute_initial_temp"> compute_initial_temp
-      </link></member>
+    <programlisting role="example">
+      <![CDATA[ 
+function f = quad ( x )
+  p = [4 3];
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+endfunction
+
+// We produce a neighbor by adding some noise to each component of a given vector
+function x_neigh = myneigh_func ( x_current, T , param)
+  nxrow = size(x_current,"r")
+  nxcol = size(x_current,"c")
+  sa_min_delta = -0.1*ones(nxrow,nxcol);
+  sa_max_delta = 0.1*ones(nxrow,nxcol);
+  x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;
+endfunction
+
+x0          = [2 2];
+Proba_start = 0.7;
+It_Pre      = 100;
+It_extern   = 50;
+It_intern   = 100;
+
+saparams = init_param();
+saparams = add_param(saparams,"neigh_func", myneigh_func);
+// or: saparams = add_param(saparams,"neigh_func", neigh_func_default);
+// or: saparams = add_param(saparams,"neigh_func", neigh_func_csa);
+// or: saparams = add_param(saparams,"neigh_func", neigh_func_fsa);
+// or: saparams = add_param(saparams,"neigh_func", neigh_func_vfsa);
+
+T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);
+Log = %f;
+// This should produce x_opt = [4 3]
+[x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams);
+ ]]>
+    </programlisting>
+  </refsection>
+
+
+  <refsection>
+    <title>Passing extra parameters</title>
+    <para>
+      In the following example, we use an objective function which requires
+      an extra parameter <literal>p</literal>. This parameter is the second 
+      input argument of the <literal>quadp</literal> function. In order to 
+      pass this parameter to the objective function, we define the objective 
+      function as <literal>list(quadp,p)</literal>. In this case,
+      the solver makes so that the calling sequence includes a second argument.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+  function f = quadp ( x , p )
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+  endfunction
+
+  x0 = [-1 -1];
+  p = [4 3];
+  T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);
+  [x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 10, 1000, T0, %f);
+ ]]>
+</programlisting>
+</refsection>
+
+  <refsection>
+    <title>Configuring an output function</title>
+    <para>
+      In the following example, we define an output function, which also 
+      provide a stopping rule. We define the function <literal>outfun</literal>
+      which takes as input arguments the data of the algorithm at the current 
+      iteration and returns the boolean <literal>stop</literal>. This function
+      prints a message into the console to inform the user about the 
+      current state of the algorithm. It also computes the boolean <literal>stop</literal>,
+      depending on the value of the function.
+      The stop variable becomes true when the function value is near zero. In order to let <literal>optim_sa</literal>
+      know about our output function, we configure the <literal>"output_func"</literal>
+      key to our <literal>outfun</literal> function and call the solver.
+      Notice that the number of external iterations is <literal>%inf</literal>, so
+      that the external loop never stops.
+      This allows to check that the output function really allows to
+      stop the algorithm.
+    </para>
 
-      <member><link linkend="neigh_func_default"> neigh_func_default
-      </link></member>
+    <programlisting role="example">
+      <![CDATA[ 
+function f = quad ( x )
+  p = [4 3];
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+endfunction
 
-      <member><link linkend="temp_law_default"> temp_law_default
-      </link></member>
+function stop = outfunc ( itExt , x_best , f_best , T , saparams )
+  [mythreshold,err] = get_param(saparams,"mythreshold",0);
+  mprintf ( "Iter = #%d, \t x_best=[%f %f], \t f_best = %e, \t T = %e\n", itExt , x_best(1),x_best(2) , f_best , T )
+  stop = ( abs(f_best) < mythreshold )
+endfunction
+
+x0 = [-1 -1];
+saparams = init_param();
+saparams = add_param(saparams,"output_func", outfunc );
+saparams = add_param(saparams,"mythreshold", 1.e-2 );
+
+rand("seed",0);
+
+T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);
+[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , It ] = optim_sa(x0, quad , %inf, 100, T0, %f, saparams);
+ ]]>
+</programlisting>
+
+    <para>
+      The previous script produces the following output. Notice that the actual 
+      output of the algorithm depends on the state of the random number generator <literal>rand</literal>:
+      if we had not initialize the seed of the uniform random number generator, 
+      we would have produced a different result.
+    </para>
+    <programlisting role="example">
+      Iter = #1,        x_best=[-1.000000 -1.000000],   f_best = 4.100000e+001,         T = 1.453537e+000
+      Iter = #2,        x_best=[-0.408041 -0.318262],   f_best = 3.044169e+001,         T = 1.308183e+000
+      Iter = #3,        x_best=[-0.231406 -0.481078],   f_best = 3.002270e+001,         T = 1.177365e+000
+      Iter = #4,        x_best=[0.661827 0.083743],     f_best = 1.964796e+001,         T = 1.059628e+000
+      Iter = #5,        x_best=[0.931415 0.820681],     f_best = 1.416565e+001,         T = 9.536654e-001
+      Iter = #6,        x_best=[1.849796 1.222178],     f_best = 7.784028e+000,         T = 8.582988e-001
+      Iter = #7,        x_best=[2.539775 1.414591],     f_best = 4.645780e+000,         T = 7.724690e-001
+      Iter = #8,        x_best=[3.206047 2.394497],     f_best = 9.969957e-001,         T = 6.952221e-001
+      Iter = #9,        x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 6.256999e-001
+      Iter = #10,       x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.631299e-001
+      Iter = #11,       x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.068169e-001
+      Iter = #12,       x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.561352e-001
+      Iter = #13,       x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.105217e-001
+      Iter = #14,       x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 3.694695e-001
+      Iter = #15,       x_best=[3.931929 3.003181],     f_best = 4.643767e-003,         T = 3.325226e-001
+    </programlisting>
+
+
+  </refsection>
+
+
+<refsection>
+    <title>See Also</title>
+
+    <simplelist type="inline">
+      <member>
+        <link linkend="compute_initial_temp">
+          compute_initial_temp
+        </link>
+      </member>
+
+      <member>
+        <link linkend="neigh_func_default">
+          neigh_func_default
+        </link>
+      </member>
+
+      <member>
+        <link linkend="temp_law_default">
+          temp_law_default
+        </link>
+      </member>
     </simplelist>
   </refsection>
 
   <refsection>
     <title>Authors</title>
 
-    <variablelist>
-      <varlistentry>
-        <term>collette</term>
+    <para>2008, Yann COLLETTE (ycollet@freesurf.fr)</para>
+    <para>Michael Baudin - Digiteo - 2010</para>
+  </refsection>
 
-        <listitem>
-          <para>Yann COLLETTE (ycollet@freesurf.fr)</para>
-        </listitem>
-      </varlistentry>
-    </variablelist>
+  <refsection>
+    <title>Bibliography</title>
+
+    <para>
+      "Simulated annealing : theory and applications", P.J.M. Laarhoven and E.H.L. Aarts, Mathematics and its applications, Dordrecht : D. Reidel, 1988
+    </para>
+    <para>
+      "Theoretical and computational aspects of simulated annealing", P.J.M. van Laarhoven, Amsterdam, Netherlands : Centrum voor Wiskunde en Informatica, 1988
+    </para>
+    <para>
+      "Genetic algorithms and simulated annealing", Lawrence Davis, London : Pitman Los Altos, Calif. Morgan Kaufmann Publishers, 1987
+    </para>
   </refsection>
 </refentry>
+
index 5115eb3..da5903c 100644 (file)
@@ -1,5 +1,6 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
+// Copyright (C) 2010 - 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
 // are also available at
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
-function [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history] = optim_sa(x0, sa_f, ItExt, ItInt, T0, Log, param)
-// Simulated annealing
-// x0         : initial solution
-// f          : objective function
-// ItExt      : number of temperature decrease
-// ItInt      : number of iterations during one temperature step
-// T0         : initial temperature
-// Log        : print some message during the run of the optimization
-// param      : a parameter list. this list contains the neighobrhood ('neigh_func') and some parameters related to this neighborhood functions (see the 
-//              related function to list the available parameters)
-
-[nargout, nargin] = argn();
-
-if ~isdef('param','local') then
-  param = [];
-end
-
-[temp_law,err]    = get_param(param,'temp_law',temp_law_default);
-[neigh_func,err]  = get_param(param,'neigh_func',neigh_func_default);
-[accept_func,err] = get_param(param,'accept_func',accept_func_default);
-
-if (~isdef('Log','local')) then
-  Log = %F;
-end
-
-if (nargout>=6) then
-  f_history_defined = %T;
-  f_history = list();
-else
-  f_history_defined = %F;
-end
-
-if (nargout>=5) then
-  temp_list_defined = %T;
-  temp_list = [];
-else
-  temp_list_defined = %F;
-end
-
-if (nargout>=7) then
-  x_history_defined = %T;
-  x_history = list();
-else
-  x_history_defined = %F;
-end
-
-if ~isdef('sa_f','local') then
-  error(gettext("optim_sa: sa_f is mandatory"));
-else
-  if typeof(sa_f)=='list' then
-    deff('y=_sa_f(x)','y=sa_f(1)(x, sa_f(2:$))');
+function [x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , iter ] = optim_sa(x0, sa_f, ItExt, ItInt, T0, Log, param)
+  // Simulated annealing
+  // x0         : initial solution
+  // f          : objective function
+  // ItExt      : number of temperature decrease
+  // ItInt      : number of iterations during one temperature step
+  // T0         : initial temperature
+  // Log        : print some message during the run of the optimization
+  // param      : a parameter list. this list contains the neighobrhood ('neigh_func') and some parameters related to this neighborhood functions (see the 
+  //              related function to list the available parameters)
+  
+  [nargout, nargin] = argn();
+  
+  if ~isdef("param","local") then
+    param = [];
+  end
+  
+  [temp_law,err]    = get_param(param,"temp_law",temp_law_default);
+  [neigh_func,err]  = get_param(param,"neigh_func",neigh_func_default);
+  [accept_func,err] = get_param(param,"accept_func",accept_func_default);
+  [output_func,err] = get_param(param,"output_func",[]);
+  
+  if (~isdef("Log","local")) then
+    Log = %F;
+  end
+  
+  if (nargout>=6) then
+    f_history_defined = %T;
+    f_history = list();
   else
-    deff('y=_sa_f(x)','y=sa_f(x)');
+    f_history_defined = %F;
   end
-end
-
-T = T0;
-
-// Some variables needed to record the behavior of the SA
-var_list  = [];
-mean_list = [];
-temp_list = [];
-
-x_current = x0;
-f_current = _sa_f(x_current);
-
-x_best = x_current;
-f_best = f_current;
-
-for i=1:ItExt
-  f_list = [];
-  x_list = list();
-  for j=1:ItInt
-    x_neigh = neigh_func(x_current,T,param);
-    f_neigh = _sa_f(x_neigh);
-    if ((f_neigh<=f_current)|(accept_func(f_current,f_neigh,T)>rand(1,1))) then
-      x_current = x_neigh;
-      f_current = f_neigh;
+  
+  if (nargout>=5) then
+    temp_list_defined = %T;
+    temp_list = [];
+  else
+    temp_list_defined = %F;
+  end
+  
+  if (nargout>=7) then
+    x_history_defined = %T;
+    x_history = list();
+  else
+    x_history_defined = %F;
+  end
+  
+  if ~isdef("sa_f","local") then
+    error(gettext("optim_sa: sa_f is mandatory"));
+  else
+    if typeof(sa_f)=="list" then
+      deff("y=_sa_f(x)","y=sa_f(1)(x, sa_f(2:$))");
+    else
+      deff("y=_sa_f(x)","y=sa_f(x)");
+    end
+  end
+  
+  T = T0;
+  
+  // Some variables needed to record the behavior of the SA
+  var_list  = [];
+  mean_list = [];
+  temp_list = [];
+  
+  x_current = x0;
+  f_current = _sa_f(x_current);
+  
+  x_best = x_current;
+  f_best = f_current;
+  
+  for iter=1:ItExt
+    if ( output_func <> [] ) then
+      stop = output_func ( iter , x_best , f_best , T , saparams );
+      if (stop) then
+        break
+      end
     end
     
-    f_list = [f_list f_current];
-    
-    if (f_best>f_current) then
-      x_best = x_current;
-      f_best = f_current;
+    f_list = [];
+    x_list = list();
+    for j=1:ItInt
+      x_neigh = neigh_func(x_current,T,param);
+      f_neigh = _sa_f(x_neigh);
+      if ((f_neigh<=f_current)|(accept_func(f_current,f_neigh,T)>rand(1,1))) then
+        x_current = x_neigh;
+        f_current = f_neigh;
+      end
+      
+      f_list = [f_list f_current];
+      
+      if (f_best>f_current) then
+        x_best = x_current;
+        f_best = f_current;
+      end
+      
+      if (x_history_defined) then
+        x_list($+1) = x_current;
+      end
     end
     
+    if (temp_list_defined) then
+      temp_list = [temp_list T];
+    end
     if (x_history_defined) then
-      x_list($+1) = x_current;
+      x_history($+1) = x_list;
     end
+    if (f_history_defined) then
+      f_history($+1) = f_list;
+    end
+    
+    // Computation of step_mean and step_var
+    step_mean = mean(f_list);
+    step_var  = stdev(f_list);
+    mean_list = [mean_list step_mean];
+    var_list  = [var_list step_var];
+    
+    if (Log) then
+      printf(gettext("%s: Temperature step %d / %d - T = %f, E(f(T)) = %f var(f(T)) = %f f_best = %f\n"), "optim_sa", iter, ItExt, T, step_mean, step_var, f_best);
+    end
+    
+    T = temp_law(T, step_mean, step_var, iter, max(size(x_current)), param);
   end
-
-  if (temp_list_defined) then
-    temp_list = [temp_list T];
-  end
-  if (x_history_defined) then
-    x_history($+1) = x_list;
-  end
-  if (f_history_defined) then
-    f_history($+1) = f_list;
-  end
-
-  // Computation of step_mean and step_var
-  step_mean = mean(f_list);
-  step_var  = stdev(f_list);
-  mean_list = [mean_list step_mean];
-  var_list  = [var_list step_var];
-  
-  if (Log) then
-    printf(gettext("%s: Temperature step %d / %d - T = %f, E(f(T)) = %f var(f(T)) = %f f_best = %f\n"), "optim_sa", i, ItExt, T, step_mean, step_var, f_best);
-  end
-
-  T = temp_law(T, step_mean, step_var, i, max(size(x_current)), param);
-end
 endfunction
 
+
index 406672c..7ec0d35 100644 (file)
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
-//
-// 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
-deff('y=test_func(x)','y=x^2');
-x0 = 10;
-T_init = compute_initial_temp(x0, test_func, 0.8, 1000, []);
-[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history] = optim_sa(x0, test_func, 10, 1000, T_init, %F, []);
-if (x_best==%nan) | (x_best==%inf) then bugmes();quit;end
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab\r
+// Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>\r
+// Copyright (C) 2010 - DIGITEO - Michael Baudin\r
+//\r
+// This file must be used under the terms of the CeCILL.\r
+// This source file is licensed as described in the file COPYING, which\r
+// you should have received as part of this distribution.  The terms\r
+// are also available at\r
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt\r
+// <-- JVM_NOT_MANDATORY -->\r
+// <-- ENGLISH IMPOSED -->\r
+//\r
+// assert_close --\r
+//   Returns 1 if the two real matrices computed and expected are close,\r
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.\r
+// Arguments\r
+//   computed, expected : the two matrices to compare\r
+//   epsilon : a small number\r
+//\r
+function flag = assert_close ( computed, expected, epsilon )\r
+  if expected==0.0 then\r
+    shift = norm(computed-expected);\r
+  else\r
+    shift = norm(computed-expected)/norm(expected);\r
+  end\r
+  if shift < epsilon then\r
+    flag = 1;\r
+  else\r
+    flag = 0;\r
+  end\r
+  if flag <> 1 then bugmes();quit;end\r
+endfunction\r
+//\r
+// assert_equal --\r
+//   Returns 1 if the two real matrices computed and expected are equal.\r
+// Arguments\r
+//   computed, expected : the two matrices to compare\r
+//   epsilon : a small number\r
+//\r
+function flag = assert_equal ( computed , expected )\r
+  if computed==expected then\r
+    flag = 1;\r
+  else\r
+    flag = 0;\r
+  end\r
+  if flag <> 1 then bugmes();quit;end\r
+endfunction\r
+///////////////////////////////////////////\r
+// Test that we can run with default values for parameters\r
+function y=test_func(x)\r
+  y=x^2\r
+endfunction\r
+x0 = 10;\r
+T_init = compute_initial_temp(x0, test_func, 0.8, 1000);\r
+[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history,iter] = optim_sa(x0, test_func, 10, 1000, T_init, %F);\r
+assert_close ( x_best , 0.0 , 1.e-2 );\r
+assert_close ( f_best , 0.0 , 1.e-2 );\r
+assert_equal ( size(mean_list) , [1 10] );\r
+assert_equal ( size(var_list) , [1 10] );\r
+assert_equal ( size(temp_list) , [1 10] );\r
+assert_equal ( size(f_history) , 10 );\r
+assert_equal ( size(x_history) , 10 );\r
+assert_equal ( iter>0 , %t );\r
+///////////////////////////////////////////\r
+// Test that we can configure our own neighbour function\r
+function f = quad ( x )\r
+  p = [4 3];\r
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2\r
+endfunction\r
+// We produce a neighbor by adding some noise to each component of a given vector\r
+function x_neigh = myneigh_func ( x_current, T , param)\r
+  nxrow = size(x_current,"r")\r
+  nxcol = size(x_current,"c")\r
+  sa_min_delta = -0.1*ones(nxrow,nxcol);\r
+  sa_max_delta = 0.1*ones(nxrow,nxcol);\r
+  x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;\r
+endfunction\r
+x0          = [2 2];\r
+Proba_start = 0.7;\r
+It_Pre      = 100;\r
+It_extern   = 50;\r
+It_intern   = 100;\r
+saparams = init_param();\r
+saparams = add_param(saparams,"neigh_func", myneigh_func);\r
+T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);\r
+Warning : redefining function: result                  . Use funcprot(0) to avoid this message\r
+\r
+Log = %f;\r
+[x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams);\r
+Warning : redefining function: result                  . Use funcprot(0) to avoid this message\r
+\r
+assert_close ( x_opt , [4 3] ,  1.e-1 );\r
+assert_close ( f_opt , 0 ,  1.e-1 );\r
+///////////////////////////////////////////\r
+// Test that an additionnal parameter can be passed to the cost function\r
+function f = quadp ( x , p )\r
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2\r
+endfunction\r
+x0 = [-1 -1];\r
+p = [4 3];\r
+T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);\r
+[x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 10, 1000, T0, %f);\r
+assert_close ( x_opt , [4 3] ,  1.e-1 );\r
+assert_close ( f_opt , 0 ,  1.e-1 );\r
+///////////////////////////////////////////\r
+// Test with a plot function, which serves also as a stop function.\r
+function f = quad ( x )\r
+  p = [4 3];\r
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2\r
+endfunction\r
+// See that the stop variable becomes true when the function value is near zero.\r
+// The threshold is rather loose.\r
+function stop = outfunc ( itExt , x_best , f_best , T , saparams )\r
+  [mythreshold,err] = get_param(saparams,"mythreshold",0);\r
+  mprintf ( "Iter = #%d, \t x_best=[%f %f], \t f_best = %e, \t T = %e\n", itExt , x_best(1),x_best(2) , f_best , T )\r
+  stop = ( abs(f_best) < mythreshold )\r
+endfunction\r
+x0 = [-1 -1];\r
+saparams = init_param();\r
+saparams = add_param(saparams,"output_func", outfunc );\r
+saparams = add_param(saparams,"mythreshold", 1.e-2 );\r
+rand("seed",0);\r
+T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);\r
+// Notice that the number of external iterations is %inf, so \r
+// that the external loop never stops.\r
+// This allows to check that the output function really allows to \r
+// stop the algorithm.\r
+[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , iter ] = optim_sa(x0, quad , %inf, 100, T0, %f, saparams);\r
+Iter = #1,      x_best=[-1.000000 -1.000000],   f_best = 4.100000e+001,         T = 1.453537e+000\r\r
+Iter = #2,      x_best=[-0.408041 -0.318262],   f_best = 3.044169e+001,         T = 1.308183e+000\r\r
+Iter = #3,      x_best=[-0.231406 -0.481078],   f_best = 3.002270e+001,         T = 1.177365e+000\r\r
+Iter = #4,      x_best=[0.661827 0.083743],     f_best = 1.964796e+001,         T = 1.059628e+000\r\r
+Iter = #5,      x_best=[0.931415 0.820681],     f_best = 1.416565e+001,         T = 9.536654e-001\r\r
+Iter = #6,      x_best=[1.849796 1.222178],     f_best = 7.784028e+000,         T = 8.582988e-001\r\r
+Iter = #7,      x_best=[2.539775 1.414591],     f_best = 4.645780e+000,         T = 7.724690e-001\r\r
+Iter = #8,      x_best=[3.206047 2.394497],     f_best = 9.969957e-001,         T = 6.952221e-001\r\r
+Iter = #9,      x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 6.256999e-001\r\r
+Iter = #10,     x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.631299e-001\r\r
+Iter = #11,     x_best=[3.164998 2.633170],     f_best = 8.317924e-001,         T = 5.068169e-001\r\r
+Iter = #12,     x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.561352e-001\r\r
+Iter = #13,     x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 4.105217e-001\r\r
+Iter = #14,     x_best=[3.961464 2.903763],     f_best = 1.074654e-002,         T = 3.694695e-001\r\r
+Iter = #15,     x_best=[3.931929 3.003181],     f_best = 4.643767e-003,         T = 3.325226e-001\r\r
+assert_close ( x_best , [4 3] ,  1.e-1 );\r
+assert_close ( f_best , 0 ,  1.e-2 );\r
+assert_equal ( iter > 0 , %t );\r
index 8e53bc6..914c457 100755 (executable)
@@ -1,6 +1,6 @@
-
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) 2008 - Yann COLLETTE <yann.collette@renault.com>
+// Copyright (C) 2010 - 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
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
 // <-- JVM_NOT_MANDATORY -->
+// <-- ENGLISH IMPOSED -->
+
 
-deff('y=test_func(x)','y=x^2');
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+
+///////////////////////////////////////////
+// Test that we can run with default values for parameters
+function y=test_func(x)
+  y=x^2
+endfunction
 
 x0 = 10;
+T_init = compute_initial_temp(x0, test_func, 0.8, 1000);
+[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history,iter] = optim_sa(x0, test_func, 10, 1000, T_init, %F);
+
+assert_close ( x_best , 0.0 , 1.e-2 );
+assert_close ( f_best , 0.0 , 1.e-2 );
+assert_equal ( size(mean_list) , [1 10] );
+assert_equal ( size(var_list) , [1 10] );
+assert_equal ( size(temp_list) , [1 10] );
+assert_equal ( size(f_history) , 10 );
+assert_equal ( size(x_history) , 10 );
+assert_equal ( iter>0 , %t );
+
+///////////////////////////////////////////
+// Test that we can configure our own neighbour function
+function f = quad ( x )
+  p = [4 3];
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+endfunction
+
+// We produce a neighbor by adding some noise to each component of a given vector
+function x_neigh = myneigh_func ( x_current, T , param)
+  nxrow = size(x_current,"r")
+  nxcol = size(x_current,"c")
+  sa_min_delta = -0.1*ones(nxrow,nxcol);
+  sa_max_delta = 0.1*ones(nxrow,nxcol);
+  x_neigh = x_current + (sa_max_delta - sa_min_delta).*rand(nxrow,nxcol) + sa_min_delta;
+endfunction
+
+x0          = [2 2];
+Proba_start = 0.7;
+It_Pre      = 100;
+It_extern   = 50;
+It_intern   = 100;
+
+saparams = init_param();
+saparams = add_param(saparams,"neigh_func", myneigh_func);
+
+T0 = compute_initial_temp(x0, quad, Proba_start, It_Pre, saparams);
+Log = %f;
+[x_opt, f_opt] = optim_sa(x0, quad, It_extern, It_intern, T0, Log,saparams);
+assert_close ( x_opt , [4 3] ,  1.e-1 );
+assert_close ( f_opt , 0 ,  1.e-1 );
+
+///////////////////////////////////////////
+// Test that an additionnal parameter can be passed to the cost function
+
+function f = quadp ( x , p )
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+endfunction
+
+x0 = [-1 -1];
+p = [4 3];
+T0 = compute_initial_temp(x0, list(quadp,p) , Proba_start, It_Pre);
+[x_opt, f_opt] = optim_sa(x0, list(quadp,p) , 10, 1000, T0, %f);
+assert_close ( x_opt , [4 3] ,  1.e-1 );
+assert_close ( f_opt , 0 ,  1.e-1 );
+
+///////////////////////////////////////////
+// Test with a plot function, which serves also as a stop function.
+
+function f = quad ( x )
+  p = [4 3];
+  f = (x(1) - p(1))^2 + (x(2) - p(2))^2
+endfunction
+
+// See that the stop variable becomes true when the function value is near zero.
+// The threshold is rather loose.
+function stop = outfunc ( itExt , x_best , f_best , T , saparams )
+  [mythreshold,err] = get_param(saparams,"mythreshold",0);
+  mprintf ( "Iter = #%d, \t x_best=[%f %f], \t f_best = %e, \t T = %e\n", itExt , x_best(1),x_best(2) , f_best , T )
+  stop = ( abs(f_best) < mythreshold )
+endfunction
+
+x0 = [-1 -1];
+saparams = init_param();
+saparams = add_param(saparams,"output_func", outfunc );
+saparams = add_param(saparams,"mythreshold", 1.e-2 );
+
+rand("seed",0);
 
-T_init = compute_initial_temp(x0, test_func, 0.8, 1000, []);
+T0 = compute_initial_temp(x0, quad , 0.7, 100, saparams);
+// Notice that the number of external iterations is %inf, so 
+// that the external loop never stops.
+// This allows to check that the output function really allows to 
+// stop the algorithm.
+[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history , iter ] = optim_sa(x0, quad , %inf, 100, T0, %f, saparams);
+assert_close ( x_best , [4 3] ,  1.e-1 );
+assert_close ( f_best , 0 ,  1.e-2 );
+assert_equal ( iter > 0 , %t );
 
-[x_best, f_best, mean_list, var_list, temp_list, f_history, x_history] = optim_sa(x0, test_func, 10, 1000, T_init, %F, []);
 
-if (x_best==%nan) | (x_best==%inf) then pause,end