1 #include <../src/tao/unconstrained/impls/neldermead/neldermead.h> 2 #include <petscvec.h> 3 4 static PetscErrorCode NelderMeadSort(TAO_NelderMead *nm) 5 { 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(PETSC_SUCCESS); 20 } 21 22 static PetscErrorCode NelderMeadReplace(TAO_NelderMead *nm, PetscInt index, Vec Xmu, PetscReal f) 23 { 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(PETSC_SUCCESS); 35 } 36 37 static PetscErrorCode TaoSetUp_NM(Tao tao) 38 { 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(PETSC_SUCCESS); 57 } 58 59 static PetscErrorCode TaoDestroy_NM(Tao tao) 60 { 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(PETSC_SUCCESS); 75 } 76 77 static PetscErrorCode TaoSetFromOptions_NM(Tao tao, PetscOptionItems PetscOptionsObject) 78 { 79 TAO_NelderMead *nm = (TAO_NelderMead *)tao->data; 80 81 PetscFunctionBegin; 82 PetscOptionsHeadBegin(PetscOptionsObject, "Nelder-Mead options"); 83 PetscCall(PetscOptionsDeprecated("-tao_nm_lamda", "-tao_nm_lambda", "3.18.4", NULL)); 84 PetscCall(PetscOptionsReal("-tao_nm_lambda", "initial step length", "", nm->lambda, &nm->lambda, NULL)); 85 PetscCall(PetscOptionsReal("-tao_nm_mu", "mu", "", nm->mu_oc, &nm->mu_oc, NULL)); 86 nm->mu_ic = -nm->mu_oc; 87 nm->mu_r = nm->mu_oc * 2.0; 88 nm->mu_e = nm->mu_oc * 4.0; 89 PetscOptionsHeadEnd(); 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 static PetscErrorCode TaoView_NM(Tao tao, PetscViewer viewer) 94 { 95 TAO_NelderMead *nm = (TAO_NelderMead *)tao->data; 96 PetscBool isascii; 97 98 PetscFunctionBegin; 99 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 100 if (isascii) { 101 PetscCall(PetscViewerASCIIPushTab(viewer)); 102 PetscCall(PetscViewerASCIIPrintf(viewer, "expansions: %" PetscInt_FMT "\n", nm->nexpand)); 103 PetscCall(PetscViewerASCIIPrintf(viewer, "reflections: %" PetscInt_FMT "\n", nm->nreflect)); 104 PetscCall(PetscViewerASCIIPrintf(viewer, "inside contractions: %" PetscInt_FMT "\n", nm->nincontract)); 105 PetscCall(PetscViewerASCIIPrintf(viewer, "outside contractionss: %" PetscInt_FMT "\n", nm->noutcontract)); 106 PetscCall(PetscViewerASCIIPrintf(viewer, "Shrink steps: %" PetscInt_FMT "\n", nm->nshrink)); 107 PetscCall(PetscViewerASCIIPopTab(viewer)); 108 } 109 PetscFunctionReturn(PETSC_SUCCESS); 110 } 111 112 static PetscErrorCode TaoSolve_NM(Tao tao) 113 { 114 TAO_NelderMead *nm = (TAO_NelderMead *)tao->data; 115 PetscReal *x; 116 PetscInt i; 117 Vec Xmur = nm->Xmur, Xmue = nm->Xmue, Xmuc = nm->Xmuc, Xbar = nm->Xbar; 118 PetscReal fr, fe, fc; 119 PetscInt shrink; 120 PetscInt low, high; 121 122 PetscFunctionBegin; 123 nm->nshrink = 0; 124 nm->nreflect = 0; 125 nm->nincontract = 0; 126 nm->noutcontract = 0; 127 nm->nexpand = 0; 128 129 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")); 130 131 PetscCall(VecCopy(tao->solution, nm->simplex[0])); 132 PetscCall(TaoComputeObjective(tao, nm->simplex[0], &nm->f_values[0])); 133 nm->indices[0] = 0; 134 for (i = 1; i < nm->N + 1; i++) { 135 PetscCall(VecCopy(tao->solution, nm->simplex[i])); 136 PetscCall(VecGetOwnershipRange(nm->simplex[i], &low, &high)); 137 if (i - 1 >= low && i - 1 < high) { 138 PetscCall(VecGetArray(nm->simplex[i], &x)); 139 x[i - 1 - low] += nm->lambda; 140 PetscCall(VecRestoreArray(nm->simplex[i], &x)); 141 } 142 143 PetscCall(TaoComputeObjective(tao, nm->simplex[i], &nm->f_values[i])); 144 nm->indices[i] = i; 145 } 146 147 /* Xbar = (Sum of all simplex vectors - worst vector)/N */ 148 PetscCall(NelderMeadSort(nm)); 149 PetscCall(VecSet(Xbar, 0.0)); 150 for (i = 0; i < nm->N; i++) PetscCall(VecAXPY(Xbar, 1.0, nm->simplex[nm->indices[i]])); 151 PetscCall(VecScale(Xbar, nm->oneOverN)); 152 tao->reason = TAO_CONTINUE_ITERATING; 153 while (1) { 154 /* Call general purpose update function */ 155 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 156 ++tao->niter; 157 shrink = 0; 158 PetscCall(VecCopy(nm->simplex[nm->indices[0]], tao->solution)); 159 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)); 160 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)); 161 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 162 if (tao->reason != TAO_CONTINUE_ITERATING) break; 163 164 /* x(mu) = (1 + mu)Xbar - mu*X_N+1 */ 165 PetscCall(VecAXPBYPCZ(Xmur, 1 + nm->mu_r, -nm->mu_r, 0, Xbar, nm->simplex[nm->indices[nm->N]])); 166 PetscCall(TaoComputeObjective(tao, Xmur, &fr)); 167 168 if (nm->f_values[nm->indices[0]] <= fr && fr < nm->f_values[nm->indices[nm->N - 1]]) { 169 /* reflect */ 170 nm->nreflect++; 171 PetscCall(PetscInfo(0, "Reflect\n")); 172 PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr)); 173 } else if (fr < nm->f_values[nm->indices[0]]) { 174 /* expand */ 175 nm->nexpand++; 176 PetscCall(PetscInfo(0, "Expand\n")); 177 PetscCall(VecAXPBYPCZ(Xmue, 1 + nm->mu_e, -nm->mu_e, 0, Xbar, nm->simplex[nm->indices[nm->N]])); 178 PetscCall(TaoComputeObjective(tao, Xmue, &fe)); 179 if (fe < fr) { 180 PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmue, fe)); 181 } else { 182 PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmur, fr)); 183 } 184 } else if (nm->f_values[nm->indices[nm->N - 1]] <= fr && fr < nm->f_values[nm->indices[nm->N]]) { 185 /* outside contraction */ 186 nm->noutcontract++; 187 PetscCall(PetscInfo(0, "Outside Contraction\n")); 188 PetscCall(VecAXPBYPCZ(Xmuc, 1 + nm->mu_oc, -nm->mu_oc, 0, Xbar, nm->simplex[nm->indices[nm->N]])); 189 190 PetscCall(TaoComputeObjective(tao, Xmuc, &fc)); 191 if (fc <= fr) 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]]) PetscCall(NelderMeadReplace(nm, nm->indices[nm->N], Xmuc, fc)); 200 else shrink = 1; 201 } 202 203 if (shrink) { 204 nm->nshrink++; 205 PetscCall(PetscInfo(0, "Shrink\n")); 206 207 for (i = 1; i < nm->N + 1; i++) { 208 PetscCall(VecAXPBY(nm->simplex[nm->indices[i]], 1.5, -0.5, nm->simplex[nm->indices[0]])); 209 PetscCall(TaoComputeObjective(tao, nm->simplex[nm->indices[i]], &nm->f_values[nm->indices[i]])); 210 } 211 PetscCall(VecAXPBY(Xbar, 1.5 * nm->oneOverN, -0.5, nm->simplex[nm->indices[0]])); 212 213 /* Add last vector's fraction of average */ 214 PetscCall(VecAXPY(Xbar, nm->oneOverN, nm->simplex[nm->indices[nm->N]])); 215 PetscCall(NelderMeadSort(nm)); 216 /* Subtract new last vector from average */ 217 PetscCall(VecAXPY(Xbar, -nm->oneOverN, nm->simplex[nm->indices[nm->N]])); 218 } 219 } 220 PetscFunctionReturn(PETSC_SUCCESS); 221 } 222 223 /*MC 224 TAONM - Nelder-Mead solver for derivative free, unconstrained minimization 225 226 Options Database Keys: 227 + -tao_nm_lambda - initial step length 228 - -tao_nm_mu - expansion/contraction factor 229 230 Level: beginner 231 M*/ 232 233 PETSC_EXTERN PetscErrorCode TaoCreate_NM(Tao tao) 234 { 235 TAO_NelderMead *nm; 236 237 PetscFunctionBegin; 238 PetscCall(PetscNew(&nm)); 239 tao->data = (void *)nm; 240 241 tao->ops->setup = TaoSetUp_NM; 242 tao->ops->solve = TaoSolve_NM; 243 tao->ops->view = TaoView_NM; 244 tao->ops->setfromoptions = TaoSetFromOptions_NM; 245 tao->ops->destroy = TaoDestroy_NM; 246 247 /* Override default settings (unless already changed) */ 248 PetscCall(TaoParametersInitialize(tao)); 249 PetscObjectParameterSetDefault(tao, max_it, 2000); 250 PetscObjectParameterSetDefault(tao, max_funcs, 4000); 251 252 nm->simplex = NULL; 253 nm->lambda = 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 PetscFunctionReturn(PETSC_SUCCESS); 260 } 261