Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature/allow mixture of bounded and unbounded fit parameters #932

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
155 changes: 6 additions & 149 deletions src/Numerics/Optimization/NonlinearMinimizerBase.cs
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
using MathNet.Numerics.LinearAlgebra;
using System;
using System.Collections.Generic;
using System.Linq;

namespace MathNet.Numerics.Optimization
Expand Down Expand Up @@ -58,20 +59,12 @@ protected void ValidateBounds(Vector<double> parameters, Vector<double> 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.");
Expand Down Expand Up @@ -161,150 +154,14 @@ protected double EvaluateFunction(IObjectiveModel objective, Vector<double> 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<double> ProjectToInternalParameters(Vector<double> 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<double> ProjectToExternalParameters(Vector<double> 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<double> parameterValues) =>
new ProjectableParameters(parameterValues, LowerBound, UpperBound, Scales);

return Pext;
}
protected Vector<double> ProjectToInternalParameters(Vector<double> Pext) => Projectable(Pext).ToInternal();

protected Vector<double> ScaleFactorsOfJacobian(Vector<double> Pint)
{
var scale = Vector<double>.Build.Dense(Pint.Count, 1.0);
protected Vector<double> ProjectToExternalParameters(Vector<double> 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<double> ScaleFactorsOfJacobian(Vector<double> Pint) => Projectable(Pint).JacobianScaleFactors();

#endregion Projection of Parameters
}
Expand Down
88 changes: 88 additions & 0 deletions src/Numerics/Optimization/ProjectableParameters.cs
Original file line number Diff line number Diff line change
@@ -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<ProjectableParameter> _parameters;

internal ProjectableParameters(IEnumerable<double> values, IList<double> lowerBounds, IList<double> upperBounds,
IList<double> scales)
{
_parameters = values.Select((value, index) =>
new ProjectableParameter(value, lowerBounds?[index], upperBounds?[index], scales?[index]));
}

internal Vector<double> ToInternal() =>
DenseVector.OfEnumerable(_parameters.Select(p => p.ToInternal()));

internal Vector<double> ToExternal() =>
DenseVector.OfEnumerable(_parameters.Select(p => p.ToExternal()));

internal Vector<double> 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;
}
}
}