1 #include <../src/tao/unconstrained/impls/neldermead/neldermead.h> 2 #include <petscvec.h> 3 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; 12 13 PetscFunctionBegin; 14 for (i=1;i<dim;i++) { 15 index = indices[i]; 16 val = values[index]; 17 for (j=i-1; j>=0 && values[indices[j]] > val; j--) { 18 indices[j+1] = indices[j]; 19 } 20 indices[j+1] = index; 21 } 22 PetscFunctionReturn(0); 23 } 24 25 /*------------------------------------------------------------*/ 26 static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f) 27 { 28 PetscFunctionBegin; 29 /* Add new vector's fraction of average */ 30 PetscCall(VecAXPY(nm->Xbar,nm->oneOverN,Xmu)); 31 PetscCall(VecCopy(Xmu,nm->simplex[index])); 32 nm->f_values[index] = f; 33 34 PetscCall(NelderMeadSort(nm)); 35 36 /* Subtract last vector from average */ 37 PetscCall(VecAXPY(nm->Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]])); 38 PetscFunctionReturn(0); 39 } 40 41 /* ---------------------------------------------------------- */ 42 static PetscErrorCode TaoSetUp_NM(Tao tao) 43 { 44 TAO_NelderMead *nm = (TAO_NelderMead *)tao->data; 45 PetscInt n; 46 47 PetscFunctionBegin; 48 PetscCall(VecGetSize(tao->solution,&n)); 49 nm->N = n; 50 nm->oneOverN = 1.0/n; 51 PetscCall(VecDuplicateVecs(tao->solution,nm->N+1,&nm->simplex)); 52 PetscCall(PetscMalloc1(nm->N+1,&nm->f_values)); 53 PetscCall(PetscMalloc1(nm->N+1,&nm->indices)); 54 PetscCall(VecDuplicate(tao->solution,&nm->Xbar)); 55 PetscCall(VecDuplicate(tao->solution,&nm->Xmur)); 56 PetscCall(VecDuplicate(tao->solution,&nm->Xmue)); 57 PetscCall(VecDuplicate(tao->solution,&nm->Xmuc)); 58 59 tao->gradient=NULL; 60 tao->step=0; 61 PetscFunctionReturn(0); 62 } 63 64 /* ---------------------------------------------------------- */ 65 static PetscErrorCode TaoDestroy_NM(Tao tao) 66 { 67 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 68 69 PetscFunctionBegin; 70 if (tao->setupcalled) { 71 PetscCall(VecDestroyVecs(nm->N+1,&nm->simplex)); 72 PetscCall(VecDestroy(&nm->Xmuc)); 73 PetscCall(VecDestroy(&nm->Xmue)); 74 PetscCall(VecDestroy(&nm->Xmur)); 75 PetscCall(VecDestroy(&nm->Xbar)); 76 } 77 PetscCall(PetscFree(nm->indices)); 78 PetscCall(PetscFree(nm->f_values)); 79 PetscCall(PetscFree(tao->data)); 80 PetscFunctionReturn(0); 81 } 82 83 /*------------------------------------------------------------*/ 84 static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao) 85 { 86 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 87 88 PetscFunctionBegin; 89 PetscOptionsHeadBegin(PetscOptionsObject,"Nelder-Mead options"); 90 PetscCall(PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL)); 91 PetscCall(PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL)); 92 nm->mu_ic = -nm->mu_oc; 93 nm->mu_r = nm->mu_oc*2.0; 94 nm->mu_e = nm->mu_oc*4.0; 95 PetscOptionsHeadEnd(); 96 PetscFunctionReturn(0); 97 } 98 99 /*------------------------------------------------------------*/ 100 static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer) 101 { 102 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 103 PetscBool isascii; 104 105 PetscFunctionBegin; 106 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 107 if (isascii) { 108 PetscCall(PetscViewerASCIIPushTab(viewer)); 109 PetscCall(PetscViewerASCIIPrintf(viewer,"expansions: %" PetscInt_FMT "\n",nm->nexpand)); 110 PetscCall(PetscViewerASCIIPrintf(viewer,"reflections: %" PetscInt_FMT "\n",nm->nreflect)); 111 PetscCall(PetscViewerASCIIPrintf(viewer,"inside contractions: %" PetscInt_FMT "\n",nm->nincontract)); 112 PetscCall(PetscViewerASCIIPrintf(viewer,"outside contractionss: %" PetscInt_FMT "\n",nm->noutcontract)); 113 PetscCall(PetscViewerASCIIPrintf(viewer,"Shrink steps: %" PetscInt_FMT "\n",nm->nshrink)); 114 PetscCall(PetscViewerASCIIPopTab(viewer)); 115 } 116 PetscFunctionReturn(0); 117 } 118 119 /*------------------------------------------------------------*/ 120 static PetscErrorCode TaoSolve_NM(Tao tao) 121 { 122 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 123 PetscReal *x; 124 PetscInt i; 125 Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar; 126 PetscReal fr,fe,fc; 127 PetscInt shrink; 128 PetscInt low,high; 129 130 PetscFunctionBegin; 131 nm->nshrink = 0; 132 nm->nreflect = 0; 133 nm->nincontract = 0; 134 nm->noutcontract = 0; 135 nm->nexpand = 0; 136 137 if (tao->XL || tao->XU || tao->ops->computebounds) { 138 PetscCall(PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n")); 139 } 140 141 PetscCall(VecCopy(tao->solution,nm->simplex[0])); 142 PetscCall(TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0])); 143 nm->indices[0]=0; 144 for (i=1;i<nm->N+1;i++) { 145 PetscCall(VecCopy(tao->solution,nm->simplex[i])); 146 PetscCall(VecGetOwnershipRange(nm->simplex[i],&low,&high)); 147 if (i-1 >= low && i-1 < high) { 148 PetscCall(VecGetArray(nm->simplex[i],&x)); 149 x[i-1-low] += nm->lamda; 150 PetscCall(VecRestoreArray(nm->simplex[i],&x)); 151 } 152 153 PetscCall(TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i])); 154 nm->indices[i] = i; 155 } 156 157 /* Xbar = (Sum of all simplex vectors - worst vector)/N */ 158 PetscCall(NelderMeadSort(nm)); 159 PetscCall(VecSet(Xbar,0.0)); 160 for (i=0;i<nm->N;i++) { 161 PetscCall(VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]])); 162 } 163 PetscCall(VecScale(Xbar,nm->oneOverN)); 164 tao->reason = TAO_CONTINUE_ITERATING; 165 while (1) { 166 /* Call general purpose update function */ 167 if (tao->ops->update) PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 168 ++tao->niter; 169 shrink = 0; 170 PetscCall(VecCopy(nm->simplex[nm->indices[0]],tao->solution)); 171 PetscCall(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)); 172 PetscCall(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)); 173 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 174 if (tao->reason != TAO_CONTINUE_ITERATING) break; 175 176 /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */ 177 PetscCall(VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]])); 178 PetscCall(TaoComputeObjective(tao,Xmur,&fr)); 179 180 if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) { 181 /* reflect */ 182 nm->nreflect++; 183 PetscCall(PetscInfo(0,"Reflect\n")); 184 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr)); 185 } else if (fr < nm->f_values[nm->indices[0]]) { 186 /* expand */ 187 nm->nexpand++; 188 PetscCall(PetscInfo(0,"Expand\n")); 189 PetscCall(VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]])); 190 PetscCall(TaoComputeObjective(tao,Xmue,&fe)); 191 if (fe < fr) { 192 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe)); 193 } else { 194 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr)); 195 } 196 } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) { 197 /* outside contraction */ 198 nm->noutcontract++; 199 PetscCall(PetscInfo(0,"Outside Contraction\n")); 200 PetscCall(VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]])); 201 202 PetscCall(TaoComputeObjective(tao,Xmuc,&fc)); 203 if (fc <= fr) { 204 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc)); 205 } else shrink=1; 206 } else { 207 /* inside contraction */ 208 nm->nincontract++; 209 PetscCall(PetscInfo(0,"Inside Contraction\n")); 210 PetscCall(VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]])); 211 PetscCall(TaoComputeObjective(tao,Xmuc,&fc)); 212 if (fc < nm->f_values[nm->indices[nm->N]]) { 213 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc)); 214 } else shrink = 1; 215 } 216 217 if (shrink) { 218 nm->nshrink++; 219 PetscCall(PetscInfo(0,"Shrink\n")); 220 221 for (i=1;i<nm->N+1;i++) { 222 PetscCall(VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]])); 223 PetscCall(TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]])); 224 } 225 PetscCall(VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]])); 226 227 /* Add last vector's fraction of average */ 228 PetscCall(VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]])); 229 PetscCall(NelderMeadSort(nm)); 230 /* Subtract new last vector from average */ 231 PetscCall(VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]])); 232 } 233 } 234 PetscFunctionReturn(0); 235 } 236 237 /* ---------------------------------------------------------- */ 238 /*MC 239 TAONM - Nelder-Mead solver for derivative free, unconstrained minimization 240 241 Options Database Keys: 242 + -tao_nm_lamda - initial step length 243 - -tao_nm_mu - expansion/contraction factor 244 245 Level: beginner 246 M*/ 247 248 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao) 249 { 250 TAO_NelderMead *nm; 251 252 PetscFunctionBegin; 253 PetscCall(PetscNewLog(tao,&nm)); 254 tao->data = (void*)nm; 255 256 tao->ops->setup = TaoSetUp_NM; 257 tao->ops->solve = TaoSolve_NM; 258 tao->ops->view = TaoView_NM; 259 tao->ops->setfromoptions = TaoSetFromOptions_NM; 260 tao->ops->destroy = TaoDestroy_NM; 261 262 /* Override default settings (unless already changed) */ 263 if (!tao->max_it_changed) tao->max_it = 2000; 264 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 265 266 nm->simplex = NULL; 267 nm->lamda = 1; 268 269 nm->mu_ic = -0.5; 270 nm->mu_oc = 0.5; 271 nm->mu_r = 1.0; 272 nm->mu_e = 2.0; 273 274 PetscFunctionReturn(0); 275 } 276