Actual source code: neldermead.c
1: #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
2: #include <petscvec.h>
4: /*------------------------------------------------------------*/
5: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
6: {
7: PetscReal *values = nm->f_values;
8: PetscInt *indices = nm->indices;
9: PetscInt dim = nm->N+1;
10: PetscInt i,j,index;
11: PetscReal val;
14: for (i=1;i<dim;i++) {
15: index = indices[i];
16: val = values[index];
17: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
18: indices[j+1] = indices[j];
19: }
20: indices[j+1] = index;
21: }
22: return(0);
23: }
25: /*------------------------------------------------------------*/
26: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
27: {
31: /* Add new vector's fraction of average */
32: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
33: VecCopy(Xmu,nm->simplex[index]);
34: nm->f_values[index] = f;
36: NelderMeadSort(nm);
38: /* Subtract last vector from average */
39: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
40: return(0);
41: }
43: /* ---------------------------------------------------------- */
44: static PetscErrorCode TaoSetUp_NM(Tao tao)
45: {
47: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
48: PetscInt n;
51: VecGetSize(tao->solution,&n);
52: nm->N = n;
53: nm->oneOverN = 1.0/n;
54: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
55: PetscMalloc1(nm->N+1,&nm->f_values);
56: PetscMalloc1(nm->N+1,&nm->indices);
57: VecDuplicate(tao->solution,&nm->Xbar);
58: VecDuplicate(tao->solution,&nm->Xmur);
59: VecDuplicate(tao->solution,&nm->Xmue);
60: VecDuplicate(tao->solution,&nm->Xmuc);
62: tao->gradient=NULL;
63: tao->step=0;
64: return(0);
65: }
67: /* ---------------------------------------------------------- */
68: static PetscErrorCode TaoDestroy_NM(Tao tao)
69: {
70: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
74: if (tao->setupcalled) {
75: VecDestroyVecs(nm->N+1,&nm->simplex);
76: VecDestroy(&nm->Xmuc);
77: VecDestroy(&nm->Xmue);
78: VecDestroy(&nm->Xmur);
79: VecDestroy(&nm->Xbar);
80: }
81: PetscFree(nm->indices);
82: PetscFree(nm->f_values);
83: PetscFree(tao->data);
84: return(0);
85: }
87: /*------------------------------------------------------------*/
88: static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao)
89: {
90: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
94: PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");
95: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL);
96: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);
97: nm->mu_ic = -nm->mu_oc;
98: nm->mu_r = nm->mu_oc*2.0;
99: nm->mu_e = nm->mu_oc*4.0;
100: PetscOptionsTail();
101: return(0);
102: }
104: /*------------------------------------------------------------*/
105: static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
106: {
107: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
108: PetscBool isascii;
112: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
113: if (isascii) {
114: PetscViewerASCIIPushTab(viewer);
115: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
116: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
117: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
118: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
119: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
120: PetscViewerASCIIPopTab(viewer);
121: }
122: return(0);
123: }
125: /*------------------------------------------------------------*/
126: static PetscErrorCode TaoSolve_NM(Tao tao)
127: {
128: PetscErrorCode ierr;
129: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
130: PetscReal *x;
131: PetscInt i;
132: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
133: PetscReal fr,fe,fc;
134: PetscInt shrink;
135: PetscInt low,high;
138: nm->nshrink = 0;
139: nm->nreflect = 0;
140: nm->nincontract = 0;
141: nm->noutcontract = 0;
142: nm->nexpand = 0;
144: if (tao->XL || tao->XU || tao->ops->computebounds) {
145: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
146: }
148: VecCopy(tao->solution,nm->simplex[0]);
149: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
150: nm->indices[0]=0;
151: for (i=1;i<nm->N+1;i++) {
152: VecCopy(tao->solution,nm->simplex[i]);
153: VecGetOwnershipRange(nm->simplex[i],&low,&high);
154: if (i-1 >= low && i-1 < high) {
155: VecGetArray(nm->simplex[i],&x);
156: x[i-1-low] += nm->lamda;
157: VecRestoreArray(nm->simplex[i],&x);
158: }
160: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
161: nm->indices[i] = i;
162: }
164: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
165: NelderMeadSort(nm);
166: VecSet(Xbar,0.0);
167: for (i=0;i<nm->N;i++) {
168: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
169: }
170: VecScale(Xbar,nm->oneOverN);
171: tao->reason = TAO_CONTINUE_ITERATING;
172: while (1) {
173: /* Call general purpose update function */
174: if (tao->ops->update) {
175: (*tao->ops->update)(tao, tao->niter, tao->user_update);
176: }
177: ++tao->niter;
178: shrink = 0;
179: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
180: TaoLogConvergenceHistory(tao, nm->f_values[nm->indices[0]], nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]], 0.0, tao->ksp_its);
181: TaoMonitor(tao,tao->niter, nm->f_values[nm->indices[0]], nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]], 0.0, 1.0);
182: (*tao->ops->convergencetest)(tao,tao->cnvP);
183: if (tao->reason != TAO_CONTINUE_ITERATING) break;
185: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
186: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
187: TaoComputeObjective(tao,Xmur,&fr);
189: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
190: /* reflect */
191: nm->nreflect++;
192: PetscInfo(0,"Reflect\n");
193: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
194: } else if (fr < nm->f_values[nm->indices[0]]) {
195: /* expand */
196: nm->nexpand++;
197: PetscInfo(0,"Expand\n");
198: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
199: TaoComputeObjective(tao,Xmue,&fe);
200: if (fe < fr) {
201: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
202: } else {
203: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
204: }
205: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
206: /* outside contraction */
207: nm->noutcontract++;
208: PetscInfo(0,"Outside Contraction\n");
209: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
211: TaoComputeObjective(tao,Xmuc,&fc);
212: if (fc <= fr) {
213: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
214: } else shrink=1;
215: } else {
216: /* inside contraction */
217: nm->nincontract++;
218: PetscInfo(0,"Inside Contraction\n");
219: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
220: TaoComputeObjective(tao,Xmuc,&fc);
221: if (fc < nm->f_values[nm->indices[nm->N]]) {
222: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
223: } else shrink = 1;
224: }
226: if (shrink) {
227: nm->nshrink++;
228: PetscInfo(0,"Shrink\n");
230: for (i=1;i<nm->N+1;i++) {
231: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
232: TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
233: }
234: VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);
236: /* Add last vector's fraction of average */
237: VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
238: NelderMeadSort(nm);
239: /* Subtract new last vector from average */
240: VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
241: }
242: }
243: return(0);
244: }
246: /* ---------------------------------------------------------- */
247: /*MC
248: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
250: Options Database Keys:
251: + -tao_nm_lamda - initial step length
252: - -tao_nm_mu - expansion/contraction factor
254: Level: beginner
255: M*/
257: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
258: {
259: TAO_NelderMead *nm;
263: PetscNewLog(tao,&nm);
264: tao->data = (void*)nm;
266: tao->ops->setup = TaoSetUp_NM;
267: tao->ops->solve = TaoSolve_NM;
268: tao->ops->view = TaoView_NM;
269: tao->ops->setfromoptions = TaoSetFromOptions_NM;
270: tao->ops->destroy = TaoDestroy_NM;
272: /* Override default settings (unless already changed) */
273: if (!tao->max_it_changed) tao->max_it = 2000;
274: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
276: nm->simplex = NULL;
277: nm->lamda = 1;
279: nm->mu_ic = -0.5;
280: nm->mu_oc = 0.5;
281: nm->mu_r = 1.0;
282: nm->mu_e = 2.0;
284: return(0);
285: }