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