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