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