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;
13: for (i=1;i<dim;i++) {
14: index = indices[i];
15: val = values[index];
16: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
17: indices[j+1] = indices[j];
18: }
19: indices[j+1] = index;
20: }
21: return 0;
22: }
24: /*------------------------------------------------------------*/
25: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
26: {
27: /* Add new vector's fraction of average */
28: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
29: VecCopy(Xmu,nm->simplex[index]);
30: nm->f_values[index] = f;
32: NelderMeadSort(nm);
34: /* Subtract last vector from average */
35: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
36: return 0;
37: }
39: /* ---------------------------------------------------------- */
40: static PetscErrorCode TaoSetUp_NM(Tao tao)
41: {
42: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
43: PetscInt n;
45: VecGetSize(tao->solution,&n);
46: nm->N = n;
47: nm->oneOverN = 1.0/n;
48: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
49: PetscMalloc1(nm->N+1,&nm->f_values);
50: PetscMalloc1(nm->N+1,&nm->indices);
51: VecDuplicate(tao->solution,&nm->Xbar);
52: VecDuplicate(tao->solution,&nm->Xmur);
53: VecDuplicate(tao->solution,&nm->Xmue);
54: VecDuplicate(tao->solution,&nm->Xmuc);
56: tao->gradient=NULL;
57: tao->step=0;
58: return 0;
59: }
61: /* ---------------------------------------------------------- */
62: static PetscErrorCode TaoDestroy_NM(Tao tao)
63: {
64: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
66: if (tao->setupcalled) {
67: VecDestroyVecs(nm->N+1,&nm->simplex);
68: VecDestroy(&nm->Xmuc);
69: VecDestroy(&nm->Xmue);
70: VecDestroy(&nm->Xmur);
71: VecDestroy(&nm->Xbar);
72: }
73: PetscFree(nm->indices);
74: PetscFree(nm->f_values);
75: PetscFree(tao->data);
76: return 0;
77: }
79: /*------------------------------------------------------------*/
80: static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao)
81: {
82: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
84: PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");
85: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL);
86: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);
87: nm->mu_ic = -nm->mu_oc;
88: nm->mu_r = nm->mu_oc*2.0;
89: nm->mu_e = nm->mu_oc*4.0;
90: PetscOptionsTail();
91: return 0;
92: }
94: /*------------------------------------------------------------*/
95: static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
96: {
97: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
98: PetscBool isascii;
100: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
101: if (isascii) {
102: PetscViewerASCIIPushTab(viewer);
103: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
104: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
105: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
106: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
107: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
108: PetscViewerASCIIPopTab(viewer);
109: }
110: return 0;
111: }
113: /*------------------------------------------------------------*/
114: static PetscErrorCode TaoSolve_NM(Tao tao)
115: {
116: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
117: PetscReal *x;
118: PetscInt i;
119: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
120: PetscReal fr,fe,fc;
121: PetscInt shrink;
122: PetscInt low,high;
124: nm->nshrink = 0;
125: nm->nreflect = 0;
126: nm->nincontract = 0;
127: nm->noutcontract = 0;
128: nm->nexpand = 0;
130: if (tao->XL || tao->XU || tao->ops->computebounds) {
131: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
132: }
134: VecCopy(tao->solution,nm->simplex[0]);
135: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
136: nm->indices[0]=0;
137: for (i=1;i<nm->N+1;i++) {
138: VecCopy(tao->solution,nm->simplex[i]);
139: VecGetOwnershipRange(nm->simplex[i],&low,&high);
140: if (i-1 >= low && i-1 < high) {
141: VecGetArray(nm->simplex[i],&x);
142: x[i-1-low] += nm->lamda;
143: VecRestoreArray(nm->simplex[i],&x);
144: }
146: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
147: nm->indices[i] = i;
148: }
150: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
151: NelderMeadSort(nm);
152: VecSet(Xbar,0.0);
153: for (i=0;i<nm->N;i++) {
154: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
155: }
156: VecScale(Xbar,nm->oneOverN);
157: tao->reason = TAO_CONTINUE_ITERATING;
158: while (1) {
159: /* Call general purpose update function */
160: if (tao->ops->update) {
161: (*tao->ops->update)(tao, tao->niter, tao->user_update);
162: }
163: ++tao->niter;
164: shrink = 0;
165: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
166: 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);
167: 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);
168: (*tao->ops->convergencetest)(tao,tao->cnvP);
169: if (tao->reason != TAO_CONTINUE_ITERATING) break;
171: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
172: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
173: TaoComputeObjective(tao,Xmur,&fr);
175: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
176: /* reflect */
177: nm->nreflect++;
178: PetscInfo(0,"Reflect\n");
179: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
180: } else if (fr < nm->f_values[nm->indices[0]]) {
181: /* expand */
182: nm->nexpand++;
183: PetscInfo(0,"Expand\n");
184: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
185: TaoComputeObjective(tao,Xmue,&fe);
186: if (fe < fr) {
187: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
188: } else {
189: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
190: }
191: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
192: /* outside contraction */
193: nm->noutcontract++;
194: PetscInfo(0,"Outside Contraction\n");
195: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
197: TaoComputeObjective(tao,Xmuc,&fc);
198: if (fc <= fr) {
199: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
200: } else shrink=1;
201: } else {
202: /* inside contraction */
203: nm->nincontract++;
204: PetscInfo(0,"Inside Contraction\n");
205: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
206: TaoComputeObjective(tao,Xmuc,&fc);
207: if (fc < nm->f_values[nm->indices[nm->N]]) {
208: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
209: } else shrink = 1;
210: }
212: if (shrink) {
213: nm->nshrink++;
214: PetscInfo(0,"Shrink\n");
216: for (i=1;i<nm->N+1;i++) {
217: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
218: TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &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: return 0;
230: }
232: /* ---------------------------------------------------------- */
233: /*MC
234: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
236: Options Database Keys:
237: + -tao_nm_lamda - initial step length
238: - -tao_nm_mu - expansion/contraction factor
240: Level: beginner
241: M*/
243: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
244: {
245: TAO_NelderMead *nm;
247: PetscNewLog(tao,&nm);
248: tao->data = (void*)nm;
250: tao->ops->setup = TaoSetUp_NM;
251: tao->ops->solve = TaoSolve_NM;
252: tao->ops->view = TaoView_NM;
253: tao->ops->setfromoptions = TaoSetFromOptions_NM;
254: tao->ops->destroy = TaoDestroy_NM;
256: /* Override default settings (unless already changed) */
257: if (!tao->max_it_changed) tao->max_it = 2000;
258: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
260: nm->simplex = NULL;
261: nm->lamda = 1;
263: nm->mu_ic = -0.5;
264: nm->mu_oc = 0.5;
265: nm->mu_r = 1.0;
266: nm->mu_e = 2.0;
268: return 0;
269: }