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: PetscFunctionBegin;
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--) indices[j + 1] = indices[j];
18: indices[j + 1] = index;
19: }
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: /*------------------------------------------------------------*/
24: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
25: {
26: PetscFunctionBegin;
27: /* Add new vector's fraction of average */
28: PetscCall(VecAXPY(nm->Xbar, nm->oneOverN, Xmu));
29: PetscCall(VecCopy(Xmu, nm->simplex[index]));
30: nm->f_values[index] = f;
32: PetscCall(NelderMeadSort(nm));
34: /* Subtract last vector from average */
35: PetscCall(VecAXPY(nm->Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: /* ---------------------------------------------------------- */
40: static PetscErrorCode TaoSetUp_NM(Tao tao)
41: {
42: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
43: PetscInt n;
45: PetscFunctionBegin;
46: PetscCall(VecGetSize(tao->solution, &n));
47: nm->N = n;
48: nm->oneOverN = 1.0 / n;
49: PetscCall(VecDuplicateVecs(tao->solution, nm->N + 1, &nm->simplex));
50: PetscCall(PetscMalloc1(nm->N + 1, &nm->f_values));
51: PetscCall(PetscMalloc1(nm->N + 1, &nm->indices));
52: PetscCall(VecDuplicate(tao->solution, &nm->Xbar));
53: PetscCall(VecDuplicate(tao->solution, &nm->Xmur));
54: PetscCall(VecDuplicate(tao->solution, &nm->Xmue));
55: PetscCall(VecDuplicate(tao->solution, &nm->Xmuc));
57: tao->gradient = NULL;
58: tao->step = 0;
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: /* ---------------------------------------------------------- */
63: static PetscErrorCode TaoDestroy_NM(Tao tao)
64: {
65: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
67: PetscFunctionBegin;
68: if (tao->setupcalled) {
69: PetscCall(VecDestroyVecs(nm->N + 1, &nm->simplex));
70: PetscCall(VecDestroy(&nm->Xmuc));
71: PetscCall(VecDestroy(&nm->Xmue));
72: PetscCall(VecDestroy(&nm->Xmur));
73: PetscCall(VecDestroy(&nm->Xbar));
74: }
75: PetscCall(PetscFree(nm->indices));
76: PetscCall(PetscFree(nm->f_values));
77: PetscCall(PetscFree(tao->data));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*------------------------------------------------------------*/
82: static PetscErrorCode TaoSetFromOptions_NM(Tao tao, PetscOptionItems *PetscOptionsObject)
83: {
84: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
86: PetscFunctionBegin;
87: PetscOptionsHeadBegin(PetscOptionsObject, "Nelder-Mead options");
88: PetscCall(PetscOptionsDeprecated("-tao_nm_lamda", "-tao_nm_lambda", "3.18.4", NULL));
89: PetscCall(PetscOptionsReal("-tao_nm_lambda", "initial step length", "", nm->lambda, &nm->lambda, NULL));
90: PetscCall(PetscOptionsReal("-tao_nm_mu", "mu", "", nm->mu_oc, &nm->mu_oc, NULL));
91: nm->mu_ic = -nm->mu_oc;
92: nm->mu_r = nm->mu_oc * 2.0;
93: nm->mu_e = nm->mu_oc * 4.0;
94: PetscOptionsHeadEnd();
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: /*------------------------------------------------------------*/
99: static PetscErrorCode TaoView_NM(Tao tao, PetscViewer viewer)
100: {
101: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
102: PetscBool isascii;
104: PetscFunctionBegin;
105: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
106: if (isascii) {
107: PetscCall(PetscViewerASCIIPushTab(viewer));
108: PetscCall(PetscViewerASCIIPrintf(viewer, "expansions: %" PetscInt_FMT "\n", nm->nexpand));
109: PetscCall(PetscViewerASCIIPrintf(viewer, "reflections: %" PetscInt_FMT "\n", nm->nreflect));
110: PetscCall(PetscViewerASCIIPrintf(viewer, "inside contractions: %" PetscInt_FMT "\n", nm->nincontract));
111: PetscCall(PetscViewerASCIIPrintf(viewer, "outside contractionss: %" PetscInt_FMT "\n", nm->noutcontract));
112: PetscCall(PetscViewerASCIIPrintf(viewer, "Shrink steps: %" PetscInt_FMT "\n", nm->nshrink));
113: PetscCall(PetscViewerASCIIPopTab(viewer));
114: }
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: /*------------------------------------------------------------*/
119: static PetscErrorCode TaoSolve_NM(Tao tao)
120: {
121: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
122: PetscReal *x;
123: PetscInt i;
124: Vec Xmur = nm->Xmur, Xmue = nm->Xmue, Xmuc = nm->Xmuc, Xbar = nm->Xbar;
125: PetscReal fr, fe, fc;
126: PetscInt shrink;
127: PetscInt low, high;
129: PetscFunctionBegin;
130: nm->nshrink = 0;
131: nm->nreflect = 0;
132: nm->nincontract = 0;
133: nm->noutcontract = 0;
134: nm->nexpand = 0;
136: if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n"));
138: PetscCall(VecCopy(tao->solution, nm->simplex[0]));
139: PetscCall(TaoComputeObjective(tao, nm->simplex[0], &nm->f_values[0]));
140: nm->indices[0] = 0;
141: for (i = 1; i < nm->N + 1; i++) {
142: PetscCall(VecCopy(tao->solution, nm->simplex[i]));
143: PetscCall(VecGetOwnershipRange(nm->simplex[i], &low, &high));
144: if (i - 1 >= low && i - 1 < high) {
145: PetscCall(VecGetArray(nm->simplex[i], &x));
146: x[i - 1 - low] += nm->lambda;
147: PetscCall(VecRestoreArray(nm->simplex[i], &x));
148: }
150: PetscCall(TaoComputeObjective(tao, nm->simplex[i], &nm->f_values[i]));
151: nm->indices[i] = i;
152: }
154: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
155: PetscCall(NelderMeadSort(nm));
156: PetscCall(VecSet(Xbar, 0.0));
157: for (i = 0; i < nm->N; i++) PetscCall(VecAXPY(Xbar, 1.0, nm->simplex[nm->indices[i]]));
158: PetscCall(VecScale(Xbar, nm->oneOverN));
159: tao->reason = TAO_CONTINUE_ITERATING;
160: while (1) {
161: /* Call general purpose update function */
162: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
163: ++tao->niter;
164: shrink = 0;
165: PetscCall(VecCopy(nm->simplex[nm->indices[0]], tao->solution));
166: PetscCall(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: PetscCall(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: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
169: if (tao->reason != TAO_CONTINUE_ITERATING) break;
171: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
172: PetscCall(VecAXPBYPCZ(Xmur, 1 + nm->mu_r, -nm->mu_r, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
173: PetscCall(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: PetscCall(PetscInfo(0, "Reflect\n"));
179: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr));
180: } else if (fr < nm->f_values[nm->indices[0]]) {
181: /* expand */
182: nm->nexpand++;
183: PetscCall(PetscInfo(0, "Expand\n"));
184: PetscCall(VecAXPBYPCZ(Xmue, 1 + nm->mu_e, -nm->mu_e, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
185: PetscCall(TaoComputeObjective(tao, Xmue, &fe));
186: if (fe < fr) {
187: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmue, fe));
188: } else {
189: PetscCall(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: PetscCall(PetscInfo(0, "Outside Contraction\n"));
195: PetscCall(VecAXPBYPCZ(Xmuc, 1 + nm->mu_oc, -nm->mu_oc, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
197: PetscCall(TaoComputeObjective(tao, Xmuc, &fc));
198: if (fc <= fr) {
199: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc));
200: } else shrink = 1;
201: } else {
202: /* inside contraction */
203: nm->nincontract++;
204: PetscCall(PetscInfo(0, "Inside Contraction\n"));
205: PetscCall(VecAXPBYPCZ(Xmuc, 1 + nm->mu_ic, -nm->mu_ic, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
206: PetscCall(TaoComputeObjective(tao, Xmuc, &fc));
207: if (fc < nm->f_values[nm->indices[nm->N]]) {
208: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc));
209: } else shrink = 1;
210: }
212: if (shrink) {
213: nm->nshrink++;
214: PetscCall(PetscInfo(0, "Shrink\n"));
216: for (i = 1; i < nm->N + 1; i++) {
217: PetscCall(VecAXPBY(nm->simplex[nm->indices[i]], 1.5, -0.5, nm->simplex[nm->indices[0]]));
218: PetscCall(TaoComputeObjective(tao, nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]));
219: }
220: PetscCall(VecAXPBY(Xbar, 1.5 * nm->oneOverN, -0.5, nm->simplex[nm->indices[0]]));
222: /* Add last vector's fraction of average */
223: PetscCall(VecAXPY(Xbar, nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
224: PetscCall(NelderMeadSort(nm));
225: /* Subtract new last vector from average */
226: PetscCall(VecAXPY(Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
227: }
228: }
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /* ---------------------------------------------------------- */
233: /*MC
234: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
236: Options Database Keys:
237: + -tao_nm_lambda - 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: PetscFunctionBegin;
248: PetscCall(PetscNew(&nm));
249: tao->data = (void *)nm;
251: tao->ops->setup = TaoSetUp_NM;
252: tao->ops->solve = TaoSolve_NM;
253: tao->ops->view = TaoView_NM;
254: tao->ops->setfromoptions = TaoSetFromOptions_NM;
255: tao->ops->destroy = TaoDestroy_NM;
257: /* Override default settings (unless already changed) */
258: PetscCall(TaoParametersInitialize(tao));
259: PetscObjectParameterSetDefault(tao, max_it, 2000);
260: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
262: nm->simplex = NULL;
263: nm->lambda = 1;
265: nm->mu_ic = -0.5;
266: nm->mu_oc = 0.5;
267: nm->mu_r = 1.0;
268: nm->mu_e = 2.0;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }