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
00034
00035
00036
00037
00038
00039
00040
00041
00042 float uniform ()
00043 {
00044 return ( (float)drand48() );
00045 }
00046
00047
00048
00049
00050
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
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
00090
00091
00092 void allocate_arrays
00093 (
00094 int dimension,
00095 float *** simplex,
00096 float ** centroid,
00097 float ** response,
00098 float ** step_size,
00099 float ** test1,
00100 float ** test2
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
00123
00124
00125 void initialize_simplex
00126 (
00127 int dimension,
00128 vfp nmodel,
00129 vfp smodel,
00130 int r,
00131 int p,
00132 int nabs,
00133 float * min_nconstr,
00134 float * max_nconstr,
00135 float * min_sconstr,
00136 float * max_sconstr,
00137 float * par_rdcd,
00138 float * parameters,
00139 float ** simplex,
00140 float * response,
00141 float * step_size,
00142 int ts_length,
00143 float ** x_array,
00144 float * ts_array
00145 )
00146
00147 {
00148 int i, j;
00149 float minval, maxval;
00150
00151
00152
00153 for (i = 0; i < dimension; i++)
00154 simplex[0][i] = parameters[i];
00155
00156
00157
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
00165 for (i = 1; i < dimension+1; i++)
00166 {
00167
00168 for (j = 0; j < r; j++)
00169 {
00170 minval = simplex[0][j] - step_size[j];
00171 if (nabs)
00172 {
00173 if (minval < min_nconstr[j])
00174 minval = min_nconstr[j];
00175 }
00176 else
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)
00184 {
00185 if (maxval > max_nconstr[j])
00186 maxval = max_nconstr[j];
00187 }
00188 else
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
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
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
00223
00224
00225
00226 void eval_vertices
00227 (
00228 int dimension,
00229 float * response,
00230 int * worst,
00231 int * next,
00232 int * best
00233 )
00234
00235 {
00236 int i;
00237
00238
00239 *worst = 0;
00240 *best = 0;
00241
00242
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
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
00266
00267
00268
00269 void restart
00270 (
00271 int dimension,
00272 vfp nmodel,
00273 vfp smodel,
00274 int r,
00275 int p,
00276 int nabs,
00277 float * min_nconstr,
00278 float * max_nconstr,
00279 float * min_sconstr,
00280 float * max_sconstr,
00281 float * par_rdcd,
00282 float ** simplex,
00283 float * response,
00284 float * step_size,
00285 int ts_length,
00286 float ** x_array,
00287 float * ts_array
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
00298 eval_vertices (dimension, response, &worst, &next, &best);
00299
00300
00301
00302 for (i = 0; i < dimension; i++)
00303 simplex[0][i] = simplex[best][i];
00304
00305
00306
00307 for (i = 0; i < dimension; i++)
00308 step_size[i] *= STEP_FACTOR;
00309
00310
00311
00312 for (i = 1; i < dimension+1; i++)
00313 {
00314
00315 for (j = 0; j < r; j++)
00316 {
00317 minval = simplex[0][j] - step_size[j];
00318 if (nabs)
00319 {
00320 if (minval < min_nconstr[j])
00321 minval = min_nconstr[j];
00322 }
00323 else
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
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
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
00370
00371
00372 void calc_centroid
00373 (
00374 int dimension,
00375 float ** simplex,
00376 int worst,
00377 float * centroid
00378 )
00379
00380 {
00381 int i, j;
00382
00383 for (i = 0; i < dimension; i++)
00384 {
00385 centroid[i] = 0.0;
00386
00387
00388 for (j = 0; j < dimension+1; j++)
00389 if (j != worst)
00390 centroid[i] += simplex[j][i];
00391 }
00392
00393
00394 for (i = 0; i < dimension; i++)
00395 centroid[i] /= dimension;
00396 }
00397
00398
00399
00400
00401
00402
00403
00404 void calc_reflection
00405 (
00406 int dimension,
00407 float ** simplex,
00408 float * centroid,
00409 int worst,
00410 float coef,
00411 float * 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
00425
00426
00427 void replace
00428 (
00429 int dimension,
00430 float ** simplex,
00431 float * response,
00432 int index,
00433 float * vertex,
00434 float resp
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
00450
00451
00452
00453 float calc_good_fit
00454 (
00455 int dimension,
00456 float * response
00457 )
00458
00459 {
00460 int i;
00461
00462 float avg, sd, tmp;
00463
00464
00465 avg = 0.0;
00466 for (i = 0; i < dimension+1; i++)
00467 avg += response[i];
00468 avg /= dimension+1;
00469
00470
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
00487
00488
00489 void deallocate_arrays
00490 (
00491 int dimension,
00492 float *** simplex,
00493 float ** centroid,
00494 float ** response,
00495 float ** step_size,
00496 float ** test1,
00497 float ** test2
00498 )
00499
00500 {
00501 int iv;
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
00524
00525
00526 void simplex_optimization
00527 (
00528 vfp nmodel,
00529 vfp smodel,
00530 int r,
00531 int p,
00532 float * min_nconstr,
00533 float * max_nconstr,
00534 float * min_sconstr,
00535 float * max_sconstr,
00536 int nabs,
00537 int ts_length,
00538 float ** x_array,
00539 float * ts_array,
00540 float * par_rdcd,
00541 float * parameters,
00542 float * sse
00543 )
00544
00545 {
00546 const int MAX_ITERATIONS = 50;
00547 const int MAX_RESTARTS = 5;
00548 const float EXPANSION_COEF = 2.0;
00549 const float REFLECTION_COEF = 1.0;
00550 const float CONTRACTION_COEF = 0.5;
00551 const float TOLERANCE = 1.0e-4;
00552
00553 float ** simplex = NULL;
00554 float * centroid = NULL;
00555 float * response = NULL;
00556 float * step_size = NULL;
00557 float * test1 = NULL;
00558 float * test2 = NULL;
00559 float resp1, resp2;
00560 int i;
00561 int worst;
00562 int best;
00563 int next;
00564 int num_iter;
00565 int num_restarts;
00566 int done;
00567 float fit;
00568 int dimension;
00569
00570
00571
00572 dimension = r + p;
00573
00574
00575 allocate_arrays (dimension, &simplex, ¢roid, &response, &step_size,
00576 &test1, &test2);
00577
00578
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
00585 num_iter = 0;
00586 num_restarts = 0;
00587 done = 0;
00588
00589 while (!done)
00590 {
00591
00592
00593 eval_vertices (dimension, response, &worst, &next, &best);
00594 calc_centroid (dimension, simplex, worst, centroid);
00595
00596
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
00605
00606 if (resp1 < response[best])
00607 {
00608
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)
00618 replace (dimension, simplex, response, worst, test2, resp2);
00619 else
00620 replace (dimension, simplex, response, worst, test1, resp1);
00621 }
00622
00623 else if (resp1 < response[next])
00624 {
00625
00626
00627 replace (dimension, simplex, response, worst, test1, resp1);
00628 }
00629 else
00630 {
00631
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
00645 if (resp2 > response[worst])
00646 {
00647
00648
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
00657 replace (dimension, simplex, response, worst, test2, resp2);
00658 }
00659
00660
00661
00662 num_iter += 1;
00663 if (num_iter >= MAX_ITERATIONS)
00664 {
00665
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
00675 if (num_restarts == MAX_RESTARTS) done = 1;
00676
00677
00678
00679 fit = calc_good_fit (dimension, response);
00680 if (fit <= TOLERANCE) done = 1;
00681
00682
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 }
00692
00693 deallocate_arrays (dimension, &simplex, ¢roid, &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