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