Actual source code: neldermead.c

  1: #include "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(TaoSolver tao)
 10: {
 12:   TAO_NelderMead *nm = (TAO_NelderMead *)tao->data;
 13:   PetscInt size;

 16:   VecGetSize(tao->solution,&size); 
 17:   nm->N = size;
 18:   nm->oneOverN = 1.0/size;
 19:   VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex); 
 20:   PetscMalloc((nm->N+1)*sizeof(PetscReal),&nm->f_values); 
 21:   PetscMalloc((nm->N+1)*sizeof(PetscInt),&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;

 30:   return(0);
 31: }

 33: /* ---------------------------------------------------------- */
 36: PetscErrorCode TaoDestroy_NM(TaoSolver tao)
 37: {
 38:   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(TaoSolver tao)
 59: {
 60:   
 61:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
 62:   PetscBool flg;
 64:   
 66:   PetscOptionsHead("Nelder-Mead options"); 

 68:   PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,&flg);  

 70:   PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,&flg); 
 71:   nm->mu_ic = -nm->mu_oc;
 72:   nm->mu_r = nm->mu_oc*2.0;
 73:   nm->mu_e = nm->mu_oc*4.0;

 75:   PetscOptionsTail(); 
 76:   return(0);
 77: }

 79: /*------------------------------------------------------------*/
 82: PetscErrorCode TaoView_NM(TaoSolver tao,PetscViewer viewer)
 83: {
 84:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
 85:   PetscBool isascii;
 87:   
 89:   PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 90:   if (isascii) {
 91:     PetscViewerASCIIPushTab(viewer); 
 92:     PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand); 
 93:     PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect); 
 94:     PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract); 
 95:     PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract); 
 96:     PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink); 
 97:     PetscViewerASCIIPopTab(viewer); 
 98:   } else {
 99:     SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO NelderMead",((PetscObject)viewer)->type_name);
100:   }

102:   return(0);
103: }

105: /*------------------------------------------------------------*/
108: PetscErrorCode TaoSolve_NM(TaoSolver tao)
109: {
111:   TAO_NelderMead *nm = (TAO_NelderMead*)tao->data;
112:   TaoSolverTerminationReason reason;
113:   PetscReal *x;
114:   PetscInt iter=0,i;
115:   Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar;
116:   PetscReal fr,fe,fc;
117:   PetscInt shrink;
118:   PetscInt low,high;
119:   
120:   
122:   nm->nshrink =      0;
123:   nm->nreflect =     0;
124:   nm->nincontract =  0;
125:   nm->noutcontract = 0;
126:   nm->nexpand =      0;
127:   
128:   if (tao->XL || tao->XU || tao->ops->computebounds) {
129:     PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n"); 
130:   }

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

144:     TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]); 
145:     nm->indices[i] = i;
146:   }
147:   
148:   /*  Xbar  = (Sum of all simplex vectors - worst vector)/N */
149:   NelderMeadSort(nm); 
150:   VecSet(Xbar,0.0); 
151:   for (i=0;i<nm->N;i++) {
152:     VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);
153:   }
154:   VecScale(Xbar,nm->oneOverN);
155:   reason = TAO_CONTINUE_ITERATING;
156:   while (1) {
157:     shrink = 0;
158:     VecCopy(nm->simplex[nm->indices[0]],tao->solution); 
159:     TaoMonitor(tao,iter++,nm->f_values[nm->indices[0]],nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]],0.0,1.0,&reason); 
160:     if (reason != TAO_CONTINUE_ITERATING) break;

162:     
163:     
164:     /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */
165:     VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]); 
166:     TaoComputeObjective(tao,Xmur,&fr); 


169:     if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
170:       /*  reflect */
171:       nm->nreflect++;
172:       PetscInfo(0,"Reflect\n"); 
173:       NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr); 
174:     }

176:     else if (fr < nm->f_values[nm->indices[0]]) {
177:       /*  expand */
178:       nm->nexpand++;
179:       PetscInfo(0,"Expand\n"); 
180:       VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]); 
181:       TaoComputeObjective(tao,Xmue,&fe); 
182:       if (fe < fr) {
183:         NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe); 
184:       } else {
185:         NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr); 
186:       }

188:     } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
189:       /* outside contraction */
190:       nm->noutcontract++;
191:       PetscInfo(0,"Outside Contraction\n"); 
192:       VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]); 

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

211:     if (shrink) {
212:       nm->nshrink++;
213:       PetscInfo(0,"Shrink\n"); 
214:       
215:       for (i=1;i<nm->N+1;i++) {
216:           VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);
217:         TaoComputeObjective(tao,nm->simplex[nm->indices[i]],
218:                                  &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:     
230:   }
231:   return(0);
232: }

234: /* ---------------------------------------------------------- */
238: PetscErrorCode TaoCreate_NM(TaoSolver tao)
239: {
240:   TAO_NelderMead *nm;

244:   PetscNewLog(tao,TAO_NelderMead,&nm); 
245:   tao->data = (void*)nm;

247:   tao->ops->setup = TaoSetUp_NM;
248:   tao->ops->solve = TaoSolve_NM;
249:   tao->ops->view = TaoView_NM;
250:   tao->ops->setfromoptions = TaoSetFromOptions_NM;
251:   tao->ops->destroy = TaoDestroy_NM;

253:   tao->max_it = 2000;
254:   tao->max_funcs = 4000;
255:   tao->fatol = 1e-8;
256:   tao->frtol = 1e-8;

258:   nm->simplex = 0;
259:   nm->lamda = 1;

261:   nm->mu_ic = -0.5;
262:   nm->mu_oc = 0.5;
263:   nm->mu_r = 1.0;
264:   nm->mu_e = 2.0;

266:   return(0);
267: }
269:     

271: /*------------------------------------------------------------*/
274: PetscErrorCode NelderMeadSort(TAO_NelderMead *nm) {
275:   PetscReal *values = nm->f_values;
276:   PetscInt *indices = nm->indices;
277:   PetscInt dim = nm->N+1;

279:   PetscInt i,j,index;
280:   PetscReal val;
282:   for (i=1;i<dim;i++) {
283:     index = indices[i];
284:     val = values[index];
285:     for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
286:       indices[j+1] = indices[j];
287:     }
288:     indices[j+1] = index;
289:   }  
290:   return(0);
291: }


294: /*------------------------------------------------------------*/
297: PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f)
298: {
301:   /*  Add new vector's fraction of average */
302:   VecAXPY(nm->Xbar,nm->oneOverN,Xmu); 
303:   VecCopy(Xmu,nm->simplex[index]); 
304:   nm->f_values[index] = f;

306:   NelderMeadSort(nm);

308:   /*  Subtract last vector from average */
309:   VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]); 
310:   return(0);
311: }