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 }

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 }

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 }

00414 {
00415   int i;
00416 
00417   for (i = 0;  i < dimension;  i++)
00418     vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
00419 }

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 }

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 }

00078 {
00079   float fval;
00080 
00081   fval = a + uniform() * (b-a);
00082 
00083   return (fval);
00084 }

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 }

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 }

00437 {
00438   int i;
00439 
00440   for (i = 0;  i < dimension;  i++)
00441     simplex[index][i] = vertex[i];
00442 
00443   response[index] = resp;
00444 }

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 }

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 }

00043 {
00044   return ( (float)drand48() );
00045 }