(* Mathematica package *) Clear[setAttr] setAttr[mod_, attr_, val_] := (Unprotect[mod]; attr[mod] ^= val; Protect[mod]); Attributes[simplexStep] = HoldAll; simplexStep[function_, ps_,verbose_:False] := Module[ {countFn, maxVal, delt, deltB, pM,dp,vPrint}, vPrint[msg_] := If[ verbose, Print[msg] ]; delt[p_] := coord[p] - coord[pM]; deltB[p_] := coord[p] - coord[pB]; shrinkToBest[p_] := setAttr[p, coord, coord[p] - deltB[p]/2]; countFn[arg_] := (cnt += 1; function[arg]); (* wrap function with a counter *) maxVal = Max[value /@ ps]; pX = First[Select[ps, value[#] == maxVal &]]; setAttr[pM, coord, Mean[coord /@ Complement[ps,{pX}]]]; vPrint[coord[pM]]; vPrint["Reflecting"]; dp = {coord[pM],coord[pX],delt[pX]}; setAttr[pX, coord, coord[pM] - delt[pX]]; setAttr[pX, value, countFn[coord[pX]]]; (* addPrims[prims,ps]; *) Which[ (* CASE: pX is now the best *) value[pX] == Min[value /@ ps], (* grow? *) vPrint["checking for grow"]; If[ countFn[coord[pX] + delt[pX]] < value[pX], vPrint["grow"]; setAttr[pX, coord, coord[pX] + delt[pX]], vPrint["no grow"] ], (* CASE: pX is still the worst *) value[pX] == Max[value /@ ps], (* shrink? *) vPrint["still the worst, going to try shrinking"]; dp = {coord[pM],coord[pX],delt[pX]}; setAttr[pX, coord, coord[pX] - delt[pX]/2]; setAttr[pX, value, countFn[coord[pX]]]; If[ (* CASE: pX is still the worst *) value[pX] == Max[value /@ ps], vPrint["still still worst, going to try reflecting back across"]; setAttr[pX, coord, coord[pX] - 2 delt[pX]]; setAttr[pX, value, countFn[coord[pX]]]; If[ (* CASE: pX is still the worst *) value[pX] == Max[value /@ ps], vPrint["still still still worst, going to shrink towards best"]; bestVal = Min[value /@ ps]; pB = First[Select[ps, value[#] == bestVal &]]; shrinkToBest /@ ps ] ]]; addPrims[prims,ps]; {Min[value /@ ps], Max[value /@ ps]} ]; Clear[rosenbrock] rosenbrock[{x1_, x2_}] := (1 - x1)^2 + 100 (x2 - x1^2)^2; Plot3D[rosenbrock[{x1, x2}], {x1, -2, 2}, {x2, -4, 8}, Mesh -> False, PlotRange -> {0, 2000}, AxesLabel -> {"x1", "x2"}]; Attributes[addPrims] = HoldAll; addPrims[var_,ps_] := (var = Join[var,Subsets[ps, {2}] // Map[coord, #, {2}] & // Line /@ # & ]; var = Join[var, {Point[Mean[coord/@ ps]]}]);