Actual source code: neldermead.c
1: #include "neldermead.h"
2: #include "petscvec.h"
4: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm);
5: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f);
6: /* ---------------------------------------------------------- */
9: static PetscErrorCode TaoSetUp_NM(TaoSolver tao)
10: {
12: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
13: PetscInt size;
16: VecGetSize(tao->solution,&size);
17: nm->N = size;
18: nm->oneOverN = 1.0/size;
19: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
20: PetscMalloc((nm->N+1)*sizeof(PetscReal),&nm->f_values);
21: PetscMalloc((nm->N+1)*sizeof(PetscInt),&nm->indices);
22: VecDuplicate(tao->solution,&nm->Xbar);
23: VecDuplicate(tao->solution,&nm->Xmur);
24: VecDuplicate(tao->solution,&nm->Xmue);
25: VecDuplicate(tao->solution,&nm->Xmuc);
27: tao->gradient=0;
28: tao->step=0;
30: return(0);
31: }
33: /* ---------------------------------------------------------- */
36: PetscErrorCode TaoDestroy_NM(TaoSolver tao)
37: {
38: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
41: if (tao->setupcalled) {
42: VecDestroyVecs(nm->N+1,&nm->simplex);
43: VecDestroy(&nm->Xmuc);
44: VecDestroy(&nm->Xmue);
45: VecDestroy(&nm->Xmur);
46: VecDestroy(&nm->Xbar);
47: }
48: PetscFree(nm->indices);
49: PetscFree(nm->f_values);
50: PetscFree(tao->data);
51: tao->data = 0;
52: return(0);
53: }
55: /*------------------------------------------------------------*/
58: PetscErrorCode TaoSetFromOptions_NM(TaoSolver tao)
59: {
60:
61: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
62: PetscBool flg;
64:
66: PetscOptionsHead("Nelder-Mead options");
68: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,&flg);
70: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,&flg);
71: nm->mu_ic = -nm->mu_oc;
72: nm->mu_r = nm->mu_oc*2.0;
73: nm->mu_e = nm->mu_oc*4.0;
75: PetscOptionsTail();
76: return(0);
77: }
79: /*------------------------------------------------------------*/
82: PetscErrorCode TaoView_NM(TaoSolver tao,PetscViewer viewer)
83: {
84: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
85: PetscBool isascii;
87:
89: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
90: if (isascii) {
91: PetscViewerASCIIPushTab(viewer);
92: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
93: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
94: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
95: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
96: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
97: PetscViewerASCIIPopTab(viewer);
98: } else {
99: SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NelderMead",((PetscObject)viewer)->type_name);
100: }
102: return(0);
103: }
105: /*------------------------------------------------------------*/
108: PetscErrorCode TaoSolve_NM(TaoSolver tao)
109: {
111: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
112: TaoSolverTerminationReason reason;
113: PetscReal *x;
114: PetscInt iter=0,i;
115: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
116: PetscReal fr,fe,fc;
117: PetscInt shrink;
118: PetscInt low,high;
119:
120:
122: nm->nshrink = 0;
123: nm->nreflect = 0;
124: nm->nincontract = 0;
125: nm->noutcontract = 0;
126: nm->nexpand = 0;
127:
128: if (tao->XL || tao->XU || tao->ops->computebounds) {
129: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
130: }
132: VecCopy(tao->solution,nm->simplex[0]);
133: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
134: nm->indices[0]=0;
135: for (i=1;i<nm->N+1;i++){
136: VecCopy(tao->solution,nm->simplex[i]);
137: VecGetOwnershipRange(nm->simplex[i],&low,&high);
138: if (i-1 >= low && i-1 < high) {
139: VecGetArray(nm->simplex[i],&x);
140: x[i-1-low] += nm->lamda;
141: VecRestoreArray(nm->simplex[i],&x);
142: }
144: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
145: nm->indices[i] = i;
146: }
147:
148: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
149: NelderMeadSort(nm);
150: VecSet(Xbar,0.0);
151: for (i=0;i<nm->N;i++) {
152: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
153: }
154: VecScale(Xbar,nm->oneOverN);
155: reason = TAO_CONTINUE_ITERATING;
156: while (1) {
157: shrink = 0;
158: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
159: TaoMonitor(tao,iter++,nm->f_values[nm->indices[0]],nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]],0.0,1.0,&reason);
160: if (reason != TAO_CONTINUE_ITERATING) break;
162:
163:
164: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
165: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
166: TaoComputeObjective(tao,Xmur,&fr);
169: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
170: /* reflect */
171: nm->nreflect++;
172: PetscInfo(0,"Reflect\n");
173: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
174: }
176: else if (fr < nm->f_values[nm->indices[0]]) {
177: /* expand */
178: nm->nexpand++;
179: PetscInfo(0,"Expand\n");
180: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
181: TaoComputeObjective(tao,Xmue,&fe);
182: if (fe < fr) {
183: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
184: } else {
185: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
186: }
188: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
189: /* outside contraction */
190: nm->noutcontract++;
191: PetscInfo(0,"Outside Contraction\n");
192: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
194: TaoComputeObjective(tao,Xmuc,&fc);
195: if (fc <= fr) {
196: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
197: } else
198: shrink=1;
199: } else {
200: /* inside contraction */
201: nm->nincontract++;
202: PetscInfo(0,"Inside Contraction\n");
203: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
204: TaoComputeObjective(tao,Xmuc,&fc);
205: if (fc < nm->f_values[nm->indices[nm->N]]) {
206: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
207: } else
208: shrink = 1;
209: }
211: if (shrink) {
212: nm->nshrink++;
213: PetscInfo(0,"Shrink\n");
214:
215: for (i=1;i<nm->N+1;i++) {
216: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
217: TaoComputeObjective(tao,nm->simplex[nm->indices[i]],
218: &nm->f_values[nm->indices[i]]);
219: }
220: VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);
222: /* Add last vector's fraction of average */
223: VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
224: NelderMeadSort(nm);
225: /* Subtract new last vector from average */
226: VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
227: }
228:
229:
230: }
231: return(0);
232: }
234: /* ---------------------------------------------------------- */
238: PetscErrorCode TaoCreate_NM(TaoSolver tao)
239: {
240: TAO_NelderMead *nm;
244: PetscNewLog(tao,TAO_NelderMead,&nm);
245: tao->data = (void*)nm;
247: tao->ops->setup = TaoSetUp_NM;
248: tao->ops->solve = TaoSolve_NM;
249: tao->ops->view = TaoView_NM;
250: tao->ops->setfromoptions = TaoSetFromOptions_NM;
251: tao->ops->destroy = TaoDestroy_NM;
253: tao->max_it = 2000;
254: tao->max_funcs = 4000;
255: tao->fatol = 1e-8;
256: tao->frtol = 1e-8;
258: nm->simplex = 0;
259: nm->lamda = 1;
261: nm->mu_ic = -0.5;
262: nm->mu_oc = 0.5;
263: nm->mu_r = 1.0;
264: nm->mu_e = 2.0;
266: return(0);
267: }
269:
271: /*------------------------------------------------------------*/
274: PetscErrorCode NelderMeadSort(TAO_NelderMead *nm) {
275: PetscReal *values = nm->f_values;
276: PetscInt *indices = nm->indices;
277: PetscInt dim = nm->N+1;
279: PetscInt i,j,index;
280: PetscReal val;
282: for (i=1;i<dim;i++) {
283: index = indices[i];
284: val = values[index];
285: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
286: indices[j+1] = indices[j];
287: }
288: indices[j+1] = index;
289: }
290: return(0);
291: }
294: /*------------------------------------------------------------*/
297: PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
298: {
301: /* Add new vector's fraction of average */
302: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
303: VecCopy(Xmu,nm->simplex[index]);
304: nm->f_values[index] = f;
306: NelderMeadSort(nm);
308: /* Subtract last vector from average */
309: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
310: return(0);
311: }