diff --git a/src/Numerics/Optimization/NonlinearMinimizerBase.cs b/src/Numerics/Optimization/NonlinearMinimizerBase.cs index 1e06410a4..cddc2c914 100644 --- a/src/Numerics/Optimization/NonlinearMinimizerBase.cs +++ b/src/Numerics/Optimization/NonlinearMinimizerBase.cs @@ -1,5 +1,6 @@ using MathNet.Numerics.LinearAlgebra; using System; +using System.Collections.Generic; using System.Linq; namespace MathNet.Numerics.Optimization @@ -58,20 +59,12 @@ protected void ValidateBounds(Vector parameters, Vector lowerBou throw new ArgumentNullException(nameof(parameters)); } - if (lowerBound != null && lowerBound.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0) - { - throw new ArgumentException("The lower bounds must be finite."); - } if (lowerBound != null && lowerBound.Count != parameters.Count) { throw new ArgumentException("The lower bounds can't have different size from the parameters."); } LowerBound = lowerBound; - if (upperBound != null && upperBound.Count(x => double.IsInfinity(x) || double.IsNaN(x)) > 0) - { - throw new ArgumentException("The upper bounds must be finite."); - } if (upperBound != null && upperBound.Count != parameters.Count) { throw new ArgumentException("The upper bounds can't have different size from the parameters."); @@ -161,150 +154,14 @@ protected double EvaluateFunction(IObjectiveModel objective, Vector Pint // Except when it is initial guess, the parameters argument is always internal parameter. // So, first map the parameters argument to the external parameters in order to calculate function values. - protected Vector ProjectToInternalParameters(Vector Pext) - { - var Pint = Pext.Clone(); - - if (LowerBound != null && UpperBound != null) - { - for (int i = 0; i < Pext.Count; i++) - { - Pint[i] = Math.Asin((2.0 * (Pext[i] - LowerBound[i]) / (UpperBound[i] - LowerBound[i])) - 1.0); - } - - return Pint; - } - - if (LowerBound != null && UpperBound == null) - { - for (int i = 0; i < Pext.Count; i++) - { - Pint[i] = (Scales == null) - ? Math.Sqrt(Math.Pow(Pext[i] - LowerBound[i] + 1.0, 2) - 1.0) - : Math.Sqrt(Math.Pow((Pext[i] - LowerBound[i]) / Scales[i] + 1.0, 2) - 1.0); - } - - return Pint; - } - - if (LowerBound == null && UpperBound != null) - { - for (int i = 0; i < Pext.Count; i++) - { - Pint[i] = (Scales == null) - ? Math.Sqrt(Math.Pow(UpperBound[i] - Pext[i] + 1.0, 2) - 1.0) - : Math.Sqrt(Math.Pow((UpperBound[i] - Pext[i]) / Scales[i] + 1.0, 2) - 1.0); - } - - return Pint; - } - - if (Scales != null) - { - for (int i = 0; i < Pext.Count; i++) - { - Pint[i] = Pext[i] / Scales[i]; - } - - return Pint; - } - - return Pint; - } - - protected Vector ProjectToExternalParameters(Vector Pint) - { - var Pext = Pint.Clone(); - - if (LowerBound != null && UpperBound != null) - { - for (int i = 0; i < Pint.Count; i++) - { - Pext[i] = LowerBound[i] + (UpperBound[i] / 2.0 - LowerBound[i] / 2.0) * (Math.Sin(Pint[i]) + 1.0); - } - - return Pext; - } - - if (LowerBound != null && UpperBound == null) - { - for (int i = 0; i < Pint.Count; i++) - { - Pext[i] = (Scales == null) - ? LowerBound[i] + Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0 - : LowerBound[i] + Scales[i] * (Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0); - } - - return Pext; - } - - if (LowerBound == null && UpperBound != null) - { - for (int i = 0; i < Pint.Count; i++) - { - Pext[i] = (Scales == null) - ? UpperBound[i] - Math.Sqrt(Pint[i] * Pint[i] + 1.0) + 1.0 - : UpperBound[i] - Scales[i] * (Math.Sqrt(Pint[i] * Pint[i] + 1.0) - 1.0); - } - - return Pext; - } - - if (Scales != null) - { - for (int i = 0; i < Pint.Count; i++) - { - Pext[i] = Pint[i] * Scales[i]; - } - - return Pext; - } + private ProjectableParameters Projectable(IEnumerable parameterValues) => + new ProjectableParameters(parameterValues, LowerBound, UpperBound, Scales); - return Pext; - } + protected Vector ProjectToInternalParameters(Vector Pext) => Projectable(Pext).ToInternal(); - protected Vector ScaleFactorsOfJacobian(Vector Pint) - { - var scale = Vector.Build.Dense(Pint.Count, 1.0); + protected Vector ProjectToExternalParameters(Vector Pint) => Projectable(Pint).ToExternal(); - if (LowerBound != null && UpperBound != null) - { - for (int i = 0; i < Pint.Count; i++) - { - scale[i] = (UpperBound[i] - LowerBound[i]) / 2.0 * Math.Cos(Pint[i]); - } - return scale; - } - - if (LowerBound != null && UpperBound == null) - { - for (int i = 0; i < Pint.Count; i++) - { - scale[i] = (Scales == null) - ? Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0) - : Scales[i] * Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0); - } - return scale; - } - - if (LowerBound == null && UpperBound != null) - { - for (int i = 0; i < Pint.Count; i++) - { - scale[i] = (Scales == null) - ? -Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0) - : -Scales[i] * Pint[i] / Math.Sqrt(Pint[i] * Pint[i] + 1.0); - } - return scale; - } - - if (Scales != null) - { - return Scales; - } - - return scale; - } + protected Vector ScaleFactorsOfJacobian(Vector Pint) => Projectable(Pint).JacobianScaleFactors(); #endregion Projection of Parameters } diff --git a/src/Numerics/Optimization/ProjectableParameters.cs b/src/Numerics/Optimization/ProjectableParameters.cs new file mode 100644 index 000000000..6383afc9d --- /dev/null +++ b/src/Numerics/Optimization/ProjectableParameters.cs @@ -0,0 +1,88 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Double; + +namespace MathNet.Numerics.Optimization +{ + internal class ProjectableParameters + { + private readonly IEnumerable _parameters; + + internal ProjectableParameters(IEnumerable values, IList lowerBounds, IList upperBounds, + IList scales) + { + _parameters = values.Select((value, index) => + new ProjectableParameter(value, lowerBounds?[index], upperBounds?[index], scales?[index])); + } + + internal Vector ToInternal() => + DenseVector.OfEnumerable(_parameters.Select(p => p.ToInternal())); + + internal Vector ToExternal() => + DenseVector.OfEnumerable(_parameters.Select(p => p.ToExternal())); + + internal Vector JacobianScaleFactors() => + DenseVector.OfEnumerable(_parameters.Select(p => p.JacobianScaleFactor())); + } + + + internal class ProjectableParameter + { + private readonly double _value; + private readonly double _lowerBound; + private readonly double _upperBound; + private readonly double _scaleFactor; + + internal ProjectableParameter(double value, double? lowerBound, double? upperBound, double? scaleFactor) + { + _value = value; + _lowerBound = lowerBound ?? double.NegativeInfinity; + _upperBound = upperBound ?? double.PositiveInfinity; + _scaleFactor = scaleFactor ?? 1; + } + + internal double ToInternal() + { + if (_lowerBound.IsFinite() && _upperBound.IsFinite()) + return Math.Asin(2 * (_value - _lowerBound) / (_upperBound - _lowerBound) - 1); + + if (_lowerBound.IsFinite()) + return Math.Sqrt(Math.Pow((_value - _lowerBound) / _scaleFactor + 1, 2) - 1); + + if (_upperBound.IsFinite()) + return Math.Sqrt(Math.Pow((_upperBound - _value) / _scaleFactor + 1, 2) - 1); + + return _value / _scaleFactor; + } + + internal double ToExternal() + { + if (_lowerBound.IsFinite() && _upperBound.IsFinite()) + return _lowerBound + (_upperBound / 2 - _lowerBound / 2) * (Math.Sin(_value) + 1); + + if (_lowerBound.IsFinite()) + return _lowerBound + _scaleFactor * (Math.Sqrt(_value * _value + 1) - 1); + + if (_upperBound.IsFinite()) + return _upperBound - _scaleFactor * (Math.Sqrt(_value * _value + 1) - 1); + + return _value * _scaleFactor; + } + + internal double JacobianScaleFactor() + { + if (_lowerBound.IsFinite() && _upperBound.IsFinite()) + return (_upperBound - _lowerBound) / 2 * Math.Cos(_value); + + if (_lowerBound.IsFinite()) + return _scaleFactor * _value / Math.Sqrt(_value * _value + 1); + + if (_upperBound.IsFinite()) + return -_scaleFactor * _value / Math.Sqrt(_value * _value + 1); + + return _scaleFactor; + } + } +}