Actual source code: neldermead.c

petsc-3.11.3 2019-06-26
Report Typos and Errors
  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: }