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=NULL; 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 PetscFunctionReturn(0); 87 } 88 89 /*------------------------------------------------------------*/ 90 static PetscErrorCode TaoSetFromOptions_NM(PetscOptionItems *PetscOptionsObject,Tao tao) 91 { 92 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 93 PetscErrorCode ierr; 94 95 PetscFunctionBegin; 96 ierr = PetscOptionsHead(PetscOptionsObject,"Nelder-Mead options");CHKERRQ(ierr); 97 ierr = PetscOptionsReal("-tao_nm_lamda","initial step length","",nm->lamda,&nm->lamda,NULL); CHKERRQ(ierr); 98 ierr = PetscOptionsReal("-tao_nm_mu","mu","",nm->mu_oc,&nm->mu_oc,NULL);CHKERRQ(ierr); 99 nm->mu_ic = -nm->mu_oc; 100 nm->mu_r = nm->mu_oc*2.0; 101 nm->mu_e = nm->mu_oc*4.0; 102 ierr = PetscOptionsTail();CHKERRQ(ierr); 103 PetscFunctionReturn(0); 104 } 105 106 /*------------------------------------------------------------*/ 107 static PetscErrorCode TaoView_NM(Tao tao,PetscViewer viewer) 108 { 109 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 110 PetscBool isascii; 111 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 115 if (isascii) { 116 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 117 ierr = PetscViewerASCIIPrintf(viewer,"expansions: %D\n",nm->nexpand);CHKERRQ(ierr); 118 ierr = PetscViewerASCIIPrintf(viewer,"reflections: %D\n",nm->nreflect);CHKERRQ(ierr); 119 ierr = PetscViewerASCIIPrintf(viewer,"inside contractions: %D\n",nm->nincontract);CHKERRQ(ierr); 120 ierr = PetscViewerASCIIPrintf(viewer,"outside contractionss: %D\n",nm->noutcontract);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPrintf(viewer,"Shrink steps: %D\n",nm->nshrink);CHKERRQ(ierr); 122 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 123 } 124 PetscFunctionReturn(0); 125 } 126 127 /*------------------------------------------------------------*/ 128 static PetscErrorCode TaoSolve_NM(Tao tao) 129 { 130 PetscErrorCode ierr; 131 TAO_NelderMead *nm = (TAO_NelderMead*)tao->data; 132 PetscReal *x; 133 PetscInt i; 134 Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar; 135 PetscReal fr,fe,fc; 136 PetscInt shrink; 137 PetscInt low,high; 138 139 PetscFunctionBegin; 140 nm->nshrink = 0; 141 nm->nreflect = 0; 142 nm->nincontract = 0; 143 nm->noutcontract = 0; 144 nm->nexpand = 0; 145 146 if (tao->XL || tao->XU || tao->ops->computebounds) { 147 ierr = PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");CHKERRQ(ierr); 148 } 149 150 ierr = VecCopy(tao->solution,nm->simplex[0]);CHKERRQ(ierr); 151 ierr = TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);CHKERRQ(ierr); 152 nm->indices[0]=0; 153 for (i=1;i<nm->N+1;i++){ 154 ierr = VecCopy(tao->solution,nm->simplex[i]);CHKERRQ(ierr); 155 ierr = VecGetOwnershipRange(nm->simplex[i],&low,&high);CHKERRQ(ierr); 156 if (i-1 >= low && i-1 < high) { 157 ierr = VecGetArray(nm->simplex[i],&x);CHKERRQ(ierr); 158 x[i-1-low] += nm->lamda; 159 ierr = VecRestoreArray(nm->simplex[i],&x);CHKERRQ(ierr); 160 } 161 162 ierr = TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);CHKERRQ(ierr); 163 nm->indices[i] = i; 164 } 165 166 /* Xbar = (Sum of all simplex vectors - worst vector)/N */ 167 ierr = NelderMeadSort(nm);CHKERRQ(ierr); 168 ierr = VecSet(Xbar,0.0);CHKERRQ(ierr); 169 for (i=0;i<nm->N;i++) { 170 ierr = VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);CHKERRQ(ierr); 171 } 172 ierr = VecScale(Xbar,nm->oneOverN);CHKERRQ(ierr); 173 tao->reason = TAO_CONTINUE_ITERATING; 174 while (1) { 175 /* Call general purpose update function */ 176 if (tao->ops->update) { 177 ierr = (*tao->ops->update)(tao, tao->niter, tao->user_update);CHKERRQ(ierr); 178 } 179 ++tao->niter; 180 shrink = 0; 181 ierr = VecCopy(nm->simplex[nm->indices[0]],tao->solution);CHKERRQ(ierr); 182 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); 183 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); 184 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 185 if (tao->reason != TAO_CONTINUE_ITERATING) break; 186 187 /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */ 188 ierr = VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 189 ierr = TaoComputeObjective(tao,Xmur,&fr);CHKERRQ(ierr); 190 191 if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) { 192 /* reflect */ 193 nm->nreflect++; 194 ierr = PetscInfo(0,"Reflect\n");CHKERRQ(ierr); 195 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr); 196 } else if (fr < nm->f_values[nm->indices[0]]) { 197 /* expand */ 198 nm->nexpand++; 199 ierr = PetscInfo(0,"Expand\n");CHKERRQ(ierr); 200 ierr = VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 201 ierr = TaoComputeObjective(tao,Xmue,&fe);CHKERRQ(ierr); 202 if (fe < fr) { 203 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);CHKERRQ(ierr); 204 } else { 205 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr); 206 } 207 } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) { 208 /* outside contraction */ 209 nm->noutcontract++; 210 ierr = PetscInfo(0,"Outside Contraction\n");CHKERRQ(ierr); 211 ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 212 213 ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr); 214 if (fc <= fr) { 215 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr); 216 } else shrink=1; 217 } else { 218 /* inside contraction */ 219 nm->nincontract++; 220 ierr = PetscInfo(0,"Inside Contraction\n");CHKERRQ(ierr); 221 ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 222 ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr); 223 if (fc < nm->f_values[nm->indices[nm->N]]) { 224 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr); 225 } else shrink = 1; 226 } 227 228 if (shrink) { 229 nm->nshrink++; 230 ierr = PetscInfo(0,"Shrink\n");CHKERRQ(ierr); 231 232 for (i=1;i<nm->N+1;i++) { 233 ierr = VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr); 234 ierr = TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);CHKERRQ(ierr); 235 } 236 ierr = VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr); 237 238 /* Add last vector's fraction of average */ 239 ierr = VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 240 ierr = NelderMeadSort(nm);CHKERRQ(ierr); 241 /* Subtract new last vector from average */ 242 ierr = VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 243 } 244 } 245 PetscFunctionReturn(0); 246 } 247 248 /* ---------------------------------------------------------- */ 249 /*MC 250 TAONM - Nelder-Mead solver for derivative free, unconstrained minimization 251 252 Options Database Keys: 253 + -tao_nm_lamda - initial step length 254 - -tao_nm_mu - expansion/contraction factor 255 256 Level: beginner 257 M*/ 258 259 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao) 260 { 261 TAO_NelderMead *nm; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = PetscNewLog(tao,&nm);CHKERRQ(ierr); 266 tao->data = (void*)nm; 267 268 tao->ops->setup = TaoSetUp_NM; 269 tao->ops->solve = TaoSolve_NM; 270 tao->ops->view = TaoView_NM; 271 tao->ops->setfromoptions = TaoSetFromOptions_NM; 272 tao->ops->destroy = TaoDestroy_NM; 273 274 /* Override default settings (unless already changed) */ 275 if (!tao->max_it_changed) tao->max_it = 2000; 276 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 277 278 nm->simplex = NULL; 279 nm->lamda = 1; 280 281 nm->mu_ic = -0.5; 282 nm->mu_oc = 0.5; 283 nm->mu_r = 1.0; 284 nm->mu_e = 2.0; 285 286 PetscFunctionReturn(0); 287 } 288 289