Actual source code: neldermead.c
petsc-3.12.5 2020-03-29
1: #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
2: #include <petscvec.h>
5: /*------------------------------------------------------------*/
6: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
7: {
8: PetscReal *values = nm->f_values;
9: PetscInt *indices = nm->indices;
10: PetscInt dim = nm->N+1;
11: PetscInt i,j,index;
12: PetscReal val;
15: for (i=1;i<dim;i++) {
16: index = indices[i];
17: val = values[index];
18: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
19: indices[j+1] = indices[j];
20: }
21: indices[j+1] = index;
22: }
23: return(0);
24: }
27: /*------------------------------------------------------------*/
28: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
29: {
33: /* Add new vector's fraction of average */
34: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
35: VecCopy(Xmu,nm->simplex[index]);
36: nm->f_values[index] = f;
38: NelderMeadSort(nm);
40: /* Subtract last vector from average */
41: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
42: return(0);
43: }
45: /* ---------------------------------------------------------- */
46: static PetscErrorCode TaoSetUp_NM(Tao tao)
47: {
49: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
50: PetscInt n;
53: VecGetSize(tao->solution,&n);
54: nm->N = n;
55: nm->oneOverN = 1.0/n;
56: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
57: PetscMalloc1(nm->N+1,&nm->f_values);
58: PetscMalloc1(nm->N+1,&nm->indices);
59: VecDuplicate(tao->solution,&nm->Xbar);
60: VecDuplicate(tao->solution,&nm->Xmur);
61: VecDuplicate(tao->solution,&nm->Xmue);
62: VecDuplicate(tao->solution,&nm->Xmuc);
64: tao->gradient=0;
65: tao->step=0;
66: return(0);
67: }
69: /* ---------------------------------------------------------- */
70: static PetscErrorCode TaoDestroy_NM(Tao tao)
71: {
72: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
76: if (tao->setupcalled) {
77: VecDestroyVecs(nm->N+1,&nm->simplex);
78: VecDestroy(&nm->Xmuc);
79: VecDestroy(&nm->Xmue);
80: VecDestroy(&nm->Xmur);
81: VecDestroy(&nm->Xbar);
82: }
83: PetscFree(nm->indices);
84: PetscFree(nm->f_values);
85: PetscFree(tao->data);
86: tao->data = 0;
87: return(0);
88: }
90: /*------------------------------------------------------------*/
91: static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao)
92: {
93: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
97: PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");
98: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL);
99: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);
100: nm->mu_ic = -nm->mu_oc;
101: nm->mu_r = nm->mu_oc*2.0;
102: nm->mu_e = nm->mu_oc*4.0;
103: PetscOptionsTail();
104: return(0);
105: }
107: /*------------------------------------------------------------*/
108: static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
109: {
110: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
111: PetscBool isascii;
115: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
116: if (isascii) {
117: PetscViewerASCIIPushTab(viewer);
118: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
119: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
120: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
121: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
122: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
123: PetscViewerASCIIPopTab(viewer);
124: }
125: return(0);
126: }
128: /*------------------------------------------------------------*/
129: static PetscErrorCode TaoSolve_NM(Tao tao)
130: {
131: PetscErrorCode ierr;
132: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
133: PetscReal *x;
134: PetscInt i;
135: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
136: PetscReal fr,fe,fc;
137: PetscInt shrink;
138: PetscInt low,high;
141: nm->nshrink = 0;
142: nm->nreflect = 0;
143: nm->nincontract = 0;
144: nm->noutcontract = 0;
145: nm->nexpand = 0;
147: if (tao->XL || tao->XU || tao->ops->computebounds) {
148: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
149: }
151: VecCopy(tao->solution,nm->simplex[0]);
152: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
153: nm->indices[0]=0;
154: for (i=1;i<nm->N+1;i++){
155: VecCopy(tao->solution,nm->simplex[i]);
156: VecGetOwnershipRange(nm->simplex[i],&low,&high);
157: if (i-1 >= low && i-1 < high) {
158: VecGetArray(nm->simplex[i],&x);
159: x[i-1-low] += nm->lamda;
160: VecRestoreArray(nm->simplex[i],&x);
161: }
163: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
164: nm->indices[i] = i;
165: }
167: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
168: NelderMeadSort(nm);
169: VecSet(Xbar,0.0);
170: for (i=0;i<nm->N;i++) {
171: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
172: }
173: VecScale(Xbar,nm->oneOverN);
174: tao->reason = TAO_CONTINUE_ITERATING;
175: while (1) {
176: /* Call general purpose update function */
177: if (tao->ops->update) {
178: (*tao->ops->update)(tao, tao->niter, tao->user_update);
179: }
180: ++tao->niter;
181: shrink = 0;
182: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
183: 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);
184: 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);
185: (*tao->ops->convergencetest)(tao,tao->cnvP);
186: if (tao->reason != TAO_CONTINUE_ITERATING) break;
188: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
189: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
190: TaoComputeObjective(tao,Xmur,&fr);
192: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
193: /* reflect */
194: nm->nreflect++;
195: PetscInfo(0,"Reflect\n");
196: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
197: } else if (fr < nm->f_values[nm->indices[0]]) {
198: /* expand */
199: nm->nexpand++;
200: PetscInfo(0,"Expand\n");
201: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
202: TaoComputeObjective(tao,Xmue,&fe);
203: if (fe < fr) {
204: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
205: } else {
206: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
207: }
208: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
209: /* outside contraction */
210: nm->noutcontract++;
211: PetscInfo(0,"Outside Contraction\n");
212: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
214: TaoComputeObjective(tao,Xmuc,&fc);
215: if (fc <= fr) {
216: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
217: } else shrink=1;
218: } else {
219: /* inside contraction */
220: nm->nincontract++;
221: PetscInfo(0,"Inside Contraction\n");
222: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
223: TaoComputeObjective(tao,Xmuc,&fc);
224: if (fc < nm->f_values[nm->indices[nm->N]]) {
225: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
226: } else shrink = 1;
227: }
229: if (shrink) {
230: nm->nshrink++;
231: PetscInfo(0,"Shrink\n");
233: for (i=1;i<nm->N+1;i++) {
234: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
235: TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
236: }
237: VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);
239: /* Add last vector's fraction of average */
240: VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
241: NelderMeadSort(nm);
242: /* Subtract new last vector from average */
243: VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
244: }
245: }
246: return(0);
247: }
249: /* ---------------------------------------------------------- */
250: /*MC
251: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
253: Options Database Keys:
254: + -tao_nm_lamda - initial step length
255: - -tao_nm_mu - expansion/contraction factor
257: Level: beginner
258: M*/
260: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
261: {
262: TAO_NelderMead *nm;
266: PetscNewLog(tao,&nm);
267: tao->data = (void*)nm;
269: tao->ops->setup = TaoSetUp_NM;
270: tao->ops->solve = TaoSolve_NM;
271: tao->ops->view = TaoView_NM;
272: tao->ops->setfromoptions = TaoSetFromOptions_NM;
273: tao->ops->destroy = TaoDestroy_NM;
275: /* Override default settings (unless already changed) */
276: if (!tao->max_it_changed) tao->max_it = 2000;
277: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
279: nm->simplex = 0;
280: nm->lamda = 1;
282: nm->mu_ic = -0.5;
283: nm->mu_oc = 0.5;
284: nm->mu_r = 1.0;
285: nm->mu_e = 2.0;
287: return(0);
288: }