00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 int count = 0;
00034 int number_restarts = 0;
00035
00036
00037
00038
00039 void allocate_arrays (float *** simplex, float ** centroid,
00040 float ** response, float ** step_size,
00041 float ** test1, float ** test2)
00042 {
00043 int i;
00044
00045 *centroid = (float *) malloc (sizeof(float) * DIMENSION);
00046 *response = (float *) malloc (sizeof(float) * (DIMENSION+1));
00047 *step_size = (float *) malloc (sizeof(float) * DIMENSION);
00048 *test1 = (float *) malloc (sizeof(float) * DIMENSION);
00049 *test2 = (float *) malloc (sizeof(float) * DIMENSION);
00050
00051 *simplex = (float **) malloc (sizeof(float *) * (DIMENSION+1));
00052
00053 for (i = 0; i < DIMENSION+1; i++)
00054 (*simplex)[i] = (float *) malloc (sizeof(float) * DIMENSION);
00055
00056 }
00057
00058
00059
00060
00061 void deallocate_arrays (float *** simplex, float ** centroid,
00062 float ** response, float ** step_size,
00063 float ** test1, float ** test2)
00064 {
00065 int i;
00066
00067 free (*centroid); *centroid = NULL;
00068 free (*response); *response = NULL;
00069 free (*step_size); *step_size = NULL;
00070 free (*test1); *test1 = NULL;
00071 free (*test2); *test2 = NULL;
00072
00073 for (i = 0; i < DIMENSION+1; i++)
00074 {
00075 free ((*simplex)[i]);
00076 (*simplex)[i] = NULL;
00077 }
00078
00079 free (*simplex); *simplex = NULL;
00080
00081 }
00082
00083
00084
00085
00086 void eval_vertices (float * response, int * worst, int * next, int * best)
00087 {
00088 int i;
00089
00090
00091 *worst = 0;
00092 *best = 0;
00093
00094
00095 for (i = 1; i < DIMENSION+1; i++)
00096 {
00097 if (response[i] > response[*worst])
00098 *worst = i;
00099 if (response[i] < response[*best])
00100 *best = i;
00101 }
00102
00103
00104 if (*worst == 0)
00105 *next = 1;
00106 else
00107 *next = 0;
00108
00109 for (i = 0; i < DIMENSION+1; i++)
00110 if ((i != *worst) && (response[i] > response[*next]))
00111 *next = i;
00112 }
00113
00114
00115
00116
00117 void restart (float ** simplex, float * response,
00118 float * step_size)
00119 {
00120 const float STEP_FACTOR = 0.9;
00121 int i, j;
00122 int worst, next, best;
00123 float minval, maxval;
00124
00125
00126
00127 eval_vertices (response, &worst, &next, &best);
00128
00129
00130 for (i = 0; i < DIMENSION; i++)
00131 simplex[0][i] = simplex[best][i];
00132
00133
00134 for (i = 0; i < DIMENSION; i++)
00135 step_size[i] *= STEP_FACTOR;
00136
00137
00138 for (i = 1; i < DIMENSION+1; i++)
00139 for (j = 0; j < DIMENSION; j++)
00140 {
00141 minval = simplex[0][j] - step_size[j];
00142 maxval = simplex[0][j] + step_size[j];
00143 simplex[i][j] = rand_uniform (minval, maxval);
00144 }
00145
00146
00147 for (i = 0; i < DIMENSION+1; i++)
00148 response[i] = calc_error (simplex[i]);
00149 }
00150
00151
00152
00153
00154
00155
00156
00157 void calc_centroid (float ** simplex, int worst, float * centroid)
00158 {
00159 int i, j;
00160
00161 for (i = 0; i < DIMENSION; i++)
00162 {
00163 centroid[i] = 0.0;
00164
00165
00166 for (j = 0; j < DIMENSION+1; j++)
00167 if (j != worst)
00168 centroid[i] += simplex[j][i];
00169 }
00170
00171
00172 for (i = 0; i < DIMENSION; i++)
00173 centroid[i] /= DIMENSION;
00174 }
00175
00176
00177
00178
00179
00180
00181
00182 void calc_reflection (float ** simplex, float * centroid,
00183 int worst, float coef, float * vertex)
00184 {
00185 int i;
00186
00187 for (i = 0; i < DIMENSION; i++)
00188 vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
00189 }
00190
00191
00192
00193
00194
00195
00196
00197 void replace (float ** simplex, float * response,
00198 int index, float * vertex, float resp)
00199 {
00200 int i;
00201
00202 for (i = 0; i < DIMENSION; i++)
00203 simplex[index][i] = vertex[i];
00204
00205 response[index] = resp;
00206 }
00207
00208
00209
00210
00211
00212
00213
00214 float calc_good_fit (float * response)
00215 {
00216 int i;
00217
00218 float avg, sd, tmp;
00219
00220
00221 avg = 0.0;
00222 for (i = 0; i < DIMENSION+1; i++)
00223 avg += response[i];
00224 avg /= DIMENSION+1;
00225
00226
00227 sd = 0.0;
00228 for (i = 0; i < DIMENSION+1; i++)
00229 {
00230 tmp = response[i] - avg;
00231 sd += tmp*tmp;
00232 }
00233 sd /= DIMENSION;
00234
00235 return (sqrt(sd));
00236 }
00237
00238
00239
00240
00241
00242
00243
00244 void simplex_initialize (float * parameters, float ** simplex,
00245 float * response, float * step_size)
00246 {
00247 int i, j;
00248 int worst, next, best;
00249 float resp;
00250 float minval, maxval;
00251
00252
00253 for (j = 0; j < DIMENSION; j++)
00254 {
00255 simplex[0][j] = parameters[j];
00256 step_size[j] = 0.5 * parameters[j];
00257 }
00258
00259 for (i = 1; i < DIMENSION+1; i++)
00260 for (j = 0; j < DIMENSION; j++)
00261 {
00262 minval = simplex[0][j] - step_size[j];
00263 maxval = simplex[0][j] + step_size[j];
00264 simplex[i][j] = rand_uniform (minval, maxval);
00265 }
00266
00267 for (i = 0; i < DIMENSION+1; i++)
00268 response[i] = calc_error (simplex[i]);
00269
00270 for (i = 1; i < 500; i++)
00271 {
00272 for (j = 0; j < DIMENSION; j++)
00273 {
00274 minval = simplex[0][j] - step_size[j];
00275 maxval = simplex[0][j] + step_size[j];
00276 parameters[j] = rand_uniform (minval, maxval);
00277 }
00278
00279 resp = calc_error (parameters);
00280 eval_vertices (response, &worst, &next, &best);
00281 if (resp < response[worst])
00282 replace (simplex, response, worst, parameters, resp);
00283 }
00284 }
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294 void simplex_optimization (float * parameters, float * sse)
00295 {
00296 const int MAX_ITERATIONS = 100;
00297 const int MAX_RESTARTS = 25;
00298 const float EXPANSION_COEF = 2.0;
00299 const float REFLECTION_COEF = 1.0;
00300 const float CONTRACTION_COEF = 0.5;
00301 const float TOLERANCE = 1.0e-10;
00302
00303 float ** simplex = NULL;
00304 float * centroid = NULL;
00305 float * response = NULL;
00306 float * step_size = NULL;
00307 float * test1 = NULL;
00308 float * test2 = NULL;
00309 float resp1, resp2;
00310 int i, worst, best, next;
00311 int num_iter, num_restarts;
00312 int done;
00313 float fit;
00314
00315
00316 allocate_arrays (&simplex, ¢roid, &response, &step_size, &test1, &test2);
00317
00318 simplex_initialize (parameters, simplex, response, step_size);
00319
00320
00321 num_iter = 0;
00322 num_restarts = 0;
00323 done = 0;
00324
00325 while (!done)
00326 {
00327
00328
00329 eval_vertices (response, &worst, &next, &best);
00330 calc_centroid (simplex, worst, centroid);
00331
00332
00333 calc_reflection (simplex, centroid, worst,
00334 REFLECTION_COEF, test1);
00335 resp1 = calc_error (test1);
00336
00337
00338
00339 if (resp1 < response[best])
00340 {
00341
00342 calc_reflection (simplex, centroid, worst, EXPANSION_COEF, test2);
00343 resp2 = calc_error (test2);
00344 if (resp2 <= resp1)
00345 replace (simplex, response, worst, test2, resp2);
00346 else
00347 replace (simplex, response, worst, test1, resp1);
00348 }
00349 else if (resp1 < response[next])
00350 {
00351
00352
00353 replace (simplex, response, worst, test1, resp1);
00354 }
00355 else
00356 {
00357
00358 if (resp1 >= response[worst])
00359 calc_reflection (simplex, centroid, worst,
00360 -CONTRACTION_COEF, test2);
00361 else
00362 calc_reflection (simplex, centroid, worst,
00363 CONTRACTION_COEF, test2);
00364 resp2 = calc_error (test2);
00365
00366
00367 if (resp2 > response[worst])
00368 {
00369
00370
00371 num_iter = 0;
00372 num_restarts += 1;
00373 restart (simplex, response, step_size);
00374 }
00375 else
00376 replace (simplex, response, worst, test2, resp2);
00377 }
00378
00379
00380
00381 num_iter += 1;
00382 if (num_iter >= MAX_ITERATIONS)
00383 {
00384
00385 num_iter = 0;
00386 num_restarts += 1;
00387 restart (simplex, response, step_size);
00388 }
00389
00390
00391 if (num_restarts == MAX_RESTARTS) done = 1;
00392
00393
00394
00395 fit = calc_good_fit (response);
00396 if (fit <= TOLERANCE) done = 1;
00397
00398
00399 if (done)
00400 {
00401 eval_vertices (response, &worst, &next, &best);
00402 for (i = 0; i < DIMENSION; i++)
00403 parameters[i] = simplex[best][i];
00404 *sse = response[best];
00405 }
00406
00407 }
00408
00409 number_restarts = num_restarts;
00410 deallocate_arrays (&simplex, ¢roid, &response, &step_size,
00411 &test1, &test2);
00412
00413 }
00414
00415
00416