| [2650] | 1 | #include "sopnamsp.h" | 
|---|
|  | 2 | #include "simplex.h" | 
|---|
|  | 3 | #include "ntuple.h" | 
|---|
|  | 4 | #include <math.h> | 
|---|
|  | 5 |  | 
|---|
|  | 6 | #include "timing.h" | 
|---|
|  | 7 |  | 
|---|
|  | 8 | //--------------------------------------------------------------- | 
|---|
|  | 9 | //-------------------  Classe   MinZFunction  ------------------- | 
|---|
|  | 10 | //--------------------------------------------------------------- | 
|---|
|  | 11 | // Interface de classe de function multivariable pour le SimplexMinmizer | 
|---|
|  | 12 |  | 
|---|
|  | 13 | MinZFunction::MinZFunction(unsigned int nvar) | 
|---|
|  | 14 | : mNVar(nvar) | 
|---|
|  | 15 | { | 
|---|
|  | 16 | } | 
|---|
|  | 17 |  | 
|---|
|  | 18 | MinZFunction::~MinZFunction() | 
|---|
|  | 19 | { | 
|---|
|  | 20 | } | 
|---|
|  | 21 |  | 
|---|
|  | 22 | //--------------------------------------------------------------- | 
|---|
|  | 23 | //-------------------  Classe   MinZFuncXi2  -------------------- | 
|---|
|  | 24 | //--------------------------------------------------------------- | 
|---|
|  | 25 | MinZFuncXi2::MinZFuncXi2(GeneralXi2* gxi2, GeneralFitData* gd) | 
|---|
|  | 26 | : mGXi2(gxi2) , mGData(gd), MinZFunction(gxi2->NPar()) | 
|---|
|  | 27 | { | 
|---|
|  | 28 | } | 
|---|
|  | 29 |  | 
|---|
|  | 30 | MinZFuncXi2::~MinZFuncXi2() | 
|---|
|  | 31 | { | 
|---|
|  | 32 | } | 
|---|
|  | 33 |  | 
|---|
|  | 34 | double MinZFuncXi2::Value(double const xp[]) | 
|---|
|  | 35 | { | 
|---|
|  | 36 | int ndataused; | 
|---|
|  | 37 | return mGXi2->Value(*mGData, const_cast<double *>(xp), ndataused); | 
|---|
|  | 38 | } | 
|---|
|  | 39 |  | 
|---|
|  | 40 | //--------------------------------------------------------------- | 
|---|
|  | 41 | //-------------------  Classe  MinZTestFunc   ------------------- | 
|---|
|  | 42 | //--------------------------------------------------------------- | 
|---|
|  | 43 | class MinZTestFunc : public  MinZFunction { | 
|---|
|  | 44 | public: | 
|---|
|  | 45 | MinZTestFunc(int sel); | 
|---|
|  | 46 | virtual double Value(double const xp[]); | 
|---|
|  | 47 | string  ToString(); | 
|---|
|  | 48 | Vector  OptParms(); | 
|---|
|  | 49 | protected: | 
|---|
|  | 50 | static int ISelToNvar(int isel); | 
|---|
|  | 51 | int mSel; | 
|---|
|  | 52 | }; | 
|---|
|  | 53 |  | 
|---|
|  | 54 | int MinZTestFunc::ISelToNvar(int isel) | 
|---|
|  | 55 | { | 
|---|
|  | 56 | if (isel == 0) return 1; | 
|---|
|  | 57 | if (isel == 1) return 1; | 
|---|
|  | 58 | else if (isel == 2) return 1; | 
|---|
|  | 59 | else if (isel == 3) return 2; | 
|---|
|  | 60 | else if (isel == 4) return 3; | 
|---|
|  | 61 | else return 1; | 
|---|
|  | 62 | } | 
|---|
|  | 63 |  | 
|---|
|  | 64 | MinZTestFunc::MinZTestFunc(int sel) | 
|---|
|  | 65 | : MinZFunction(ISelToNvar(sel)) | 
|---|
|  | 66 | { | 
|---|
|  | 67 | if ((sel < 0) || (sel > 4)) sel = 0; | 
|---|
|  | 68 | mSel = sel; | 
|---|
|  | 69 | } | 
|---|
|  | 70 |  | 
|---|
|  | 71 | string MinZTestFunc::ToString() | 
|---|
|  | 72 | { | 
|---|
|  | 73 | string rs; | 
|---|
|  | 74 | if (mSel == 0) { | 
|---|
|  | 75 | rs = "-x+(x-2)^2"; | 
|---|
|  | 76 | } | 
|---|
|  | 77 | else if (mSel == 1) { | 
|---|
|  | 78 | rs = "0.1*x^2-3exp(-(x-2)^2)-5*exp(-0.5*(x+3)^2)"; | 
|---|
|  | 79 | } | 
|---|
|  | 80 | else if (mSel == 2) { | 
|---|
|  | 81 | rs = "0.1*x^2-3exp(-(x-2)^2)+5*exp(-0.5*(x+3)^2)"; | 
|---|
|  | 82 | } | 
|---|
|  | 83 | else if (mSel == 3) { | 
|---|
|  | 84 | rs = "1.3*(x-50.35)^2+25*(y+3.14)^2"; | 
|---|
|  | 85 | } | 
|---|
|  | 86 | else if (mSel == 4) { | 
|---|
|  | 87 | rs = "(x-2.2)^2+2.*(y+3.6)^2+3.*(z-1.1)^2"; | 
|---|
|  | 88 | } | 
|---|
|  | 89 | else rs = "????"; | 
|---|
|  | 90 | return  rs; | 
|---|
|  | 91 | } | 
|---|
|  | 92 |  | 
|---|
|  | 93 | Vector MinZTestFunc::OptParms() | 
|---|
|  | 94 | { | 
|---|
|  | 95 | Vector xx; | 
|---|
|  | 96 | if (mSel == 0) { | 
|---|
|  | 97 | Vector rv(1); | 
|---|
|  | 98 | rv = 2.5; | 
|---|
|  | 99 | return rv; | 
|---|
|  | 100 | } | 
|---|
|  | 101 | else if (mSel == 1) { | 
|---|
|  | 102 | Vector rv(1); | 
|---|
|  | 103 | rv = -2.883; | 
|---|
|  | 104 | return rv; | 
|---|
|  | 105 | } | 
|---|
|  | 106 | else if (mSel == 2) { | 
|---|
|  | 107 | Vector rv(1); | 
|---|
|  | 108 | rv = 1.812; | 
|---|
|  | 109 | return rv; | 
|---|
|  | 110 | } | 
|---|
|  | 111 | else if (mSel == 3) { | 
|---|
|  | 112 | Vector rv(2); | 
|---|
|  | 113 | rv(0) = 50.35; | 
|---|
|  | 114 | rv(1) = -3.14; | 
|---|
|  | 115 | return rv; | 
|---|
|  | 116 | } | 
|---|
|  | 117 | else if (mSel == 4) { | 
|---|
|  | 118 | Vector rv(3); | 
|---|
|  | 119 | rv(0) = 2.2; | 
|---|
|  | 120 | rv(1) = -3.6; | 
|---|
|  | 121 | rv(2) = 1.1; | 
|---|
|  | 122 | return rv; | 
|---|
|  | 123 | } | 
|---|
|  | 124 | else xx = 0.; | 
|---|
|  | 125 | return xx ; | 
|---|
|  | 126 | } | 
|---|
|  | 127 |  | 
|---|
|  | 128 |  | 
|---|
|  | 129 | double MinZTestFunc::Value(double const xp[]) | 
|---|
|  | 130 | { | 
|---|
|  | 131 | double retval = 0; | 
|---|
|  | 132 | if (mSel == 0) { | 
|---|
|  | 133 | double x = xp[0]; | 
|---|
|  | 134 | retval = -x+(x-2.)*(x-2.); | 
|---|
|  | 135 | } | 
|---|
|  | 136 | else if ((mSel == 1) || (mSel == 2)) { | 
|---|
|  | 137 | double x = xp[0]; | 
|---|
|  | 138 | retval = 0.1*x*x; | 
|---|
|  | 139 | x = xp[0]-2.; | 
|---|
|  | 140 | x = x*x; | 
|---|
|  | 141 | retval -= 3*exp(-x); | 
|---|
|  | 142 | x = xp[0]+3.; | 
|---|
|  | 143 | x = 0.5*x*x; | 
|---|
|  | 144 | if (mSel == 1) retval -= 5*exp(-x); | 
|---|
|  | 145 | else retval += 5*exp(-x); | 
|---|
|  | 146 | } | 
|---|
|  | 147 | else if (mSel == 3) { | 
|---|
|  | 148 | double x = xp[0]-50.35; | 
|---|
|  | 149 | double y = xp[1]+3.14; | 
|---|
|  | 150 | retval = 1.3*x*x+25.*y*y; | 
|---|
|  | 151 | } | 
|---|
|  | 152 | else if (mSel == 4) { | 
|---|
|  | 153 | double x = xp[0]-2.2; | 
|---|
|  | 154 | double y = xp[1]+3.6; | 
|---|
|  | 155 | double z = xp[2]-1.1; | 
|---|
|  | 156 | retval = x*x+2.*y*y+3.*z*z; | 
|---|
|  | 157 | } | 
|---|
|  | 158 | else retval = 0.; | 
|---|
|  | 159 | return retval; | 
|---|
|  | 160 | } | 
|---|
|  | 161 |  | 
|---|
|  | 162 | //--------------------------------------------------------------- | 
|---|
|  | 163 | //-------------------  Classe   MinZSimplex  -------------------- | 
|---|
|  | 164 | //--------------------------------------------------------------- | 
|---|
|  | 165 | string __Vec2Str4MinZ_AutoTest(Vector& xx) | 
|---|
|  | 166 | { | 
|---|
|  | 167 | string rs; | 
|---|
|  | 168 | char buff[32]; | 
|---|
|  | 169 | for(int i=0; i<xx.Size(); i++) { | 
|---|
|  | 170 | sprintf(buff," %g " , xx(i)); | 
|---|
|  | 171 | rs += buff; | 
|---|
|  | 172 | } | 
|---|
|  | 173 | return rs; | 
|---|
|  | 174 | } | 
|---|
|  | 175 |  | 
|---|
|  | 176 | int MinZSimplex::AutoTest(int tsel, int prtlev) | 
|---|
|  | 177 | { | 
|---|
|  | 178 | int rc = 0; | 
|---|
|  | 179 | cout << " --- MinZSimplex::AutoTest() --- TSel= " << tsel << " PrtLev=" << prtlev << endl; | 
|---|
|  | 180 | for(int i=0; i<5; i++) { | 
|---|
|  | 181 | if ((tsel >= 0) && (tsel != i))  continue; | 
|---|
|  | 182 | cout << " ======= Test avec ISel= " << i; | 
|---|
|  | 183 | Vector xx; | 
|---|
|  | 184 | MinZTestFunc mzf(i); | 
|---|
|  | 185 | cout << " - Func= " << mzf.ToString() << endl; | 
|---|
|  | 186 | Vector rv = mzf.OptParms(); | 
|---|
|  | 187 | xx = rv; | 
|---|
|  | 188 | for(int j=0; j<2; j++) { | 
|---|
|  | 189 | double vi = 50.*(j-0.5); | 
|---|
|  | 190 | for(int k=0; k<2; k++) { | 
|---|
|  | 191 | double vs = (k == 0) ? 1. : 10. ; | 
|---|
|  | 192 | cout << "--[" << j << "," << k | 
|---|
|  | 193 | << "] Initialisation avec IniPoint= " << vi << " IniStep= " << vs << endl; | 
|---|
|  | 194 | MinZSimplex simplex(&mzf); | 
|---|
|  | 195 | xx = vi; | 
|---|
|  | 196 | simplex.SetInitialPoint(xx); | 
|---|
|  | 197 | xx = vs; | 
|---|
|  | 198 | simplex.SetInitialStep(xx); | 
|---|
|  | 199 | simplex.SetPrtLevel(prtlev); | 
|---|
|  | 200 | int rcs = simplex.Minimize(xx); | 
|---|
|  | 201 | Vector diff = rv-xx; | 
|---|
|  | 202 | double d2 = diff.Norm2(); | 
|---|
|  | 203 | cout << " Rc(simplex.Minimize() = " << rc << " NIter= " | 
|---|
|  | 204 | << simplex.NbIter() << " ===> Distance^2= " << d2 | 
|---|
|  | 205 | << "\nConverged to " <<  __Vec2Str4MinZ_AutoTest(xx) | 
|---|
|  | 206 | << "  Best Value= " << __Vec2Str4MinZ_AutoTest(rv) | 
|---|
|  | 207 | << "  Diff = " << __Vec2Str4MinZ_AutoTest(diff) << endl; | 
|---|
|  | 208 | if ((rcs > 5) || (d2 > 0.5))  rc ++; | 
|---|
|  | 209 | } | 
|---|
|  | 210 | } | 
|---|
|  | 211 | } | 
|---|
|  | 212 | cout << " --- MinZSimplex::AutoTest() --- Rc=" << rc << " -- END ----- " << endl; | 
|---|
|  | 213 | return rc; | 
|---|
|  | 214 | } | 
|---|
|  | 215 |  | 
|---|
|  | 216 | MinZSimplex::MinZSimplex(MinZFunction *mzf) | 
|---|
|  | 217 | : mZF(mzf) , mPoint0(mZF->NVar()) , mStep0(mZF->NVar()) | 
|---|
|  | 218 | { | 
|---|
|  | 219 | SetMaxIter(); | 
|---|
|  | 220 | SetControls(); | 
|---|
|  | 221 | Vector xx(NDim()); | 
|---|
|  | 222 | xx = 0.; | 
|---|
|  | 223 | SetInitialPoint(xx); | 
|---|
|  | 224 | xx = 1.0; | 
|---|
|  | 225 | SetInitialStep(xx); | 
|---|
|  | 226 | SetStopTolerance(); | 
|---|
|  | 227 | mIter = -1; | 
|---|
|  | 228 | mStop = -1; | 
|---|
|  | 229 | SetPrtLevel(); | 
|---|
|  | 230 | } | 
|---|
|  | 231 |  | 
|---|
|  | 232 | MinZSimplex::~MinZSimplex() | 
|---|
|  | 233 | { | 
|---|
|  | 234 | } | 
|---|
|  | 235 |  | 
|---|
|  | 236 | int MinZSimplex::Minimize(Vector& fpoint) | 
|---|
|  | 237 | { | 
|---|
|  | 238 | // vector< TVector<r_8> > splx; | 
|---|
|  | 239 | Vector splx[100]; | 
|---|
|  | 240 | Vector Y(NDim()+1); | 
|---|
|  | 241 | // On calcule le simplex initial | 
|---|
|  | 242 | // N = NDim, N+1 points (pp) ds l'espace a N dimensions | 
|---|
|  | 243 | // Point0, Point0 + Step0(i) e_i | 
|---|
|  | 244 | Vector pp,ppc; | 
|---|
|  | 245 | pp = mPoint0; | 
|---|
|  | 246 | //ppc = pp; | 
|---|
|  | 247 | //splx.push_back(ppc); | 
|---|
|  | 248 | splx[0] = pp; | 
|---|
|  | 249 | int i,j,k; | 
|---|
|  | 250 | for(i=0; i<NDim(); i++) { | 
|---|
|  | 251 | Vector pps; | 
|---|
|  | 252 | pps = mPoint0; | 
|---|
|  | 253 | pps(i) += mStep0(i); | 
|---|
|  | 254 | //splx.push_back(pps); | 
|---|
|  | 255 | splx[i+1] = pps; | 
|---|
|  | 256 | } | 
|---|
|  | 257 | int mpts = NDim()+1; | 
|---|
|  | 258 | // calcul des valeurs de la fonction sur les sommets | 
|---|
|  | 259 | for(i=0; i<mpts; i++) | 
|---|
|  | 260 | Y(i) = Value(splx[i]); | 
|---|
|  | 261 |  | 
|---|
|  | 262 | int iter = 0; | 
|---|
|  | 263 | mIter = iter; | 
|---|
|  | 264 | mStop = 0; | 
|---|
|  | 265 |  | 
|---|
|  | 266 | int nbugrtol2 = 0; | 
|---|
|  | 267 | bool stop = false, stop0=false; | 
|---|
|  | 268 | int rc = 0; | 
|---|
|  | 269 | int ilo, ihi, inhi; | 
|---|
|  | 270 | int move = 0; | 
|---|
|  | 271 | char* smov[6] = { "None", "Reflection", "ReflecExpand", "ContractHigh", "ContractLow", "ExpandHigh" }; | 
|---|
|  | 272 | int movcnt[6] = {0,0,0,0,0,0}; | 
|---|
|  | 273 |  | 
|---|
|  | 274 | int nrep1=0, nrep2=0; | 
|---|
|  | 275 | FindMinMax12(Y, ilo, ihi, inhi); | 
|---|
|  | 276 | double yhilast = Y(ihi); | 
|---|
|  | 277 | yhilast += fabs(yhilast); | 
|---|
|  | 278 |  | 
|---|
|  | 279 | while (!stop) {  // | 
|---|
|  | 280 | FindMinMax12(Y, ilo, ihi, inhi); | 
|---|
|  | 281 | double ymean = (fabs(Y(ihi))+fabs(Y(ilo))); | 
|---|
|  | 282 | if (ymean < mTol0) { stop0 = true; ymean = mTol0; } | 
|---|
|  | 283 | double rtol1 = 2.*fabs(Y(ihi)-Y(ilo))/ymean; | 
|---|
|  | 284 | double ym2 = (fabs(yhilast)+fabs(Y(ihi))); | 
|---|
|  | 285 | if (ym2 < mTol0) ym2 = mTol0; | 
|---|
|  | 286 | double rtol2 = 2.*(yhilast-Y(ihi))/ym2; | 
|---|
|  | 287 | yhilast = Y(ihi); | 
|---|
|  | 288 | if (rtol2 < 0.) { | 
|---|
|  | 289 | if (move != 40) { | 
|---|
|  | 290 | cout << " !!!! MinZSimplex::Minimize() BUG RTol2< 0. --> Chs " << endl; | 
|---|
|  | 291 | nbugrtol2++; | 
|---|
|  | 292 | } | 
|---|
|  | 293 | else nrep2 = 0; | 
|---|
|  | 294 | rtol2 = -rtol2; | 
|---|
|  | 295 | } | 
|---|
|  | 296 | if (PrtLevel() > 1) | 
|---|
|  | 297 | cout << "--MinZSimplex::Minimize() - Iter=" << iter | 
|---|
|  | 298 | << " Move= " << move << " (" <<  smov[move/10] << ")" << endl; | 
|---|
|  | 299 | if (PrtLevel() > 2) | 
|---|
|  | 300 | cout << "..ILO=" << ilo << " IHI=" << ihi << " INHI=" << inhi | 
|---|
|  | 301 | << " Y(ILO)=" << Y(ilo) << " Y(IHI)=" << Y(ihi) << "\n" | 
|---|
|  | 302 | << "...YMean_Abs=" << ymean <<  " RTOL1=" << rtol1 << " RTOL2=" << rtol2 <<  endl; | 
|---|
|  | 303 | if (PrtLevel() > 3) { | 
|---|
|  | 304 | for(i=0; i<mpts; i++) { | 
|---|
|  | 305 | cout << "....Simplex[" << i << "]= "; | 
|---|
|  | 306 | for(j=0; j<NDim(); j++) cout << splx[i](j) << " , "; | 
|---|
|  | 307 | cout << " Y=Value= " << Y(i) << endl; | 
|---|
|  | 308 | } | 
|---|
|  | 309 | } | 
|---|
|  | 310 | if (rtol1 < mTol1) nrep1++; | 
|---|
|  | 311 | else nrep1 = 0; | 
|---|
|  | 312 | if (rtol2 < mTol2) nrep2++; | 
|---|
|  | 313 | else nrep2 = 0; | 
|---|
|  | 314 |  | 
|---|
|  | 315 | if (stop0) { mStop = 1; rc = 0; stop = true; break; } | 
|---|
|  | 316 | if (nrep1 > mRep1) { mStop = 2; rc = 0; stop = true; break; } | 
|---|
|  | 317 | if (nrep2 > mRep2) { mStop = 3; rc = 0; stop = true; break; } | 
|---|
|  | 318 | if (iter > MaxIter() ) { mStop = 0, rc = iter;  break; } | 
|---|
|  | 319 | iter++; | 
|---|
|  | 320 | if (iter > 0)  movcnt[move/10]++; | 
|---|
|  | 321 |  | 
|---|
|  | 322 | // Next iteration, on modifie le simplex | 
|---|
|  | 323 | // Calcul du centre de gravite su simplex, hors le point le + haut | 
|---|
|  | 324 | Vector pbar(NDim()); | 
|---|
|  | 325 | pbar = 0.; | 
|---|
|  | 326 | for(i=0; i<mpts; i++) { | 
|---|
|  | 327 | if (i == ihi)  continue; | 
|---|
|  | 328 | pbar += splx[i]; | 
|---|
|  | 329 | } | 
|---|
|  | 330 | pbar /= (double)NDim(); | 
|---|
|  | 331 | // On calcule le sommet oppose a point IHI (le + haut) | 
|---|
|  | 332 | Vector pr, prr; | 
|---|
|  | 333 | double YPR, YPRR; | 
|---|
|  | 334 | pr = (1.+Alpha())*pbar-Alpha()*splx[ihi]; | 
|---|
|  | 335 | YPR = Value(pr); | 
|---|
|  | 336 | if (YPR < Y(ilo)) {   // Amelioaration par rapport au meilleur point, | 
|---|
|  | 337 | // on va plus loin d'un facteur gamma | 
|---|
|  | 338 | prr = Gamma()*pr+(1.-Gamma())*pbar; | 
|---|
|  | 339 | YPRR = Value(prr); | 
|---|
|  | 340 | if (YPRR < Y(ilo)) {  // On remplace le IHI par YPRR | 
|---|
|  | 341 | splx[ihi] = prr; | 
|---|
|  | 342 | Y(ihi) = YPRR; | 
|---|
|  | 343 | move = 20; | 
|---|
|  | 344 | } | 
|---|
|  | 345 | else {  // sinon, on remplace par YPR | 
|---|
|  | 346 | splx[ihi] = pr; | 
|---|
|  | 347 | Y(ihi) = YPR; | 
|---|
|  | 348 | move = 10; | 
|---|
|  | 349 | } | 
|---|
|  | 350 | } | 
|---|
|  | 351 | else {  // Moins bon que le meilleur point .. | 
|---|
|  | 352 | if (YPR > Y(inhi)) {  // Plus mauvais que le second plus haut (INHI) | 
|---|
|  | 353 | if (YPR < Y(ihi)) {   // Mais meilleur que le plus haut (IHI) | 
|---|
|  | 354 | splx[ihi] = pr;     // On remplace donc le plus haut | 
|---|
|  | 355 | Y(ihi) = YPR; | 
|---|
|  | 356 | move = 11; | 
|---|
|  | 357 | } | 
|---|
|  | 358 | else { // Plus mauvais que le plus mauvais IHI | 
|---|
|  | 359 | // on tente avec un point intermediaire | 
|---|
|  | 360 | prr = Beta()*splx[ihi]+(1.-Beta())*pbar; | 
|---|
|  | 361 | YPRR = Value(prr); | 
|---|
|  | 362 | if (YPRR < Y(ihi)) {   // Le point intermediaire ameliore les choses | 
|---|
|  | 363 | splx[ihi] = prr;     // On remplace donc le point le + haut | 
|---|
|  | 364 | Y(ihi) = YPRR; | 
|---|
|  | 365 | move = 30; | 
|---|
|  | 366 | } | 
|---|
|  | 367 | else { | 
|---|
|  | 368 | // On tente aussi de rester du meme cote, mais aller plus loin | 
|---|
|  | 369 | prr = Gamma2()*splx[ihi]+(1.-Gamma2())*pbar; | 
|---|
|  | 370 | YPRR = Value(prr); | 
|---|
|  | 371 | if (YPRR < Y(ihi)) {   // Le point intermediaire ameliore les choses | 
|---|
|  | 372 | splx[ihi] = prr;     // On remplace donc le point le + haut | 
|---|
|  | 373 | Y(ihi) = YPRR; | 
|---|
|  | 374 | move = 50; | 
|---|
|  | 375 | } | 
|---|
|  | 376 | else { | 
|---|
|  | 377 | // Rien n'y fait, on contracte autour du meilleur point | 
|---|
|  | 378 | for(i=0; i<mpts; i++) { | 
|---|
|  | 379 | if (i == ilo)  continue; | 
|---|
|  | 380 | splx[i] = Beta2()*splx[i]+(1.-Beta())*splx[ilo]; | 
|---|
|  | 381 | Y(i) = Value(splx[i]); | 
|---|
|  | 382 | move = 40; | 
|---|
|  | 383 | } | 
|---|
|  | 384 | } | 
|---|
|  | 385 | } | 
|---|
|  | 386 | } | 
|---|
|  | 387 | } | 
|---|
|  | 388 | else {  // Meilleur que le IHI et le INHI | 
|---|
|  | 389 | splx[ihi] = pr;     // On remplace le plus haut | 
|---|
|  | 390 | Y(ihi) = YPR; | 
|---|
|  | 391 | move = 12; | 
|---|
|  | 392 | } | 
|---|
|  | 393 | } | 
|---|
|  | 394 | }   // Fin de la boucle while principale | 
|---|
|  | 395 |  | 
|---|
|  | 396 | fpoint = splx[ilo]; | 
|---|
|  | 397 | mIter = iter; | 
|---|
|  | 398 |  | 
|---|
|  | 399 | if (PrtLevel() > 0) { | 
|---|
|  | 400 | string sr; | 
|---|
|  | 401 | StopReason(sr); | 
|---|
|  | 402 | cout << "-----MinZSimplex::Minimize()/Ended - NIter=" << iter | 
|---|
|  | 403 | << " Moves[0..5]= " << movcnt[0] << "," << movcnt[1] << "," | 
|---|
|  | 404 | << movcnt[2] << "," << movcnt[3] << "," | 
|---|
|  | 405 | << movcnt[4] << "," << movcnt[5] | 
|---|
|  | 406 | << "\n..MinZSimplex Stop=" << StopReason() << " -> " << sr << endl; | 
|---|
|  | 407 |  | 
|---|
|  | 408 | if (nbugrtol2 > 0)  cout << "MinZSimplex::Minimize()/Warning - nbugrtol2= " << nbugrtol2 << endl; | 
|---|
|  | 409 | } | 
|---|
|  | 410 | return rc; | 
|---|
|  | 411 | } | 
|---|
|  | 412 |  | 
|---|
|  | 413 | int MinZSimplex::StopReason(string& s) | 
|---|
|  | 414 | { | 
|---|
|  | 415 | char* sr[5] = { "NoConverg, MaxIterReached", "OK, fm<Tol0", "OK, Df/f<Tol1", | 
|---|
|  | 416 | "OK, [Df/f max]Iter<Tol2" "Error - Wrong StopReason" }; | 
|---|
|  | 417 | int stop = mStop; | 
|---|
|  | 418 | if ((stop < 0) || (stop > 3)) stop = 4; | 
|---|
|  | 419 | s = sr[stop]; | 
|---|
|  | 420 | return mStop; | 
|---|
|  | 421 | } | 
|---|
|  | 422 |  | 
|---|
|  | 423 | int MinZSimplex::FindMinMax12(Vector& fval, int& ilo, int& ihi, int& inhi) | 
|---|
|  | 424 | { | 
|---|
|  | 425 | ilo = 0; | 
|---|
|  | 426 | if (fval(0) > fval(1)) { ihi = 0;  inhi = 1; } | 
|---|
|  | 427 | else { ihi = 1;  inhi = 0; } | 
|---|
|  | 428 |  | 
|---|
|  | 429 | for(int k=0; k<fval.Size(); k++) { | 
|---|
|  | 430 | if (fval(k) < fval(ilo))  ilo = k; | 
|---|
|  | 431 | if (fval(k) > fval(ihi)) { | 
|---|
|  | 432 | inhi = ihi; | 
|---|
|  | 433 | ihi = k; | 
|---|
|  | 434 | } | 
|---|
|  | 435 | else if (fval(k) > fval(inhi)) { | 
|---|
|  | 436 | if (k != ihi)  inhi = k;  // ce test n'est peut-etre pas necessaire ??? | 
|---|
|  | 437 | } | 
|---|
|  | 438 | } | 
|---|
|  | 439 | return ilo; | 
|---|
|  | 440 | } | 
|---|