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