00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008    This file implements the Simplex algorithm for non-linear optimization.
00009    Some routines are adapted from:  A Jump Start Course in C++ Programming
00010 
00011   File:     simplex.c
00012   Author:   B. Douglas Ward
00013   Date:     23 May 1997
00014 
00015   Mod:      Changed random number generator function from rand to drand48.
00016             29 August 1997
00017 */
00018 
00019 
00020 /*---------------------------------------------------------------------------*/
00021 /*
00022   This software is Copyright 1997 by
00023 
00024             Medical College of Wisconsin
00025             8701 Watertown Plank Road
00026             Milwaukee, WI 53226
00027 
00028   License is granted to use this program for nonprofit research purposes only.
00029   It is specifically against the license to use this program for any clinical
00030   application. The Medical College of Wisconsin makes no warranty of usefulness
00031   of this program for any particular purpose.  The redistribution of this
00032   program for a fee, or the derivation of for-profit works from this program
00033   is not allowed.
00034 */
00035 
00036 
00037 /*---------------------------------------------------------------------------*/
00038 /*
00039   Routine to generate a uniform U(0,1) random variate.
00040 */
00041 
00042 float uniform ()
00043 {
00044   return ( (float)drand48() );
00045 }
00046 
00047 
00048 /*---------------------------------------------------------------------------*/
00049 /*
00050   Routine to generate a normal N(0,1) random variate.
00051 */
00052 
00053 void normal (float * n1, float * n2)
00054 {
00055   float u1, u2;
00056   float r;
00057 
00058 
00059   u1 = 0.0;
00060   while (u1 <= 0.0)
00061     {
00062       u1 = uniform();
00063     }
00064   u2 = uniform();
00065 
00066   r   = sqrt(-2.0*log(u1));
00067   *n1 = r * cos(2.0*PI*u2);
00068   *n2 = r * sin(2.0*PI*u2);
00069 }
00070 
00071 
00072 /*---------------------------------------------------------------------------*/
00073 /*
00074   Routine to generate a U(a,b) random variate.
00075 */
00076 
00077 float get_random_value(float a, float b)
00078 {
00079   float fval;
00080 
00081   fval = a + uniform() * (b-a);
00082 
00083   return (fval);
00084 }
00085 
00086 
00087 /*---------------------------------------------------------------------------*/
00088 /*
00089   Allocate memory for simplex algorithm.
00090 */
00091 
00092 void allocate_arrays 
00093 (
00094   int dimension,              /* dimension of parameter space */
00095   float *** simplex,          /* the simplex itself */
00096   float ** centroid,          /* center of mass of the simplex */
00097   float ** response,          /* error sum of squares at each vertex */
00098   float ** step_size,         /* controls random placement of new vertex */
00099   float ** test1,             /* test vertex */
00100   float ** test2              /* test vertex */
00101 )
00102 
00103 {
00104   int i;
00105 
00106   *centroid = (float *) malloc (sizeof(float) * dimension);
00107   *step_size = (float *) malloc (sizeof(float) * dimension);
00108   *test1 = (float *) malloc (sizeof(float) * dimension);
00109   *test2 = (float *) malloc (sizeof(float) * dimension);
00110 
00111   *response = (float *) malloc (sizeof(float) * (dimension+1));   
00112   *simplex = (float **) malloc (sizeof(float *) * (dimension+1));
00113 
00114   for (i = 0;  i < dimension+1;  i++)
00115     (*simplex)[i] = (float *) malloc (sizeof(float) * dimension);
00116        
00117 }
00118 
00119 
00120 /*---------------------------------------------------------------------------*/
00121 /*
00122   Set up initial values for the simplex vertices.
00123 */
00124 
00125 void initialize_simplex 
00126 (
00127   int dimension,          /* dimension of the full model */
00128   vfp nmodel,             /* pointer to noise model */
00129   vfp smodel,             /* pointer to signal model */
00130   int r,                  /* number of parameters in the noise model */
00131   int p,                  /* number of parameters in the signal model */
00132   int nabs,               /* use absolute constraints for noise parameters */
00133   float * min_nconstr,    /* minimum parameter constraints for noise model */
00134   float * max_nconstr,    /* maximum parameter constraints for noise model */
00135   float * min_sconstr,    /* minimum parameter constraints for signal model */
00136   float * max_sconstr,    /* maximum parameter constraints for signal model */
00137   float * par_rdcd,       /* estimated parameters for the reduced model */
00138   float * parameters,     /* starting point */
00139   float ** simplex,       /* the simplex itself */
00140   float * response,       /* sse at each vertex of the simplex */
00141   float * step_size,      /* amount of allowed variation at each parameter */
00142   int ts_length,          /* length of time series array */
00143   float ** x_array,       /* independent variable matrix */
00144   float * ts_array        /* observed time series */
00145 )
00146 
00147 {
00148   int i, j;
00149   float minval, maxval;
00150 
00151 
00152   /*----- copy parameter vector into first vertex of simplex -----*/
00153   for (i = 0;  i < dimension;  i++)    
00154     simplex[0][i] = parameters[i];
00155 
00156     
00157   /*----- set up initial step sizes -----*/
00158   for (i = 0;  i < r;  i++)
00159     step_size[i] = 0.1 * (max_nconstr[i] - min_nconstr[i]);
00160   for (i = r;  i < dimension;  i++)
00161     step_size[i] = 0.1 * (max_sconstr[i-r] - min_sconstr[i-r]);
00162 
00163     
00164   /*----- choose random vectors for remaining vertices -----*/
00165   for (i = 1;  i < dimension+1;  i++)
00166     {
00167       /*----- choose noise model parameters -----*/
00168       for (j = 0;  j < r;  j++)
00169         {
00170           minval = simplex[0][j] - step_size[j];
00171           if (nabs)   /*--- absolute noise parameter constraints ---*/
00172             {
00173               if (minval < min_nconstr[j])  
00174                 minval = min_nconstr[j];
00175             }
00176           else        /*--- relative noise parameter constraints ---*/
00177             {
00178               if (minval < min_nconstr[j] + par_rdcd[j])  
00179                 minval = min_nconstr[j] + par_rdcd[j];
00180             }
00181 
00182           maxval = simplex[0][j] + step_size[j];
00183           if (nabs)   /*--- absolute noise parameter constraints ---*/
00184             {
00185               if (maxval > max_nconstr[j])  
00186                 maxval = max_nconstr[j];
00187             }
00188           else        /*--- relative noise parameter constraints ---*/
00189             {
00190               if (maxval > max_nconstr[j] + par_rdcd[j])  
00191                 maxval = max_nconstr[j] + par_rdcd[j];
00192             }
00193           simplex[i][j] = get_random_value (minval, maxval);
00194         }
00195 
00196 
00197       /*----- choose signal model parameters -----*/
00198       for (j = r;  j < dimension;  j++)
00199         {
00200           minval = simplex[0][j] - step_size[j];
00201           if (minval < min_sconstr[j-r])  
00202             minval = min_sconstr[j-r];
00203           maxval = simplex[0][j] + step_size[j];
00204           if (maxval > max_sconstr[j-r])  
00205             maxval = max_sconstr[j-r];
00206           simplex[i][j] = get_random_value (minval, maxval);
00207         }
00208     }
00209 
00210 
00211   /*----- calculate and save sse for each vertex of simplex -----*/
00212   for (i = 0;  i < dimension+1;  i++)
00213     response[i] = calc_sse(nmodel, smodel, r, p, nabs, 
00214                            min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00215                            par_rdcd, simplex[i], ts_length, x_array, ts_array);
00216 
00217 }
00218 
00219 
00220 /*---------------------------------------------------------------------------*/
00221 /*
00222   Evaluate the vertices of the simplex.  Find indices of the best, worst, and 
00223   next-to-worst vertices.
00224 */
00225 
00226 void eval_vertices 
00227 (
00228  int dimension,            /* dimension of parameter space */
00229  float * response,         /* error sum of squares at each vertex */
00230  int * worst,              /* index of worst vertex in simplex */
00231  int * next,               /* index of next-to-worst vertex in simplex */
00232  int * best                /* index of best vertex in simplex */
00233 )
00234 
00235 {
00236   int i;
00237 
00238   /* initialize values */
00239   *worst = 0;
00240   *best = 0;
00241 
00242   /* find the best and worst */
00243   for (i = 1;  i < dimension+1;  i++)
00244     {
00245       if (response[i] > response[*worst])
00246         *worst = i;
00247       if (response[i] < response[*best])
00248         *best = i;
00249     }
00250 
00251   /* find the next worst index */
00252   if (*worst == 0)
00253     *next = 1;
00254   else
00255     *next = 0;
00256   
00257   for (i = 0;  i < dimension+1;  i++)
00258     if ((i != *worst) && (response[i] > response[*next]))
00259       *next = i;
00260 }
00261 
00262 
00263 /*---------------------------------------------------------------------------*/
00264 /*
00265   Routine to restart the simplex algorithm by putting random vectors into
00266   the simplex (and keeping the best previous vertex).
00267 */
00268 
00269 void restart 
00270 (
00271   int dimension,          /* dimension of the full model */
00272   vfp nmodel,             /* pointer to noise model */
00273   vfp smodel,             /* pointer to signal model */
00274   int r,                  /* number of parameters in the noise model */
00275   int p,                  /* number of parameters in the signal model */
00276   int nabs,               /* use absolute constraints for noise parameters */
00277   float * min_nconstr,    /* minimum parameter constraints for noise model */
00278   float * max_nconstr,    /* maximum parameter constraints for noise model */
00279   float * min_sconstr,    /* minimum parameter constraints for signal model */
00280   float * max_sconstr,    /* maximum parameter constraints for signal model */
00281   float * par_rdcd,       /* estimated parameters for the reduced model */
00282   float ** simplex,       /* the simplex itself */
00283   float * response,       /* sse at each vertex of the simplex */
00284   float * step_size,      /* amount of allowed variation at each parameter */
00285   int ts_length,          /* length of time series array */
00286   float ** x_array,       /* independent variable matrix */
00287   float * ts_array        /* observed time series */
00288 )
00289 
00290 {
00291   const float STEP_FACTOR = 0.9;
00292   int i, j;
00293   int worst, next, best;
00294   float minval, maxval;
00295 
00296 
00297   /* find the current best vertex */
00298   eval_vertices (dimension, response, &worst, &next, &best); 
00299 
00300 
00301   /* set the first vertex to the current best */
00302   for (i = 0; i < dimension;  i++)
00303     simplex[0][i] = simplex[best][i];
00304 
00305 
00306   /* decrease step size */
00307   for (i = 0;  i < dimension;  i++)
00308     step_size[i] *= STEP_FACTOR;
00309 
00310 
00311   /* set up remaining vertices of simplex using new step size */
00312   for (i = 1;  i < dimension+1;  i++)
00313     {
00314       /*----- choose noise model parameters -----*/
00315       for (j = 0;  j < r;  j++)
00316         {
00317           minval = simplex[0][j] - step_size[j];
00318           if (nabs)   /*--- absolute noise parameter constraints ---*/
00319             {
00320               if (minval < min_nconstr[j])  
00321                 minval = min_nconstr[j];
00322             }
00323           else        /*--- relative noise parameter constraints ---*/
00324             {
00325               if (minval < min_nconstr[j] + par_rdcd[j])  
00326                 minval = min_nconstr[j] + par_rdcd[j];
00327             }
00328 
00329           maxval = simplex[0][j] + step_size[j];
00330           if (nabs)
00331             {
00332               if (maxval > max_nconstr[j])  
00333                 maxval = max_nconstr[j];
00334             }
00335           else
00336             {
00337               if (maxval > max_nconstr[j] + par_rdcd[j])  
00338                 maxval = max_nconstr[j] + par_rdcd[j];
00339             }
00340 
00341           simplex[i][j] = get_random_value (minval, maxval);
00342         }
00343 
00344 
00345       /*----- choose signal model parameters -----*/
00346       for (j = r;  j < dimension;  j++)
00347         {
00348           minval = simplex[0][j] - step_size[j];
00349           if (minval < min_sconstr[j-r])  
00350             minval = min_sconstr[j-r];
00351           maxval = simplex[0][j] + step_size[j];
00352           if (maxval > max_sconstr[j-r])  
00353             maxval = max_sconstr[j-r];
00354           simplex[i][j] = get_random_value (minval, maxval);
00355         }
00356     }
00357 
00358   
00359   /* initialize response for each vector */
00360   for (i = 0;  i < dimension+1;  i++)
00361     response[i] = calc_sse (nmodel, smodel, r, p, nabs, 
00362                            min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00363                            par_rdcd, simplex[i], ts_length, x_array, ts_array);
00364 }
00365 
00366 
00367 /*---------------------------------------------------------------------------*/
00368 /*
00369   Calculate the centroid of the simplex, ignoring the worst vertex.
00370 */
00371 
00372 void calc_centroid 
00373 (
00374   int dimension,         /* dimension of parameter space */
00375   float ** simplex,      /* the simplex itself */
00376   int worst,             /* index of worst vertex in simplex */
00377   float * centroid       /* center of mass of the simplex minus worst vertex */
00378 )
00379 
00380 {
00381   int i, j;
00382 
00383   for (i = 0;  i < dimension;  i++)
00384     {
00385       centroid[i] = 0.0;
00386 
00387       /* add each vertex, except the worst */
00388       for (j = 0;  j < dimension+1;  j++)
00389         if (j != worst)
00390           centroid[i] += simplex[j][i];
00391     }
00392 
00393   /* divide by the number of vertices */
00394   for (i = 0;  i < dimension;  i++)
00395     centroid[i] /= dimension;
00396 }
00397 
00398 
00399 /*---------------------------------------------------------------------------*/
00400 /*
00401   Calculate the reflection of the worst vertex about the centroid.
00402 */
00403 
00404 void calc_reflection 
00405 (
00406   int dimension,               /* dimension of parameter space */
00407   float ** simplex,            /* the simplex itself */
00408   float * centroid,            /* center of mass of the simplex */
00409   int worst,                   /* index of worst vertex in simplex */
00410   float coef,                  /* expansion or contraction factor */
00411   float * vertex               /* new vertex */
00412 )
00413 
00414 {
00415   int i;
00416 
00417   for (i = 0;  i < dimension;  i++)
00418     vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
00419 }
00420 
00421 
00422 /*---------------------------------------------------------------------------*/
00423 /*
00424   Replace a vertex of the simplex.
00425 */
00426 
00427 void replace 
00428 (
00429   int dimension,              /* dimension of parameter space */
00430   float ** simplex,           /* the simplex itself */
00431   float * response,           /* error sum of squares at each vertex */
00432   int index,                  /* index of vertex to be replaced */
00433   float * vertex,             /* new vertex */
00434   float resp                  /* error sum of squares at new vertex */
00435 )
00436 
00437 {
00438   int i;
00439 
00440   for (i = 0;  i < dimension;  i++)
00441     simplex[index][i] = vertex[i];
00442 
00443   response[index] = resp;
00444 }
00445 
00446 
00447 /*---------------------------------------------------------------------------*/
00448 /*
00449   Calculate the goodness of fit.  This is measured by the variation in 
00450   responses at the different vertices relative to the average response.
00451 */
00452 
00453 float calc_good_fit 
00454 (
00455   int dimension,                   /* dimension of parameter space */
00456   float * response                 /* error sum of squares at each vertex */
00457 )
00458 
00459 {
00460   int i;
00461 
00462   float avg, sd, tmp;
00463 
00464   /*----- average the response values -----*/
00465   avg = 0.0;
00466   for (i = 0;  i < dimension+1;  i++)
00467     avg += response[i];
00468   avg /= dimension+1;
00469 
00470   /*----- compute standard deviation of response -----*/
00471   sd = 0.0;
00472   for (i = 0;  i < dimension+1;  i++)
00473     {
00474       tmp = response[i] - avg;
00475       sd += tmp*tmp;
00476     }
00477 
00478   sd /= dimension;
00479 
00480   return (sqrt(sd) / avg);
00481 }
00482 
00483 
00484 /*---------------------------------------------------------------------------*/
00485 /*
00486   Release memory required for simplex optimization.
00487 */
00488 
00489 void deallocate_arrays 
00490 (
00491   int dimension,              /* dimension of parameter space */
00492   float *** simplex,          /* the simplex itself */
00493   float ** centroid,          /* center of mass of the simplex */
00494   float ** response,          /* error sum of squares at each vertex */
00495   float ** step_size,         /* controls random placement of new vertex */
00496   float ** test1,             /* test vertex */
00497   float ** test2              /* test vertex */
00498 )
00499 
00500 {
00501   int iv;                     /* vertex index */   
00502 
00503 
00504   free (*centroid);    *centroid = NULL;
00505   free (*response);    *response = NULL;
00506   free (*step_size);   *step_size = NULL;
00507   free (*test1);       *test1 = NULL;
00508   free (*test2);       *test2 = NULL;
00509 
00510   for (iv = 0;  iv < dimension+1;  iv++)
00511     {
00512       free ((*simplex)[iv]);
00513       (*simplex)[iv] = NULL;
00514     }
00515 
00516   free (*simplex);     *simplex = NULL;
00517        
00518 }
00519 
00520 
00521 /*---------------------------------------------------------------------------*/
00522 /*
00523   Implementation of the (non-linear) simplex optimization algorithm.
00524 */
00525 
00526 void simplex_optimization
00527 ( 
00528   vfp nmodel,             /* pointer to noise model */
00529   vfp smodel,             /* pointer to signal model */
00530   int r,                  /* number of parameters in the noise model */
00531   int p,                  /* number of parameters in the signal model */
00532   float * min_nconstr,    /* minimum parameter constraints for noise model */
00533   float * max_nconstr,    /* maximum parameter constraints for noise model */
00534   float * min_sconstr,    /* minimum parameter constraints for signal model */
00535   float * max_sconstr,    /* maximum parameter constraints for signal model */
00536   int nabs,               /* use absolute constraints for noise parameters */
00537   int ts_length,          /* length of time series array */
00538   float ** x_array,       /* independent variable matrix */
00539   float * ts_array,       /* observed time series */
00540   float * par_rdcd,       /* estimated parameters for the reduced model */
00541   float * parameters,     /* estimated parameters */
00542   float * sse             /* error sum of squares */
00543 )
00544 
00545 {
00546   const int MAX_ITERATIONS = 50;         /* maximum number of iterations */
00547   const int MAX_RESTARTS = 5;            /* maximum number of restarts */
00548   const float EXPANSION_COEF = 2.0;      /* expansion coefficient */
00549   const float REFLECTION_COEF = 1.0;     /* reflection coefficient */
00550   const float CONTRACTION_COEF = 0.5;    /* contraction coefficient */
00551   const float TOLERANCE = 1.0e-4;        /* solution convergence tolerance */
00552 
00553   float ** simplex   = NULL;    /* the simplex itself */
00554   float * centroid   = NULL;    /* center of mass of the simplex */
00555   float * response   = NULL;    /* error sum of squares at each vertex */
00556   float * step_size  = NULL;    /* controls random placement of new vertex */
00557   float * test1      = NULL;    /* test vertex */
00558   float * test2      = NULL;    /* test vertex */
00559   float resp1, resp2;           /* error sum of squares for test vertex */
00560   int i;                        /* vertex index */ 
00561   int worst;                    /* index of worst vertex in simplex */
00562   int best;                     /* index of best vertex in simplex */
00563   int next;                     /* index of next-to-worst vertex in simplex */ 
00564   int num_iter;                 /* number of simplex algorithm iterations */
00565   int num_restarts;             /* number of restarts of simplex algorithm */
00566   int done;                     /* boolean for search finished */
00567   float fit;                    /* array of fitted time series values */
00568   int dimension;                /* dimension of parameter space */
00569 
00570 
00571   /*----- dimension of parameter space -----*/
00572   dimension = r + p;
00573 
00574   /*----- allocate memory -----*/
00575   allocate_arrays (dimension, &simplex, &centroid, &response, &step_size,
00576                    &test1, &test2);
00577   
00578   /*----- initialization for simplex algorithm -----*/
00579   initialize_simplex (dimension, nmodel, smodel, r, p, nabs,
00580                       min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00581                       par_rdcd, parameters, simplex, response, step_size, 
00582                       ts_length, x_array, ts_array);
00583 
00584   /* start loop to do simplex optimization */
00585   num_iter = 0;
00586   num_restarts = 0;
00587   done = 0;
00588   
00589   while (!done)
00590     {
00591       /*----- find the worst vertex and compute centroid of remaining simplex, 
00592         discarding the worst vertex -----*/
00593       eval_vertices (dimension, response, &worst, &next, &best);
00594       calc_centroid (dimension, simplex, worst, centroid);
00595       
00596       /*----- reflect the worst point through the centroid -----*/
00597       calc_reflection (dimension, simplex, centroid, worst, 
00598                        REFLECTION_COEF, test1);
00599       resp1 = calc_sse (nmodel, smodel, r, p, nabs, min_nconstr, max_nconstr, 
00600                         min_sconstr, max_sconstr, par_rdcd, test1, 
00601                         ts_length, x_array, ts_array);
00602 
00603 
00604       /*----- test the reflection against the best vertex and expand it if the
00605         reflection is better.  if not, keep the reflection -----*/
00606       if (resp1 < response[best])
00607         {
00608           /*----- try expanding -----*/
00609           calc_reflection (dimension, simplex, centroid, worst, 
00610                            EXPANSION_COEF, test2);
00611 
00612           resp2 = calc_sse (nmodel, smodel, r, p, nabs, 
00613                             min_nconstr, max_nconstr, 
00614                             min_sconstr, max_sconstr, par_rdcd, test2, 
00615                             ts_length, x_array, ts_array);
00616 
00617           if (resp2 <= resp1)    /* keep expansion */     
00618             replace (dimension, simplex, response, worst, test2, resp2);
00619           else                   /* keep reflection */
00620             replace (dimension, simplex, response, worst, test1, resp1);
00621         }
00622 
00623       else if (resp1 < response[next])
00624         {
00625           /*----- new response is between the best and next worst 
00626             so keep reflection -----*/
00627           replace (dimension, simplex, response, worst, test1, resp1); 
00628         }
00629           else
00630         {
00631           /*----- try contraction -----*/
00632           if (resp1 >= response[worst])
00633             calc_reflection (dimension, simplex, centroid, worst, 
00634                              -CONTRACTION_COEF, test2);
00635           else
00636             calc_reflection (dimension, simplex, centroid, worst, 
00637                              CONTRACTION_COEF, test2);
00638 
00639           resp2 = calc_sse (nmodel, smodel, r, p, nabs, 
00640                             min_nconstr, max_nconstr, 
00641                             min_sconstr, max_sconstr, par_rdcd, test2, 
00642                             ts_length, x_array, ts_array);
00643           
00644           /*---- test the contracted response against the worst response ----*/
00645           if (resp2 > response[worst])
00646             {
00647               /*----- new contracted response is worse, so decrease step size
00648                 and restart -----*/
00649               num_iter = 0;
00650               num_restarts += 1;
00651               restart (dimension, nmodel, smodel, r, p, nabs,
00652                        min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00653                        par_rdcd, simplex, response, step_size, 
00654                        ts_length, x_array, ts_array);
00655             }
00656           else       /*----- keep contraction -----*/
00657             replace (dimension, simplex, response, worst, test2, resp2);
00658         }
00659 
00660       /*----- test to determine when to stop.  
00661         first, check the number of iterations -----*/
00662       num_iter += 1;    /*----- increment iteration counter -----*/
00663       if (num_iter >= MAX_ITERATIONS)
00664         {
00665           /*----- restart with smaller steps -----*/
00666           num_iter = 0;
00667           num_restarts += 1;
00668           restart (dimension, nmodel, smodel, r, p, nabs,
00669                    min_nconstr, max_nconstr, min_sconstr, max_sconstr, 
00670                    par_rdcd, simplex, response, step_size, 
00671                    ts_length, x_array, ts_array);
00672         }
00673 
00674       /*----- limit the number of restarts -----*/
00675       if (num_restarts == MAX_RESTARTS)  done = 1;
00676 
00677       /*----- compare relative standard deviation of vertex responses 
00678         against a defined tolerance limit -----*/
00679       fit = calc_good_fit (dimension, response);
00680       if (fit <= TOLERANCE)  done = 1;
00681 
00682       /*----- if done, copy the best solution to the output array -----*/
00683       if (done) 
00684         {
00685           eval_vertices (dimension, response, &worst, &next, &best);
00686           for (i = 0;  i < dimension;  i++)
00687             parameters[i] = simplex[best][i];
00688           *sse = response[best];
00689         }
00690 
00691     }  /*----- while (!done) -----*/
00692  
00693   deallocate_arrays (dimension, &simplex, &centroid, &response, &step_size,
00694                      &test1, &test2);
00695 
00696 }
00697 
00698 
00699 
00700 
00701 
00702 
00703 
00704 
00705 
00706 
00707 
00708 
00709 
00710 
00711 
00712 
00713 
00714 
00715