NelderMeadMinimize`Dump`CompiledNelderMead

gistfile1.c

(* Produces compiled code for the Nelder-Mead algorithm with the objective function inlined. *)
(* The objective function takes the form F[parametersToOptimize..,constantParameters] *)
NelderMeadMinimize`Dump`CompiledNelderMead[
   objectiveFunction_Function | objectiveFunction_CompiledFunction, vars : {__Symbol}, const: {__Symbol},
   opts : OptionsPattern[NelderMeadMinimize`Dump`CompiledNelderMead]
  ] :=
  NelderMeadMinimize`Dump`CompiledNelderMead[
    objectiveFunction, vars, const
    opts
   ] =
   With[{
     (* Inlined option values *)
     historyLength = If[# === Automatic, 10 Length[vars], #] & @ OptionValue["HistoryLength"],
     reflectRatio = OptionValue["ReflectRatio"], expandRatio = OptionValue["ExpandRatio"],
     contractRatio = OptionValue["ContractRatio"], shrinkRatio = OptionValue["ShrinkRatio"],
     (* Other inlined values *)
     origin = ConstantArray[0., Length[vars]],
     infinity = $MaxMachineNumber,
     epsilon = $MachineEpsilon,
     (* Inlined functions *)
     f = apply[objectiveFunction, Evaluate[vars~Join~const]],
     diffs = cumulativeAbsoluteDifferences,
     (* Options to be passed to Compile *)
     compileopts = Sequence @@ If[$VersionNumber >= 8, {
         (* Mathematica 8 and above offer improved behaviour using these options *)
         RuntimeOptions -> {"Speed", "CompareWithTolerance" -> True, "EvaluateSymbolically" -> False},
         CompilationTarget -> OptionValue[CompilationTarget],
         CompilationOptions->{"ExpressionOptimization"->True, "InlineCompiledFunctions"->Automatic, "InlineExternalDefinitions"->True}
        }, {
         (* Ordering is an external call in Mathematica 7 and so needs type information *)
         {{_Ordering, _Integer, 1}}
        }
       ]
    },
    Compile[{{pts, _Real, 2}, {cst, _Real, 1},{tol, _Real, 0}, {maxit, _Integer, 0}},
     Block[{
       (* Housekeeping *)
       history = Table[infinity, {historyLength}], iteration = maxit,
       (* Basic quantities *)
       simplex = pts, vals = f[#~Join~cst]& /@ pts, ordering,
       (* Calculated points and function values *)
       centroid = origin,
       reflectedPoint = origin, reflectedValue = infinity,
       expandedPoint = origin, expandedValue = infinity,
       contractedPoint = origin, contractedValue = infinity,
       (* More readable indices into the simplex array *)
       best = 1, worst = -1, rest = Rest@Range@Length[pts],
       (* Operation counts (for debugging purposes) *)
       evaluations = Length[pts],
       reflections = 0, expansions = 0, contractions = 0, shrinkages = 0
      },
      While[
       (* Order simplex points by function value *)
       ordering = Ordering[vals];
       vals = vals[[ordering]]; simplex = simplex[[ordering]];
       (* Decrement and test iterator *)
       (iteration--) != 0,
       (* Check for convergence *)
       history[[1]] = vals[[best]]; history = RotateLeft[history];
       If[diffs[history] <= tol + epsilon diffs[history],
        Break[]
       ];
       (* Find centroid of first (N - 1) points *)
       centroid = Mean@Most[simplex];
       (* Reflect *)
       reflectedPoint = centroid + reflectRatio (centroid - simplex[[worst]]);
       reflectedValue = f[reflectedPoint~Join~cst]; ++evaluations;
       If[vals[[best]] <= reflectedValue < vals[[-2]],
        vals[[worst]] = reflectedValue; simplex[[worst]] = reflectedPoint;
        ++reflections; Continue[]
       ];
       (* Expand *)
       If[reflectedValue < vals[[best]],
        expandedPoint = centroid + expandRatio (reflectedPoint - centroid);
        expandedValue = f[expandedPoint~Join~cst]; ++evaluations;
        If[expandedValue < reflectedValue,
         vals[[worst]] = expandedValue; simplex[[worst]] = expandedPoint;
         ++expansions; Continue[],
         vals[[worst]] = reflectedValue; simplex[[worst]] = reflectedPoint;
         ++reflections; Continue[]
        ];
       ];
       (* Contract *)
       If[reflectedValue < vals[[worst]],
        (* Outside contraction *)
        contractedPoint = centroid + contractRatio (reflectedPoint - centroid);
        contractedValue = f[contractedPoint~Join~cst]; ++evaluations;
        If[contractedValue <= reflectedValue,
         vals[[worst]] = contractedValue; simplex[[worst]] = contractedPoint;
         ++contractions; Continue[]
        ];,
        (* Inside contraction *)
        contractedPoint = centroid - contractRatio (centroid - simplex[[worst]]);
        contractedValue = f[contractedPoint~Join~cst]; ++evaluations;
        If[contractedValue < vals[[worst]],
         vals[[worst]] = contractedValue; simplex[[worst]] = contractedPoint;
         ++contractions; Continue[]
        ];
       ];
       (* Shrink *)
       simplex[[rest]] = simplex[[best]] + shrinkRatio (simplex[[rest]] - simplex[[best]]);
       vals[[rest]] = f /@ simplex[[rest]]; 
       evaluations += Length[rest] - 1;
       ++shrinkages;
      ];
      (* A call out of the VM is necessary to return the results *)
     (* results = {vals, simplex, {evaluations, reflections, expansions, contractions, shrinkages}};*)
      First[simplex]~Join~{evaluations, reflections, expansions, contractions, shrinkages}
      (*{evaluations, reflections, expansions, contractions, shrinkages}*)
     ], compileopts
    ]
   ];