Clarified constructor of optimsimplex. Added cost function checker in optimbase
Michaël Baudin [Mon, 28 Sep 2009 15:10:42 +0000 (17:10 +0200)]
43 files changed:
scilab/modules/optimization/demos/neldermead/neldermead_outputcmd.sce [new file with mode: 0644]
scilab/modules/optimization/help/en_US/optimbase/optimbase.xml
scilab/modules/optimization/help/en_US/optimsimplex/optimsimplex.xml
scilab/modules/optimization/macros/neldermead/fminsearch.sci
scilab/modules/optimization/macros/neldermead/neldermead_cget.sci
scilab/modules/optimization/macros/neldermead/neldermead_configure.sci
scilab/modules/optimization/macros/neldermead/neldermead_new.sci
scilab/modules/optimization/macros/neldermead/neldermead_search.sci
scilab/modules/optimization/macros/optimbase/optimbase_checkcostfun.sci [new file with mode: 0644]
scilab/modules/optimization/macros/optimbase/optimbase_function.sci
scilab/modules/optimization/macros/optimbase/optimbase_isinbounds.sci [new file with mode: 0644]
scilab/modules/optimization/macros/optimbase/optimbase_isinnonlincons.sci [new file with mode: 0644]
scilab/modules/optimization/macros/optimsimplex/optimsimplex_axes.sci [deleted file]
scilab/modules/optimization/macros/optimsimplex/optimsimplex_new.sci
scilab/modules/optimization/macros/optimsimplex/optimsimplex_pfeffer.sci [deleted file]
scilab/modules/optimization/macros/optimsimplex/optimsimplex_randbounds.sci [deleted file]
scilab/modules/optimization/macros/optimsimplex/optimsimplex_spendley.sci [deleted file]
scilab/modules/optimization/macros/optimsimplex/optimsimplex_xbar.sci
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.tst [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkx0.dia.ref
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.tst [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isfeasible.dia.ref
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.tst [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.tst [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_proj2bnds.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_axes.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_axes.tst
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_getset.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_pfeffer.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_pfeffer.tst
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_print.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_randbounds.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_randbounds.tst
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_spendley.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_spendley.tst
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_tostring.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_tostring.tst
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_xbar.dia.ref
scilab/modules/optimization/tests/unit_tests/optimsimplex/optimsimplex_xbar.tst

diff --git a/scilab/modules/optimization/demos/neldermead/neldermead_outputcmd.sce b/scilab/modules/optimization/demos/neldermead/neldermead_outputcmd.sce
new file mode 100644 (file)
index 0000000..78334d7
--- /dev/null
@@ -0,0 +1,68 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+function y = rosenbrock (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+
+
+//
+// myoutputcmd --
+//  This command is called back by the Nelder-Mead
+//  algorithm.
+// Arguments
+//  state : the current state of the algorithm
+//    "init", "iter", "done"
+//  data : the data at the current state
+//    This is a tlist with the following entries:
+//    * x : the optimal vector of parameters
+//    * fval : the minimum function value
+//    * simplex : the simplex, as a simplex object
+//    * iteration : the number of iterations performed
+//    * funccount : the number of function evaluations
+//
+function myoutputcmd ( state , data )
+  iter = data.iteration
+  if ( state == "init" ) then
+    mprintf ( "=================================\n");
+    mprintf ( "Initialization\n");
+  elseif ( state == "done" ) then
+    mprintf ( "=================================\n");
+    mprintf ( "End of Optimization\n");
+  end
+  fc = data.funccount
+  fval = data.fval
+  x = data.x
+  simplex = data.simplex
+  // Simplex is a data structure, which can be managed
+  // by the optimsimplex class.
+  ssize = optimsimplex_size ( simplex )
+  mprintf ( "Iteration #%d, Feval #%d, Fval = %e, x = %s, Size = %e\n", iter, fc, fval, strcat(string(x)," "), ssize);
+endfunction
+
+
+nm = neldermead_new ();
+nm = neldermead_configure(nm,"-numberofvariables",2);
+nm = neldermead_configure(nm,"-function",rosenbrock);
+nm = neldermead_configure(nm,"-x0",[-1.2 1.0]');
+nm = neldermead_configure(nm,"-maxiter",200);
+nm = neldermead_configure(nm,"-maxfunevals",300);
+nm = neldermead_configure(nm,"-tolfunrelative",10*%eps);
+nm = neldermead_configure(nm,"-tolxrelative",10*%eps);
+nm = neldermead_configure(nm,"-simplex0method","axes");
+nm = neldermead_configure(nm,"-simplex0length",1.0);
+nm = neldermead_configure(nm,"-method","variable");
+nm = neldermead_configure(nm,"-verbose",0);
+nm = neldermead_configure(nm,"-verbosetermination",0);
+nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
+nm = neldermead_search(nm);
+nm = neldermead_destroy(nm);
+
+
index fb74698..5bf1cf6 100644 (file)
@@ -45,45 +45,48 @@ data = optimbase_outstruct ( this )
 [ this , p ] = optimbase_proj2bnds ( this ,  x )
 this = optimbase_stoplog ( this , msg )
 [this , terminate , status] = optimbase_terminate (this , previousfopt , currentfopt , previousxopt , currentxopt )
+this = optimbase_checkcostfun ( this )
+[ this , isfeasible ] = optimbase_isinbounds ( this , x )
+[ this , isfeasible ] = optimbase_isinnonlinconst ( this , x )
 </synopsis>
   </refsynopsisdiv>
 
   <refsection>
     <title>Purpose</title>
 
-    <para>
-      The goal of this component is to provide a building block for
-      optimization methods. The goal is to provide a building block for a large
-      class of specialized optimization methods. This component manages</para>
-      <itemizedlist>
-        <listitem>
-          <para>the number of variables,</para>
-        </listitem>
+    <para>The goal of this component is to provide a building block for
+    optimization methods. The goal is to provide a building block for a large
+    class of specialized optimization methods. This component manages</para>
 
-        <listitem>
-          <para>the minimum and maximum bounds,</para>
-        </listitem>
+    <itemizedlist>
+      <listitem>
+        <para>the number of variables,</para>
+      </listitem>
 
-        <listitem>
-          <para>the number of non linear inequality constraints,</para>
-        </listitem>
+      <listitem>
+        <para>the minimum and maximum bounds,</para>
+      </listitem>
 
-        <listitem>
-          <para>the cost function,</para>
-        </listitem>
+      <listitem>
+        <para>the number of non linear inequality constraints,</para>
+      </listitem>
 
-        <listitem>
-          <para>the logging system,</para>
-        </listitem>
+      <listitem>
+        <para>the cost function,</para>
+      </listitem>
 
-        <listitem>
-          <para>various termination criteria,</para>
-        </listitem>
+      <listitem>
+        <para>the logging system,</para>
+      </listitem>
 
-        <listitem>
-          <para>etc...</para>
-        </listitem>
-      </itemizedlist>
+      <listitem>
+        <para>various termination criteria,</para>
+      </listitem>
+
+      <listitem>
+        <para>etc...</para>
+      </listitem>
+    </itemizedlist>
   </refsection>
 
   <refsection>
@@ -115,23 +118,24 @@ this = optimbase_stoplog ( this , msg )
       </listitem>
 
       <listitem>
-        <para>Manage various termination criteria, including </para>
-       <itemizedlist>
-            <listitem>
-              <para>maximum number of iterations,</para>
-            </listitem>
+        <para>Manage various termination criteria, including</para>
 
-            <listitem>
-              <para>tolerance on function value (relative or absolute),</para>
-            </listitem>
+        <itemizedlist>
+          <listitem>
+            <para>maximum number of iterations,</para>
+          </listitem>
 
-            <listitem>
-              <para>tolerance on x (relative or absolute),</para>
-            </listitem>
+          <listitem>
+            <para>tolerance on function value (relative or absolute),</para>
+          </listitem>
 
-            <listitem>
-              <para>maximum number of evaluations of cost function,</para>
-            </listitem>
+          <listitem>
+            <para>tolerance on x (relative or absolute),</para>
+          </listitem>
+
+          <listitem>
+            <para>maximum number of evaluations of cost function,</para>
+          </listitem>
         </itemizedlist>
       </listitem>
 
@@ -197,11 +201,11 @@ this = optimbase_stoplog ( this , msg )
 
     <para>The optimization problem to solve is the following</para>
 
-    <programlisting role="example"><![CDATA[ 
+    <programlisting role="example"> 
 min f(x)
-l_i <= x_i <= h_i, i = 1,n
-g_i(x) <= 0, i = 1,nbineq
- ]]></programlisting>
+l_i &lt;= x_i &lt;= h_i, i = 1,n
+g_i(x) &lt;= 0, i = 1,nbineq
+ </programlisting>
 
     <para>where</para>
 
@@ -1167,6 +1171,88 @@ g_i(x) <= 0, i = 1,nbineq
           </variablelist>
         </listitem>
       </varlistentry>
+
+      <varlistentry>
+        <term>this = optimbase_checkcostfun ( this )</term>
+
+        <listitem>
+          <para>Check that the cost function is correctly connected. Generate
+          an error if there is one. Takes into account for the cost function
+          at the initial guess x0 only. If there are nonlinear constraints,
+          check that the index flag is correctly used for index=1,2 or 3.
+          </para>
+
+          <para>This function requires at least one call to the cost function
+          to make the necessary checks.</para>
+
+          <variablelist>
+            <varlistentry>
+              <term>this</term>
+
+              <listitem>
+                <para>The current object.</para>
+              </listitem>
+            </varlistentry>
+          </variablelist>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>[ this , isfeasible ] = optimbase_isinbounds ( this , x )</term>
+
+        <listitem>
+          <para>Returns isfeasible = %t if the given point satisfies bounds
+          constraints. Returns isfeasible = %f if the given point is not in
+          the bounds.</para>
+
+          <variablelist>
+            <varlistentry>
+              <term>this</term>
+
+              <listitem>
+                <para>The current object.</para>
+              </listitem>
+            </varlistentry>
+
+            <varlistentry>
+              <term>isfeasible</term>
+
+              <listitem>
+                <para>a boolean</para>
+              </listitem>
+            </varlistentry>
+          </variablelist>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>[ this , isfeasible ] = optimbase_isinnonlinconst ( this , x
+        )</term>
+
+        <listitem>
+          <para>Returns isfeasible = %t if the given point satisfies the
+          nonlinear constraints. Returns isfeasible = %f if the given point
+          does not satisfy the nonlinear constraints.</para>
+
+          <variablelist>
+            <varlistentry>
+              <term>this</term>
+
+              <listitem>
+                <para>The current object.</para>
+              </listitem>
+            </varlistentry>
+
+            <varlistentry>
+              <term>isfeasible</term>
+
+              <listitem>
+                <para>a boolean</para>
+              </listitem>
+            </varlistentry>
+          </variablelist>
+        </listitem>
+      </varlistentry>
     </variablelist>
   </refsection>
 
@@ -1180,9 +1266,9 @@ g_i(x) <= 0, i = 1,nbineq
     <para>In the more general case, the cost function is expected to have the
     following header</para>
 
-    <programlisting role="example"><![CDATA[ 
+    <programlisting role="example"> 
 function y = myfunction(x, index, data)
- ]]></programlisting>
+ </programlisting>
 
     <para>where</para>
 
@@ -1257,9 +1343,9 @@ function y = myfunction(x, index, data)
     <para>In the most simple case, the cost function is expected to have the
     following header</para>
 
-    <programlisting role="example"><![CDATA[ 
+    <programlisting role="example"> 
 function y = myfunction(x)
- ]]></programlisting>
+ </programlisting>
 
     <para>where x is the current point and y is the value of the cost. This
     case is associated with an unconstrained problem without any additionnal
@@ -1275,9 +1361,9 @@ function y = myfunction(x)
 
     <para>The output function must have the following header</para>
 
-    <programlisting role="example"><![CDATA[ 
+    <programlisting role="example"> 
 function outputcmd(state, data, myobj)
- ]]></programlisting>
+ </programlisting>
 
     <para>where</para>
 
@@ -1382,9 +1468,9 @@ function outputcmd(state, data, myobj)
         <para>The number of iterations is examined and compared to the
         -maxiter option : if the following condition</para>
 
-        <programlisting role="example"><![CDATA[ 
-iterations >= maxiter
- ]]></programlisting>
+        <programlisting role="example"> 
+iterations &gt;= maxiter
+ </programlisting>
 
         <para>is true, then the status is set to "maxiter" and terminate is
         set to 1.</para>
@@ -1394,9 +1480,9 @@ iterations >= maxiter
         <para>The number of function evaluations and compared to the
         -maxfunevals option is examined : if the following condition</para>
 
-        <programlisting role="example"><![CDATA[ 
-funevals >= maxfunevals
- ]]></programlisting>
+        <programlisting role="example"> 
+funevals &gt;= maxfunevals
+ </programlisting>
 
         <para>is true, then the status is set to "maxfuneval" and terminate is
         set to 1.</para>
@@ -1421,9 +1507,9 @@ funevals >= maxfunevals
             <listitem>
               <para>if the following condition</para>
 
-              <programlisting role="example"><![CDATA[ 
-abs(currentfopt) < tolfunrelative * abs(previousfopt) + tolfunabsolute
- ]]></programlisting>
+              <programlisting role="example"> 
+abs(currentfopt) &lt; tolfunrelative * abs(previousfopt) + tolfunabsolute
+ </programlisting>
 
               <para>is true, then the status is set to "tolf" and terminate is
               set to 1.</para>
@@ -1459,9 +1545,9 @@ abs(currentfopt) < tolfunrelative * abs(previousfopt) + tolfunabsolute
             <listitem>
               <para>if the following condition</para>
 
-              <programlisting role="example"><![CDATA[ 
-norm(currentxopt - previousxopt) < tolxrelative * norm(currentxopt) + tolxabsolute
- ]]></programlisting>
+              <programlisting role="example"> 
+norm(currentxopt - previousxopt) &lt; tolxrelative * norm(currentxopt) + tolxabsolute
+ </programlisting>
 
               <para>is true, then the status is set to "tolx" and terminate is
               set to 1.</para>
@@ -1493,7 +1579,7 @@ norm(currentxopt - previousxopt) < tolxrelative * norm(currentxopt) + tolxabsolu
     the algorithm. After the optimization is performed, the optimum is
     retrieved with quiery features.</para>
 
-    <programlisting role="example"><![CDATA[ 
+    <programlisting role="example"> 
 function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
@@ -1506,7 +1592,7 @@ nbvar = optimbase_cget(opt,"-numberofvariables");
 opt = optimbase_configure(opt,"-function",rosenbrock);
 [this,f] = optimbase_function(opt, [0.0 0.0]);
 opt = optimbase_destroy(opt);
- ]]></programlisting>
+ </programlisting>
   </refsection>
 
   <refsection>
@@ -1519,26 +1605,26 @@ opt = optimbase_destroy(opt);
     <title>TODO</title>
 
     <itemizedlist>
-        <listitem>
-          <para>manage equality constraints</para>
-        </listitem>
+      <listitem>
+        <para>manage equality constraints</para>
+      </listitem>
 
-        <listitem>
-          <para>manage linear constraints</para>
-        </listitem>
+      <listitem>
+        <para>manage linear constraints</para>
+      </listitem>
 
-        <listitem>
-          <para>manage quadratic objective</para>
-        </listitem>
+      <listitem>
+        <para>manage quadratic objective</para>
+      </listitem>
 
-        <listitem>
-          <para>manage linear objective</para>
-        </listitem>
+      <listitem>
+        <para>manage linear objective</para>
+      </listitem>
 
-        <listitem>
-          <para>methods with derivatives : add a flag to manage for the
-          derivatives of the cost function</para>
-        </listitem>
+      <listitem>
+        <para>methods with derivatives : add a flag to manage for the
+        derivatives of the cost function</para>
+      </listitem>
     </itemizedlist>
   </refsection>
 </refentry>
index 11b5336..5ca237c 100644 (file)
     <title>SYNOPSIS</title>
 
     <synopsis>[ newobj , data ] = optimsimplex_new ( coords , fun , data )
+[ newobj , data ] = optimsimplex_new ( "axes" , x0 , fun , len , data )
+[ newobj , data ] = optimsimplex_new ( "spendley" , x0 , fun , len , data )
+[ newobj , data ] = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual , deltazero , data )
+[ newobj , data ] = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax , nbve  , data )
 this = optimsimplex_destroy (this)
-[ this , data ] = optimsimplex_axes ( this , x0 , fun , len , data )
-[ this , data ] = optimsimplex_pfeffer ( this , x0 , fun , deltausual , deltazero , data )
-[ this , data ] = optimsimplex_randbounds ( this , x0 , fun , boundsmin , boundsmax , nbpoints  , data )
-[ this , data ] = optimsimplex_spendley ( this , x0 , fun , len , data )
 this = optimsimplex_setall ( this , simplex )
 this = optimsimplex_setallfv ( this , fv )
 this = optimsimplex_setallx ( this , x )
@@ -217,6 +217,12 @@ cen = optimsimplex_xbar ( this , iexcl )</synopsis>
 
     <variablelist>
       <varlistentry>
+        <term>newobj = optimsimplex_new ( )</term>
+
+        <term>newobj = optimsimplex_new ( coords )</term>
+
+        <term>newobj = optimsimplex_new ( coords , fun )</term>
+
         <term>[ newobj , data ] = optimsimplex_new ( coords , fun , data
         )</term>
 
@@ -290,29 +296,17 @@ function [ y , data ] = myfunction ( x , data )
       </varlistentry>
 
       <varlistentry>
-        <term>this = optimsimplex_destroy (this)</term>
+        <term>this = optimsimplex_new ( "axes" , x0 )</term>
 
-        <listitem>
-          <para>Destroy the given object.</para>
+        <term>this = optimsimplex_new ( "axes" , x0 , fun )</term>
 
-          <variablelist>
-            <varlistentry>
-              <term>this</term>
+        <term>this = optimsimplex_new ( "axes" , x0 , fun , len )</term>
 
-              <listitem>
-                <para>The current simplex object.</para>
-              </listitem>
-            </varlistentry>
-          </variablelist>
-        </listitem>
-      </varlistentry>
-
-      <varlistentry>
-        <term>[ this , data ] = optimsimplex_axes ( this , x0 , fun , len ,
+        <term>[ this , data ] = optimsimplex_new ( "axes" , x0 , fun , len ,
         data )</term>
 
         <listitem>
-          <para>Configure the current simplex so that it is computed axis by
+          <para>Creates a new simplex object so that it is computed axis by
           axis, with the given length.</para>
 
           <variablelist>
@@ -386,13 +380,23 @@ function [ y , data ] = myfunction ( x , data )
       </varlistentry>
 
       <varlistentry>
-        <term>[ this , data ] = optimsimplex_pfeffer ( this , x0 , fun ,
+        <term>this = optimsimplex_new ( "pfeffer" , x0 )</term>
+
+        <term>this = optimsimplex_new ( "pfeffer" , x0 , fun )</term>
+
+        <term>this = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual
+        )</term>
+
+        <term>this = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual ,
+        deltazero )</term>
+
+        <term>[ this , data ] = optimsimplex_new ( "pfeffer" , x0 , fun ,
         deltausual , deltazero , data )</term>
 
         <listitem>
-          <para>Configure the current simplex so that it is computed from
-          Pfeffer's method, i.e. a relative delta for non-zero values and an
-          absolute delta for zero values.</para>
+          <para>Creates a new simplex so that it is computed from Pfeffer's
+          method, i.e. a relative delta for non-zero values and an absolute
+          delta for zero values.</para>
 
           <variablelist>
             <varlistentry>
@@ -472,13 +476,16 @@ function [ y , data ] = myfunction ( x , data )
       </varlistentry>
 
       <varlistentry>
-        <term>[ this , data ] = optimsimplex_randbounds ( this , x0 , fun ,
+        <term>this = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin ,
+        boundsmax , nbpoints )</term>
+
+        <term>[ this , data ] = optimsimplex_new ( "randbounds" , x0 , fun ,
         boundsmin , boundsmax , nbpoints , data )</term>
 
         <listitem>
-          <para>Configure the current simplex so that it is computed by taking
-          the bounds into account with random scaling. The number of vertices
-          in the simplex is arbitrary.</para>
+          <para>Creates a new simplex so that it is computed by taking the
+          bounds into account with random scaling. The number of vertices in
+          the simplex is arbitrary.</para>
 
           <variablelist>
             <varlistentry>
@@ -572,12 +579,18 @@ function [ y , data ] = myfunction ( x , data )
       </varlistentry>
 
       <varlistentry>
-        <term>[ this , data ] = optimsimplex_spendley ( this , x0 , fun , len
+        <term>this = optimsimplex_new ( "spendley" , x0 )</term>
+
+        <term>this = optimsimplex_new ( "spendley" , x0 , fun )</term>
+
+        <term>this = optimsimplex_new ( "spendley" , x0 , fun , len )</term>
+
+        <term>[ this , data ] = optimsimplex_new ( "spendley" , x0 , fun , len
         , data )</term>
 
         <listitem>
-          <para>Configure the current simplex so that it is computed from
-          Spendley's et al. method, i.e. a regular simplex made of nbve = n+1
+          <para>Creates a new simplex so that it is computed from Spendley's
+          et al. method, i.e. a regular simplex made of nbve = n+1
           vertices.</para>
 
           <variablelist>
@@ -649,6 +662,24 @@ function [ y , data ] = myfunction ( x , data )
       </varlistentry>
 
       <varlistentry>
+        <term>this = optimsimplex_destroy (this)</term>
+
+        <listitem>
+          <para>Destroy the given object.</para>
+
+          <variablelist>
+            <varlistentry>
+              <term>this</term>
+
+              <listitem>
+                <para>The current simplex object.</para>
+              </listitem>
+            </varlistentry>
+          </variablelist>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
         <term>this = optimsimplex_setall ( this , simplex )</term>
 
         <listitem>
@@ -1621,7 +1652,7 @@ function [ y , data ] = myfunction ( x , data )
               <term>method</term>
 
               <listitem>
-                <para>optional, the method to use to compute the size. </para>
+                <para>optional, the method to use to compute the size.</para>
 
                 <para>The available methods are the following :</para>
 
index 827d8ab..f9c728d 100644 (file)
@@ -53,8 +53,6 @@ function [x,fval,exitflag,output] = fminsearch ( varargin )
   end
   TolFun = options.TolFun;
   TolX = options.TolX;
-  // Get options from the options struct
-  options.MaxFunEvals
   // Perform Optimization
   nm = neldermead_new ();
   nm = neldermead_configure(nm,"-x0",x0.');
@@ -72,6 +70,7 @@ function [x,fval,exitflag,output] = fminsearch ( varargin )
   nm = neldermead_configure(nm,"-tolsimplexizemethod","disabled");
   nm = neldermead_configure(nm,"-toldeltafv",TolFun);
   nm = neldermead_configure(nm,"-tolsimplexizeabsolute",TolX);
+  nm = neldermead_configure(nm,"-checkcostfunction",0);
   nm = neldermead_search(nm);
   x = neldermead_get(nm,"-xopt").';
   fval = neldermead_get(nm,"-fopt");
index 1ef16a9..60d8803 100644 (file)
@@ -73,9 +73,12 @@ function value = neldermead_cget (this,key)
     value = this.nbineqloops;
   case "-ineqscaling" then
     value = this.ineqscaling;
+  case "-checkcostfunction" then
+    value = this.checkcostfunction;
+  case "-scalingmethod" then
+    value = this.scalingmethod;
   else
     // Delegate to the optimization object
     value = optimbase_cget ( this.optbase , key );
   end
 endfunction
-
index f84c1fd..2a968a2 100644 (file)
@@ -119,9 +119,20 @@ function this = neldermead_configure (this,key,value)
     this.nbineqloops = value;
   case "-ineqscaling" then
     this.ineqscaling = value;
+  case "-checkcostfunction" then
+    select value
+    case 0 then
+      this.checkcostfunction = value;
+    case 1 then
+      this.checkcostfunction = value;
+    else
+      errmsg = msprintf(gettext("%s: Unknown value %s for -checkcostfunction option"),"neldermead_configure", value);
+      error(errmsg);
+    end
+  case "-scalingmethod" then
+    this.scalingmethod = value;
   else
     // Delegate to the optimization object
     this.optbase = optimbase_configure ( this.optbase , key , value );
   end
 endfunction
-
index 023c417..cbacc96 100644 (file)
@@ -28,7 +28,8 @@ function newobj = neldermead_new ()
     "kelleynormalizationflag","kelleystagnationalpha0", ...
     "kelleyalpha","restartnb","restartflag","restartdetection" , ...
     "startupflag" , ...
-    "boxnbpoints" , "boxnbpointseff" , "nbineqloops" , "ineqscaling" ]);
+    "boxnbpoints" , "boxnbpointseff" , "nbineqloops" , "ineqscaling" , ...
+    "checkcostfunction" , "scalingmethod" ]);
   newobj.optbase = optimbase_new();
   // Possible values "variable", "fixed".
   newobj.method = "variable";
@@ -113,5 +114,8 @@ function newobj = neldermead_new ()
   // The scaling coefficient in nonlinear inequality constraints
   // in Box method, in (0,1) range
   newobj.ineqscaling = 0.5
+  // Set to 0 to disable the checking of the connection of the cost function
+  newobj.checkcostfunction = 1;
+  // The scaling algorithm : "tox0", "tocentroid"
+  newobj.scalingmethod = "tox0";
 endfunction
-
index c5e7516..07eba3b 100644 (file)
@@ -666,20 +666,13 @@ endfunction
 //   Computes the initial simplex, depending on the -simplex0method.
 //
 function this = neldermead_startup (this)
-  // 5. Store initial data into the base optimization component
+  // 0. Check that the cost function is correctly connected
   // Note: this call to the cost function is not used, but helps the
   // user while he is tuning his object.
-  x0 = optimbase_cget ( this.optbase , "-x0" );
-  cmd = "[ this , fx0 ] = neldermead_function ( this , x0 )";
-  ierr=execstr(cmd,"errcatch");
-  if ierr <> 0 then
-    errmsg = msprintf ( gettext ( "%s: Cannot evaluate cost function (use neldermead_function to check your configuration)." ) , "neldermead_startup" )
-    error ( errmsg );
+  checkfun = this.checkcostfunction;
+  if checkfun == 1 then
+    optimbase_checkcostfun ( this.optbase );
   end
-  this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
-  this.optbase = optimbase_set ( this.optbase , "-xopt" , x0.' );
-  this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
-  this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
   // 1. If the problem has bounds, check that they are consistent
   [ this.optbase , hasbounds ] = optimbase_hasbounds ( this.optbase );
   if ( hasbounds ) then
@@ -689,25 +682,22 @@ function this = neldermead_startup (this)
     end
   end
   // 2. Get the initial guess and compute the initial simplex
+  x0 = optimbase_cget ( this.optbase , "-x0" );
   select this.simplex0method
   case "given" then
     [ simplex0 , this ] = optimsimplex_new ( this.coords0 , ...
       neldermead_costf , this );
   case "axes" then
-    simplex0 = optimsimplex_new ( );
-    [ simplex0 , this ] = optimsimplex_axes ( simplex0 , ...
+    [ simplex0 , this ] = optimsimplex_new ( "axes" , ...
       x0.' , neldermead_costf , this.simplex0length , this );
   case "spendley" then
-    simplex0 = optimsimplex_new ( );
-    [ simplex0 , this ] = optimsimplex_spendley ( simplex0 , ...
+    [ simplex0 , this ] = optimsimplex_new ( "spendley" , ...
       x0.' , neldermead_costf , this.simplex0length , this );
   case "pfeffer" then
-    simplex0 = optimsimplex_new ( );
-    [ simplex0 , this ] = optimsimplex_pfeffer ( simplex0 , ...
+    [ simplex0 , this ] = optimsimplex_new ( "pfeffer" , ...
       x0.' , neldermead_costf , this.simplex0deltausual , ...
       this.simplex0deltazero , this );
   case "randbounds" then
-    simplex0 = optimsimplex_new ( );
     if ( this.boxnbpoints == "2n" ) then
       this.boxnbpointseff = 2 * this.optbase.numberofvariables;
     else
@@ -716,7 +706,7 @@ function this = neldermead_startup (this)
     if ( ~hasbounds ) then
       error ( msprintf(gettext("%s: Randomized bounds initial simplex is not available without bounds." ), "neldermead_startup"))
     end
-    [ simplex0 , this ] = optimsimplex_randbounds ( simplex0 , x0.' , ...
+    [ simplex0 , this ] = optimsimplex_new ( "randbounds" , x0.' , ...
       neldermead_costf , this.optbase.boundsmin , this.optbase.boundsmax , ...
       this.boxnbpointseff  , this );
   else
@@ -728,15 +718,91 @@ function this = neldermead_startup (this)
   //
   if ( hasbounds | this.optbase.nbineqconst > 0 ) then
     this = neldermead_log (this,sprintf("Scaling initial simplex into nonlinear inequality constraints..."));
+    select this.scalingmethod
+    case "tox0" then
+      [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 );
+    case "tocenter" then
+      [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 );
+    else
+      errmsg = msprintf(gettext("%s: Unknown value %s for -scalingmethod option"),"neldermead_startup", this.scalingmethod );
+      error(errmsg);
+    end
+  end
+  //
+  // 4. Store the simplex
+  //
+  this.simplex0 = optimsimplex_destroy ( this.simplex0 );
+  this.simplex0 = simplex0;
+  this.simplexsize0 = optimsimplex_size ( simplex0 );
+  // 5. Store initial data into the base optimization component
+  fx0 = optimsimplex_getfv ( this.simplex0 , 1 );
+  this.optbase = optimbase_set ( this.optbase , "-fx0" , fx0 );
+  this.optbase = optimbase_set ( this.optbase , "-xopt" , x0.' );
+  this.optbase = optimbase_set ( this.optbase , "-fopt" , fx0 );
+  this.optbase = optimbase_set ( this.optbase , "-iterations" , 0 );
+  // 6. If Kelley's stagnation is enabled, initialize Kelley's stagnation detection system.
+  if ( this.kelleystagnationflag == 1 ) then
+    this = neldermead_kelleystag ( this );
+  end
+endfunction
+//
+// neldermead_scaletox0 --
+//   Scale the simplex into the bounds and the 
+//   nonlinear inequality constraints, if any.
+//   Scale toward x0, which is feasible.
+// Arguments
+//   
+//
+function [ this , simplex0 ] = neldermead_scaletox0 ( this , simplex0 )
     nbve = optimsimplex_getnbve ( simplex0 );
-    for ive = 1 : nbve
+    xref = optimsimplex_getx ( simplex0 , 1 );
+    for ive = 2 : nbve
       // optimsimplex returns a row vector
       x = optimsimplex_getx ( simplex0 , ive );
       this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
         strcat(string(x)," ")+"]... " , ...
         ive , nbve ));
-      // Transpose x, because x0 is a column vector
-      [ this , status , xp ] = _scaleinconstraints ( this , x.' , x0 );
+      // Transpose x into a row vector
+      [ this , status , xp ] = _scaleinconstraints ( this , x.' , xref );
+      if ( ~status ) then
+        errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
+          "neldermead_startup", ive , nbve , strcat(string(x)," "));
+        error(errmsg);
+      end
+      if ( or ( x <> xp ) ) then
+        [ this , fv ] = neldermead_function ( this , xp );
+        // Transpose xp, which is a column vector
+        simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
+      end
+    end
+endfunction
+//
+// neldermead_scaletocenter --
+//   Scale the simplex into the bounds and the 
+//   nonlinear inequality constraints, if any.
+//   Scale to the centroid of the points
+//   which satisfy the constraints.
+// Notes
+//   This is Box's method for scaling.
+//   It is unsure, since the centroid of the points
+//   which satisfy the constraints may not be feasible.
+// TODO : test !
+// TODO : insert a note for that specific point
+// Arguments
+//   
+//
+function [ this , simplex0 ] = neldermead_scaletocenter ( this , simplex0 , x0 )
+    nbve = optimsimplex_getnbve ( simplex0 );
+    xref = optimsimplex_getx ( simplex0 , 1 );
+    for ive = 2 : nbve
+      xref = optimsimplex_xbar ( simplex0 , ive:nbve );
+      // optimsimplex returns a row vector
+      x = optimsimplex_getx ( simplex0 , ive );
+      this = neldermead_log (this,sprintf("Scaling vertex #%d/%d at ["+...
+        strcat(string(x)," ")+"]... " , ...
+        ive , nbve ));
+      // Transpose x into a row vector
+      [ this , status , xp ] = _scaleinconstraints ( this , x.' , xref );
       if ( ~status ) then
         errmsg = msprintf(gettext("%s: Impossible to scale the vertex #%d/%d at [%s] into inequality constraints"), ...
           "neldermead_startup", ive , nbve , strcat(string(x)," "));
@@ -748,19 +814,7 @@ function this = neldermead_startup (this)
         simplex0 = optimsimplex_setve ( simplex0 , ive , fv , xp.' );
       end
     end
-  end
-  //
-  // 4. Store the simplex
-  //
-  this.simplex0 = optimsimplex_destroy ( this.simplex0 );
-  this.simplex0 = simplex0;
-  this.simplexsize0 = optimsimplex_size ( simplex0 );
-  // 6. If Kelley's stagnation is enabled, initialize Kelley's stagnation detection system.
-  if ( this.kelleystagnationflag == 1 ) then
-    this = neldermead_kelleystag ( this );
-  end
 endfunction
-
 //
 // neldermead_kelleystag --
 //   Initialize Kelley's stagnation detection system when normalization is required,
@@ -1055,5 +1109,4 @@ function this = neldermead_box ( this )
   this.optbase = optimbase_set ( this.optbase , "-fopt" , flow );
   this.optbase = optimbase_set ( this.optbase , "-status" , status );
   this.simplexopt = simplex;
-endfunction
-
+endfunction
\ No newline at end of file
diff --git a/scilab/modules/optimization/macros/optimbase/optimbase_checkcostfun.sci b/scilab/modules/optimization/macros/optimbase/optimbase_checkcostfun.sci
new file mode 100644 (file)
index 0000000..998a215
--- /dev/null
@@ -0,0 +1,76 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+//
+// optimbase_checkcostfun --
+//   Check that the cost function is correctly connected.
+//   Generate an error if there is one.
+//
+function this = optimbase_checkcostfun ( this )
+  //
+  // Check that f(x0) works and is a scalar
+  //
+  cmd = "[ this , fx0 ] = optimbase_function ( this , this.x0 )";
+  ierr=execstr(cmd,"errcatch");
+  if ierr <> 0 then
+    errmsg = msprintf ( gettext ( "%s: Cannot evaluate cost function at x0." ) , "optimbase_checkcostfun" )
+    error ( errmsg );
+  end
+  if ( size(fx0,1) <> 1 ) | ( size(fx0,2) <> 1 ) then
+    errmsg = msprintf ( gettext ( "%s: Call to cost function with x0 is not a scalar, but a %d x %d matrix." ) , "optimbase_checkcostfun" , size(fx0,1) , size(fx0,2) )
+    error ( errmsg );
+  end
+  //
+  // If there are nonlinear constraints, check that the index is correctly managed.
+  //
+  if this.nbineqconst > 0 then
+    // index = 1
+    index = 1;
+    cmd = "[ this , fx0 ] = optimbase_function ( this , this.x0 , index )";
+    ierr=execstr(cmd,"errcatch");
+    if ierr <> 0 then
+      errmsg = msprintf ( gettext ( "%s: Cannot evaluate cost function at (x0,index=1)." ) , "optimbase_checkcostfun" )
+      error ( errmsg );
+    end
+    // index = 2
+    index = 2;
+    cmd = "[ this , result ] = optimbase_function ( this , this.x0 , index )";
+    ierr=execstr(cmd,"errcatch");
+    if ierr <> 0 then
+      errmsg = msprintf ( gettext ( "%s: Cannot evaluate cost function at (x0,index=2)." ) , "optimbase_checkcostfun" )
+      error ( errmsg );
+    end
+    if size(result,1) <> 1 then
+      errmsg = msprintf ( gettext ( "%s: The result of the cost function at (x0,index=2) has %d rows, instead of only 1." ) , "optimbase_checkcostfun" , size(result,1))
+      error ( errmsg );
+    end
+    if ( size(result,2) <> this.nbineqconst ) then
+      errmsg = msprintf ( gettext ( "%s: The result of the cost function at (x0,index=2) has %d columns, instead of the number of constraints %d." ) , "optimbase_checkcostfun" , size(result,2) , this.nbineqconst )
+      error ( errmsg );
+    end
+    // index = 3
+    index = 3;
+    cmd = "[ this , result ] = optimbase_function ( this , this.x0 , index )";
+    ierr=execstr(cmd,"errcatch");
+    if ierr <> 0 then
+      errmsg = msprintf ( gettext ( "%s: Cannot evaluate cost function at (x0,index=3)." ) , "optimbase_checkcostfun" )
+      error ( errmsg );
+    end
+    if size(result,1) <> 1 then
+      errmsg = msprintf ( gettext ( "%s: The result of the cost function at (x0,index=3) has %d rows, instead of only 1." ) , "optimbase_checkcostfun" , size(result,1))
+      error ( errmsg );
+    end
+    if ( size(result,2) <> this.nbineqconst + 1 ) then
+      errmsg = msprintf ( gettext ( "%s: The result of the cost function at (x0,index=3) has %d columns, instead of the number of constraints %d + 1." ) , "optimbase_checkcostfun" , size(result,2) , this.nbineqconst )
+      error ( errmsg );
+    end
+  end    
+endfunction
+
index 4d76705..18b732c 100644 (file)
@@ -39,15 +39,28 @@ function [ this , result ] = optimbase_function ( this , x , index )
     error(errmsg)
   end
   if typeof(this.costfargument)=="string" then
+    // There is no additionnal argument for the cost function
     if (~isdef('index','local')) then
       result = this.fun ( x );
     else
       result = this.fun ( x , index );
     end
   else
+    // There IS one additionnal argument for the cost function
     if (~isdef('index','local')) then
-      result = this.fun ( x , this.costfargument );
+      // The caller did not provide the value of index
+      if ( this.nbineqconst == 0 ) then
+        // There is one additionnal argument, but no nonlinear constraints,
+        // therefore, there is no need for a index value.
+        result = this.fun ( x , this.costfargument );
+      else
+        // Set the index, so that, if an additionnal cost function argument is provided,
+        // it can be appended at the end.
+        index = 1;
+        result = this.fun ( x , index , this.costfargument );
+      end
     else
+      // There is one additionnal argument, and the caller provided the value of index.
       result = this.fun ( x , index , this.costfargument );
     end
   end
@@ -59,3 +72,4 @@ function [ this , result ] = optimbase_function ( this , x , index )
   end
 endfunction
 
+
diff --git a/scilab/modules/optimization/macros/optimbase/optimbase_isinbounds.sci b/scilab/modules/optimization/macros/optimbase/optimbase_isinbounds.sci
new file mode 100644 (file)
index 0000000..0769708
--- /dev/null
@@ -0,0 +1,41 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - INRIA - 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
+
+//
+// optimbase_isinbounds --
+//   Returns %t if the given point satisfies bounds constraints.
+//   Returns %f if the given point is not in the bounds.
+// Arguments
+//   x : the point to analyse
+//   isfeasible : = %t, %f
+//
+function [ this , isfeasible ] = optimbase_isinbounds ( this , x )
+    isfeasible = %t
+      [ this , hasbounds ] = optimbase_hasbounds ( this );
+      if ( hasbounds ) then
+        for ix = 1 : this.numberofvariables
+          xmin = this.boundsmin ( ix )
+          xmax = this.boundsmax ( ix )
+          xix = x ( ix )
+          if ( xix < xmin ) then
+            isfeasible = %f
+            optimbase_log ( this , sprintf ( "Component #%d/%d of x is lower than min bound %e", ...
+              ix , this.numberofvariables , xmin ) )
+            break
+          end
+          if (xix > xmax) then
+            isfeasible = %f
+            optimbase_log ( this , sprintf ( "Component #%d/%d of x is greater than max bound %e", ...
+              ix , this.numberofvariables , xmax ) )
+            break
+          end
+        end
+      end
+endfunction
+
diff --git a/scilab/modules/optimization/macros/optimbase/optimbase_isinnonlincons.sci b/scilab/modules/optimization/macros/optimbase/optimbase_isinnonlincons.sci
new file mode 100644 (file)
index 0000000..f0bd793
--- /dev/null
@@ -0,0 +1,33 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2009 - INRIA - 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
+
+//
+// optimbase_isinnonlincons --
+//   Returns %t if the given point satisfies inequality constraints.
+//   Returns %f if the given point does not satisfies inequality constraints.
+// Arguments
+//   x : the point to analyse
+//   isfeasible : = %t or %f
+//
+function [ this , isfeasible ] = optimbase_isinnonlincons ( this , x )
+    isfeasible = %t
+      if ( this.nbineqconst > 0) then
+        [ this , const ] = optimbase_function ( this , x , 2 );
+        index = 0
+        for ic = 1 : this.nbineqconst
+          if ( const ( ic ) < 0.0 ) then
+            optimbase_log ( this , sprintf ( "Inequality constraint #%d/%d is not satisfied for x", ...
+              ic , this.nbineqconst ) )
+            isfeasible = %f
+            break
+          end
+        end
+      end
+endfunction
+
diff --git a/scilab/modules/optimization/macros/optimsimplex/optimsimplex_axes.sci b/scilab/modules/optimization/macros/optimsimplex/optimsimplex_axes.sci
deleted file mode 100644 (file)
index 391fdfe..0000000
+++ /dev/null
@@ -1,57 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - 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
-
-//
-// optimsimplex_axes --
-//   Configure the current simplex so that it is computed from the axes and a length.
-// Arguments
-//   x0 : the initial point, as a row vector
-//   fun : name of the function
-//   length : the length of the simplex
-//     If length is a value, that unique length
-//     is used in all directions.
-//     If length is a vector with n values, each 
-//     length is used with the corresponding 
-//     direction.
-//   data : user-defined data
-//
-function [ this , data ] = optimsimplex_axes ( this , x0 , fun , len , data )
-  if (~isdef('len','local')) then
-    len = 1.0;
-  end
-  n = length(x0);
-  this.n = n;
-  this.nbve = n + 1;
-  nl=length(len)
-  if nl==1 then
-    xlen = len * ones(n,1)
-  else
-    xlen = len
-  end
-  this.x = zeros ( this.nbve , n )
-  this.fv = zeros ( this.nbve , 1 )
-  //
-  // Set 1st point
-  //
-  this.x ( 1 , 1:n ) = x0 (1:n)
-  //
-  // Set points #2 to #n+1
-  //
-  for j = 2 : this.nbve
-    this.x ( j , 1:n ) = x0 (1:n)
-    this.x ( j , j-1 ) = this.x ( j , j-1 ) + xlen(j-1)
-  end
-  // Compute Function Value
-  if (~isdef('data','local')) then
-    this = optimsimplex_computefv ( this , fun )
-  else
-    [ this , data ] = optimsimplex_computefv ( this , fun , data )
-  end
-endfunction
-
index 24907f6..2e36aa6 100644 (file)
 //   coords : list of point coordinates in the simplex
 //   fun : the function to compute at vertices
 //   data : user-defined data, passed to the function
+// Uses :
+//   newobj = optimsimplex_new ( )
+//   newobj = optimsimplex_new ( coords )
+//   newobj = optimsimplex_new ( coords , fun )
+//   [ newobj , data ] = optimsimplex_new ( coords , fun , data )
+//   newobj = optimsimplex_new ( "axes" , x0 )
+//   newobj = optimsimplex_new ( "axes" , x0 , fun )
+//   newobj = optimsimplex_new ( "axes" , x0 , fun , len )
+//   [ newobj , data ] = optimsimplex_new ( "axes" , x0 , fun , len , data )
+//   newobj = optimsimplex_new ( "spendley" , x0 )
+//   newobj = optimsimplex_new ( "spendley" , x0 , fun )
+//   newobj = optimsimplex_new ( "spendley" , x0 , fun , len )
+//   [ newobj , data ] = optimsimplex_new ( "spendley" , x0 , fun , len , data )
+//   newobj = optimsimplex_new ( "pfeffer" , x0 )
+//   newobj = optimsimplex_new ( "pfeffer" , x0 , fun )
+//   newobj = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual )
+//   newobj = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual , deltazero )
+//   [ newobj , data ] = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual , deltazero , data )
+//   newobj = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax , nbve )
+//   [ newobj , data ] = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax , nbve  , data )
+// 
 //
-function [ newobj , data ] = optimsimplex_new ( coords , fun , data )
+function [ newobj , data ] = optimsimplex_new ( varargin )
+  [lhs,rhs]=argn();
+  if rhs>7 then
+    errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 0 to 7 are expected."), "optimsimplex_new", rhs);
+    error(errmsg)
+  end
+  if rhs == 0 then 
+    //   newobj = optimsimplex_new ( )
+    newobj = optimsimplex_coords ( )
+    return;
+  end
+  var1 = varargin(1);
+  if typeof(var1) == "string" then
+    stype = varargin(1);
+    select stype
+    case "axes" then
+      //   newobj = optimsimplex_new ( "axes" , x0 )
+      //   newobj = optimsimplex_new ( "axes" , x0 , fun )
+      //   newobj = optimsimplex_new ( "axes" , x0 , fun , len )
+      //   [ newobj , data ] = optimsimplex_new ( "axes" , x0 , fun , len , data )
+      if rhs<2 | rhs > 5 then
+        errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 2 to 5 are expected."), "optimsimplex_new", rhs);
+        error(errmsg)
+      end
+      x0   = varargin(2);
+      if rhs==2 then
+        newobj = optimsimplex_axes ( x0 )
+      elseif rhs==3 then
+        fun  = varargin(3);
+        newobj = optimsimplex_axes ( x0 , fun )
+      elseif rhs==4 then
+        fun  = varargin(3);
+        len  = varargin(4);
+        newobj = optimsimplex_axes ( x0 , fun , len )
+      elseif rhs==5 then
+        fun  = varargin(3);
+        len  = varargin(4);
+        data = varargin(5);
+        [ newobj , data ] = optimsimplex_axes ( x0 , fun , len , data )
+      end
+    case "spendley" then
+      //   newobj = optimsimplex_new ( "spendley" , x0 )
+      //   newobj = optimsimplex_new ( "spendley" , x0 , fun )
+      //   newobj = optimsimplex_new ( "spendley" , x0 , fun , len )
+      //   [ newobj , data ] = optimsimplex_new ( "spendley" , x0 , fun , len , data )
+      if rhs<2 | rhs > 5 then
+        errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 2 to 5 are expected."), "optimsimplex_new", rhs);
+        error(errmsg)
+      end
+      if rhs==2 then
+        x0   = varargin(2);
+        newobj = optimsimplex_spendley ( x0 )
+      elseif rhs==3 then
+        x0   = varargin(2);
+        fun  = varargin(3);
+        newobj = optimsimplex_spendley ( x0 , fun )
+      elseif rhs==4 then
+        x0   = varargin(2);
+        fun  = varargin(3);
+        len  = varargin(4);
+        newobj = optimsimplex_spendley ( x0 , fun , len )
+      elseif rhs==5 then
+        x0   = varargin(2);
+        fun  = varargin(3);
+        len  = varargin(4);
+        data = varargin(5);
+        [ newobj , data ] = optimsimplex_spendley ( x0 , fun , len , data )
+      end
+    case "pfeffer" then
+      //   newobj = optimsimplex_new ( "pfeffer" , x0 )
+      //   newobj = optimsimplex_new ( "pfeffer" , x0 , fun )
+      //   newobj = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual )
+      //   newobj = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual , deltazero )
+      //   [ newobj , data ] = optimsimplex_new ( "pfeffer" , x0 , fun , deltausual , deltazero , data )
+      if rhs<2 | rhs > 6 then
+        errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 2 to 6 are expected."), "optimsimplex_new", rhs);
+        error(errmsg)
+      end
+      x0   = varargin(2);
+      if rhs==2 then
+        newobj = optimsimplex_pfeffer ( x0 )
+      elseif rhs==3 then
+        fun  = varargin(3);
+        newobj = optimsimplex_pfeffer ( x0 , fun )
+      elseif rhs==4 then
+        fun         = varargin(3);
+        deltausual  = varargin(4);
+        newobj = optimsimplex_pfeffer ( x0 , fun , deltausual )
+      elseif rhs==5 then
+        fun        = varargin(3);
+        deltausual = varargin(4);
+        deltazero  = varargin(5);
+        newobj = optimsimplex_pfeffer ( x0 , fun , deltausual , deltazero )
+      elseif rhs==6 then
+        fun        = varargin(3);
+        deltausual = varargin(4);
+        deltazero  = varargin(5);
+        data       = varargin(6);
+        [ newobj , data ] = optimsimplex_pfeffer ( x0 , fun , deltausual , deltazero , data )
+      end
+    case "randbounds" then
+      //   newobj = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax )
+      //   newobj = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax , nbve )
+      //   [ newobj , data ] = optimsimplex_new ( "randbounds" , x0 , fun , boundsmin , boundsmax , nbve  , data )
+      if rhs<5 | rhs > 7 then
+        errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 6 to 7 are expected."), "optimsimplex_new", rhs);
+        error(errmsg)
+      end
+      x0   = varargin(2);
+      fun  = varargin(3);
+      boundsmin = varargin(4);
+      boundsmax = varargin(5);
+      if rhs==5 then
+        newobj = optimsimplex_randbounds ( x0 , fun , boundsmin , boundsmax )
+      elseif rhs==6 then
+        nbve = varargin(6);
+        newobj = optimsimplex_randbounds ( x0 , fun , boundsmin , boundsmax , nbve )
+      elseif rhs==7 then
+        nbve = varargin(6);
+        data       = varargin(7);
+        [ newobj , data ] = optimsimplex_randbounds ( x0 , fun , boundsmin , boundsmax , nbve , data )
+      end
+    else
+      errmsg = msprintf(gettext("%s: Unknown key %s"),"optimsimplex_new",key)
+      error(errmsg)
+    end
+  else
+    //   newobj = optimsimplex_new ( coords )
+    //   newobj = optimsimplex_new ( coords , fun )
+    //   [ newobj , data ] = optimsimplex_new ( coords , fun , data )
+    if rhs < 1 | rhs > 3 then
+      errmsg = msprintf(gettext("%s: Unexpected number of input arguments : %d provided while 1 to 3 are expected."), "optimsimplex_new", rhs);
+      error(errmsg)
+    end
+    if rhs == 1
+      coords = varargin(1);
+      newobj = optimsimplex_coords ( coords )
+    elseif rhs == 2
+      coords = varargin(1);
+      fun = varargin(2);
+      newobj = optimsimplex_coords ( coords , fun )
+    elseif rhs == 3
+      coords = varargin(1);
+      fun = varargin(2);
+      data = varargin(3);
+      [ newobj , data ] = optimsimplex_coords ( coords , fun , data )
+    end
+  end    
+endfunction
+
+function [ newobj , data ] = optimsimplex_coords ( coords , fun , data )
   if (~isdef('coords','local')) then
     coords = [];
   end
-   if (~isdef('fun','local')) then
+  if (~isdef('fun','local')) then
     fun = [];
   end
   newobj = tlist(["T_SIMPLEX",...
@@ -67,4 +238,232 @@ function [ newobj , data ] = optimsimplex_new ( coords , fun , data )
     end
   end
 endfunction
+//
+// optimsimplex_axes --
+//   Configure the current simplex so that it is computed from the axes and a length.
+// Arguments
+//   x0 : the initial point, as a row vector
+//   fun : name of the function
+//   length : the length of the simplex
+//     If length is a value, that unique length
+//     is used in all directions.
+//     If length is a vector with n values, each 
+//     length is used with the corresponding 
+//     direction.
+//   data : user-defined data
+//
+function [ newobj , data ] = optimsimplex_axes ( x0 , fun , len , data )
+  newobj = optimsimplex_coords ( )
+  if size(x0,1)<>1 then
+    errmsg = msprintf(gettext("%s: The x0 vector is expected to be a row matrix, but current shape is %d x %d"),"optimsimplex_axes",size(x0,1),size(x0,2));
+    error(errmsg);
+  end
+  if (~isdef('fun','local')) then
+    fun = [];
+  end
+  if (~isdef('len','local')) then
+    len = 1.0;
+  end
+  n = length(x0);
+  newobj.n = n;
+  newobj.nbve = n + 1;
+  nl=length(len)
+  if nl==1 then
+    xlen = len * ones(n,1)
+  else
+    xlen = len
+  end
+  newobj.x = zeros ( newobj.nbve , n )
+  newobj.fv = zeros ( newobj.nbve , 1 )
+  //
+  // Set 1st point
+  //
+  newobj.x ( 1 , 1:n ) = x0 (1:n)
+  //
+  // Set points #2 to #n+1
+  //
+  for j = 2 : newobj.nbve
+    newobj.x ( j , 1:n ) = x0 (1:n)
+    newobj.x ( j , j-1 ) = newobj.x ( j , j-1 ) + xlen(j-1)
+  end
+  // Compute Function Value
+  if fun <> [] then
+    if (~isdef('data','local')) then
+      newobj = optimsimplex_computefv ( newobj , fun )
+    else
+      [ newobj , data ] = optimsimplex_computefv ( newobj , fun , data )
+    end
+  end
+endfunction
+
+//
+// optimsimplex_spendley --
+//   Configure the current simplex so that it is computed from Spendley's method, 
+//   i.e. a regular simplex made of k = n+1 vertices.
+// Arguments
+//   x0 : the initial point, as a row vector
+//   fun : name of the function
+//   len : the length of the simplex
+//   data : user-defined data
+//
+function [ newobj , data ] = optimsimplex_spendley ( x0 , fun , len , data )
+  newobj = optimsimplex_coords ( )
+  if size(x0,1)<>1 then
+    errmsg = msprintf(gettext("%s: The x0 vector is expected to be a row matrix, but current shape is %d x %d"),"optimsimplex_spendley",size(x0,1),size(x0,2));
+    error(errmsg);
+  end
+  if (~isdef('fun','local')) then
+    fun = [];
+  end
+  if (~isdef('len','local')) then
+   len = 1.0;
+  end
+  n = length(x0);
+  newobj.n = n;
+  newobj.nbve = n + 1;
+  newobj.x = zeros ( n+1  , n )
+  newobj.fv = zeros ( n+1 , 1 )
+  //
+  // Compute p (diagonal term) , q (off-diagonal term)
+  //
+  p  = (n - 1.0 + sqrt(n + 1))/(n * sqrt(2.0))
+  q = (sqrt(n + 1) - 1.0)/(n * sqrt(2.0))
+  //
+  // Set 1st point
+  //
+  newobj.x ( 1 , 1:n ) = x0 (1:n)
+  //
+  // Set points #2 to #n+1
+  //
+  for j = 2 : newobj.n+1
+    // Note : Vectorize when possible
+    // In order to vectorize, add q everywhere, then substract p just for one diagonal term j-1
+    newobj.x ( j , 1:n ) = x0 (1:n) + len * q * ones(1,n)
+    newobj.x ( j , j-1 ) = newobj.x ( j , j-1 ) - len * q + len * p
+  end
+  // Compute Function Value
+  if fun <> [] then
+    if (~isdef('data','local')) then
+      newobj = optimsimplex_computefv ( newobj , fun )
+    else
+      [ newobj , data ] = optimsimplex_computefv ( newobj , fun , data )
+    end
+  end
+endfunction
+
+//
+// optimsimplex_pfeffer --
+//   Configure the current simplex so that it is computed from Pfeffer's method, 
+//   i.e. a relative delta for non-zero values and an absolute delta 
+//   for zero values.
+// Arguments
+//   x0 : the initial point, as a row vector
+//   fun : name of the function
+//   deltausual : the absolute delta for non-zero values
+//   deltazero : the absolute delta for zero values
+//   data : user-defined data
+// References
+//   "Global Optimization Of Lennard-Jones Atomic Clusters"
+//   Ellen Fan, Thesis, February 26, 2002, McMaster University
+//   Method due to L. Pfeffer at Stanford
+//
+function [ newobj , data ] = optimsimplex_pfeffer ( x0 , fun , deltausual , deltazero , data )
+  newobj = optimsimplex_coords ( )
+  if size(x0,1)<>1 then
+    errmsg = msprintf(gettext("%s: The x0 vector is expected to be a row matrix, but current shape is %d x %d"),"optimsimplex_pfeffer",size(x0,1),size(x0,2));
+    error(errmsg);
+  end
+  if (~isdef('fun','local')) then
+    fun = [];
+  end
+   if (~isdef('deltausual','local')) then
+    deltausual = 0.05
+  end
+   if (~isdef('deltazero','local')) then
+    deltazero = 0.0075
+  end
+  if size(x0,1)<>1 then
+    errmsg = msprintf(gettext("%s: The x0 vector is expected to be a row matrix, but current shape is %d x %d"),"optimsimplex_pfeffer",size(x0,1),size(x0,2));
+    error(errmsg);
+  end
+  n = length(x0);
+  newobj.n = n;
+  newobj.nbve = n + 1;
+  newobj.x = zeros ( n+1 , n )
+  newobj.fv = zeros ( n+1 , 1 )
+  //
+  // Set 1st point
+  //
+  newobj.x ( 1 , 1:n ) = x0 (1,1:n)
+  //
+  // Set points #2 to #n+1
+  //
+  for j = 2 : newobj.n+1
+    newobj.x ( j,1:n ) = x0 (1:n)
+    if x0( j-1 ) == 0.0 then
+      newobj.x ( j , j-1 ) = deltazero
+    else
+      newobj.x ( j , j-1 ) = newobj.x ( j , j-1 ) + deltausual * x0( j-1 )
+    end
+  end
+  // Compute Function Value
+  if fun <> [] then
+    if (~isdef('data','local')) then
+      newobj = optimsimplex_computefv ( newobj , fun )
+    else
+      [ newobj , data ] = optimsimplex_computefv ( newobj , fun , data )
+    end
+  end
+endfunction
+
+//
+// optimsimplex_randbounds --
+//   Configure the current simplex so that it is computed by taking the bounds
+//   into account with random scaling.
+// Arguments
+//   x0 : the initial point
+//   nbpts : the number of vertices in the simplex (+ the initial point)
+//   fun : name of the function
+//   boundsmin : array of minimum bounds
+//   boundsmax : array of maximum bounds
+//   nbve : total number of vertices in the simplex
+//   data : user-defined data
+//
+function [ newobj , data ] = optimsimplex_randbounds ( x0 , fun , boundsmin , boundsmax , nbve  , data )
+  newobj = optimsimplex_coords ( )
+  if size(x0,1)<>1 then
+    errmsg = msprintf(gettext("%s: The x0 vector is expected to be a row matrix, but current shape is %d x %d"),"optimsimplex_randbounds",size(x0,1),size(x0,2));
+    error(errmsg);
+  end
+    n = length ( x0 )
+    if (~isdef('nbve','local')) then
+      nbve = n + 1;
+    end
+    newobj.n = n;
+    newobj.nbve = nbve;
+    newobj.x = zeros ( nbve , n )
+    newobj.fv = zeros ( nbve , 1 )
+    //
+    // Set 1st point
+    //
+    newobj.x ( 1 , 1:n ) = x0 (1:n)
+    //
+    // Set points #2 to #nbve, by randomizing the bounds
+    //
+    for ive = 2 : nbve 
+      //
+      // Compute vertex coordinates, as a random number 
+      // between min and max bounds.
+      //
+      for ix  = 1 : n
+        newobj.x ( ive , ix ) = boundsmin( ix ) + rand() * (boundsmax( ix ) - boundsmin( ix ))
+      end
+    end
+  // Compute Function Value
+    if (~isdef('data','local')) then
+      newobj = optimsimplex_computefv ( newobj , fun )
+    else
+      [ newobj , data ] = optimsimplex_computefv ( newobj , fun , data )
+    end
+endfunction
 
diff --git a/scilab/modules/optimization/macros/optimsimplex/optimsimplex_pfeffer.sci b/scilab/modules/optimization/macros/optimsimplex/optimsimplex_pfeffer.sci
deleted file mode 100644 (file)
index 4852d5d..0000000
+++ /dev/null
@@ -1,60 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - 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
-
-//
-// optimsimplex_pfeffer --
-//   Configure the current simplex so that it is computed from Pfeffer's method, 
-//   i.e. a relative delta for non-zero values and an absolute delta 
-//   for zero values.
-// Arguments
-//   x0 : the initial point, as a row vector
-//   fun : name of the function
-//   deltausual : the absolute delta for non-zero values
-//   deltazero : the absolute delta for zero values
-//   data : user-defined data
-// References
-//   "Global Optimization Of Lennard-Jones Atomic Clusters"
-//   Ellen Fan, Thesis, February 26, 2002, McMaster University
-//   Method due to L. Pfeffer at Stanford
-//
-function [ this , data ] = optimsimplex_pfeffer ( this , x0 , fun , deltausual , deltazero , data )
-   if (~isdef('deltausual','local')) then
-    deltausual = 0.05
-  end
-   if (~isdef('deltazero','local')) then
-    deltazero = 0.0075
-  end
-  n = length(x0);
-  this.n = n;
-  this.nbve = n + 1;
-  this.x = zeros ( n+1 , n )
-  this.fv = zeros ( n+1 , 1 )
-  //
-  // Set 1st point
-  //
-  this.x ( 1 , 1:n ) = x0 (1:n)
-  //
-  // Set points #2 to #n+1
-  //
-  for j = 2 : this.n+1
-    this.x ( j,1:n ) = x0 (1:n)
-    if x0( j-1 ) == 0.0 then
-      this.x ( j , j-1 ) = deltazero
-    else
-      this.x ( j , j-1 ) = this.x ( j , j-1 ) + deltausual * x0( j-1 )
-    end
-  end
-  // Compute Function Value
-  if (~isdef('data','local')) then
-    this = optimsimplex_computefv ( this , fun )
-  else
-    [ this , data ] = optimsimplex_computefv ( this , fun , data )
-  end
-endfunction
-
diff --git a/scilab/modules/optimization/macros/optimsimplex/optimsimplex_randbounds.sci b/scilab/modules/optimization/macros/optimsimplex/optimsimplex_randbounds.sci
deleted file mode 100644 (file)
index 58c6a7d..0000000
+++ /dev/null
@@ -1,53 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - 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
-
-
-//
-// optimsimplex_randbounds --
-//   Configure the current simplex so that it is computed by taking the bounds
-//   into account with random scaling.
-// Arguments
-//   x0 : the initial point
-//   nbpts : the number of vertices in the simplex (+ the initial point)
-//   fun : name of the function
-//   boundsmin : array of minimum bounds
-//   boundsmax : array of maximum bounds
-//   nbve : total number of vertices in the simplex
-//   data : user-defined data
-//
-function [ this , data ] = optimsimplex_randbounds ( this , x0 , fun , boundsmin , boundsmax , nbve  , data )
-    n = length ( x0 )
-    this.n = n;
-    this.nbve = nbve;
-    this.x = zeros ( nbve , n )
-    this.fv = zeros ( nbve , 1 )
-    //
-    // Set 1st point
-    //
-    this.x ( 1 , 1:n ) = x0 (1:n)
-    //
-    // Set points #2 to #nbve, by randomizing the bounds
-    //
-    for ive = 2 : nbve 
-      //
-      // Compute vertex coordinates, as a random number 
-      // between min and max bounds.
-      //
-      for ix  = 1 : n
-        this.x ( ive , ix ) = boundsmin( ix ) + rand() * (boundsmax( ix ) - boundsmin( ix ))
-      end
-    end
-  // Compute Function Value
-  if (~isdef('data','local')) then
-    this = optimsimplex_computefv ( this , fun )
-  else
-    [ this , data ] = optimsimplex_computefv ( this , fun , data )
-  end
-endfunction
-
diff --git a/scilab/modules/optimization/macros/optimsimplex/optimsimplex_spendley.sci b/scilab/modules/optimization/macros/optimsimplex/optimsimplex_spendley.sci
deleted file mode 100644 (file)
index b85d7de..0000000
+++ /dev/null
@@ -1,54 +0,0 @@
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
-// Copyright (C) 2008-2009 - INRIA - 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
-
-//
-// optimsimplex_spendley --
-//   Configure the current simplex so that it is computed from Spendley's method, 
-//   i.e. a regular simplex made of k = n+1 vertices.
-// Arguments
-//   x0 : the initial point, as a row vector
-//   fun : name of the function
-//   len : the length of the simplex
-//   data : user-defined data
-//
-function [ this , data ] = optimsimplex_spendley ( this , x0 , fun , len , data )
-  if (~isdef('len','local')) then
-   len = 1.0;
-  end
-  n = length(x0);
-  this.n = n;
-  this.nbve = n + 1;
-  this.x = zeros ( n+1  , n )
-  this.fv = zeros ( n+1 , 1 )
-  //
-  // Compute p (diagonal term) , q (off-diagonal term)
-  //
-  p  = (n - 1.0 + sqrt(n + 1))/(n * sqrt(2.0))
-  q = (sqrt(n + 1) - 1.0)/(n * sqrt(2.0))
-  //
-  // Set 1st point
-  //
-  this.x ( 1 , 1:n ) = x0 (1:n)
-  //
-  // Set points #2 to #n+1
-  //
-  for j = 2 : this.n+1
-    // Note : Vectorize when possible
-    // In order to vectorize, add q everywhere, then substract p just for one diagonal term j-1
-    this.x ( j , 1:n ) = x0 (1:n) + len * q * ones(1,n)
-    this.x ( j , j-1 ) = this.x ( j , j-1 ) - len * q + len * p
-  end
-  // Compute Function Value
-  if (~isdef('data','local')) then
-    this = optimsimplex_computefv ( this , fun )
-  else
-    [ this , data ] = optimsimplex_computefv ( this , fun , data )
-  end
-endfunction
-
index 16dffd3..4e440a7 100644 (file)
@@ -22,9 +22,14 @@ function cen = optimsimplex_xbar ( this , iexcl )
    if (~isdef('iexcl','local')) then
     iexcl = this.nbve;
   end
+  if ( size(iexcl,1)<>1 ) then
+    errmsg = msprintf(gettext("%s: The exclusion index vector has %d instead of 1."), ...
+      "optimsimplex_xbar", size(iexcl,1) );
+    error(errmsg);
+  end
   // Vectorize by making the sum of all, substracting only one vector
   cen = sum(this.x(1:this.nbve,1:this.n),'r')
-  cen = cen - this.x(iexcl,1:this.n)
-  cen = cen / ( this.nbve - 1)
+  cen = cen - sum(this.x(iexcl,1:this.n),'r')
+  nexcl = size(iexcl,2)
+  cen = cen / ( this.nbve - nexcl ) 
 endfunction
-
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.dia.ref b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.dia.ref
new file mode 100644 (file)
index 0000000..5beb040
--- /dev/null
@@ -0,0 +1,297 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+//
+// 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 bugmes();quit;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 bugmes();quit;end
+endfunction
+//
+// Here, the cost function is OK
+//
+function y = rosenbrockOk (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrockOk);
+nm = optimbase_checkcostfun(nm);
+nm = optimbase_destroy(nm);
+//
+// Here, the cost function is not callable
+//
+function y = rosenbrock2 (x)
+  y = fdsmklqfjdsf;
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrock2);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: Cannot evaluate cost function at x0.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+//
+// Here, the cost function is callable, but returns a matrix,
+// instead of a scalar.
+//
+function y = rosenbrock2 (x)
+  y = ones(10,10);
+endfunction
+Attention: redéfinition de la fonction: rosenbrock2             . Utilisez funcprot(0) pour éviter ce message
+
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrock2);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: Call to cost function with x0 is not a scalar, but a 10 x 10 matrix.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+//
+// Test with good non linear constraints
+//
+function result = optimtestcase ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+nm = optimbase_checkcostfun(nm);
+nm = optimbase_destroy(nm);
+//
+// Test with wrong  non linear constraints f(x0,2) is not a row vector
+//
+function result = optimtestcase2 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3].'
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase2);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=2) has 3 rows, instead of only 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+//
+// Test with wrong  non linear constraints f(x0,2) is a row vector with 5 components instead of 3
+//
+function result = optimtestcase3 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3 0.0 0.0]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase3);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=2) has 5 columns, instead of the number of constraints 3.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+//
+// Test with wrong  non linear constraints f(x0,3) is a column vector
+//
+function result = optimtestcase4 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3].'
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase4);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=3) has 4 rows, instead of only 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+//
+// Test with wrong  non linear constraints f(x0,3) is a row vector with 5 columns instead of 4
+//
+function result = optimtestcase5 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3 0.0]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase5);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=3) has 5 columns, instead of the number of constraints 3 + 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.tst b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_checkcostfun.tst
new file mode 100644 (file)
index 0000000..0d31725
--- /dev/null
@@ -0,0 +1,306 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+//
+// 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
+
+//
+// Here, the cost function is OK
+//
+function y = rosenbrockOk (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrockOk);
+nm = optimbase_checkcostfun(nm);
+nm = optimbase_destroy(nm);
+
+//
+// Here, the cost function is not callable
+//
+function y = rosenbrock2 (x)
+  y = fdsmklqfjdsf;
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrock2);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: Cannot evaluate cost function at x0.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
+//
+// Here, the cost function is callable, but returns a matrix,
+// instead of a scalar.
+//
+function y = rosenbrock2 (x)
+  y = ones(10,10);
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",2);
+nm = optimbase_configure(nm,"-x0",[1.1 1.1]');
+nm = optimbase_configure(nm,"-function",rosenbrock2);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: Call to cost function with x0 is not a scalar, but a 10 x 10 matrix.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
+//
+// Test with good non linear constraints
+//
+function result = optimtestcase ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+nm = optimbase_checkcostfun(nm);
+nm = optimbase_destroy(nm);
+
+//
+// Test with wrong  non linear constraints f(x0,2) is not a row vector
+//
+function result = optimtestcase2 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3].'
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase2);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=2) has 3 rows, instead of only 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
+//
+// Test with wrong  non linear constraints f(x0,2) is a row vector with 5 components instead of 3
+//
+function result = optimtestcase3 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3 0.0 0.0]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase3);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=2) has 5 columns, instead of the number of constraints 3.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
+//
+// Test with wrong  non linear constraints f(x0,3) is a column vector
+//
+function result = optimtestcase4 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3].'
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase4);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=3) has 4 rows, instead of only 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
+//
+// Test with wrong  non linear constraints f(x0,3) is a row vector with 5 columns instead of 4
+//
+function result = optimtestcase5 ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3 0.0]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+nm = optimbase_new ();
+nm = optimbase_configure(nm,"-numberofvariables",4);
+nm = optimbase_configure(nm,"-function",optimtestcase5);
+nm = optimbase_configure(nm,"-x0",[0.0 0.0 0.0 0.0]');
+nm = optimbase_configure(nm,"-nbineqconst",3);
+cmd = "nm = optimbase_checkcostfun(nm);";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_checkcostfun: The result of the cost function at (x0,index=3) has 5 columns, instead of the number of constraints 3 + 1.";
+assert_equal ( computed , expected );
+nm = optimbase_destroy(nm);
+
index 4198737..82d22e8 100644 (file)
@@ -48,13 +48,13 @@ endfunction
 // Arguments
 //    x : the point where to compute the cost
 //    index : a flag which states what is to compute
-//    * if index=1, or no index, returns the value of the cost
+//    * if index=1, or no index, returns the value of the cost 
 //      function (default case)
-//    * if index=2, returns the value of the nonlinear inequality
+//    * if index=2, returns the value of the nonlinear inequality 
 //      constraints, as a row array
 //    * if index=3, returns an array which contains
-//      at index #0, the value of the cost function
-//      at index #1 to the end is the list of the values of the nonlinear
+//      at index #0, the value of the cost function  
+//      at index #1 to the end is the list of the values of the nonlinear 
 //      constraints
 //  Discussion:
 //    The problem is to minimize a cost function with 4 non linear constraints.
@@ -112,7 +112,9 @@ opt = optimbase_configure ( opt , "-numberofvariables",2);
 opt = optimbase_configure ( opt , "-verbose",1);
 [ opt , isok ] = optimbase_checkx0 ( opt );
 Checking initial guess...
+
 ... initial guess is feasible.
+
 assert_equal ( isok , %T );
 opt = optimbase_destroy(opt);
 //
@@ -125,13 +127,18 @@ opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
 opt = optimbase_configure ( opt , "-x0", [1.0 1.0]' );
 [ opt , isok ] = optimbase_checkx0 ( opt );
 Checking initial guess...
+
 ... initial guess is feasible.
+
 assert_equal ( isok , %T );
 opt = optimbase_configure ( opt , "-x0",[-6.0 1.0]');
 [ opt , isok ] = optimbase_checkx0 ( opt );
 Checking initial guess...
+
 Component #1/2 of x is lower than min bound -5.000000e+000
+
 ... initial guess is not feasible.
+
 assert_equal ( isok , %F );
 opt = optimbase_destroy(opt);
 //
@@ -144,16 +151,25 @@ opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
 opt = optimbase_configure ( opt , "-x0" , [ 14.0950013 , 0.8429636 ]');
 [ opt , isok ] = optimbase_checkx0 ( opt );
 Checking initial guess...
+
 Computed constraints = 1.095001e+000 2.779266e-007 2.322073e-006 8.429636e-001
+
 Function Evaluation #1 is [1.0950013 0.0000003 0.0000023 0.8429636] at [14.095001 0.8429636]
+
 ... initial guess is feasible.
+
 assert_equal ( isok , %T );
 opt = optimbase_configure ( opt , "-x0" , [ 14.0950013 , 0.0 ]');
 [ opt , isok ] = optimbase_checkx0 ( opt );
 Checking initial guess...
+
 Computed constraints = 1.095001e+000 7.719049e+000 -7.719046e+000 0.000000e+000
+
 Function Evaluation #1 is [1.0950013 7.7190486 -7.719046 0] at [14.095001 0]
+
 Inequality constraint #3/4 is not satisfied for x
+
 ... initial guess is not feasible.
+
 assert_equal ( isok , %F );
 opt = optimbase_destroy(opt);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.dia.ref b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.dia.ref
new file mode 100644 (file)
index 0000000..3596aeb
--- /dev/null
@@ -0,0 +1,272 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+//
+// 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 bugmes();quit;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 bugmes();quit;end
+endfunction
+//
+// Test the call to the cost function
+//
+// Test when there is no cost function
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+cmd = "[this,f] = optimbase_function ( opt , [0.0 0.0] );";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_function: Empty function (use -function option).";
+assert_equal ( computed , expected );
+opt = optimbase_destroy(opt);
+// Test simple case
+function y = rosenbrock (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-function",rosenbrock);
+[this,f] = optimbase_function ( opt , [0.0 0.0] );
+assert_close ( f , 1.0 , %eps );
+opt = optimbase_destroy(opt);
+//
+// Test with an additional argument
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function y = rosenbrock2 ( x , mydata )
+  a = mydata.a
+  y = 100*(x(2)-x(1)^2)^2 + ( a - x(1))^2;
+endfunction
+mystuff = tlist(["T_MYSTUFF","a"]);
+mystuff.a = 12.0;
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-function",rosenbrock2);
+opt = optimbase_configure(opt,"-costfargument",mystuff);
+[this,f] = optimbase_function ( opt , [0.0 0.0] );
+assert_close ( f , 144. , %eps );
+opt = optimbase_destroy(opt);
+//
+// Test with non linear constraints : there is an index
+// optimtestcase --
+//   Non linear inequality constraints are positive.
+//    
+// Arguments
+//   x: the point where to compute the function
+//   index : the stuff to compute
+// Note
+//  The following protocol is used
+//  * if index=1, or no index, returns the value of the cost 
+//    function (default case)
+//  * if index=2, returns the value of the nonlinear inequality 
+//    constraints, as a row array
+//  * if index=3, returns an array which contains
+//    at index #0, the value of the cost function  
+//    at index #1 to the end is the list of the values of the nonlinear 
+//    constraints
+//  The inequality constraints are expected to be positive.
+//
+function result = optimtestcase ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+x0 = [0.0 0.0 0.0 0.0].';
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",4);
+opt= optimbase_configure(opt,"-function",optimtestcase);
+[this,f] = optimbase_function ( opt , x0 );
+assert_close ( f , 0. , %eps );
+[this,f] = optimbase_function ( opt , x0 , 1 );
+assert_close ( f , 0. , %eps );
+[this,f] = optimbase_function ( opt , x0 , 2 );
+assert_close ( f , [8.    10.    5.] , %eps );
+[this,f] = optimbase_function ( opt , x0 , 3 );
+assert_close ( f , [0.    8.    10.    5.] , %eps );
+opt = optimbase_destroy(opt);
+//
+// Test when there are both nonlinear constraints and a customized data
+// boxproblemA --
+//   Computes the Box problem A cost function and 
+//   inequality constraints.
+//
+// Arguments
+//   x: the point where to compute the function
+//   index : the stuff to compute
+//   data : the parameters of Box cost function
+// Note
+//  The following protocol is used
+//  * if index=1, or no index, returns the value of the cost 
+//    function (default case)
+//  * if index=2, returns the value of the nonlinear inequality 
+//    constraints, as a row array
+//  * if index=3, returns an array which contains
+//    at index #0, the value of the cost function  
+//    at index #1 to the end is the list of the values of the nonlinear 
+//    constraints
+//  The inequality constraints are expected to be positive.
+//
+function result = boxproblemA ( x , index , data )
+  b = x(2) + 0.01 * x(3)
+  x6 = (data.k1 + data.k2 * x(2) ...
+            + data.k3 * x(3) + data.k4 * x(4) + data.k5 * x(5)) * x(1)
+  y1 = data.k6 + data.k7 * x(2) + data.k8 * x(3) ...
+            + data.k9 * x(4) + data.k10 * x(5)
+  y2 = data.k11 + data.k12 * x(2) + data.k13 * x(3) ...
+            + data.k14 * x(4) + data.k15 * x(5)
+  y3 = data.k16 + data.k17 * x(2) + data.k18 * x(3) ...
+            + data.k19 * x(4) + data.k20 * x(5)
+  y4 = data.k21 + data.k22 * x(2) + data.k23 * x(3) ...
+            + data.k24 * x(4) + data.k25 * x(5)
+  x7 = ( y1 + y2 + y3 ) * x(1)
+  x8 = (data.k26 + data.k27 * x(2) + data.k28 * x(3) ...
+            + data.k29 * x(4) ...
+            + data.k30 * x(5) ) * x(1) + x6 + x7
+  if ( index==1 | index==3 ) then
+    f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...
+             + 7840 * data.a6 - 100000 * data.a0 ...
+             - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
+             + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...
+             - 24345 + data.a1 * x6
+    f = -f
+  end
+  if ( index==2 | index==3 ) then
+      c1 = x6
+      c2 = 294000 - x6
+      c3 = x7
+      c4 = 294000 - x7
+      c5 = x8
+      c6 = 277200 - x8
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3 c4 c5 c6]
+  case 3 then
+    result = [f c1 c2 c3 c4 c5 c6]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+boxparams = struct();
+boxparams.a0 = 9;
+boxparams.a1 = 15;
+boxparams.a2 = 50;
+boxparams.a3 = 9.583;
+boxparams.a4 = 20;
+boxparams.a5 = 15;
+boxparams.a6 = 6;
+boxparams.a7 = 0.75;
+boxparams.k1 = -145421.402;
+boxparams.k2 = 2931.1506;
+boxparams.k3 = -40.427932;
+boxparams.k4 = 5106.192;
+boxparams.k5 = 15711.36;
+boxparams.k6 = -161622.577;
+boxparams.k7 = 4176.15328;
+boxparams.k8 = 2.8260078;
+boxparams.k9 = 9200.476;
+boxparams.k10 = 13160.295;
+boxparams.k11 = -21686.9194;
+boxparams.k12 = 123.56928;
+boxparams.k13 = -21.1188894;
+boxparams.k14 = 706.834;
+boxparams.k15 = 2898.573;
+boxparams.k16 = 28298.388;
+boxparams.k17 = 60.81096;
+boxparams.k18 = 31.242116;
+boxparams.k19 = 329.574;
+boxparams.k20 = -2882.082;
+boxparams.k21 = 74095.3845;
+boxparams.k22 = -306.262544;
+boxparams.k23 = 16.243649;
+boxparams.k24 = -3094.252;
+boxparams.k25 = -5566.2628;
+boxparams.k26 = -26237.0;
+boxparams.k27 = 99.0;
+boxparams.k28 = -0.42;
+boxparams.k29 = 1300.0;
+boxparams.k30 = 2100.0;
+boxparams.k31 = 925548.252;
+boxparams.k32 = -61968.8432;
+boxparams.k33 = 23.3088196;
+boxparams.k34 = -27097.648;
+boxparams.k35 = -50843.766;
+x0 = [2.52 2.0 37.5 9.25 6.8].';
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",5);
+opt= optimbase_configure(opt,"-function",boxproblemA);
+opt = optimbase_configure(opt,"-costfargument",boxparams);
+opt = optimbase_configure(opt,"-nbineqconst",6);
+[this,f] = optimbase_function ( opt , x0 );
+assert_close ( f , -2351243.5  , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 1 );
+assert_close ( f , -2351243.5  , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 2 );
+assert_close ( f , [32745.827    261254.17    96991.969    197008.03    130368.43    146831.57] , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 3 );
+assert_close ( f , [-2351243.5    32745.827    261254.17    96991.969    197008.03    130368.43    146831.57] , 1.e-7 );
+opt = optimbase_destroy(opt);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.tst b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_function.tst
new file mode 100644 (file)
index 0000000..b3faabc
--- /dev/null
@@ -0,0 +1,290 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+//
+// 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 the call to the cost function
+//
+
+
+// Test when there is no cost function
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+cmd = "[this,f] = optimbase_function ( opt , [0.0 0.0] );";
+execstr(cmd,"errcatch");
+computed = lasterror();
+expected = "optimbase_function: Empty function (use -function option).";
+assert_equal ( computed , expected );
+opt = optimbase_destroy(opt);
+
+// Test simple case
+function y = rosenbrock (x)
+  y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
+endfunction
+
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-function",rosenbrock);
+[this,f] = optimbase_function ( opt , [0.0 0.0] );
+assert_close ( f , 1.0 , %eps );
+opt = optimbase_destroy(opt);
+
+
+//
+// Test with an additional argument
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function y = rosenbrock2 ( x , mydata )
+  a = mydata.a
+  y = 100*(x(2)-x(1)^2)^2 + ( a - x(1))^2;
+endfunction
+
+mystuff = tlist(["T_MYSTUFF","a"]);
+mystuff.a = 12.0;
+
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-function",rosenbrock2);
+opt = optimbase_configure(opt,"-costfargument",mystuff);
+[this,f] = optimbase_function ( opt , [0.0 0.0] );
+assert_close ( f , 144. , %eps );
+opt = optimbase_destroy(opt);
+
+//
+// Test with non linear constraints : there is an index
+// optimtestcase --
+//   Non linear inequality constraints are positive.
+//    
+// Arguments
+//   x: the point where to compute the function
+//   index : the stuff to compute
+// Note
+//  The following protocol is used
+//  * if index=1, or no index, returns the value of the cost 
+//    function (default case)
+//  * if index=2, returns the value of the nonlinear inequality 
+//    constraints, as a row array
+//  * if index=3, returns an array which contains
+//    at index #0, the value of the cost function  
+//    at index #1 to the end is the list of the values of the nonlinear 
+//    constraints
+//  The inequality constraints are expected to be positive.
+//
+function result = optimtestcase ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index == 1 | index == 3 ) then
+    f = x(1)^2 + x(2)^2 + 2.0 * x(3)^2 + x(4)^2 ...
+      - 5.0 * x(1) - 5.0 * x(2) - 21.0 * x(3) + 7.0 * x(4)
+  end
+  if ( index == 2 | index == 3 ) then
+    c1 = - x(1)^2 - x(2)^2 - x(3)^2 - x(4)^2 ...
+              - x(1) + x(2) - x(3) + x(4) + 8
+    c2 = - x(1)^2 - 2.0 * x(2)^2 - x(3)^2 - 2.0 * x(4)^2 ...
+              + x(1) + x(4) + 10.0
+    c3 = - 2.0 * x(1)^2 - x(2)^2 - x(3)^2 - 2.0 * x(1) ...
+              + x(2) + x(4) + 5.0
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3]
+  case 3 then
+    result = [f c1 c2 c3]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+x0 = [0.0 0.0 0.0 0.0].';
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",4);
+opt= optimbase_configure(opt,"-function",optimtestcase);
+[this,f] = optimbase_function ( opt , x0 );
+assert_close ( f , 0. , %eps );
+[this,f] = optimbase_function ( opt , x0 , 1 );
+assert_close ( f , 0. , %eps );
+[this,f] = optimbase_function ( opt , x0 , 2 );
+assert_close ( f , [8.    10.    5.] , %eps );
+[this,f] = optimbase_function ( opt , x0 , 3 );
+assert_close ( f , [0.    8.    10.    5.] , %eps );
+opt = optimbase_destroy(opt);
+
+
+//
+// Test when there are both nonlinear constraints and a customized data
+// boxproblemA --
+//   Computes the Box problem A cost function and 
+//   inequality constraints.
+//
+// Arguments
+//   x: the point where to compute the function
+//   index : the stuff to compute
+//   data : the parameters of Box cost function
+// Note
+//  The following protocol is used
+//  * if index=1, or no index, returns the value of the cost 
+//    function (default case)
+//  * if index=2, returns the value of the nonlinear inequality 
+//    constraints, as a row array
+//  * if index=3, returns an array which contains
+//    at index #0, the value of the cost function  
+//    at index #1 to the end is the list of the values of the nonlinear 
+//    constraints
+//  The inequality constraints are expected to be positive.
+//
+function result = boxproblemA ( x , index , data )
+  b = x(2) + 0.01 * x(3)
+  x6 = (data.k1 + data.k2 * x(2) ...
+            + data.k3 * x(3) + data.k4 * x(4) + data.k5 * x(5)) * x(1)
+  y1 = data.k6 + data.k7 * x(2) + data.k8 * x(3) ...
+            + data.k9 * x(4) + data.k10 * x(5)
+  y2 = data.k11 + data.k12 * x(2) + data.k13 * x(3) ...
+            + data.k14 * x(4) + data.k15 * x(5)
+  y3 = data.k16 + data.k17 * x(2) + data.k18 * x(3) ...
+            + data.k19 * x(4) + data.k20 * x(5)
+  y4 = data.k21 + data.k22 * x(2) + data.k23 * x(3) ...
+            + data.k24 * x(4) + data.k25 * x(5)
+  x7 = ( y1 + y2 + y3 ) * x(1)
+  x8 = (data.k26 + data.k27 * x(2) + data.k28 * x(3) ...
+            + data.k29 * x(4) ...
+            + data.k30 * x(5) ) * x(1) + x6 + x7
+  if ( index==1 | index==3 ) then
+    f = (data.a2 * y1 + data.a3 * y2 + data.a4 * y3 + data.a5 * y4 ...
+             + 7840 * data.a6 - 100000 * data.a0 ...
+             - 50800 * b * data.a7 + data.k31 + data.k32 * x(2) + data.k33 * x(3) ...
+             + data.k34 * x(4) + data.k35 * x(5)) * x(1) ...
+             - 24345 + data.a1 * x6
+    f = -f
+  end
+  if ( index==2 | index==3 ) then
+      c1 = x6
+      c2 = 294000 - x6
+      c3 = x7
+      c4 = 294000 - x7
+      c5 = x8
+      c6 = 277200 - x8
+  end
+  select index
+  case 1 then
+    result = f
+  case 2 then
+    result = [c1 c2 c3 c4 c5 c6]
+  case 3 then
+    result = [f c1 c2 c3 c4 c5 c6]
+  else
+    errmsg = sprintf("Unexpected index %d" , index);
+    error(errmsg);
+  end
+endfunction
+
+boxparams = struct();
+boxparams.a0 = 9;
+boxparams.a1 = 15;
+boxparams.a2 = 50;
+boxparams.a3 = 9.583;
+boxparams.a4 = 20;
+boxparams.a5 = 15;
+boxparams.a6 = 6;
+boxparams.a7 = 0.75;
+boxparams.k1 = -145421.402;
+boxparams.k2 = 2931.1506;
+boxparams.k3 = -40.427932;
+boxparams.k4 = 5106.192;
+boxparams.k5 = 15711.36;
+boxparams.k6 = -161622.577;
+boxparams.k7 = 4176.15328;
+boxparams.k8 = 2.8260078;
+boxparams.k9 = 9200.476;
+boxparams.k10 = 13160.295;
+boxparams.k11 = -21686.9194;
+boxparams.k12 = 123.56928;
+boxparams.k13 = -21.1188894;
+boxparams.k14 = 706.834;
+boxparams.k15 = 2898.573;
+boxparams.k16 = 28298.388;
+boxparams.k17 = 60.81096;
+boxparams.k18 = 31.242116;
+boxparams.k19 = 329.574;
+boxparams.k20 = -2882.082;
+boxparams.k21 = 74095.3845;
+boxparams.k22 = -306.262544;
+boxparams.k23 = 16.243649;
+boxparams.k24 = -3094.252;
+boxparams.k25 = -5566.2628;
+boxparams.k26 = -26237.0;
+boxparams.k27 = 99.0;
+boxparams.k28 = -0.42;
+boxparams.k29 = 1300.0;
+boxparams.k30 = 2100.0;
+boxparams.k31 = 925548.252;
+boxparams.k32 = -61968.8432;
+boxparams.k33 = 23.3088196;
+boxparams.k34 = -27097.648;
+boxparams.k35 = -50843.766;
+x0 = [2.52 2.0 37.5 9.25 6.8].';
+
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",5);
+opt= optimbase_configure(opt,"-function",boxproblemA);
+opt = optimbase_configure(opt,"-costfargument",boxparams);
+opt = optimbase_configure(opt,"-nbineqconst",6);
+[this,f] = optimbase_function ( opt , x0 );
+assert_close ( f , -2351243.5  , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 1 );
+assert_close ( f , -2351243.5  , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 2 );
+assert_close ( f , [32745.827    261254.17    96991.969    197008.03    130368.43    146831.57] , 1.e-7 );
+[this,f] = optimbase_function ( opt , x0 , 3 );
+assert_close ( f , [-2351243.5    32745.827    261254.17    96991.969    197008.03    130368.43    146831.57] , 1.e-7 );
+opt = optimbase_destroy(opt);
+
index 17866a6..47a9a54 100644 (file)
@@ -48,13 +48,13 @@ endfunction
 // Arguments
 //    x : the point where to compute the cost
 //    index : a flag which states what is to compute
-//    * if index=1, or no index, returns the value of the cost
+//    * if index=1, or no index, returns the value of the cost 
 //      function (default case)
-//    * if index=2, returns the value of the nonlinear inequality
+//    * if index=2, returns the value of the nonlinear inequality 
 //      constraints, as a row array
 //    * if index=3, returns an array which contains
-//      at index #1, the value of the cost function
-//      at index #2 to the end is the list of the values of the nonlinear
+//      at index #1, the value of the cost function  
+//      at index #2 to the end is the list of the values of the nonlinear 
 //      inequality constraints
 //  Discussion:
 //    The problem is to minimize a cost function with 4 non linear constraints.
@@ -104,8 +104,8 @@ function result = gouldnonconvex ( x , index )
   end
 endfunction
 //
-// The same cost function as before, with an
-// additionnal argument, which contains parameters of the
+// The same cost function as before, with an 
+// additionnal argument, which contains parameters of the 
 // cost function and constraints.
 // In this case, the mydata variable is passed
 // explicitely by the optimization class.
@@ -153,9 +153,11 @@ opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
 assert_equal ( isfeasible , 1 );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [-6.0 0.0] );
 Component #1/2 of x is lower than min bound -5.000000e+000
+
 assert_equal ( isfeasible , 0 );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [0.0 6.0] );
 Component #2/2 of x is greater than max bound 5.000000e+000
+
 assert_equal ( isfeasible , 0 );
 opt = optimbase_destroy(opt);
 //
@@ -167,12 +169,17 @@ opt = optimbase_configure(opt,"-nbineqconst",4);
 opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [ 14.0950013 , 0.8429636 ] );
 Computed constraints = 1.095001e+000 2.779266e-007 2.322073e-006 8.429636e-001
+
 Function Evaluation #1 is [1.0950013 0.0000003 0.0000023 0.8429636] at [14.095001 0.8429636]
+
 assert_equal ( isfeasible , 1 );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [ 14.0950013 , 0.0 ] );
 Computed constraints = 1.095001e+000 7.719049e+000 -7.719046e+000 0.000000e+000
+
 Function Evaluation #2 is [1.0950013 7.7190486 -7.719046 0] at [14.095001 0]
+
 Inequality constraint #3/4 is not satisfied for x
+
 assert_equal ( isfeasible , -1 );
 opt = optimbase_destroy(opt);
 //
@@ -219,11 +226,16 @@ opt = optimbase_configure ( opt , "-function" , gouldnonconvex2 );
 opt = optimbase_configure ( opt , "-costfargument" , mystuff );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [ 14.0950013 , 0.8429636 ] );
 Computed constraints = 1.095001e+000 2.779266e-007 2.322073e-006 8.429636e-001
+
 Function Evaluation #1 is [1.0950013 0.0000003 0.0000023 0.8429636] at [14.095001 0.8429636]
+
 assert_equal ( isfeasible , 1 );
 [ opt , isfeasible ] = optimbase_isfeasible ( opt ,  [ 14.0950013 , 0.0 ] );
 Computed constraints = 1.095001e+000 7.719049e+000 -7.719046e+000 0.000000e+000
+
 Function Evaluation #2 is [1.0950013 7.7190486 -7.719046 0] at [14.095001 0]
+
 Inequality constraint #3/4 is not satisfied for x
+
 assert_equal ( isfeasible , -1 );
 opt = optimbase_destroy(opt);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.dia.ref b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.dia.ref
new file mode 100644 (file)
index 0000000..80097a4
--- /dev/null
@@ -0,0 +1,193 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+//
+// 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 bugmes();quit;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 bugmes();quit;end
+endfunction
+//
+// gould.nonconvex --
+//   The Gould test case with additionnal inequality constraints.
+// Arguments
+//    x : the point where to compute the cost
+//    index : a flag which states what is to compute
+//    * if index=1, or no index, returns the value of the cost 
+//      function (default case)
+//    * if index=2, returns the value of the nonlinear inequality 
+//      constraints, as a row array
+//    * if index=3, returns an array which contains
+//      at index #1, the value of the cost function  
+//      at index #2 to the end is the list of the values of the nonlinear 
+//      inequality constraints
+//  Discussion:
+//    The problem is to minimize a cost function with 4 non linear constraints.
+//    This is Problem 4.1 in Subrahmanyam, extracted from Gould.
+//    Non convex.
+//    The constraint region is a narrow winding (half-moon shaped) valley.
+//    Solution showed with tolerance 1.e-8.
+//
+//  Reference:
+//    An extension of the simplex method to constrained
+//    nonlinear optimization
+//    M.B. Subrahmanyam
+//    Journal of optimization theory and applications
+//    Vol. 62, August 1989
+//
+//    Gould F.J.
+//    Nonlinear Tolerance Programming
+//    Numerical methods for Nonlinear optimization
+//    Edited by F.A. Lootsma, pp 349-366, 1972
+//
+function result = gouldnonconvex ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - 10.0 )^3 + ( x(2) - 20.0 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - 13.0
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - 100.0
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + 82.81
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// The same cost function as before, with an 
+// additionnal argument, which contains parameters of the 
+// cost function and constraints.
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function result = gouldnonconvex2 ( x , index , mydata )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - mydata.f1 )^3 + ( x(2) - mydata.f2 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - mydata.a1
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - mydata.a2
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + mydata.a3
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// Test optimbase_isfeasible method
+//
+// Test with bounds
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure ( opt , "-boundsmin" , [-5.0 -5.0] );
+opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [0.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [-6.0 0.0] );
+Component #1/2 of x is lower than min bound -5.000000e+000
+
+assert_equal ( isfeasible , %f );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [0.0 6.0] );
+Component #2/2 of x is greater than max bound 5.000000e+000
+
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure(opt,"-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints and additionnal argument in cost function
+mystuff = struct();
+mystuff.f1 = 10.0;
+mystuff.f2 = 20.0;
+mystuff.a1 = 13.0;
+mystuff.a2 = 100.0;
+mystuff.a3 = 82.81;
+opt = optimbase_new ();
+opt = optimbase_configure ( opt , "-numberofvariables",2);
+opt = optimbase_configure ( opt , "-verbose",1);
+opt = optimbase_configure ( opt , "-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex2 );
+opt = optimbase_configure ( opt , "-costfargument" , mystuff );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.tst b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinbounds.tst
new file mode 100644 (file)
index 0000000..a101a3a
--- /dev/null
@@ -0,0 +1,192 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+//
+// 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
+//
+// gould.nonconvex --
+//   The Gould test case with additionnal inequality constraints.
+// Arguments
+//    x : the point where to compute the cost
+//    index : a flag which states what is to compute
+//    * if index=1, or no index, returns the value of the cost 
+//      function (default case)
+//    * if index=2, returns the value of the nonlinear inequality 
+//      constraints, as a row array
+//    * if index=3, returns an array which contains
+//      at index #1, the value of the cost function  
+//      at index #2 to the end is the list of the values of the nonlinear 
+//      inequality constraints
+//  Discussion:
+//    The problem is to minimize a cost function with 4 non linear constraints.
+//    This is Problem 4.1 in Subrahmanyam, extracted from Gould.
+//    Non convex.
+//    The constraint region is a narrow winding (half-moon shaped) valley.
+//    Solution showed with tolerance 1.e-8.
+//
+//  Reference:
+//    An extension of the simplex method to constrained
+//    nonlinear optimization
+//    M.B. Subrahmanyam
+//    Journal of optimization theory and applications
+//    Vol. 62, August 1989
+//
+//    Gould F.J.
+//    Nonlinear Tolerance Programming
+//    Numerical methods for Nonlinear optimization
+//    Edited by F.A. Lootsma, pp 349-366, 1972
+//
+function result = gouldnonconvex ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - 10.0 )^3 + ( x(2) - 20.0 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - 13.0
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - 100.0
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + 82.81
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// The same cost function as before, with an 
+// additionnal argument, which contains parameters of the 
+// cost function and constraints.
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function result = gouldnonconvex2 ( x , index , mydata )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - mydata.f1 )^3 + ( x(2) - mydata.f2 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - mydata.a1
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - mydata.a2
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + mydata.a3
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// Test optimbase_isfeasible method
+//
+// Test with bounds
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure ( opt , "-boundsmin" , [-5.0 -5.0] );
+opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [0.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [-6.0 0.0] );
+assert_equal ( isfeasible , %f );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [0.0 6.0] );
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure(opt,"-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints and additionnal argument in cost function
+mystuff = struct();
+mystuff.f1 = 10.0;
+mystuff.f2 = 20.0;
+mystuff.a1 = 13.0;
+mystuff.a2 = 100.0;
+mystuff.a3 = 82.81;
+opt = optimbase_new ();
+opt = optimbase_configure ( opt , "-numberofvariables",2);
+opt = optimbase_configure ( opt , "-verbose",1);
+opt = optimbase_configure ( opt , "-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex2 );
+opt = optimbase_configure ( opt , "-costfargument" , mystuff );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinbounds ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
+
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.dia.ref b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.dia.ref
new file mode 100644 (file)
index 0000000..a9db2b9
--- /dev/null
@@ -0,0 +1,209 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+//
+// 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 bugmes();quit;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 bugmes();quit;end
+endfunction
+//
+// gould.nonconvex --
+//   The Gould test case with additionnal inequality constraints.
+// Arguments
+//    x : the point where to compute the cost
+//    index : a flag which states what is to compute
+//    * if index=1, or no index, returns the value of the cost 
+//      function (default case)
+//    * if index=2, returns the value of the nonlinear inequality 
+//      constraints, as a row array
+//    * if index=3, returns an array which contains
+//      at index #1, the value of the cost function  
+//      at index #2 to the end is the list of the values of the nonlinear 
+//      inequality constraints
+//  Discussion:
+//    The problem is to minimize a cost function with 4 non linear constraints.
+//    This is Problem 4.1 in Subrahmanyam, extracted from Gould.
+//    Non convex.
+//    The constraint region is a narrow winding (half-moon shaped) valley.
+//    Solution showed with tolerance 1.e-8.
+//
+//  Reference:
+//    An extension of the simplex method to constrained
+//    nonlinear optimization
+//    M.B. Subrahmanyam
+//    Journal of optimization theory and applications
+//    Vol. 62, August 1989
+//
+//    Gould F.J.
+//    Nonlinear Tolerance Programming
+//    Numerical methods for Nonlinear optimization
+//    Edited by F.A. Lootsma, pp 349-366, 1972
+//
+function result = gouldnonconvex ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - 10.0 )^3 + ( x(2) - 20.0 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - 13.0
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - 100.0
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + 82.81
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// The same cost function as before, with an 
+// additionnal argument, which contains parameters of the 
+// cost function and constraints.
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function result = gouldnonconvex2 ( x , index , mydata )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - mydata.f1 )^3 + ( x(2) - mydata.f2 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - mydata.a1
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - mydata.a2
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + mydata.a3
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// Test optimbase_isfeasible method
+//
+// Test with bounds
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure ( opt , "-boundsmin" , [-5.0 -5.0] );
+opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [0.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [-6.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [0.0 6.0] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure(opt,"-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.8429636 ] );
+Computed constraints = 1.095001e+000 2.779266e-007 2.322073e-006 8.429636e-001
+
+Function Evaluation #1 is [1.0950013 0.0000003 0.0000023 0.8429636] at [14.095001 0.8429636]
+
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.0 ] );
+Computed constraints = 1.095001e+000 7.719049e+000 -7.719046e+000 0.000000e+000
+
+Function Evaluation #2 is [1.0950013 7.7190486 -7.719046 0] at [14.095001 0]
+
+Inequality constraint #3/4 is not satisfied for x
+
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints and additionnal argument in cost function
+mystuff = struct();
+mystuff.f1 = 10.0;
+mystuff.f2 = 20.0;
+mystuff.a1 = 13.0;
+mystuff.a2 = 100.0;
+mystuff.a3 = 82.81;
+opt = optimbase_new ();
+opt = optimbase_configure ( opt , "-numberofvariables",2);
+opt = optimbase_configure ( opt , "-verbose",1);
+opt = optimbase_configure ( opt , "-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex2 );
+opt = optimbase_configure ( opt , "-costfargument" , mystuff );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.8429636 ] );
+Computed constraints = 1.095001e+000 2.779266e-007 2.322073e-006 8.429636e-001
+
+Function Evaluation #1 is [1.0950013 0.0000003 0.0000023 0.8429636] at [14.095001 0.8429636]
+
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.0 ] );
+Computed constraints = 1.095001e+000 7.719049e+000 -7.719046e+000 0.000000e+000
+
+Function Evaluation #2 is [1.0950013 7.7190486 -7.719046 0] at [14.095001 0]
+
+Inequality constraint #3/4 is not satisfied for x
+
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
diff --git a/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.tst b/scilab/modules/optimization/tests/unit_tests/optimbase/optimbase_isinnonlincons.tst
new file mode 100644 (file)
index 0000000..1e81af3
--- /dev/null
@@ -0,0 +1,192 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008-2009 - INRIA - 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
+
+
+//
+// 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
+//
+// gould.nonconvex --
+//   The Gould test case with additionnal inequality constraints.
+// Arguments
+//    x : the point where to compute the cost
+//    index : a flag which states what is to compute
+//    * if index=1, or no index, returns the value of the cost 
+//      function (default case)
+//    * if index=2, returns the value of the nonlinear inequality 
+//      constraints, as a row array
+//    * if index=3, returns an array which contains
+//      at index #1, the value of the cost function  
+//      at index #2 to the end is the list of the values of the nonlinear 
+//      inequality constraints
+//  Discussion:
+//    The problem is to minimize a cost function with 4 non linear constraints.
+//    This is Problem 4.1 in Subrahmanyam, extracted from Gould.
+//    Non convex.
+//    The constraint region is a narrow winding (half-moon shaped) valley.
+//    Solution showed with tolerance 1.e-8.
+//
+//  Reference:
+//    An extension of the simplex method to constrained
+//    nonlinear optimization
+//    M.B. Subrahmanyam
+//    Journal of optimization theory and applications
+//    Vol. 62, August 1989
+//
+//    Gould F.J.
+//    Nonlinear Tolerance Programming
+//    Numerical methods for Nonlinear optimization
+//    Edited by F.A. Lootsma, pp 349-366, 1972
+//
+function result = gouldnonconvex ( x , index )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - 10.0 )^3 + ( x(2) - 20.0 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - 13.0
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - 100.0
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + 82.81
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// The same cost function as before, with an 
+// additionnal argument, which contains parameters of the 
+// cost function and constraints.
+// In this case, the mydata variable is passed
+// explicitely by the optimization class.
+// So the actual name "mydata" does not matter
+// and whatever variable name can be used.
+//
+function result = gouldnonconvex2 ( x , index , mydata )
+  if (~isdef('index','local')) then
+    index = 1
+  end
+  if ( index==1 | index==3 ) then
+    f = (x(1) - mydata.f1 )^3 + ( x(2) - mydata.f2 ) ^ 3
+  end
+  if ( index==2 | index==3 ) then
+    c1 = x(1) - mydata.a1
+    c2 = ( x(1) - 5.0 )^2  + (x(2) - 5.0 )^2 - mydata.a2
+    c3 = -( x(1) - 6.0 )^2 - (x(2) - 5.0 )^2 + mydata.a3
+    c4 = x(2)
+  end
+  select index
+  case 1 then
+      result = f
+      mprintf( "Computed f = %e\n", f);
+  case 2
+      result = [c1 c2 c3 c4]
+      mprintf( "Computed constraints = %e %e %e %e\n", c1 , c2 , c3 , c4);
+  case 3
+      result = [f c1 c2 c3 c4]
+      mprintf( "Computed f = %e and constraints = %e %e %e %e\n", f , c1 , c2 , c3 , c4);
+  else
+    errmsg = sprintf("Unknown index %d", index )
+    error(errmsg)
+  end
+endfunction
+//
+// Test optimbase_isfeasible method
+//
+// Test with bounds
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure ( opt , "-boundsmin" , [-5.0 -5.0] );
+opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [0.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [-6.0 0.0] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [0.0 6.0] );
+assert_equal ( isfeasible , %t );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints
+opt = optimbase_new ();
+opt = optimbase_configure(opt,"-numberofvariables",2);
+opt = optimbase_configure(opt,"-verbose",1);
+opt = optimbase_configure(opt,"-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
+//
+// Test with nonlinear inequality constraints and additionnal argument in cost function
+mystuff = struct();
+mystuff.f1 = 10.0;
+mystuff.f2 = 20.0;
+mystuff.a1 = 13.0;
+mystuff.a2 = 100.0;
+mystuff.a3 = 82.81;
+opt = optimbase_new ();
+opt = optimbase_configure ( opt , "-numberofvariables",2);
+opt = optimbase_configure ( opt , "-verbose",1);
+opt = optimbase_configure ( opt , "-nbineqconst",4);
+opt = optimbase_configure ( opt , "-function" , gouldnonconvex2 );
+opt = optimbase_configure ( opt , "-costfargument" , mystuff );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.8429636 ] );
+assert_equal ( isfeasible , %t );
+[ opt , isfeasible ] = optimbase_isinnonlincons ( opt ,  [ 14.0950013 , 0.0 ] );
+assert_equal ( isfeasible , %f );
+opt = optimbase_destroy(opt);
+
index 65af4b6..ffa0b14 100644 (file)
@@ -57,7 +57,9 @@ opt = optimbase_configure ( opt , "-boundsmax" , [5.0 5.0] );
 assert_equal ( p , [0.0 0.0] );
 [ opt , p ] = optimbase_proj2bnds ( opt ,  [-6.0 6.0] );
 Projecting p(1) = -6.000000e+000 on min bound -5.000000e+000
+
 Projecting p(2) = 6.000000e+000 on max bound 5.000000e+000
+
 assert_equal ( p , [-5.0 5.0] );
 opt = optimbase_destroy(opt);
 //
index 111362f..b469660 100644 (file)
@@ -46,16 +46,63 @@ function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
 //
+// Test with x0 only
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.  -1.2    1.  
+    0.  -0.2    1.  
+    0.  -1.2    2.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Test with function only
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    24.2  -1.2    1.  
+    93.6  -0.2    1.  
+    36.2  -1.2    2.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+//
 // Test with a scalar length
 //
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock , 2.0 );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    24.2   -1.2    1.  
+    13.      0.8    1.  
+    248.2  -1.2    3.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Test with a vector length
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock , [1.0 2.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+24.2 -1.2 1.
+93.6 -0.2 1.
+248.2 -1.2 3.
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Test with a scalar length and an additionnal object
+//
 myobj = tlist(["T_MYSTUFF","nb"]);
 myobj.nb = 0;
 function [ y , myobj ] = mycostf ( x , myobj )
   y = rosenbrock(x);
   myobj.nb = myobj.nb + 1
 endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_axes ( s1 , x0 = [-1.2 1.0] , fun = mycostf , len = 1.0, data=myobj );
+[ s1 , myobj ] = optimsimplex_new ( "axes" , [-1.2 1.0] , mycostf , 1.0, myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
 24.2 -1.2 1.0
@@ -67,16 +114,3 @@ assert_equal ( myobj.nb , 3 );
 nbve = optimsimplex_getnbve ( s1 );
 assert_equal ( nbve , 3 );
 s1 = optimsimplex_destroy ( s1 );
-//
-// Test with a vector length
-//
-s1 = optimsimplex_new ();
-s1 = optimsimplex_axes ( s1 , x0 = [-1.2 1.0] , fun = rosenbrock , len = [1.0 2.0] );
-computed = optimsimplex_getall ( s1 );
-expected = [
-24.2 -1.2 1.
-93.6 -0.2 1.
-248.2 -1.2 3.
-];
-assert_close ( computed , expected , %eps );
-s1 = optimsimplex_destroy ( s1 );
index 72a658c..c47cc21 100644 (file)
@@ -47,18 +47,68 @@ function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
 
+//
+// Test with x0 only
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.  -1.2    1.  
+    0.  -0.2    1.  
+    0.  -1.2    2.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Test with function only
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    24.2  -1.2    1.  
+    93.6  -0.2    1.  
+    36.2  -1.2    2.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
 
 //
 // Test with a scalar length
 //
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock , 2.0 );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    24.2   -1.2    1.  
+    13.      0.8    1.  
+    248.2  -1.2    3.  
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Test with a vector length
+//
+s1 = optimsimplex_new ( "axes" , [-1.2 1.0] , rosenbrock , [1.0 2.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+24.2 -1.2 1.
+93.6 -0.2 1.
+248.2 -1.2 3.
+];
+assert_close ( computed , expected , %eps );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Test with a scalar length and an additionnal object
+//
 myobj = tlist(["T_MYSTUFF","nb"]);
 myobj.nb = 0;
 function [ y , myobj ] = mycostf ( x , myobj )
   y = rosenbrock(x);
   myobj.nb = myobj.nb + 1
 endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_axes ( s1 , x0 = [-1.2 1.0] , fun = mycostf , len = 1.0, data=myobj );
+[ s1 , myobj ] = optimsimplex_new ( "axes" , [-1.2 1.0] , mycostf , 1.0, myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
 24.2 -1.2 1.0
@@ -70,17 +120,4 @@ assert_equal ( myobj.nb , 3 );
 nbve = optimsimplex_getnbve ( s1 );
 assert_equal ( nbve , 3 );
 s1 = optimsimplex_destroy ( s1 );
-//
-// Test with a vector length
-//
-s1 = optimsimplex_new ();
-s1 = optimsimplex_axes ( s1 , x0 = [-1.2 1.0] , fun = rosenbrock , len = [1.0 2.0] );
-computed = optimsimplex_getall ( s1 );
-expected = [
-24.2 -1.2 1.
-93.6 -0.2 1.
-248.2 -1.2 3.
-];
-assert_close ( computed , expected , %eps );
-s1 = optimsimplex_destroy ( s1 );
 
index 9d66bbd..6c9e3b0 100644 (file)
@@ -79,11 +79,11 @@ Dimension : 2
 
 Number of vertices : 3
 
-Vertex #1/3 : fv=12.000000, x=0.000000 0.000000
+Vertex #1/3 : fv=1.200000e+001, x=0.000000e+000 0.000000e+000
 
-Vertex #2/3 : fv=13.000000, x=1.000000 0.000000
+Vertex #2/3 : fv=1.300000e+001, x=1.000000e+000 0.000000e+000
 
-Vertex #3/3 : fv=14.000000, x=0.000000 2.000000
+Vertex #3/3 : fv=1.400000e+001, x=0.000000e+000 2.000000e+000
 
 // We are done !
 s1 = optimsimplex_destroy(s1);
index 9eff9f8..5b031de 100644 (file)
@@ -45,14 +45,52 @@ endfunction
 function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
+//
+// Pfeffer, case with x0 only
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.0      -1.2     0.      
+    0.0  -1.26    0.      
+    0.0  -1.2     0.0075  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Pfeffer, case with x0 and function
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] , rosenbrock );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    212.2      -1.2     0.      
+    257.15498  -1.26    0.      
+    210.04562  -1.2     0.0075  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Pfeffer, case with specified deltas
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] , rosenbrock , 0.1 , 0.01 );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    212.2      -1.2     0.    
+    308.97818  -1.32    0.    
+    209.33     -1.2     0.01  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Case with additionnal object
+//
 myobj = tlist(["T_MYSTUFF","nb"]);
 myobj.nb = 0;
 function [ y , myobj ] = mycostf ( x , myobj )
   y = rosenbrock(x);
   myobj.nb = myobj.nb + 1
 endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_pfeffer ( s1 , x0 = [-1.2 1.0] , fun = mycostf , data=myobj );
+[ s1 , myobj ] = optimsimplex_new ( "pfeffer" , [-1.2 1.0] , mycostf , 0.05 , 0.0075 , myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2       -1.2     1.    
@@ -62,16 +100,3 @@ expected = [
 assert_close ( computed , expected , 10 * %eps );
 assert_equal ( myobj.nb , 3 );
 s1 = optimsimplex_destroy ( s1 );
-//
-// Pfeffer, case 0.0
-//
-s1 = optimsimplex_new ();
-s1 = optimsimplex_pfeffer ( s1 , x0 = [-1.2 0.0] , fun = rosenbrock );
-computed = optimsimplex_getall ( s1 );
-expected = [
-    212.2      -1.2     0.      
-    257.15498  -1.26    0.      
-    210.04562  -1.2     0.0075  
-];
-assert_close ( computed , expected , 1.e-6 );
-s1 = optimsimplex_destroy ( s1 );
index ec26e74..d5117de 100644 (file)
@@ -47,15 +47,55 @@ function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
 
+//
+// Pfeffer, case with x0 only
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.0      -1.2     0.      
+    0.0  -1.26    0.      
+    0.0  -1.2     0.0075  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Pfeffer, case with x0 and function
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] , rosenbrock );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    212.2      -1.2     0.      
+    257.15498  -1.26    0.      
+    210.04562  -1.2     0.0075  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Pfeffer, case with specified deltas
+//
+s1 = optimsimplex_new ( "pfeffer" , [-1.2 0.0] , rosenbrock , 0.1 , 0.01 );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    212.2      -1.2     0.    
+    308.97818  -1.32    0.    
+    209.33     -1.2     0.01  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
 
+//
+// Case with additionnal object
+//
 myobj = tlist(["T_MYSTUFF","nb"]);
 myobj.nb = 0;
 function [ y , myobj ] = mycostf ( x , myobj )
   y = rosenbrock(x);
   myobj.nb = myobj.nb + 1
 endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_pfeffer ( s1 , x0 = [-1.2 1.0] , fun = mycostf , data=myobj );
+[ s1 , myobj ] = optimsimplex_new ( "pfeffer" , [-1.2 1.0] , mycostf , 0.05 , 0.0075 , myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2       -1.2     1.    
@@ -66,17 +106,4 @@ assert_close ( computed , expected , 10 * %eps );
 assert_equal ( myobj.nb , 3 );
 s1 = optimsimplex_destroy ( s1 );
 
-//
-// Pfeffer, case 0.0
-//
-s1 = optimsimplex_new ();
-s1 = optimsimplex_pfeffer ( s1 , x0 = [-1.2 0.0] , fun = rosenbrock );
-computed = optimsimplex_getall ( s1 );
-expected = [
-    212.2      -1.2     0.      
-    257.15498  -1.26    0.      
-    210.04562  -1.2     0.0075  
-];
-assert_close ( computed , expected , 1.e-6 );
-s1 = optimsimplex_destroy ( s1 );
 
index 3ab08bc..fac6d51 100644 (file)
@@ -63,11 +63,11 @@ Dimension : 2
 
 Number of vertices : 3
 
-Vertex #1/3 : fv=12.000000, x=0.000000 0.000000
+Vertex #1/3 : fv=1.200000e+001, x=0.000000e+000 0.000000e+000
 
-Vertex #2/3 : fv=13.000000, x=1.000000 0.000000
+Vertex #2/3 : fv=1.300000e+001, x=1.000000e+000 0.000000e+000
 
-Vertex #3/3 : fv=14.000000, x=0.000000 2.000000
+Vertex #3/3 : fv=1.400000e+001, x=0.000000e+000 2.000000e+000
 
 s1 = optimsimplex_destroy ( s1 );
 //
index aff7d3a..40000e9 100644 (file)
@@ -46,6 +46,40 @@ function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
 //
+// Test randbounds with default number of vertices
+//
+s1 = optimsimplex_new ( "randbounds" , [-1.2 1.0], rosenbrock, ...
+  [-5.0 -5.0] , [5.0 5.0] );
+computed = optimsimplex_getall ( s1 );
+assert_equal ( size(computed,1) , 3 );
+assert_equal ( size(computed,2) , 3 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:3,2:3) < 5.0 , expected );
+assert_equal ( computed(1:3,2:3) > -5.0 , expected );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Test randbounds with 5 vertices
+//
+s1 = optimsimplex_new ( "randbounds" , [-1.2 1.0], rosenbrock, ...
+  [-5.0 -5.0] , [5.0 5.0], 5 );
+computed = optimsimplex_getall ( s1 );
+assert_equal ( size(computed,1) , 5 );
+assert_equal ( size(computed,2) , 3 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:5,2:3) < 5.0 , expected );
+assert_equal ( computed(1:5,2:3) > -5.0 , expected );
+s1 = optimsimplex_destroy ( s1 );
+//
 // Test optimsimplex_randbounds
 //
 function [ y , myobj ] = mycostf ( x , myobj )
@@ -53,15 +87,24 @@ function [ y , myobj ] = mycostf ( x , myobj )
   myobj.nb = myobj.nb + 1
 endfunction
 //
-// Test randbounds
+// Test randbounds with additionnal object
 //
 mydude = tlist(["T_MYSTUFF","nb"]);
 mydude.nb = 0;
 s1 = optimsimplex_new ();
-[ s1 , mydude ] = optimsimplex_randbounds ( s1 , x0 = [-1.2 1.0], fun = mycostf, ...
-  boundsmin = [-5.0 -5.0] , boundsmax = [5.0 5.0], nbve=5 , data = mydude );
+[ s1 , mydude ] = optimsimplex_new ( "randbounds" , [-1.2 1.0], mycostf, ...
+  [-5.0 -5.0] , [5.0 5.0], 5 , mydude );
 computed = optimsimplex_getall ( s1 );
 assert_equal ( size(computed,1) , 5 );
 assert_equal ( size(computed,2) , 3 );
 assert_equal ( mydude.nb , 5 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:5,2:3) < 5.0 , expected );
+assert_equal ( computed(1:5,2:3) > -5.0 , expected );
 s1 = optimsimplex_destroy ( s1 );
index d24bc05..c7f971c 100644 (file)
@@ -48,6 +48,42 @@ function y = rosenbrock (x)
 endfunction
 
 //
+// Test randbounds with default number of vertices
+//
+s1 = optimsimplex_new ( "randbounds" , [-1.2 1.0], rosenbrock, ...
+  [-5.0 -5.0] , [5.0 5.0] );
+computed = optimsimplex_getall ( s1 );
+assert_equal ( size(computed,1) , 3 );
+assert_equal ( size(computed,2) , 3 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:3,2:3) < 5.0 , expected );
+assert_equal ( computed(1:3,2:3) > -5.0 , expected );
+s1 = optimsimplex_destroy ( s1 );
+
+//
+// Test randbounds with 5 vertices
+//
+s1 = optimsimplex_new ( "randbounds" , [-1.2 1.0], rosenbrock, ...
+  [-5.0 -5.0] , [5.0 5.0], 5 );
+computed = optimsimplex_getall ( s1 );
+assert_equal ( size(computed,1) , 5 );
+assert_equal ( size(computed,2) , 3 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:5,2:3) < 5.0 , expected );
+assert_equal ( computed(1:5,2:3) > -5.0 , expected );
+s1 = optimsimplex_destroy ( s1 );
+
+//
 // Test optimsimplex_randbounds
 //
 function [ y , myobj ] = mycostf ( x , myobj )
@@ -56,16 +92,25 @@ function [ y , myobj ] = mycostf ( x , myobj )
 endfunction
 
 //
-// Test randbounds
+// Test randbounds with additionnal object
 //
 mydude = tlist(["T_MYSTUFF","nb"]);
 mydude.nb = 0;
 s1 = optimsimplex_new ();
-[ s1 , mydude ] = optimsimplex_randbounds ( s1 , x0 = [-1.2 1.0], fun = mycostf, ...
-  boundsmin = [-5.0 -5.0] , boundsmax = [5.0 5.0], nbve=5 , data = mydude );
+[ s1 , mydude ] = optimsimplex_new ( "randbounds" , [-1.2 1.0], mycostf, ...
+  [-5.0 -5.0] , [5.0 5.0], 5 , mydude );
 computed = optimsimplex_getall ( s1 );
 assert_equal ( size(computed,1) , 5 );
 assert_equal ( size(computed,2) , 3 );
 assert_equal ( mydude.nb , 5 );
+expected = [
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+  %T %T  
+];
+assert_equal ( computed(1:5,2:3) < 5.0 , expected );
+assert_equal ( computed(1:5,2:3) > -5.0 , expected );
 s1 = optimsimplex_destroy ( s1 );
 
index 39431d8..62abd9e 100644 (file)
@@ -45,33 +45,47 @@ endfunction
 function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
-myobj = tlist(["T_MYSTUFF","nb"]);
-myobj.nb = 0;
-function [ y , myobj ] = mycostf ( x , myobj )
-  y = rosenbrock(x);
-  myobj.nb = myobj.nb + 1
-endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_spendley ( s1 , x0 = [-1.2 1.0] , fun = mycostf , data=myobj );
+//
+// Spendley basic case
+//
+s1 = optimsimplex_new ( "spendley" , [-1.2 1.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.0   -1.2    1.                        
+    0.0   -0.2340741737109317543997    1.2588190451025207394764  
+    0.0  -0.9411809548974792161147    1.9659258262890682011914  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
+//
+// Spendley basic case
+//
+s1 = optimsimplex_new ( "spendley" , [-1.2 1.0] , rosenbrock );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2   -1.2    1.                        
     146.4913601204771680386   -0.2340741737109317543997    1.2588190451025207394764  
     120.43069965448485447723  -0.9411809548974792161147    1.9659258262890682011914  
 ];
-assert_close ( computed , expected , 10 * %eps );
-assert_equal ( myobj.nb , 3 );
+assert_close ( computed , expected , 1.e-6 );
 s1 = optimsimplex_destroy ( s1 );
 //
-// Pfeffer, case 0.0
+// Spendley, case with object
 //
+myobj = tlist(["T_MYSTUFF","nb"]);
+myobj.nb = 0;
+function [ y , myobj ] = mycostf ( x , myobj )
+  y = rosenbrock(x);
+  myobj.nb = myobj.nb + 1
+endfunction
 s1 = optimsimplex_new ();
-s1 = optimsimplex_spendley ( s1 , x0 = [-1.2 1.0] , fun = rosenbrock );
+[ s1 , myobj ] = optimsimplex_new ( "spendley" , [-1.2 1.0] , mycostf , 1.0 , myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2   -1.2    1.                        
     146.4913601204771680386   -0.2340741737109317543997    1.2588190451025207394764  
     120.43069965448485447723  -0.9411809548974792161147    1.9659258262890682011914  
 ];
-assert_close ( computed , expected , 1.e-6 );
+assert_close ( computed , expected , 10 * %eps );
+assert_equal ( myobj.nb , 3 );
 s1 = optimsimplex_destroy ( s1 );
index 93433ac..22ba645 100644 (file)
@@ -47,36 +47,52 @@ function y = rosenbrock (x)
   y = 100*(x(2)-x(1)^2)^2 + (1-x(1))^2;
 endfunction
 
+//
+// Spendley basic case
+//
+s1 = optimsimplex_new ( "spendley" , [-1.2 1.0] );
+computed = optimsimplex_getall ( s1 );
+expected = [
+    0.0   -1.2    1.                        
+    0.0   -0.2340741737109317543997    1.2588190451025207394764  
+    0.0  -0.9411809548974792161147    1.9659258262890682011914  
+];
+assert_close ( computed , expected , 1.e-6 );
+s1 = optimsimplex_destroy ( s1 );
 
-myobj = tlist(["T_MYSTUFF","nb"]);
-myobj.nb = 0;
-function [ y , myobj ] = mycostf ( x , myobj )
-  y = rosenbrock(x);
-  myobj.nb = myobj.nb + 1
-endfunction
-s1 = optimsimplex_new ();
-[ s1 , myobj ] = optimsimplex_spendley ( s1 , x0 = [-1.2 1.0] , fun = mycostf , data=myobj );
+//
+// Spendley basic case
+//
+s1 = optimsimplex_new ( "spendley" , [-1.2 1.0] , rosenbrock );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2   -1.2    1.                        
     146.4913601204771680386   -0.2340741737109317543997    1.2588190451025207394764  
     120.43069965448485447723  -0.9411809548974792161147    1.9659258262890682011914  
 ];
-assert_close ( computed , expected , 10 * %eps );
-assert_equal ( myobj.nb , 3 );
+assert_close ( computed , expected , 1.e-6 );
 s1 = optimsimplex_destroy ( s1 );
 
 //
-// Pfeffer, case 0.0
+// Spendley, case with object
 //
+
+myobj = tlist(["T_MYSTUFF","nb"]);
+myobj.nb = 0;
+function [ y , myobj ] = mycostf ( x , myobj )
+  y = rosenbrock(x);
+  myobj.nb = myobj.nb + 1
+endfunction
 s1 = optimsimplex_new ();
-s1 = optimsimplex_spendley ( s1 , x0 = [-1.2 1.0] , fun = rosenbrock );
+[ s1 , myobj ] = optimsimplex_new ( "spendley" , [-1.2 1.0] , mycostf , 1.0 , myobj );
 computed = optimsimplex_getall ( s1 );
 expected = [
     24.2   -1.2    1.                        
     146.4913601204771680386   -0.2340741737109317543997    1.2588190451025207394764  
     120.43069965448485447723  -0.9411809548974792161147    1.9659258262890682011914  
 ];
-assert_close ( computed , expected , 1.e-6 );
+assert_close ( computed , expected , 10 * %eps );
+assert_equal ( myobj.nb , 3 );
 s1 = optimsimplex_destroy ( s1 );
 
+
index 1dddb19..a981ffc 100644 (file)
@@ -42,6 +42,7 @@ function flag = assert_equal ( computed , expected )
   end
   if flag <> 1 then bugmes();quit;end
 endfunction
+// Test with 3 vertices
 s1 = optimsimplex_new ();
 simplex = [
 24. -2.0 1.0
@@ -50,7 +51,44 @@ simplex = [
 ];
 s1 = optimsimplex_setall ( s1 , simplex );
 str = optimsimplex_tostring ( s1 );
-assert_equal ( str(1) , "Vertex #1/3 : fv=24.000000, x=-2.000000 1.000000" );
-assert_equal ( str(2) , "Vertex #2/3 : fv=93.000000, x=-1.000000 3.000000" );
-assert_equal ( str(3) , "Vertex #3/3 : fv=36.000000, x=-3.000000 2.000000" );
+expected = [
+"Vertex #1/3 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000"
+"Vertex #2/3 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000"
+"Vertex #3/3 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000"
+]
+ expected  =
+!Vertex #1/3 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000  !
+!                                                                !
+!Vertex #2/3 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000  !
+!                                                                !
+!Vertex #3/3 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000  !
+assert_equal ( str , expected );
+s1 = optimsimplex_destroy ( s1 );
+// Test with 4 vertices
+s1 = optimsimplex_new ();
+simplex = [
+24. -2.0 1.0
+93. -1.0 3.0
+36. -3.0 2.0
+36. -3.0 2.0
+];
+s1 = optimsimplex_setall ( s1 , simplex );
+str = optimsimplex_tostring ( s1 );
+expected = [
+"Vertex #1/4 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000" 
+"Vertex #2/4 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000" 
+"Vertex #3/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000" 
+"Vertex #4/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000" 
+]
+ expected  =
+!Vertex #1/4 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000  !
+!                                                                !
+!Vertex #2/4 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000  !
+!                                                                !
+!Vertex #3/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000  !
+!                                                                !
+!Vertex #4/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000  !
+assert_equal ( str , expected );
 s1 = optimsimplex_destroy ( s1 );
index 6f89e1d..ec3ab6d 100644 (file)
@@ -53,9 +53,12 @@ simplex = [
 ];
 s1 = optimsimplex_setall ( s1 , simplex );
 str = optimsimplex_tostring ( s1 );
-assert_equal ( str(1) , "Vertex #1/3 : fv=24.000000, x=-2.000000 1.000000" );
-assert_equal ( str(2) , "Vertex #2/3 : fv=93.000000, x=-1.000000 3.000000" );
-assert_equal ( str(3) , "Vertex #3/3 : fv=36.000000, x=-3.000000 2.000000" );
+expected = [
+"Vertex #1/3 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000"
+"Vertex #2/3 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000"
+"Vertex #3/3 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000"
+]
+assert_equal ( str , expected );
 s1 = optimsimplex_destroy ( s1 );
 
 // Test with 4 vertices
@@ -68,9 +71,12 @@ simplex = [
 ];
 s1 = optimsimplex_setall ( s1 , simplex );
 str = optimsimplex_tostring ( s1 );
-assert_equal ( str(1) , "Vertex #1/4 : fv=24.000000, x=-2.000000 1.000000" );
-assert_equal ( str(2) , "Vertex #2/4 : fv=93.000000, x=-1.000000 3.000000" );
-assert_equal ( str(3) , "Vertex #3/4 : fv=36.000000, x=-3.000000 2.000000" );
-assert_equal ( str(4) , "Vertex #4/4 : fv=36.000000, x=-3.000000 2.000000" );
+expected = [
+"Vertex #1/4 : fv=2.400000e+001, x=-2.000000e+000 1.000000e+000" 
+"Vertex #2/4 : fv=9.300000e+001, x=-1.000000e+000 3.000000e+000" 
+"Vertex #3/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000" 
+"Vertex #4/4 : fv=3.600000e+001, x=-3.000000e+000 2.000000e+000" 
+]
+assert_equal ( str , expected );
 s1 = optimsimplex_destroy ( s1 );
 
index 83bcc4e..f029378 100644 (file)
@@ -53,7 +53,7 @@ s1 = optimsimplex_setall ( s1 , simplex );
 cen = optimsimplex_xbar ( s1 );
 assert_close ( cen , [-1.2 1.5], %eps );
 s1 = optimsimplex_destroy ( s1 );
-// iexcl =2
+// Case iexcl = 2
 s1 = optimsimplex_new ();
 simplex = [
 24.2 -1.2 1.0
@@ -64,8 +64,19 @@ s1 = optimsimplex_setall ( s1 , simplex );
 cen = optimsimplex_xbar ( s1 , 2 );
 assert_close ( cen , [-0.7 1.0], %eps );
 s1 = optimsimplex_destroy ( s1 );
+// Case iexcl = 2:3
+s1 = optimsimplex_new ();
+simplex = [
+24.2 -1.2 1.0
+36.2 -1.2 2.0
+93.6 -0.2 1.0
+];
+s1 = optimsimplex_setall ( s1 , simplex );
+cen = optimsimplex_xbar ( s1 , 2:3 );
+assert_close ( cen , [-1.2 1.0], %eps );
+s1 = optimsimplex_destroy ( s1 );
 //
-// Test with 5 vertices
+// Test with 5 vertices and default exclusion
 //
 s1 = optimsimplex_new ();
 simplex = [
@@ -79,3 +90,18 @@ s1 = optimsimplex_setall ( s1 , simplex );
 cen = optimsimplex_xbar ( s1 );
 assert_close ( cen , [-0.65 1.0], %eps );
 s1 = optimsimplex_destroy ( s1 );
+//
+// Test with 5 vertices and several exclusions
+//
+s1 = optimsimplex_new ();
+simplex = [
+24.2 -1.2 1.0
+36.2 -1.2 2.0
+93.6 -0.2 1.0
+93.6 0.0 0.0
+93.6 10.0 10.0
+];
+s1 = optimsimplex_setall ( s1 , simplex );
+cen = optimsimplex_xbar ( s1 , [1 3 5]);
+assert_close ( cen , [-0.6    1.0], 10 * %eps );
+s1 = optimsimplex_destroy ( s1 );
index 0bbd90e..375b6f1 100644 (file)
@@ -54,7 +54,8 @@ s1 = optimsimplex_setall ( s1 , simplex );
 cen = optimsimplex_xbar ( s1 );
 assert_close ( cen , [-1.2 1.5], %eps );
 s1 = optimsimplex_destroy ( s1 );
-// iexcl =2
+
+// Case iexcl = 2
 s1 = optimsimplex_new ();
 simplex = [
 24.2 -1.2 1.0
@@ -65,8 +66,20 @@ s1 = optimsimplex_setall ( s1 , simplex );
 cen = optimsimplex_xbar ( s1 , 2 );
 assert_close ( cen , [-0.7 1.0], %eps );
 s1 = optimsimplex_destroy ( s1 );
+
+// Case iexcl = 2:3
+s1 = optimsimplex_new ();
+simplex = [
+24.2 -1.2 1.0
+36.2 -1.2 2.0
+93.6 -0.2 1.0
+];
+s1 = optimsimplex_setall ( s1 , simplex );
+cen = optimsimplex_xbar ( s1 , 2:3 );
+assert_close ( cen , [-1.2 1.0], %eps );
+s1 = optimsimplex_destroy ( s1 );
 //
-// Test with 5 vertices
+// Test with 5 vertices and default exclusion
 //
 s1 = optimsimplex_new ();
 simplex = [
@@ -81,3 +94,19 @@ cen = optimsimplex_xbar ( s1 );
 assert_close ( cen , [-0.65 1.0], %eps );
 s1 = optimsimplex_destroy ( s1 );
 
+//
+// Test with 5 vertices and several exclusions
+//
+s1 = optimsimplex_new ();
+simplex = [
+24.2 -1.2 1.0
+36.2 -1.2 2.0
+93.6 -0.2 1.0
+93.6 0.0 0.0
+93.6 10.0 10.0
+];
+s1 = optimsimplex_setall ( s1 , simplex );
+cen = optimsimplex_xbar ( s1 , [1 3 5]);
+assert_close ( cen , [-0.6    1.0], 10 * %eps );
+s1 = optimsimplex_destroy ( s1 );
+