Actual source code: neldermead.c
petsc-3.5.4 2015-05-23
1: #include <../src/tao/unconstrained/impls/neldermead/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(Tao tao)
10: {
12: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
13: PetscInt n;
16: VecGetSize(tao->solution,&n);
17: nm->N = n;
18: nm->oneOverN = 1.0/n;
19: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
20: PetscMalloc1((nm->N+1),&nm->f_values);
21: PetscMalloc1((nm->N+1),&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;
29: return(0);
30: }
32: /* ---------------------------------------------------------- */
35: PetscErrorCode TaoDestroy_NM(Tao tao)
36: {
37: 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(Tao tao)
59: {
60: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
61: PetscBool flg;
65: PetscOptionsHead("Nelder-Mead options");
66: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,&flg);
67: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,&flg);
68: nm->mu_ic = -nm->mu_oc;
69: nm->mu_r = nm->mu_oc*2.0;
70: nm->mu_e = nm->mu_oc*4.0;
71: PetscOptionsTail();
72: return(0);
73: }
75: /*------------------------------------------------------------*/
78: PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
79: {
80: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
81: PetscBool isascii;
85: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
86: if (isascii) {
87: PetscViewerASCIIPushTab(viewer);
88: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
89: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
90: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
91: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
92: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
93: PetscViewerASCIIPopTab(viewer);
94: }
95: return(0);
96: }
98: /*------------------------------------------------------------*/
101: PetscErrorCode TaoSolve_NM(Tao tao)
102: {
103: PetscErrorCode ierr;
104: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
105: TaoConvergedReason reason;
106: PetscReal *x;
107: PetscInt iter=0,i;
108: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
109: PetscReal fr,fe,fc;
110: PetscInt shrink;
111: PetscInt low,high;
114: nm->nshrink = 0;
115: nm->nreflect = 0;
116: nm->nincontract = 0;
117: nm->noutcontract = 0;
118: nm->nexpand = 0;
120: if (tao->XL || tao->XU || tao->ops->computebounds) {
121: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
122: }
124: VecCopy(tao->solution,nm->simplex[0]);
125: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
126: nm->indices[0]=0;
127: for (i=1;i<nm->N+1;i++){
128: VecCopy(tao->solution,nm->simplex[i]);
129: VecGetOwnershipRange(nm->simplex[i],&low,&high);
130: if (i-1 >= low && i-1 < high) {
131: VecGetArray(nm->simplex[i],&x);
132: x[i-1-low] += nm->lamda;
133: VecRestoreArray(nm->simplex[i],&x);
134: }
136: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
137: nm->indices[i] = i;
138: }
140: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
141: NelderMeadSort(nm);
142: VecSet(Xbar,0.0);
143: for (i=0;i<nm->N;i++) {
144: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
145: }
146: VecScale(Xbar,nm->oneOverN);
147: reason = TAO_CONTINUE_ITERATING;
148: while (1) {
149: shrink = 0;
150: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
151: 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);
152: if (reason != TAO_CONTINUE_ITERATING) break;
154: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
155: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
156: TaoComputeObjective(tao,Xmur,&fr);
158: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
159: /* reflect */
160: nm->nreflect++;
161: PetscInfo(0,"Reflect\n");
162: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
163: } else if (fr < nm->f_values[nm->indices[0]]) {
164: /* expand */
165: nm->nexpand++;
166: PetscInfo(0,"Expand\n");
167: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
168: TaoComputeObjective(tao,Xmue,&fe);
169: if (fe < fr) {
170: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
171: } else {
172: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
173: }
174: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
175: /* outside contraction */
176: nm->noutcontract++;
177: PetscInfo(0,"Outside Contraction\n");
178: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
180: TaoComputeObjective(tao,Xmuc,&fc);
181: if (fc <= fr) {
182: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
183: } else shrink=1;
184: } else {
185: /* inside contraction */
186: nm->nincontract++;
187: PetscInfo(0,"Inside Contraction\n");
188: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
189: TaoComputeObjective(tao,Xmuc,&fc);
190: if (fc < nm->f_values[nm->indices[nm->N]]) {
191: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
192: } else shrink = 1;
193: }
195: if (shrink) {
196: nm->nshrink++;
197: PetscInfo(0,"Shrink\n");
199: for (i=1;i<nm->N+1;i++) {
200: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
201: TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
202: }
203: VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);
205: /* Add last vector's fraction of average */
206: VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
207: NelderMeadSort(nm);
208: /* Subtract new last vector from average */
209: VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
210: }
211: }
212: return(0);
213: }
215: /* ---------------------------------------------------------- */
216: /*MC
217: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
219: Options Database Keys:
220: + -tao_nm_lamda - initial step length
221: . -tao_nm_mu - expansion/contraction factor
223: Level: beginner
224: M*/
226: EXTERN_C_BEGIN
229: PetscErrorCode TaoCreate_NM(Tao tao)
230: {
231: TAO_NelderMead *nm;
235: PetscNewLog(tao,&nm);
236: tao->data = (void*)nm;
238: tao->ops->setup = TaoSetUp_NM;
239: tao->ops->solve = TaoSolve_NM;
240: tao->ops->view = TaoView_NM;
241: tao->ops->setfromoptions = TaoSetFromOptions_NM;
242: tao->ops->destroy = TaoDestroy_NM;
244: tao->max_it = 2000;
245: tao->max_funcs = 4000;
246: #if defined(PETSC_USE_REAL_SINGLE)
247: tao->fatol = 1e-5;
248: tao->frtol = 1e-5;
249: #else
250: tao->fatol = 1e-8;
251: tao->frtol = 1e-8;
252: #endif
254: nm->simplex = 0;
255: nm->lamda = 1;
257: nm->mu_ic = -0.5;
258: nm->mu_oc = 0.5;
259: nm->mu_r = 1.0;
260: nm->mu_e = 2.0;
262: return(0);
263: }
264: EXTERN_C_END
267: /*------------------------------------------------------------*/
270: PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
271: {
272: PetscReal *values = nm->f_values;
273: PetscInt *indices = nm->indices;
274: PetscInt dim = nm->N+1;
275: PetscInt i,j,index;
276: PetscReal val;
279: for (i=1;i<dim;i++) {
280: index = indices[i];
281: val = values[index];
282: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
283: indices[j+1] = indices[j];
284: }
285: indices[j+1] = index;
286: }
287: return(0);
288: }
291: /*------------------------------------------------------------*/
294: PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
295: {
299: /* Add new vector's fraction of average */
300: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
301: VecCopy(Xmu,nm->simplex[index]);
302: nm->f_values[index] = f;
304: NelderMeadSort(nm);
306: /* Subtract last vector from average */
307: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
308: return(0);
309: }