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: }