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--) {
 17:       indices[j+1] = indices[j];
 18:     }
 19:     indices[j+1] = index;
 20:   }
 21:   return 0;
 22: }

 24: /*------------------------------------------------------------*/
 25: static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
 26: {
 27:   /*  Add new vector's fraction of average */
 28:   VecAXPY(nm->Xbar,nm->oneOverN,Xmu);
 29:   VecCopy(Xmu,nm->simplex[index]);
 30:   nm->f_values[index] = f;

 32:   NelderMeadSort(nm);

 34:   /*  Subtract last vector from average */
 35:   VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
 36:   return 0;
 37: }

 39: /* ---------------------------------------------------------- */
 40: static PetscErrorCode TaoSetUp_NM(Tao tao)
 41: {
 42:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
 43:   PetscInt       n;

 45:   VecGetSize(tao->solution,&n);
 46:   nm->N = n;
 47:   nm->oneOverN = 1.0/n;
 48:   VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex);
 49:   PetscMalloc1(nm->N+1,&nm->f_values);
 50:   PetscMalloc1(nm->N+1,&nm->indices);
 51:   VecDuplicate(tao->solution,&nm->Xbar);
 52:   VecDuplicate(tao->solution,&nm->Xmur);
 53:   VecDuplicate(tao->solution,&nm->Xmue);
 54:   VecDuplicate(tao->solution,&nm->Xmuc);

 56:   tao->gradient=NULL;
 57:   tao->step=0;
 58:   return 0;
 59: }

 61: /* ---------------------------------------------------------- */
 62: static PetscErrorCode TaoDestroy_NM(Tao tao)
 63: {
 64:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;

 66:   if (tao->setupcalled) {
 67:     VecDestroyVecs(nm->N+1,&nm->simplex);
 68:     VecDestroy(&nm->Xmuc);
 69:     VecDestroy(&nm->Xmue);
 70:     VecDestroy(&nm->Xmur);
 71:     VecDestroy(&nm->Xbar);
 72:   }
 73:   PetscFree(nm->indices);
 74:   PetscFree(nm->f_values);
 75:   PetscFree(tao->data);
 76:   return 0;
 77: }

 79: /*------------------------------------------------------------*/
 80: static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao)
 81: {
 82:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;

 84:   PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");
 85:   PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL);
 86:   PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);
 87:   nm->mu_ic = -nm->mu_oc;
 88:   nm->mu_r = nm->mu_oc*2.0;
 89:   nm->mu_e = nm->mu_oc*4.0;
 90:   PetscOptionsTail();
 91:   return 0;
 92: }

 94: /*------------------------------------------------------------*/
 95: static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer)
 96: {
 97:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
 98:   PetscBool      isascii;

100:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
101:   if (isascii) {
102:     PetscViewerASCIIPushTab(viewer);
103:     PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);
104:     PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);
105:     PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);
106:     PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);
107:     PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);
108:     PetscViewerASCIIPopTab(viewer);
109:   }
110:   return 0;
111: }

113: /*------------------------------------------------------------*/
114: static PetscErrorCode TaoSolve_NM(Tao tao)
115: {
116:   TAO_NelderMead     *nm = (TAO_NelderMead*)tao->data;
117:   PetscReal          *x;
118:   PetscInt           i;
119:   Vec                Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
120:   PetscReal          fr,fe,fc;
121:   PetscInt           shrink;
122:   PetscInt           low,high;

124:   nm->nshrink =      0;
125:   nm->nreflect =     0;
126:   nm->nincontract =  0;
127:   nm->noutcontract = 0;
128:   nm->nexpand =      0;

130:   if (tao->XL || tao->XU || tao->ops->computebounds) {
131:     PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");
132:   }

134:   VecCopy(tao->solution,nm->simplex[0]);
135:   TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);
136:   nm->indices[0]=0;
137:   for (i=1;i<nm->N+1;i++) {
138:     VecCopy(tao->solution,nm->simplex[i]);
139:     VecGetOwnershipRange(nm->simplex[i],&low,&high);
140:     if (i-1 >= low && i-1 < high) {
141:       VecGetArray(nm->simplex[i],&x);
142:       x[i-1-low] += nm->lamda;
143:       VecRestoreArray(nm->simplex[i],&x);
144:     }

146:     TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);
147:     nm->indices[i] = i;
148:   }

150:   /*  Xbar  = (Sum of all simplex vectors - worst vector)/N */
151:   NelderMeadSort(nm);
152:   VecSet(Xbar,0.0);
153:   for (i=0;i<nm->N;i++) {
154:     VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
155:   }
156:   VecScale(Xbar,nm->oneOverN);
157:   tao->reason = TAO_CONTINUE_ITERATING;
158:   while (1) {
159:     /* Call general purpose update function */
160:     if (tao->ops->update) {
161:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
162:     }
163:     ++tao->niter;
164:     shrink = 0;
165:     VecCopy(nm->simplex[nm->indices[0]],tao->solution);
166:     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:     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:     (*tao->ops->convergencetest)(tao,tao->cnvP);
169:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

171:     /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
172:     VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);
173:     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:       PetscInfo(0,"Reflect\n");
179:       NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);
180:     } else if (fr < nm->f_values[nm->indices[0]]) {
181:       /*  expand */
182:       nm->nexpand++;
183:       PetscInfo(0,"Expand\n");
184:       VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);
185:       TaoComputeObjective(tao,Xmue,&fe);
186:       if (fe < fr) {
187:         NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);
188:       } else {
189:         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:       PetscInfo(0,"Outside Contraction\n");
195:       VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);

197:       TaoComputeObjective(tao,Xmuc,&fc);
198:       if (fc <= fr) {
199:         NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
200:       } else shrink=1;
201:     } else {
202:       /* inside contraction */
203:       nm->nincontract++;
204:       PetscInfo(0,"Inside Contraction\n");
205:       VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);
206:       TaoComputeObjective(tao,Xmuc,&fc);
207:       if (fc < nm->f_values[nm->indices[nm->N]]) {
208:         NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);
209:       } else shrink = 1;
210:     }

212:     if (shrink) {
213:       nm->nshrink++;
214:       PetscInfo(0,"Shrink\n");

216:       for (i=1;i<nm->N+1;i++) {
217:         VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
218:         TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);
219:       }
220:       VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);

222:       /*  Add last vector's fraction of average */
223:       VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
224:       NelderMeadSort(nm);
225:       /*  Subtract new last vector from average */
226:       VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);
227:     }
228:   }
229:   return 0;
230: }

232: /* ---------------------------------------------------------- */
233: /*MC
234:  TAONM - Nelder-Mead solver for derivative free, unconstrained minimization

236:  Options Database Keys:
237: + -tao_nm_lamda - 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:   PetscNewLog(tao,&nm);
248:   tao->data = (void*)nm;

250:   tao->ops->setup = TaoSetUp_NM;
251:   tao->ops->solve = TaoSolve_NM;
252:   tao->ops->view = TaoView_NM;
253:   tao->ops->setfromoptions = TaoSetFromOptions_NM;
254:   tao->ops->destroy = TaoDestroy_NM;

256:   /* Override default settings (unless already changed) */
257:   if (!tao->max_it_changed) tao->max_it = 2000;
258:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

260:   nm->simplex = NULL;
261:   nm->lamda = 1;

263:   nm->mu_ic = -0.5;
264:   nm->mu_oc = 0.5;
265:   nm->mu_r = 1.0;
266:   nm->mu_e = 2.0;

268:   return 0;
269: }