1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 This file contains routines for implementing the Simplex algorithm for
10 non-linear optimization to estimate the unknown parameters.
11
12 The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
13
14 File: Simplex.c
15 Author: B. Douglas Ward
16 Date: 28 January 2000
17
18
19 Note: The calling program should include the following statements
20 prior to "#include "Simplex.c"
21
22 #define DIMENSION p where p = number of parameters to be estimated
23 #include "randgen.c"
24 float calc_error (float * vertex);
25 */
26
27
28 /*---------------------------------------------------------------------------*/
29 /*
30 Global variables and constants.
31 */
32
33 int count = 0; /* count of error evaluations */
34 int number_restarts = 0; /* count of simplex algorithm restarts */
35
36
37 /*---------------------------------------------------------------------------*/
38
allocate_arrays(float *** simplex,float ** centroid,float ** response,float ** step_size,float ** test1,float ** test2)39 void allocate_arrays (float *** simplex, float ** centroid,
40 float ** response, float ** step_size,
41 float ** test1, float ** test2)
42 {
43 int i;
44
45 *centroid = (float *) malloc (sizeof(float) * DIMENSION);
46 *response = (float *) malloc (sizeof(float) * (DIMENSION+1));
47 *step_size = (float *) malloc (sizeof(float) * DIMENSION);
48 *test1 = (float *) malloc (sizeof(float) * DIMENSION);
49 *test2 = (float *) malloc (sizeof(float) * DIMENSION);
50
51 *simplex = (float **) malloc (sizeof(float *) * (DIMENSION+1));
52
53 for (i = 0; i < DIMENSION+1; i++)
54 (*simplex)[i] = (float *) malloc (sizeof(float) * DIMENSION);
55
56 }
57
58
59 /*---------------------------------------------------------------------------*/
60
deallocate_arrays(float *** simplex,float ** centroid,float ** response,float ** step_size,float ** test1,float ** test2)61 void deallocate_arrays (float *** simplex, float ** centroid,
62 float ** response, float ** step_size,
63 float ** test1, float ** test2)
64 {
65 int i;
66
67 free (*centroid); *centroid = NULL;
68 free (*response); *response = NULL;
69 free (*step_size); *step_size = NULL;
70 free (*test1); *test1 = NULL;
71 free (*test2); *test2 = NULL;
72
73 for (i = 0; i < DIMENSION+1; i++)
74 {
75 free ((*simplex)[i]);
76 (*simplex)[i] = NULL;
77 }
78
79 free (*simplex); *simplex = NULL;
80
81 }
82
83
84 /*---------------------------------------------------------------------------*/
85
eval_vertices(float * response,int * worst,int * next,int * best)86 void eval_vertices (float * response, int * worst, int * next, int * best)
87 {
88 int i;
89
90 /* initialize values */
91 *worst = 0;
92 *best = 0;
93
94 /* find the best and worst */
95 for (i = 1; i < DIMENSION+1; i++)
96 {
97 if (response[i] > response[*worst])
98 *worst = i;
99 if (response[i] < response[*best])
100 *best = i;
101 }
102
103 /* find the next worst index */
104 if (*worst == 0)
105 *next = 1;
106 else
107 *next = 0;
108
109 for (i = 0; i < DIMENSION+1; i++)
110 if ((i != *worst) && (response[i] > response[*next]))
111 *next = i;
112 }
113
114
115 /*---------------------------------------------------------------------------*/
116
restart(float ** simplex,float * response,float * step_size)117 void restart (float ** simplex, float * response,
118 float * step_size)
119 {
120 const float STEP_FACTOR = 0.9;
121 int i, j;
122 int worst, next, best;
123 float minval, maxval;
124
125
126 /* find the current best vertex */
127 eval_vertices (response, &worst, &next, &best);
128
129 /* set the first vertex to the current best */
130 for (i = 0; i < DIMENSION; i++)
131 simplex[0][i] = simplex[best][i];
132
133 /* decrease step size */
134 for (i = 0; i < DIMENSION; i++)
135 step_size[i] *= STEP_FACTOR;
136
137 /* set up remaining vertices of simplex using new step size */
138 for (i = 1; i < DIMENSION+1; i++)
139 for (j = 0; j < DIMENSION; j++)
140 {
141 minval = simplex[0][j] - step_size[j];
142 maxval = simplex[0][j] + step_size[j];
143 simplex[i][j] = rand_uniform (minval, maxval);
144 }
145
146 /* initialize response for each vector */
147 for (i = 0; i < DIMENSION+1; i++)
148 response[i] = calc_error (simplex[i]);
149 }
150
151
152 /*---------------------------------------------------------------------------*/
153 /*
154 Calculate the centroid of the simplex, ignoring the worst vertex.
155 */
156
calc_centroid(float ** simplex,int worst,float * centroid)157 void calc_centroid (float ** simplex, int worst, float * centroid)
158 {
159 int i, j;
160
161 for (i = 0; i < DIMENSION; i++)
162 {
163 centroid[i] = 0.0;
164
165 /* add each vertex, except the worst */
166 for (j = 0; j < DIMENSION+1; j++)
167 if (j != worst)
168 centroid[i] += simplex[j][i];
169 }
170
171 /* divide by the number of vertices */
172 for (i = 0; i < DIMENSION; i++)
173 centroid[i] /= DIMENSION;
174 }
175
176
177 /*---------------------------------------------------------------------------*/
178 /*
179 Calculate the reflection of the worst vertex about the centroid.
180 */
181
calc_reflection(float ** simplex,float * centroid,int worst,float coef,float * vertex)182 void calc_reflection (float ** simplex, float * centroid,
183 int worst, float coef, float * vertex)
184 {
185 int i;
186
187 for (i = 0; i < DIMENSION; i++)
188 vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
189 }
190
191
192 /*---------------------------------------------------------------------------*/
193 /*
194 Replace a vertex of the simplex.
195 */
196
replace(float ** simplex,float * response,int index,float * vertex,float resp)197 void replace (float ** simplex, float * response,
198 int index, float * vertex, float resp)
199 {
200 int i;
201
202 for (i = 0; i < DIMENSION; i++)
203 simplex[index][i] = vertex[i];
204
205 response[index] = resp;
206 }
207
208
209 /*---------------------------------------------------------------------------*/
210 /*
211 Calculate goodness of fit.
212 */
213
calc_good_fit(float * response)214 float calc_good_fit (float * response)
215 {
216 int i;
217
218 float avg, sd, tmp;
219
220 /* average the response values */
221 avg = 0.0;
222 for (i = 0; i < DIMENSION+1; i++)
223 avg += response[i];
224 avg /= DIMENSION+1;
225
226 /* compute standard deviation of response */
227 sd = 0.0;
228 for (i = 0; i < DIMENSION+1; i++)
229 {
230 tmp = response[i] - avg;
231 sd += tmp*tmp;
232 }
233 sd /= DIMENSION;
234
235 return (sqrt(sd));
236 }
237
238
239 /*---------------------------------------------------------------------------*/
240 /*
241 Perform initialization for the Simplex algorithm.
242 */
243
simplex_initialize(float * parameters,float ** simplex,float * response,float * step_size)244 void simplex_initialize (float * parameters, float ** simplex,
245 float * response, float * step_size)
246 {
247 int i, j;
248 int worst, next, best;
249 float resp;
250 float minval, maxval;
251
252
253 for (j = 0; j < DIMENSION; j++)
254 {
255 simplex[0][j] = parameters[j];
256 step_size[j] = 0.5 * parameters[j];
257 }
258
259 for (i = 1; i < DIMENSION+1; i++)
260 for (j = 0; j < DIMENSION; j++)
261 {
262 minval = simplex[0][j] - step_size[j];
263 maxval = simplex[0][j] + step_size[j];
264 simplex[i][j] = rand_uniform (minval, maxval);
265 }
266
267 for (i = 0; i < DIMENSION+1; i++)
268 response[i] = calc_error (simplex[i]);
269
270 for (i = 1; i < 500; i++)
271 {
272 for (j = 0; j < DIMENSION; j++)
273 {
274 minval = simplex[0][j] - step_size[j];
275 maxval = simplex[0][j] + step_size[j];
276 parameters[j] = rand_uniform (minval, maxval);
277 }
278
279 resp = calc_error (parameters);
280 eval_vertices (response, &worst, &next, &best);
281 if (resp < response[worst])
282 replace (simplex, response, worst, parameters, resp);
283 }
284 }
285
286
287 /*---------------------------------------------------------------------------*/
288 /*
289 Use Simplex algorithm to estimate the voxel intensity distribution.
290
291 The Simplex algorithm is adapted from: A Jump Start Course in C++ Programming
292 */
293
simplex_optimization(float * parameters,float * sse)294 void simplex_optimization (float * parameters, float * sse)
295 {
296 const int MAX_ITERATIONS = 100;
297 const int MAX_RESTARTS = 25;
298 const float EXPANSION_COEF = 2.0;
299 const float REFLECTION_COEF = 1.0;
300 const float CONTRACTION_COEF = 0.5;
301 const float TOLERANCE = 1.0e-10;
302
303 float ** simplex = NULL;
304 float * centroid = NULL;
305 float * response = NULL;
306 float * step_size = NULL;
307 float * test1 = NULL;
308 float * test2 = NULL;
309 float resp1, resp2;
310 int i, worst, best, next;
311 int num_iter, num_restarts;
312 int done;
313 float fit;
314
315
316 allocate_arrays (&simplex, ¢roid, &response, &step_size, &test1, &test2);
317
318 simplex_initialize (parameters, simplex, response, step_size);
319
320 /* start loop to do simplex optimization */
321 num_iter = 0;
322 num_restarts = 0;
323 done = 0;
324
325 while (!done)
326 {
327 /* find the worst vertex and compute centroid of remaining simplex,
328 discarding the worst vertex */
329 eval_vertices (response, &worst, &next, &best);
330 calc_centroid (simplex, worst, centroid);
331
332 /* reflect the worst point through the centroid */
333 calc_reflection (simplex, centroid, worst,
334 REFLECTION_COEF, test1);
335 resp1 = calc_error (test1);
336
337 /* test the reflection against the best vertex and expand it if the
338 reflection is better. if not, keep the reflection */
339 if (resp1 < response[best])
340 {
341 /* try expanding */
342 calc_reflection (simplex, centroid, worst, EXPANSION_COEF, test2);
343 resp2 = calc_error (test2);
344 if (resp2 <= resp1) /* keep expansion */
345 replace (simplex, response, worst, test2, resp2);
346 else /* keep reflection */
347 replace (simplex, response, worst, test1, resp1);
348 }
349 else if (resp1 < response[next])
350 {
351 /* new response is between the best and next worst
352 so keep reflection */
353 replace (simplex, response, worst, test1, resp1);
354 }
355 else
356 {
357 /* try contraction */
358 if (resp1 >= response[worst])
359 calc_reflection (simplex, centroid, worst,
360 -CONTRACTION_COEF, test2);
361 else
362 calc_reflection (simplex, centroid, worst,
363 CONTRACTION_COEF, test2);
364 resp2 = calc_error (test2);
365
366 /* test the contracted response against the worst response */
367 if (resp2 > response[worst])
368 {
369 /* new contracted response is worse, so decrease step size
370 and restart */
371 num_iter = 0;
372 num_restarts += 1;
373 restart (simplex, response, step_size);
374 }
375 else /* keep contraction */
376 replace (simplex, response, worst, test2, resp2);
377 }
378
379 /* test to determine when to stop.
380 first, check the number of iterations */
381 num_iter += 1; /* increment iteration counter */
382 if (num_iter >= MAX_ITERATIONS)
383 {
384 /* restart with smaller steps */
385 num_iter = 0;
386 num_restarts += 1;
387 restart (simplex, response, step_size);
388 }
389
390 /* limit the number of restarts */
391 if (num_restarts == MAX_RESTARTS) done = 1;
392
393 /* compare relative standard deviation of vertex responses
394 against a defined tolerance limit */
395 fit = calc_good_fit (response);
396 if (fit <= TOLERANCE) done = 1;
397
398 /* if done, copy the best solution to the output array */
399 if (done)
400 {
401 eval_vertices (response, &worst, &next, &best);
402 for (i = 0; i < DIMENSION; i++)
403 parameters[i] = simplex[best][i];
404 *sse = response[best];
405 }
406
407 } /* while (!done) */
408
409 number_restarts = num_restarts;
410 deallocate_arrays (&simplex, ¢roid, &response, &step_size,
411 &test1, &test2);
412
413 }
414
415
416 /*---------------------------------------------------------------------------*/
417