Actual source code: neldermead.c
1: #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
2: #include <petscvec.h>
4: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
5: {
6: PetscReal *values = nm->f_values;
7: PetscInt *indices = nm->indices;
8: PetscInt dim = nm->N + 1;
9: PetscInt i, j, index;
10: PetscReal val;
12: PetscFunctionBegin;
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: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
23: {
24: PetscFunctionBegin;
25: /* Add new vector's fraction of average */
26: PetscCall(VecAXPY(nm->Xbar, nm->oneOverN, Xmu));
27: PetscCall(VecCopy(Xmu, nm->simplex[index]));
28: nm->f_values[index] = f;
30: PetscCall(NelderMeadSort(nm));
32: /* Subtract last vector from average */
33: PetscCall(VecAXPY(nm->Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode TaoSetUp_NM(Tao tao)
38: {
39: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
40: PetscInt n;
42: PetscFunctionBegin;
43: PetscCall(VecGetSize(tao->solution, &n));
44: nm->N = n;
45: nm->oneOverN = 1.0 / n;
46: PetscCall(VecDuplicateVecs(tao->solution, nm->N + 1, &nm->simplex));
47: PetscCall(PetscMalloc1(nm->N + 1, &nm->f_values));
48: PetscCall(PetscMalloc1(nm->N + 1, &nm->indices));
49: PetscCall(VecDuplicate(tao->solution, &nm->Xbar));
50: PetscCall(VecDuplicate(tao->solution, &nm->Xmur));
51: PetscCall(VecDuplicate(tao->solution, &nm->Xmue));
52: PetscCall(VecDuplicate(tao->solution, &nm->Xmuc));
54: tao->gradient = NULL;
55: tao->step = 0;
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode TaoDestroy_NM(Tao tao)
60: {
61: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
63: PetscFunctionBegin;
64: if (tao->setupcalled) {
65: PetscCall(VecDestroyVecs(nm->N + 1, &nm->simplex));
66: PetscCall(VecDestroy(&nm->Xmuc));
67: PetscCall(VecDestroy(&nm->Xmue));
68: PetscCall(VecDestroy(&nm->Xmur));
69: PetscCall(VecDestroy(&nm->Xbar));
70: }
71: PetscCall(PetscFree(nm->indices));
72: PetscCall(PetscFree(nm->f_values));
73: PetscCall(PetscFree(tao->data));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode TaoSetFromOptions_NM(Tao tao, PetscOptionItems PetscOptionsObject)
78: {
79: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
81: PetscFunctionBegin;
82: PetscOptionsHeadBegin(PetscOptionsObject, "Nelder-Mead options");
83: PetscCall(PetscOptionsDeprecated("-tao_nm_lamda", "-tao_nm_lambda", "3.18.4", NULL));
84: PetscCall(PetscOptionsReal("-tao_nm_lambda", "initial step length", "", nm->lambda, &nm->lambda, NULL));
85: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode TaoView_NM(Tao tao, PetscViewer viewer)
94: {
95: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
96: PetscBool isascii;
98: PetscFunctionBegin;
99: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
100: if (isascii) {
101: PetscCall(PetscViewerASCIIPushTab(viewer));
102: PetscCall(PetscViewerASCIIPrintf(viewer, "expansions: %" PetscInt_FMT "\n", nm->nexpand));
103: PetscCall(PetscViewerASCIIPrintf(viewer, "reflections: %" PetscInt_FMT "\n", nm->nreflect));
104: PetscCall(PetscViewerASCIIPrintf(viewer, "inside contractions: %" PetscInt_FMT "\n", nm->nincontract));
105: PetscCall(PetscViewerASCIIPrintf(viewer, "outside contractionss: %" PetscInt_FMT "\n", nm->noutcontract));
106: PetscCall(PetscViewerASCIIPrintf(viewer, "Shrink steps: %" PetscInt_FMT "\n", nm->nshrink));
107: PetscCall(PetscViewerASCIIPopTab(viewer));
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
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: PetscFunctionBegin;
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) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n"));
131: PetscCall(VecCopy(tao->solution, nm->simplex[0]));
132: PetscCall(TaoComputeObjective(tao, nm->simplex[0], &nm->f_values[0]));
133: nm->indices[0] = 0;
134: for (i = 1; i < nm->N + 1; i++) {
135: PetscCall(VecCopy(tao->solution, nm->simplex[i]));
136: PetscCall(VecGetOwnershipRange(nm->simplex[i], &low, &high));
137: if (i - 1 >= low && i - 1 < high) {
138: PetscCall(VecGetArray(nm->simplex[i], &x));
139: x[i - 1 - low] += nm->lambda;
140: PetscCall(VecRestoreArray(nm->simplex[i], &x));
141: }
143: PetscCall(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: PetscCall(NelderMeadSort(nm));
149: PetscCall(VecSet(Xbar, 0.0));
150: for (i = 0; i < nm->N; i++) PetscCall(VecAXPY(Xbar, 1.0, nm->simplex[nm->indices[i]]));
151: PetscCall(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: PetscCall(VecCopy(nm->simplex[nm->indices[0]], tao->solution));
159: 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));
160: 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));
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: PetscCall(VecAXPBYPCZ(Xmur, 1 + nm->mu_r, -nm->mu_r, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
166: PetscCall(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: PetscCall(PetscInfo(0, "Reflect\n"));
172: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr));
173: } else if (fr < nm->f_values[nm->indices[0]]) {
174: /* expand */
175: nm->nexpand++;
176: PetscCall(PetscInfo(0, "Expand\n"));
177: PetscCall(VecAXPBYPCZ(Xmue, 1 + nm->mu_e, -nm->mu_e, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
178: PetscCall(TaoComputeObjective(tao, Xmue, &fe));
179: if (fe < fr) {
180: PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmue, fe));
181: } else {
182: PetscCall(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: PetscCall(PetscInfo(0, "Outside Contraction\n"));
188: PetscCall(VecAXPBYPCZ(Xmuc, 1 + nm->mu_oc, -nm->mu_oc, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
190: PetscCall(TaoComputeObjective(tao, Xmuc, &fc));
191: if (fc <= fr) PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc));
192: else shrink = 1;
193: } else {
194: /* inside contraction */
195: nm->nincontract++;
196: PetscCall(PetscInfo(0, "Inside Contraction\n"));
197: PetscCall(VecAXPBYPCZ(Xmuc, 1 + nm->mu_ic, -nm->mu_ic, 0, Xbar, nm->simplex[nm->indices[nm->N]]));
198: PetscCall(TaoComputeObjective(tao, Xmuc, &fc));
199: if (fc < nm->f_values[nm->indices[nm->N]]) PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc));
200: else shrink = 1;
201: }
203: if (shrink) {
204: nm->nshrink++;
205: PetscCall(PetscInfo(0, "Shrink\n"));
207: for (i = 1; i < nm->N + 1; i++) {
208: PetscCall(VecAXPBY(nm->simplex[nm->indices[i]], 1.5, -0.5, nm->simplex[nm->indices[0]]));
209: PetscCall(TaoComputeObjective(tao, nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]));
210: }
211: PetscCall(VecAXPBY(Xbar, 1.5 * nm->oneOverN, -0.5, nm->simplex[nm->indices[0]]));
213: /* Add last vector's fraction of average */
214: PetscCall(VecAXPY(Xbar, nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
215: PetscCall(NelderMeadSort(nm));
216: /* Subtract new last vector from average */
217: PetscCall(VecAXPY(Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]]));
218: }
219: }
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*MC
224: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
226: Options Database Keys:
227: + -tao_nm_lambda - initial step length
228: - -tao_nm_mu - expansion/contraction factor
230: Level: beginner
231: M*/
233: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
234: {
235: TAO_NelderMead *nm;
237: PetscFunctionBegin;
238: PetscCall(PetscNew(&nm));
239: tao->data = (void *)nm;
241: tao->ops->setup = TaoSetUp_NM;
242: tao->ops->solve = TaoSolve_NM;
243: tao->ops->view = TaoView_NM;
244: tao->ops->setfromoptions = TaoSetFromOptions_NM;
245: tao->ops->destroy = TaoDestroy_NM;
247: /* Override default settings (unless already changed) */
248: PetscCall(TaoParametersInitialize(tao));
249: PetscObjectParameterSetDefault(tao, max_it, 2000);
250: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
252: nm->simplex = NULL;
253: nm->lambda = 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;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }