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