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 TaoConvergedReason reason; 134 PetscReal *x; 135 PetscInt i; 136 Vec Xmur=nm->Xmur, Xmue=nm->Xmue, Xmuc=nm->Xmuc, Xbar=nm->Xbar; 137 PetscReal fr,fe,fc; 138 PetscInt shrink; 139 PetscInt low,high; 140 141 PetscFunctionBegin; 142 nm->nshrink = 0; 143 nm->nreflect = 0; 144 nm->nincontract = 0; 145 nm->noutcontract = 0; 146 nm->nexpand = 0; 147 148 if (tao->XL || tao->XU || tao->ops->computebounds) { 149 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by NelderMead algorithm\n");CHKERRQ(ierr); 150 } 151 152 ierr = VecCopy(tao->solution,nm->simplex[0]);CHKERRQ(ierr); 153 ierr = TaoComputeObjective(tao,nm->simplex[0],&nm->f_values[0]);CHKERRQ(ierr); 154 nm->indices[0]=0; 155 for (i=1;i<nm->N+1;i++){ 156 ierr = VecCopy(tao->solution,nm->simplex[i]);CHKERRQ(ierr); 157 ierr = VecGetOwnershipRange(nm->simplex[i],&low,&high);CHKERRQ(ierr); 158 if (i-1 >= low && i-1 < high) { 159 ierr = VecGetArray(nm->simplex[i],&x);CHKERRQ(ierr); 160 x[i-1-low] += nm->lamda; 161 ierr = VecRestoreArray(nm->simplex[i],&x);CHKERRQ(ierr); 162 } 163 164 ierr = TaoComputeObjective(tao,nm->simplex[i],&nm->f_values[i]);CHKERRQ(ierr); 165 nm->indices[i] = i; 166 } 167 168 /* Xbar = (Sum of all simplex vectors - worst vector)/N */ 169 ierr = NelderMeadSort(nm);CHKERRQ(ierr); 170 ierr = VecSet(Xbar,0.0);CHKERRQ(ierr); 171 for (i=0;i<nm->N;i++) { 172 ierr = VecAXPY(Xbar,1.0,nm->simplex[nm->indices[i]]);CHKERRQ(ierr); 173 } 174 ierr = VecScale(Xbar,nm->oneOverN);CHKERRQ(ierr); 175 reason = TAO_CONTINUE_ITERATING; 176 while (1) { 177 shrink = 0; 178 ierr = VecCopy(nm->simplex[nm->indices[0]],tao->solution);CHKERRQ(ierr); 179 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); 180 if (reason != TAO_CONTINUE_ITERATING) break; 181 182 /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */ 183 ierr = VecAXPBYPCZ(Xmur,1+nm->mu_r,-nm->mu_r,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 184 ierr = TaoComputeObjective(tao,Xmur,&fr);CHKERRQ(ierr); 185 186 if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N-1]]) { 187 /* reflect */ 188 nm->nreflect++; 189 ierr = PetscInfo(0,"Reflect\n");CHKERRQ(ierr); 190 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr); 191 } else if (fr < nm->f_values[nm->indices[0]]) { 192 /* expand */ 193 nm->nexpand++; 194 ierr = PetscInfo(0,"Expand\n");CHKERRQ(ierr); 195 ierr = VecAXPBYPCZ(Xmue,1+nm->mu_e,-nm->mu_e,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 196 ierr = TaoComputeObjective(tao,Xmue,&fe);CHKERRQ(ierr); 197 if (fe < fr) { 198 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmue,fe);CHKERRQ(ierr); 199 } else { 200 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmur,fr);CHKERRQ(ierr); 201 } 202 } else if (nm->f_values[nm->indices[nm->N-1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) { 203 /* outside contraction */ 204 nm->noutcontract++; 205 ierr = PetscInfo(0,"Outside Contraction\n");CHKERRQ(ierr); 206 ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_oc,-nm->mu_oc,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 207 208 ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr); 209 if (fc <= fr) { 210 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr); 211 } else shrink=1; 212 } else { 213 /* inside contraction */ 214 nm->nincontract++; 215 ierr = PetscInfo(0,"Inside Contraction\n");CHKERRQ(ierr); 216 ierr = VecAXPBYPCZ(Xmuc,1+nm->mu_ic,-nm->mu_ic,0,Xbar,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 217 ierr = TaoComputeObjective(tao,Xmuc,&fc);CHKERRQ(ierr); 218 if (fc < nm->f_values[nm->indices[nm->N]]) { 219 ierr = NelderMeadReplace(nm,nm->indices[nm->N],Xmuc,fc);CHKERRQ(ierr); 220 } else shrink = 1; 221 } 222 223 if (shrink) { 224 nm->nshrink++; 225 ierr = PetscInfo(0,"Shrink\n");CHKERRQ(ierr); 226 227 for (i=1;i<nm->N+1;i++) { 228 ierr = VecAXPBY(nm->simplex[nm->indices[i]],1.5,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr); 229 ierr = TaoComputeObjective(tao,nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]]);CHKERRQ(ierr); 230 } 231 ierr = VecAXPBY(Xbar,1.5*nm->oneOverN,-0.5,nm->simplex[nm->indices[0]]);CHKERRQ(ierr); 232 233 /* Add last vector's fraction of average */ 234 ierr = VecAXPY(Xbar,nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 235 ierr = NelderMeadSort(nm);CHKERRQ(ierr); 236 /* Subtract new last vector from average */ 237 ierr = VecAXPY(Xbar,-nm->oneOverN,nm->simplex[nm->indices[nm->N]]);CHKERRQ(ierr); 238 } 239 } 240 PetscFunctionReturn(0); 241 } 242 243 /* ---------------------------------------------------------- */ 244 /*MC 245 TAONM - Nelder-Mead solver for derivative free, unconstrained minimization 246 247 Options Database Keys: 248 + -tao_nm_lamda - initial step length 249 . -tao_nm_mu - expansion/contraction factor 250 251 Level: beginner 252 M*/ 253 254 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao) 255 { 256 TAO_NelderMead *nm; 257 PetscErrorCode ierr; 258 259 PetscFunctionBegin; 260 ierr = PetscNewLog(tao,&nm);CHKERRQ(ierr); 261 tao->data = (void*)nm; 262 263 tao->ops->setup = TaoSetUp_NM; 264 tao->ops->solve = TaoSolve_NM; 265 tao->ops->view = TaoView_NM; 266 tao->ops->setfromoptions = TaoSetFromOptions_NM; 267 tao->ops->destroy = TaoDestroy_NM; 268 269 /* Override default settings (unless already changed) */ 270 if (!tao->max_it_changed) tao->max_it = 2000; 271 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 272 273 nm->simplex = 0; 274 nm->lamda = 1; 275 276 nm->mu_ic = -0.5; 277 nm->mu_oc = 0.5; 278 nm->mu_r = 1.0; 279 nm->mu_e = 2.0; 280 281 PetscFunctionReturn(0); 282 } 283 284