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: PetscOptionsDeprecated("-tao_nm_lamda", "-tao_nm_lambda", "3.18.4", NULL);
84: PetscOptionsReal("-tao_nm_lambda", "initial step length", "", nm->lambda, &nm->lambda, NULL);
85: PetscOptionsReal("-tao_nm_mu", "mu", "", nm->mu_oc, &nm->mu_oc, NULL);
86: nm->mu_ic = -nm->mu_oc;
87: nm->mu_r = nm->mu_oc * 2.0;
88: nm->mu_e = nm->mu_oc * 4.0;
89: PetscOptionsHeadEnd();
90: return 0;
91: }
93: /*------------------------------------------------------------*/
94: static PetscErrorCode TaoView_NM(Tao tao, PetscViewer viewer)
95: {
96: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
97: PetscBool isascii;
99: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
100: if (isascii) {
101: PetscViewerASCIIPushTab(viewer);
102: PetscViewerASCIIPrintf(viewer, "expansions: %" PetscInt_FMT "\n", nm->nexpand);
103: PetscViewerASCIIPrintf(viewer, "reflections: %" PetscInt_FMT "\n", nm->nreflect);
104: PetscViewerASCIIPrintf(viewer, "inside contractions: %" PetscInt_FMT "\n", nm->nincontract);
105: PetscViewerASCIIPrintf(viewer, "outside contractionss: %" PetscInt_FMT "\n", nm->noutcontract);
106: PetscViewerASCIIPrintf(viewer, "Shrink steps: %" PetscInt_FMT "\n", nm->nshrink);
107: PetscViewerASCIIPopTab(viewer);
108: }
109: return 0;
110: }
112: /*------------------------------------------------------------*/
113: static PetscErrorCode TaoSolve_NM(Tao tao)
114: {
115: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
116: PetscReal *x;
117: PetscInt i;
118: Vec Xmur = nm->Xmur, Xmue = nm->Xmue, Xmuc = nm->Xmuc, Xbar = nm->Xbar;
119: PetscReal fr, fe, fc;
120: PetscInt shrink;
121: PetscInt low, high;
123: nm->nshrink = 0;
124: nm->nreflect = 0;
125: nm->nincontract = 0;
126: nm->noutcontract = 0;
127: nm->nexpand = 0;
129: if (tao->XL || tao->XU || tao->ops->computebounds) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
131: VecCopy(tao->solution, nm->simplex[0]);
132: TaoComputeObjective(tao, nm->simplex[0], &nm->f_values[0]);
133: nm->indices[0] = 0;
134: for (i = 1; i < nm->N + 1; i++) {
135: VecCopy(tao->solution, nm->simplex[i]);
136: VecGetOwnershipRange(nm->simplex[i], &low, &high);
137: if (i - 1 >= low && i - 1 < high) {
138: VecGetArray(nm->simplex[i], &x);
139: x[i - 1 - low] += nm->lambda;
140: VecRestoreArray(nm->simplex[i], &x);
141: }
143: TaoComputeObjective(tao, nm->simplex[i], &nm->f_values[i]);
144: nm->indices[i] = i;
145: }
147: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
148: NelderMeadSort(nm);
149: VecSet(Xbar, 0.0);
150: for (i = 0; i < nm->N; i++) VecAXPY(Xbar, 1.0, nm->simplex[nm->indices[i]]);
151: VecScale(Xbar, nm->oneOverN);
152: tao->reason = TAO_CONTINUE_ITERATING;
153: while (1) {
154: /* Call general purpose update function */
155: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
156: ++tao->niter;
157: shrink = 0;
158: VecCopy(nm->simplex[nm->indices[0]], tao->solution);
159: 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);
160: 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);
161: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
162: if (tao->reason != TAO_CONTINUE_ITERATING) break;
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);
168: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N - 1]]) {
169: /* reflect */
170: nm->nreflect++;
171: PetscInfo(0, "Reflect\n");
172: NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr);
173: } else if (fr < nm->f_values[nm->indices[0]]) {
174: /* expand */
175: nm->nexpand++;
176: PetscInfo(0, "Expand\n");
177: VecAXPBYPCZ(Xmue, 1 + nm->mu_e, -nm->mu_e, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
178: TaoComputeObjective(tao, Xmue, &fe);
179: if (fe < fr) {
180: NelderMeadReplace(nm, nm->indices[nm->N], Xmue, fe);
181: } else {
182: NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr);
183: }
184: } else if (nm->f_values[nm->indices[nm->N - 1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
185: /* outside contraction */
186: nm->noutcontract++;
187: PetscInfo(0, "Outside Contraction\n");
188: VecAXPBYPCZ(Xmuc, 1 + nm->mu_oc, -nm->mu_oc, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
190: TaoComputeObjective(tao, Xmuc, &fc);
191: if (fc <= fr) {
192: NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc);
193: } else shrink = 1;
194: } else {
195: /* inside contraction */
196: nm->nincontract++;
197: PetscInfo(0, "Inside Contraction\n");
198: VecAXPBYPCZ(Xmuc, 1 + nm->mu_ic, -nm->mu_ic, 0, Xbar, nm->simplex[nm->indices[nm->N]]);
199: TaoComputeObjective(tao, Xmuc, &fc);
200: if (fc < nm->f_values[nm->indices[nm->N]]) {
201: NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc);
202: } else shrink = 1;
203: }
205: if (shrink) {
206: nm->nshrink++;
207: PetscInfo(0, "Shrink\n");
209: for (i = 1; i < nm->N + 1; i++) {
210: VecAXPBY(nm->simplex[nm->indices[i]], 1.5, -0.5, nm->simplex[nm->indices[0]]);
211: TaoComputeObjective(tao, nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
212: }
213: VecAXPBY(Xbar, 1.5 * nm->oneOverN, -0.5, nm->simplex[nm->indices[0]]);
215: /* Add last vector's fraction of average */
216: VecAXPY(Xbar, nm->oneOverN, nm->simplex[nm->indices[nm->N]]);
217: NelderMeadSort(nm);
218: /* Subtract new last vector from average */
219: VecAXPY(Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]);
220: }
221: }
222: return 0;
223: }
225: /* ---------------------------------------------------------- */
226: /*MC
227: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
229: Options Database Keys:
230: + -tao_nm_lambda - initial step length
231: - -tao_nm_mu - expansion/contraction factor
233: Level: beginner
234: M*/
236: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
237: {
238: TAO_NelderMead *nm;
240: PetscNew(&nm);
241: tao->data = (void *)nm;
243: tao->ops->setup = TaoSetUp_NM;
244: tao->ops->solve = TaoSolve_NM;
245: tao->ops->view = TaoView_NM;
246: tao->ops->setfromoptions = TaoSetFromOptions_NM;
247: tao->ops->destroy = TaoDestroy_NM;
249: /* Override default settings (unless already changed) */
250: if (!tao->max_it_changed) tao->max_it = 2000;
251: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
253: nm->simplex = NULL;
254: nm->lambda = 1;
256: nm->mu_ic = -0.5;
257: nm->mu_oc = 0.5;
258: nm->mu_r = 1.0;
259: nm->mu_e = 2.0;
261: return 0;
262: }