using System; using System.Collections.Generic; using System.Text; namespace Petronode.CommonData { /// /// Based on the 1965 research paper, "A Simplex Method for Function Minimization", /// by J.A. Nelder and R. Mead /// public class Nelder_Mead_Optimization { static Random m_random = new Random(1); public int Dimension; // Problem dimension public Nelder_Mead_Solution[] Solutions; // Potential solutions public double[] minX; // Solution constraints public double[] maxX; public double alpha = 1.0; // Reflection public double beta = 0.5; // Contraction public double gamma = 2.0; // Expansion public double Tolerance = 0.001; // The main solving loop stops upon reaching this target public int maxLoop; // Limits the main solving loop public delegate double ObjectiveFunctionDelegate(double[] vector); public ObjectiveFunctionDelegate onObjectiveFunction = null; public delegate void ProgressReportDelegate(int percent, object info); public ProgressReportDelegate onReportProgress = null; /// /// Constructor, generates the class /// /// size Number of Amoeba "legs" /// minX Constraints minimums /// maxX Constraints maximums /// maxLoop Maximum number of iterations to prevent locking public Nelder_Mead_Optimization(int size, double[] minX, double[] maxX, int maxLoop) { Dimension = minX.Length; this.minX = new double[Dimension]; this.maxX = new double[Dimension]; Array.Copy(minX, this.minX, this.minX.Length); Array.Copy(maxX, this.maxX, this.maxX.Length); this.maxLoop = maxLoop; Solutions = new Nelder_Mead_Solution[size]; for (int i = 0; i < Solutions.Length; ++i) Solutions[i] = NMS_Factory_Random(); Array.Sort(Solutions); } /// /// Constructor, generates the class /// /// size Number of Amoeba "legs" /// initialX First point given /// minX Constraints minimums /// maxX Constraints maximums /// maxLoop Maximum number of iterations to prevent locking public Nelder_Mead_Optimization(int size, double[] initialX, double[] minX, double[] maxX, int maxLoop) { Dimension = minX.Length; this.minX = new double[Dimension]; this.maxX = new double[Dimension]; Array.Copy(minX, this.minX, this.minX.Length); Array.Copy(maxX, this.maxX, this.maxX.Length); this.maxLoop = maxLoop; Solutions = new Nelder_Mead_Solution[size]; Solutions[0] = NMS_Factory_Deterministic(initialX); for (int i = 1; i < Solutions.Length; ++i) Solutions[i] = NMS_Factory_Random(); Array.Sort(Solutions); } /// /// Performs solution initialization (either after the random or manual) /// public void InitSolutions() { if (onObjectiveFunction == null) return; foreach (Nelder_Mead_Solution nms in Solutions) nms.Value = onObjectiveFunction(nms.Vector); Array.Sort(Solutions); } /// /// Computes the target function by the current delegate /// If delegate is not set, result is 0 /// public double ComputeTarget(double[] v) { if (onObjectiveFunction == null) return 0.0; return onObjectiveFunction(v); } /// /// Computes the target function by the current delegate /// If delegate is not set, result is 0. nms is not modified /// public double ComputeTarget(Nelder_Mead_Solution nms) { if (onObjectiveFunction == null) return 0.0; return onObjectiveFunction(nms.Vector); } /// /// Runs the solution /// /// the best solution found public Nelder_Mead_Solution LocateTarget() { for (int t = 0; t < maxLoop; t++) { if (BestSolution.Value < Tolerance) break; // solution in tolerance // progress report if ((t % 10 == 0) && onReportProgress != null) onReportProgress(t * 100 / maxLoop, (object)BestSolution.ToString()); Nelder_Mead_Solution centroid = GetCentroid(); Nelder_Mead_Solution reflected = GetReflected(centroid); if (reflected.Value < BestSolution.Value) { Nelder_Mead_Solution expanded = GetExpanded(reflected, centroid); ReplaceTheWorstWith((expanded.Value < BestSolution.Value) ? expanded : reflected); continue; } if (!IsWorseThanAllButTheWorst(reflected)) { ReplaceTheWorstWith(reflected); continue; } if (reflected.Value <= WorstSolution.Value) ReplaceTheWorstWith(reflected); Nelder_Mead_Solution contracted = GetContracted(centroid); if (contracted.Value > WorstSolution.Value) ShrinkTowardsTheBestSolution(); else ReplaceTheWorstWith(contracted); } return BestSolution; } /// /// Retrieves the best solution /// public Nelder_Mead_Solution BestSolution { get { return Solutions[0]; } } /// /// Retrieves the worst solution /// public Nelder_Mead_Solution WorstSolution { get { return Solutions[Solutions.Length - 1]; } } /// /// Converts the class to string for debug /// public override string ToString() { StringBuilder sb = new StringBuilder(); for (int i = 0; i < Solutions.Length; i++) { Nelder_Mead_Solution s = Solutions[i]; sb.Append(i.ToString("Vector[0]: ")); sb.Append(Solutions[i].ToString()); sb.Append("\r\n"); } return sb.ToString(); } #region Private methods private Nelder_Mead_Solution GetCentroid() { Nelder_Mead_Solution nms = new Nelder_Mead_Solution(BestSolution); for (int i = 1; i < Solutions.Length - 1; i++) nms.Append(Solutions[i]); nms.Scale(1.0 / Convert.ToDouble(Solutions.Length - 1)); return NMS_Factory_Deterministic(nms); } private Nelder_Mead_Solution GetReflected(Nelder_Mead_Solution centroid) { Nelder_Mead_Solution nms = new Nelder_Mead_Solution(WorstSolution, -alpha, centroid, 1.0 + alpha); return NMS_Factory_Deterministic(nms); } private Nelder_Mead_Solution GetExpanded(Nelder_Mead_Solution reflected, Nelder_Mead_Solution centroid) { Nelder_Mead_Solution nms = new Nelder_Mead_Solution(reflected, gamma, centroid, 1.0 - gamma); return NMS_Factory_Deterministic(nms); } private Nelder_Mead_Solution GetContracted(Nelder_Mead_Solution centroid) { Nelder_Mead_Solution nms = new Nelder_Mead_Solution(WorstSolution, beta, centroid, 1.0 - beta); return NMS_Factory_Deterministic(nms); } private void ShrinkTowardsTheBestSolution() { for (int i = 1; i < Solutions.Length; i++) { Solutions[i].Append(BestSolution); Solutions[i].Scale(beta); NMS_Factory_Deterministic(Solutions[i]); } Array.Sort(Solutions); } private void ReplaceTheWorstWith(Nelder_Mead_Solution newSolution) { WorstSolution.ReplaceWith(newSolution); Array.Sort(Solutions); } private bool IsWorseThanAllButTheWorst(Nelder_Mead_Solution reflected) { for (int i = 0; i < Solutions.Length - 1; i++) { if (reflected.Value <= Solutions[i].Value) return false; } return true; } private Nelder_Mead_Solution NMS_Factory_Random() { double[] d = new double[Dimension]; for (int i = 0; i < Dimension; i++) d[i] = (maxX[i] - minX[i]) * m_random.NextDouble() + minX[i]; Nelder_Mead_Solution nms = new Nelder_Mead_Solution(d); nms.Value = ComputeTarget(nms.Vector); return nms; } private Nelder_Mead_Solution NMS_Factory_Deterministic(Nelder_Mead_Solution nms) { for (int i = 0; i < nms.Vector.Length; i++) { //if (nms.Vector[i] < minX[i]) // nms.Vector[i] = (maxX[i] - minX[i]) * 0.1 * m_random.NextDouble() + minX[i]; //if (nms.Vector[i] > maxX[i]) // nms.Vector[i] = maxX[i] - (maxX[i] - minX[i]) * 0.1 * m_random.NextDouble(); if (nms.Vector[i] < minX[i]) nms.Vector[i] = minX[i]; if (nms.Vector[i] > maxX[i]) nms.Vector[i] = maxX[i]; } nms.Value = ComputeTarget(nms.Vector); return nms; } private Nelder_Mead_Solution NMS_Factory_Deterministic(double[] vc) { Nelder_Mead_Solution nms = new Nelder_Mead_Solution(vc); nms.Value = ComputeTarget(nms.Vector); return nms; } #endregion } /// /// Describes a single Nelder-Mead solution as a point in n-dimensional space /// with the associated target function of value /// public class Nelder_Mead_Solution : IComparable { public double[] Vector; public double Value = 0.0; /// /// Constructor, defines the solution based on vector /// /// vector vector points public Nelder_Mead_Solution(double[] vector) { this.Vector = new double[vector.Length]; Array.Copy(vector, this.Vector, this.Vector.Length); } /// /// Constructor, defines the solution based on another solution /// /// other external solution public Nelder_Mead_Solution(Nelder_Mead_Solution other) { this.Vector = new double[other.Vector.Length]; Array.Copy(other.Vector, this.Vector, this.Vector.Length); } /// /// Constructor, defines the solution based two other solutions /// other1 * c1 + other2 * c2 /// public Nelder_Mead_Solution(Nelder_Mead_Solution other1, double c1, Nelder_Mead_Solution other2, double c2) { this.Vector = new double[other1.Vector.Length]; for (int i = 0; i < this.Vector.Length; i++) this.Vector[i] = other1.Vector[i] * c1 + other2.Vector[i] * c2; } /// /// Performs value comparison /// public int CompareTo(Nelder_Mead_Solution other) { if (this.Value < other.Value) return -1; if (this.Value > other.Value) return 1; return 0; } /// /// Appends solution vector from other /// public void Append(Nelder_Mead_Solution other) { for (int i = 0; i < Vector.Length; i++) Vector[i] += other.Vector[i]; } /// /// Scales solution vector by c /// public void Scale(double c) { for (int i = 0; i < Vector.Length; i++) Vector[i] *= c; } /// /// Replaces this to an external solution /// public void ReplaceWith(Nelder_Mead_Solution other) { Array.Copy(other.Vector, this.Vector, this.Vector.Length); this.Value = other.Value; } /// /// Converts the solution to string for debug /// public override string ToString() { StringBuilder sb = new StringBuilder("F( "); for (int i = 0; i < Vector.Length; i++) { if (Vector[i] >= 0.0) sb.Append(" "); sb.Append(Vector[i].ToString("0.000 ")); } sb.Append(Value.ToString(") = 0.000")); return sb.ToString(); }