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) { 168 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 169 } 170 ++tao->niter; 171 shrink = 0; 172 PetscCall(VecCopy(nm->simplex[nm->indices[0]],tao->solution)); 173 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)); 174 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)); 175 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 176 if (tao->reason != TAO_CONTINUE_ITERATING) break; 177 178 /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */ 179 PetscCall(VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]])); 180 PetscCall(TaoComputeObjective(tao,Xmur,&fr)); 181 182 if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) { 183 /* reflect */ 184 nm->nreflect++; 185 PetscCall(PetscInfo(0,"Reflect\n")); 186 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr)); 187 } else if (fr < nm->f_values[nm->indices[0]]) { 188 /* expand */ 189 nm->nexpand++; 190 PetscCall(PetscInfo(0,"Expand\n")); 191 PetscCall(VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]])); 192 PetscCall(TaoComputeObjective(tao,Xmue,&fe)); 193 if (fe < fr) { 194 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe)); 195 } else { 196 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr)); 197 } 198 } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) { 199 /* outside contraction */ 200 nm->noutcontract++; 201 PetscCall(PetscInfo(0,"Outside Contraction\n")); 202 PetscCall(VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]])); 203 204 PetscCall(TaoComputeObjective(tao,Xmuc,&fc)); 205 if (fc <= fr) { 206 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc)); 207 } else shrink=1; 208 } else { 209 /* inside contraction */ 210 nm->nincontract++; 211 PetscCall(PetscInfo(0,"Inside Contraction\n")); 212 PetscCall(VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]])); 213 PetscCall(TaoComputeObjective(tao,Xmuc,&fc)); 214 if (fc < nm->f_values[nm->indices[nm->N]]) { 215 PetscCall(NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc)); 216 } else shrink = 1; 217 } 218 219 if (shrink) { 220 nm->nshrink++; 221 PetscCall(PetscInfo(0,"Shrink\n")); 222 223 for (i=1;i<nm->N+1;i++) { 224 PetscCall(VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]])); 225 PetscCall(TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]])); 226 } 227 PetscCall(VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]])); 228 229 /* Add last vector's fraction of average */ 230 PetscCall(VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]])); 231 PetscCall(NelderMeadSort(nm)); 232 /* Subtract new last vector from average */ 233 PetscCall(VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]])); 234 } 235 } 236 PetscFunctionReturn(0); 237 } 238 239 /* ---------------------------------------------------------- */ 240 /*MC 241 TAONM - Nelder-Mead solver for derivative free, unconstrained minimization 242 243 Options Database Keys: 244 + -tao_nm_lamda - initial step length 245 - -tao_nm_mu - expansion/contraction factor 246 247 Level: beginner 248 M*/ 249 250 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao) 251 { 252 TAO_NelderMead *nm; 253 254 PetscFunctionBegin; 255 PetscCall(PetscNewLog(tao,&nm)); 256 tao->data = (void*)nm; 257 258 tao->ops->setup = TaoSetUp_NM; 259 tao->ops->solve = TaoSolve_NM; 260 tao->ops->view = TaoView_NM; 261 tao->ops->setfromoptions = TaoSetFromOptions_NM; 262 tao->ops->destroy = TaoDestroy_NM; 263 264 /* Override default settings (unless already changed) */ 265 if (!tao->max_it_changed) tao->max_it = 2000; 266 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 267 268 nm->simplex = NULL; 269 nm->lamda = 1; 270 271 nm->mu_ic = -0.5; 272 nm->mu_oc = 0.5; 273 nm->mu_r = 1.0; 274 nm->mu_e = 2.0; 275 276 PetscFunctionReturn(0); 277 } 278