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