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