Actual source code: neldermead.c
petsc-3.7.5 2017-01-01
1: #include <../src/tao/unconstrained/impls/neldermead/neldermead.h>
2: #include <petscvec.h>
4: static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm);
5: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f);
6: /* ---------------------------------------------------------- */
9: static PetscErrorCode TaoSetUp_NM(Tao tao)
10: {
12: TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
13: PetscInt n;
16: VecGetSize(tao->solution,&n);
17: nm->N = n;
18: nm->oneOverN = 1.0/n;
19: VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
20: PetscMalloc1(nm->N+1,&nm->f_values);
21: PetscMalloc1(nm->N+1,&nm->indices);
22: VecDuplicate(tao->solution,&nm->Xbar);
23: VecDuplicate(tao->solution,&nm->Xmur);
24: VecDuplicate(tao->solution,&nm->Xmue);
25: VecDuplicate(tao->solution,&nm->Xmuc);
27: tao->gradient=0;
28: tao->step=0;
29: return(0);
30: }
32: /* ---------------------------------------------------------- */
35: PetscErrorCode TaoDestroy_NM(Tao tao)
36: {
37: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
41: if (tao->setupcalled) {
42: VecDestroyVecs(nm->N+1,&nm->simplex);
43: VecDestroy(&nm->Xmuc);
44: VecDestroy(&nm->Xmue);
45: VecDestroy(&nm->Xmur);
46: VecDestroy(&nm->Xbar);
47: }
48: PetscFree(nm->indices);
49: PetscFree(nm->f_values);
50: PetscFree(tao->data);
51: tao->data = 0;
52: return(0);
53: }
55: /*------------------------------------------------------------*/
58: PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao)
59: {
60: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
64: PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");
65: PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL);
66: PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);
67: nm->mu_ic = -nm->mu_oc;
68: nm->mu_r = nm->mu_oc*2.0;
69: nm->mu_e = nm->mu_oc*4.0;
70: PetscOptionsTail();
71: return(0);
72: }
74: /*------------------------------------------------------------*/
77: PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
78: {
79: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
80: PetscBool isascii;
84: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
85: if (isascii) {
86: PetscViewerASCIIPushTab(viewer);
87: PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
88: PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
89: PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
90: PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
91: PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
92: PetscViewerASCIIPopTab(viewer);
93: }
94: return(0);
95: }
97: /*------------------------------------------------------------*/
100: PetscErrorCode TaoSolve_NM(Tao tao)
101: {
102: PetscErrorCode ierr;
103: TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
104: TaoConvergedReason reason;
105: PetscReal *x;
106: PetscInt i;
107: Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
108: PetscReal fr,fe,fc;
109: PetscInt shrink;
110: PetscInt low,high;
113: nm->nshrink = 0;
114: nm->nreflect = 0;
115: nm->nincontract = 0;
116: nm->noutcontract = 0;
117: nm->nexpand = 0;
119: if (tao->XL || tao->XU || tao->ops->computebounds) {
120: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
121: }
123: VecCopy(tao->solution,nm->simplex[0]);
124: TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
125: nm->indices[0]=0;
126: for (i=1;i<nm->N+1;i++){
127: VecCopy(tao->solution,nm->simplex[i]);
128: VecGetOwnershipRange(nm->simplex[i],&low,&high);
129: if (i-1 >= low && i-1 < high) {
130: VecGetArray(nm->simplex[i],&x);
131: x[i-1-low] += nm->lamda;
132: VecRestoreArray(nm->simplex[i],&x);
133: }
135: TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
136: nm->indices[i] = i;
137: }
139: /* Xbar = (Sum of all simplex vectors - worst vector)/N */
140: NelderMeadSort(nm);
141: VecSet(Xbar,0.0);
142: for (i=0;i<nm->N;i++) {
143: VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
144: }
145: VecScale(Xbar,nm->oneOverN);
146: reason = TAO_CONTINUE_ITERATING;
147: while (1) {
148: shrink = 0;
149: VecCopy(nm->simplex[nm->indices[0]],tao->solution);
150: 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,&reason);
151: if (reason != TAO_CONTINUE_ITERATING) break;
153: /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
154: VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
155: TaoComputeObjective(tao,Xmur,&fr);
157: if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
158: /* reflect */
159: nm->nreflect++;
160: PetscInfo(0,"Reflect\n");
161: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
162: } else if (fr < nm->f_values[nm->indices[0]]) {
163: /* expand */
164: nm->nexpand++;
165: PetscInfo(0,"Expand\n");
166: VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
167: TaoComputeObjective(tao,Xmue,&fe);
168: if (fe < fr) {
169: NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
170: } else {
171: NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
172: }
173: } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
174: /* outside contraction */
175: nm->noutcontract++;
176: PetscInfo(0,"Outside Contraction\n");
177: VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);
179: TaoComputeObjective(tao,Xmuc,&fc);
180: if (fc <= fr) {
181: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
182: } else shrink=1;
183: } else {
184: /* inside contraction */
185: nm->nincontract++;
186: PetscInfo(0,"Inside Contraction\n");
187: VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
188: TaoComputeObjective(tao,Xmuc,&fc);
189: if (fc < nm->f_values[nm->indices[nm->N]]) {
190: NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
191: } else shrink = 1;
192: }
194: if (shrink) {
195: nm->nshrink++;
196: PetscInfo(0,"Shrink\n");
198: for (i=1;i<nm->N+1;i++) {
199: VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
200: TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
201: }
202: VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);
204: /* Add last vector's fraction of average */
205: VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
206: NelderMeadSort(nm);
207: /* Subtract new last vector from average */
208: VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
209: }
210: }
211: return(0);
212: }
214: /* ---------------------------------------------------------- */
215: /*MC
216: TAONM - Nelder-Mead solver for derivative free, unconstrained minimization
218: Options Database Keys:
219: + -tao_nm_lamda - initial step length
220: . -tao_nm_mu - expansion/contraction factor
222: Level: beginner
223: M*/
227: PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao)
228: {
229: TAO_NelderMead *nm;
233: PetscNewLog(tao,&nm);
234: tao->data = (void*)nm;
236: tao->ops->setup = TaoSetUp_NM;
237: tao->ops->solve = TaoSolve_NM;
238: tao->ops->view = TaoView_NM;
239: tao->ops->setfromoptions = TaoSetFromOptions_NM;
240: tao->ops->destroy = TaoDestroy_NM;
242: /* Override default settings (unless already changed) */
243: if (!tao->max_it_changed) tao->max_it = 2000;
244: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
246: nm->simplex = 0;
247: nm->lamda = 1;
249: nm->mu_ic = -0.5;
250: nm->mu_oc = 0.5;
251: nm->mu_r = 1.0;
252: nm->mu_e = 2.0;
254: return(0);
255: }
258: /*------------------------------------------------------------*/
261: PetscErrorCode NelderMeadSort(TAO_NelderMead *nm)
262: {
263: PetscReal *values = nm->f_values;
264: PetscInt *indices = nm->indices;
265: PetscInt dim = nm->N+1;
266: PetscInt i,j,index;
267: PetscReal val;
270: for (i=1;i<dim;i++) {
271: index = indices[i];
272: val = values[index];
273: for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
274: indices[j+1] = indices[j];
275: }
276: indices[j+1] = index;
277: }
278: return(0);
279: }
282: /*------------------------------------------------------------*/
285: PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
286: {
290: /* Add new vector's fraction of average */
291: VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
292: VecCopy(Xmu,nm->simplex[index]);
293: nm->f_values[index] = f;
295: NelderMeadSort(nm);
297: /* Subtract last vector from average */
298: VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
299: return(0);
300: }