Actual source code: neldermead.c

  1: #include "neldermead.h"
  2: #ifdef TAO_USE_PETSC
  3: #include "src/petsctao/vector/taovec_petsc.h"
  4: #endif

  6: int NelderMeadSort(TAO_NelderMead *nm);
  7: int NelderMeadReplace(TAO_NelderMead *nm, int index, TaoVec *Xmu, double f);
  8: /* ---------------------------------------------------------- */
 11: int TaoSetUp_NelderMead(TAO_SOLVER tao, void *solver)
 12: {
 13:   TAO_NelderMead *nm = (TAO_NelderMead *)solver;
 14:   int info;
 15:   TaoInt size;
 16:   TaoVec *X;

 18:   TaoFunctionBegin;

 20:   
 21:   info = TaoGetSolution(tao,&X); CHKERRQ(info);
 22:   info = X->GetDimension(&size); CHKERRQ(info);
 23:   nm->N = size;
 24:   nm->oneOverN = 1.0/size;
 25:   info = X->CloneVecs(nm->N+1,&nm->simplex);
 26:   nm->f_values = new double[nm->N+1];
 27:   nm->indices =new int[nm->N+1];
 28:   info = X->Clone(&nm->Xbar); CHKERRQ(info);
 29:   info = X->Clone(&nm->Xmur); CHKERRQ(info);
 30:   info = X->Clone(&nm->Xmue); CHKERRQ(info);
 31:   info = X->Clone(&nm->Xmuc); CHKERRQ(info);
 32:   info = TaoSetLagrangianGradientVector(tao,0);CHKERRQ(info);
 33:   info = TaoSetStepDirectionVector(tao,0);CHKERRQ(info);

 35:   TaoFunctionReturn(0);
 36: }

 38: /* ---------------------------------------------------------- */
 41: int TaoSetDown_NelderMead(TAO_SOLVER tao, void *solver)
 42: {
 43:   TAO_NelderMead *nm = (TAO_NelderMead*)solver;
 44:   int info;
 45:   int i;
 46:   TaoFunctionBegin;
 47:   for (i=0;i<nm->N+1;i++) {
 48:     TaoVecDestroy( nm->simplex[i]);
 49:   }
 50:   delete [] nm->simplex;

 52:   info = TaoVecDestroy(nm->Xmuc); CHKERRQ(info);
 53:   info = TaoVecDestroy(nm->Xmue); CHKERRQ(info);
 54:   info = TaoVecDestroy(nm->Xmur); CHKERRQ(info);
 55:   info = TaoVecDestroy(nm->Xbar); CHKERRQ(info);
 56:   
 57:   delete [] nm->indices;
 58:   delete [] nm->f_values;
 59:   TaoFunctionReturn(0);
 60: }

 62: /*------------------------------------------------------------*/
 65: int TaoSetOptions_NelderMead(TAO_SOLVER tao, void *solver)
 66: {
 67:   
 68:   TAO_NelderMead *nm = (TAO_NelderMead*)solver;
 69:   TaoTruth flg;
 70:   int info;
 71:   
 72:   TaoFunctionBegin;
 73:   info = TaoOptionsHead("Nelder-Mead options"); CHKERRQ(info);

 75:   info = TaoOptionDouble("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,&flg);

 77:   CHKERRQ(info);

 79:   info = TaoOptionDouble("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,&flg); CHKERRQ(info);
 80:   nm->mu_ic = -nm->mu_oc;
 81:   nm->mu_r = nm->mu_oc*2.0;
 82:   nm->mu_e = nm->mu_oc*4.0;

 84:   info = TaoOptionsTail(); CHKERRQ(info);
 85:   TaoFunctionReturn(0);
 86: }

 88: /*------------------------------------------------------------*/
 91: int TaoView_NelderMead(TAO_SOLVER tao, void *solver)
 92: {
 93:   //int info;

 95:   TaoFunctionBegin;
 96:   TaoFunctionReturn(0);
 97: }

 99: /*------------------------------------------------------------*/
102: int TaoSolve_NelderMead(TAO_SOLVER tao, void *solver)
103: {
104:   int info;
105:   TAO_NelderMead *nm = (TAO_NelderMead*)solver;
106:   TaoTerminateReason reason;
107:   TaoVec *xx;
108:   double *x;
109:   double step=0.0;
110:   int iter=0,i;
111:   TaoVec *Xmur=nm->Xmur, *Xmue=nm->Xmue, *Xmuc=nm->Xmuc, *Xbar=nm->Xbar;
112:   double fr,fe,fc;
113:   int shrink;
114:   
115:   
116:   TaoFunctionBegin;
117:   info = TaoGetSolution(tao,&xx); CHKERRQ(info);
118:   info = nm->simplex[0]->CopyFrom(xx); CHKERRQ(info);
119:   info = TaoComputeMeritFunction(tao,nm->simplex[0],&nm->f_values[0]); CHKERRQ(info);
120:   nm->indices[0] = 0;
121:   for (i=1;i<nm->N+1;i++){
122:     info = nm->simplex[i]->CopyFrom(xx); CHKERRQ(info);
123: #ifdef TAO_USE_PETSC
124:     PetscInt low,high;
125:     TaoVecPetsc *tvxx  = dynamic_cast<TaoVecPetsc*>(nm->simplex[i]);
126:     if (!tvxx) {SETERRQ(1,"Could not cast TaoVec to TaoVecPetsc");}
127:     Vec px = tvxx->GetVec();
128:     info = VecGetOwnershipRange(px,&low,&high); CHKERRQ(info);
129:     if (i-1 >= low && i-1 < high) {
130:       info = VecGetArray(px,&x); CHKERRQ(info);
131:       x[i-1-low] += nm->lamda;
132:       info = VecRestoreArray(px,&x); CHKERRQ(info);
133:     }


136: #else
137:     TaoInt dim;
138:     if (i>0) {
139:       
140:       info = nm->simplex[i]->GetArray(&x,&dim); CHKERRQ(info);
141:       x[i-1] += nm->lamda;
142:       info = nm->simplex[i]->RestoreArray(&x,&dim); CHKERRQ(info);
143:     }
144: #endif
145:     info = TaoComputeMeritFunction(tao,nm->simplex[i],&nm->f_values[i]); CHKERRQ(info);
146:     nm->indices[i] = i;
147:   }

149:   // Xbar  = (Sum of all simplex vectors - worst vector)/N
150:   info = NelderMeadSort(nm); CHKERRQ(info);
151:   info = Xbar->SetToZero(); CHKERRQ(info);
152:   for (i=0;i<nm->N;i++) {
153:     info = Xbar->Axpy(1.0,nm->simplex[nm->indices[i]]);
154:   }
155:   info = Xbar->Scale(nm->oneOverN);

157:   reason = TAO_CONTINUE_ITERATING;
158:   while (1) {
159:     shrink = 0;

161:     info = xx->CopyFrom(nm->simplex[nm->indices[0]]); CHKERRQ(info);
162:     info = TaoMonitor(tao,iter++,nm->f_values[nm->indices[0]],nm->f_values[nm->indices[nm->N]]-nm->f_values[nm->indices[0]],0,step,&reason); CHKERRQ(info);
163:     if (reason != TAO_CONTINUE_ITERATING) break;

165:     
166:     
167:     //x(mu) = (1 + mu)Xbar - mu*X_N+1
168:     info = Xmur->Waxpby(1+nm->mu_r,Xbar,-nm->mu_r,
169:                        nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
170:     info = TaoComputeMeritFunction(tao,Xmur,&fr); CHKERRQ(info);


173:     if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) {
174:       // reflect
175:       info = PetscInfo(0,"Reflect\n"); CHKERRQ(info);
176:       info = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr); CHKERRQ(info);
177:     }

179:     else if (fr < nm->f_values[nm->indices[0]]) {
180:       // expand
181:       info = PetscInfo(0,"Expand\n"); CHKERRQ(info);
182:       info = Xmue->Waxpby(1+nm->mu_e,Xbar,-nm->mu_e,
183:                           nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
184:       info = TaoComputeMeritFunction(tao,Xmue,&fe); CHKERRQ(info);
185:       if (fe < fr) {
186:         info = NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe); CHKERRQ(info);
187:       } else {
188:         info = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr); CHKERRQ(info);
189:       }

191:     } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) {
192:       //outside contraction
193:       info = PetscInfo(0,"Outside Contraction\n"); CHKERRQ(info);
194:       info = Xmuc->Waxpby(1+nm->mu_oc,Xbar,-nm->mu_oc,
195:                           nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
196:       info = TaoComputeMeritFunction(tao,Xmuc,&fc); CHKERRQ(info);
197:       if (fc <= fr) {
198:         info = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc); CHKERRQ(info);
199:       }        else 
200:         shrink=1;
201:     } else {
202:       //inside contraction
203:       info = PetscInfo(0,"Inside Contraction\n"); CHKERRQ(info);
204:       info = Xmuc->Waxpby(1+nm->mu_ic,Xbar,-nm->mu_ic,
205:                           nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
206:       info = TaoComputeMeritFunction(tao,Xmuc,&fc); CHKERRQ(info);
207:       if (fc < nm->f_values[nm->indices[nm->N]]) {
208:         info = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc); CHKERRQ(info);
209:       } else
210:         shrink = 1;
211:     }

213:     if (shrink) {
214:       info = PetscInfo(0,"Shrink\n"); CHKERRQ(info);
215:       printf("shrink!\n");
216:       
217:       for (i=1;i<nm->N+1;i++) {
218:         info = nm->simplex[nm->indices[i]]->Axpby(1.5,nm->simplex[nm->indices[0]],-.5); CHKERRQ(info);
219:         info = TaoComputeMeritFunction(tao,nm->simplex[nm->indices[i]],
220:                                  &nm->f_values[nm->indices[i]]); CHKERRQ(info);
221:       }

223:       info = Xbar->Axpby(1.5*nm->oneOverN,nm->simplex[nm->indices[0]],-0.5); CHKERRQ(info);
224:       // Add last vector's fraction of average
225:       info = nm->Xbar->Axpy(nm->oneOverN,
226:                             nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
227:       info = NelderMeadSort(nm);
228:       // Subtract new last vector from average
229:       info = nm->Xbar->Axpy(-nm->oneOverN,
230:                             nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
231:     }
232:     
233:     
234:   }
235:   TaoFunctionReturn(0);
236: }

238: /* ---------------------------------------------------------- */
242: int TaoCreate_NelderMead(TAO_SOLVER tao)
243: {
244:   TAO_NelderMead *nm;
245:   int info;

247:   TaoFunctionBegin;
248:   info = TaoNew(TAO_NelderMead,&nm); CHKERRQ(info);
249:   tao->data = (void*)nm;

251:   info = PetscLogObjectMemory(tao,sizeof(TAO_NelderMead)); CHKERRQ(info);

253:   info = TaoSetTaoSolveRoutine(tao,TaoSolve_NelderMead,(void*)nm); CHKERRQ(info);
254:   info = TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_NelderMead, 
255:                                     TaoSetDown_NelderMead); CHKERRQ(info);
256:   info = TaoSetTaoOptionsRoutine(tao,TaoSetOptions_NelderMead); CHKERRQ(info);
257:   info = TaoSetTaoViewRoutine(tao, TaoView_NelderMead); CHKERRQ(info);
258:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
259:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
260:   info = TaoSetTolerances(tao,1e-8,1e-8,0,0); CHKERRQ(info);


263:   nm->simplex = 0;
264:   nm->lamda = 1;

266:   nm->mu_ic = -0.5;
267:   nm->mu_oc = 0.5;
268:   nm->mu_r = 1.0;
269:   nm->mu_e = 2.0;

271:   TaoFunctionReturn(0);
272: }
274:     

276: /*------------------------------------------------------------*/
279: int NelderMeadSort(TAO_NelderMead *nm) {
280:   double *values = nm->f_values;
281:   int *indices = nm->indices;
282:   int dim = nm->N+1;

284:   int i,j,index;
285:   double val;
286:   TaoFunctionBegin;  
287:   for (i=1;i<dim;i++) {
288:     index = indices[i];
289:     val = values[index];
290:     for (j=i-1; j>=0 && values[indices[j]] > val; j--) {
291:       indices[j+1] = indices[j];
292:     }
293:     indices[j+1] = index;
294:   }  
295:   TaoFunctionReturn(0);
296: }


299: /*------------------------------------------------------------*/
302: int NelderMeadReplace(TAO_NelderMead *nm, int index, TaoVec *Xmu, double f)
303: {
304:   int info;
305:   TaoFunctionBegin;
306:   // Add new vector's fraction of average
307:   info = nm->Xbar->Axpy(nm->oneOverN,Xmu); CHKERRQ(info);
308:   info = nm->simplex[index]->CopyFrom(Xmu); CHKERRQ(info);
309:   nm->f_values[index] = f;

311:   info = NelderMeadSort(nm);
312:   // Subtract last vector from average
313:   info = nm->Xbar->Axpy(-nm->oneOverN,
314:                         nm->simplex[nm->indices[nm->N]]); CHKERRQ(info);
315:   TaoFunctionReturn(0);
316: }