|
|
|
@ -8,10 +8,30 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
{ |
|
|
|
internal class FittingObjectiveModel : IObjectiveModel |
|
|
|
{ |
|
|
|
#region Private Variables
|
|
|
|
|
|
|
|
readonly Func<Vector<double>, double, double> userFunction; // (p, x) => f(x; p)
|
|
|
|
readonly Func<Vector<double>, double, Vector<double>> userDerivatives; // (p, x) => df(x; p)/dp
|
|
|
|
readonly int accuracyOrder; // the desired accuracy order to evaluate the jacobian by numerical approximaiton.
|
|
|
|
|
|
|
|
Vector<double> coefficients; |
|
|
|
Vector<double> Pint; // internal(unbounded) coefficients
|
|
|
|
public Vector<double> Pext; // external(bounded) coefficients
|
|
|
|
|
|
|
|
bool hasFunctionValue; |
|
|
|
double functionValue; // the residual sum of squares. Residuals * Residuals
|
|
|
|
Vector<double> residuals; // the error values
|
|
|
|
|
|
|
|
bool hasJacobianValue; |
|
|
|
Matrix<double> jacobianValue; // the Jacobian matrix.
|
|
|
|
Vector<double> gradientValue; // the Gradient vector.
|
|
|
|
Matrix<double> hessianValue; // the Hessian matrix.
|
|
|
|
|
|
|
|
bool isBounded; |
|
|
|
|
|
|
|
#endregion Private Variables
|
|
|
|
|
|
|
|
#region Public Variables
|
|
|
|
#region Public Variables - Observed Data
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the independent variable.
|
|
|
|
@ -25,178 +45,183 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the weights for the observations.
|
|
|
|
/// inverse of the standard measurement errors
|
|
|
|
/// If null, unity weighting is used.
|
|
|
|
/// </summary>
|
|
|
|
public Matrix<double> Weights { get; private set; } |
|
|
|
// W = LL'
|
|
|
|
private Vector<double> L; |
|
|
|
private Vector<double> L; // Weights = LL'
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the parameters.
|
|
|
|
/// Get the number of observations.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> Parameters { get; private set; } |
|
|
|
public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public List<bool> IsFixed { get; set; } |
|
|
|
#endregion Public Variables - Observed Data
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> LowerBound { get; set; } |
|
|
|
#region Public Variables - Bounds of Parameter
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the values of the parameters.
|
|
|
|
/// Get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> UpperBound { get; set; } |
|
|
|
public List<bool> IsFixed { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the scale factor of the parameters.
|
|
|
|
/// Get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> Scales { get; set; } |
|
|
|
public Vector<double> LowerBound { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set of get whether or not the parameters are bounded.
|
|
|
|
/// Get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public bool IsBounded { get; set; } |
|
|
|
public Vector<double> UpperBound { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the y-values of the fitted model that correspond to the independent values.
|
|
|
|
/// Get the scale factor of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> Values { get; private set; } |
|
|
|
public Vector<double> Scales { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the error values, R(x; p) = L * (y - f(x; p)) where L = sqrt(W)
|
|
|
|
/// Get the number of unknown parameters.
|
|
|
|
/// </summary>
|
|
|
|
private Vector<double> Residuals; |
|
|
|
public int NumberOfParameters { get { return (Point == null) ? 0 : Point.Count; } } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the residual sum of squares, R.DotProduct(R)
|
|
|
|
/// </summary>
|
|
|
|
public double Residue { get; private set; } |
|
|
|
#endregion Public Variables - Bounds of Parameter
|
|
|
|
|
|
|
|
#region Public Variables - Others
|
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the Jacobian matrix of x and p, J(x; p).
|
|
|
|
/// Get the number of calls to function.
|
|
|
|
/// </summary>
|
|
|
|
public Matrix<double> Jacobian { get; private set; } |
|
|
|
|
|
|
|
public int FunctionEvaluations { get; set; } |
|
|
|
/// <summary>
|
|
|
|
/// Get the Gradient vector of x and p, J'WR
|
|
|
|
/// Get the number of calls to jacobian.
|
|
|
|
/// </summary>
|
|
|
|
public Vector<double> Gradient { get; private set; } |
|
|
|
public int JacobianEvaluations { get; set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the Hessian matrix of x and p, J'WJ
|
|
|
|
/// </summary>
|
|
|
|
public Matrix<double> Hessian { get; private set; } |
|
|
|
#endregion Public Variables - Others
|
|
|
|
|
|
|
|
public FittingObjectiveModel(Func<Vector<double>, double, double>function, Func<Vector<double>, double, Vector<double>> derivatives = null, int accuracyOrder = 2) |
|
|
|
{ |
|
|
|
this.userFunction = function; |
|
|
|
this.userDerivatives = derivatives; |
|
|
|
this.accuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder)); |
|
|
|
|
|
|
|
IsFinished = false; |
|
|
|
} |
|
|
|
|
|
|
|
public IObjectiveModel Fork() |
|
|
|
{ |
|
|
|
return new FittingObjectiveModel(userFunction, userDerivatives, accuracyOrder) |
|
|
|
{ |
|
|
|
ObservedX = ObservedX, |
|
|
|
ObservedY = ObservedY, |
|
|
|
Weights = Weights, |
|
|
|
|
|
|
|
coefficients = coefficients, |
|
|
|
Pint = Pint, |
|
|
|
Pext = Pext, |
|
|
|
|
|
|
|
hasFunctionValue = hasFunctionValue, |
|
|
|
functionValue = functionValue, |
|
|
|
|
|
|
|
hasJacobianValue = hasJacobianValue, |
|
|
|
jacobianValue = jacobianValue, |
|
|
|
gradientValue = gradientValue, |
|
|
|
hessianValue = hessianValue |
|
|
|
}; |
|
|
|
} |
|
|
|
|
|
|
|
public IObjectiveModel CreateNew() |
|
|
|
{ |
|
|
|
return new FittingObjectiveModel(userFunction, userDerivatives, accuracyOrder); |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the number of observations.
|
|
|
|
/// Set or get the values of the parameters.
|
|
|
|
/// </summary>
|
|
|
|
public int NumberOfObservations { get { return (ObservedY == null) ? 0 : ObservedY.Count; } } |
|
|
|
public Vector<double> Point { get { return coefficients; } } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the number of unknown parameters.
|
|
|
|
/// Get the y-values of the fitted model that correspond to the independent values.
|
|
|
|
/// </summary>
|
|
|
|
public int NumberOfParameters { get { return (Parameters == null) ? 0 : Parameters.Count; } } |
|
|
|
public Vector<double> ModelValues { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the degree of freedom
|
|
|
|
/// Get the residual sum of squares.
|
|
|
|
/// </summary>
|
|
|
|
public int DegreeOfFreedom |
|
|
|
public double Value |
|
|
|
{ |
|
|
|
get |
|
|
|
{ |
|
|
|
var dof = NumberOfObservations - NumberOfParameters; |
|
|
|
if (IsFixed != null) |
|
|
|
if (!hasFunctionValue) |
|
|
|
{ |
|
|
|
dof = dof + IsFixed.Count(p => p == true); |
|
|
|
EvaluateFunction(); |
|
|
|
hasFunctionValue = true; |
|
|
|
} |
|
|
|
return dof; |
|
|
|
return functionValue; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the covariance matrix.
|
|
|
|
/// </summary>
|
|
|
|
public Matrix<double> Covariance { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the correlation matrix.
|
|
|
|
/// </summary>
|
|
|
|
public Matrix<double> Correlation { get; private set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get the number of calls to function.
|
|
|
|
/// </summary>
|
|
|
|
public int FunctionEvaluations { get; set; } |
|
|
|
/// <summary>
|
|
|
|
/// Get the number of calls to jacobian.
|
|
|
|
/// </summary>
|
|
|
|
public int JacobianEvaluations { get; set; } |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set or get the desired accuracy order of the numerical jacobian.
|
|
|
|
/// Get the Gradient vector of x and p.
|
|
|
|
/// </summary>
|
|
|
|
public int AccuracyOrder { get; set; } |
|
|
|
public Vector<double> Gradient |
|
|
|
{ |
|
|
|
get |
|
|
|
{ |
|
|
|
if (!hasJacobianValue) |
|
|
|
{ |
|
|
|
EvaluateJacobian(); |
|
|
|
hasJacobianValue = true; |
|
|
|
} |
|
|
|
return gradientValue; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Get whether or not the analytical jacobian is supported.
|
|
|
|
/// Get the Hessian matrix of x and p, J'WJ
|
|
|
|
/// </summary>
|
|
|
|
public bool IsJacobianSupported { get { return userDerivatives != null; } } |
|
|
|
|
|
|
|
#endregion Public Variables
|
|
|
|
|
|
|
|
public FittingObjectiveModel(Func<Vector<double>, double, double>function, Func<Vector<double>, double, Vector<double>> derivatives, int accuracyOrder = 2) |
|
|
|
public Matrix<double> Hessian |
|
|
|
{ |
|
|
|
userFunction = function; |
|
|
|
userDerivatives = derivatives; |
|
|
|
AccuracyOrder = Math.Min(6, Math.Max(1, accuracyOrder)); |
|
|
|
get |
|
|
|
{ |
|
|
|
if (!hasJacobianValue) |
|
|
|
{ |
|
|
|
EvaluateJacobian(); |
|
|
|
hasJacobianValue = true; |
|
|
|
} |
|
|
|
return hessianValue; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
public IObjectiveModel Fork() |
|
|
|
/// <summary>
|
|
|
|
/// Get the degree of freedom
|
|
|
|
/// </summary>
|
|
|
|
public int DegreeOfFreedom |
|
|
|
{ |
|
|
|
return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder) |
|
|
|
get |
|
|
|
{ |
|
|
|
ObservedX = ObservedX, |
|
|
|
ObservedY = ObservedY, |
|
|
|
Weights = Weights, |
|
|
|
|
|
|
|
Parameters = Parameters, |
|
|
|
LowerBound = LowerBound, |
|
|
|
UpperBound = UpperBound, |
|
|
|
IsFixed = IsFixed, |
|
|
|
Scales = Scales, |
|
|
|
IsBounded = IsBounded, |
|
|
|
|
|
|
|
Residue = Residue, |
|
|
|
Jacobian = Jacobian |
|
|
|
}; |
|
|
|
var df = NumberOfObservations - NumberOfParameters; |
|
|
|
if (IsFixed != null) |
|
|
|
{ |
|
|
|
df = df + IsFixed.Count(p => p == true); |
|
|
|
} |
|
|
|
return df; |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
public IObjectiveModel CreateNew() |
|
|
|
{ |
|
|
|
return new FittingObjectiveModel(userFunction, userDerivatives, AccuracyOrder); |
|
|
|
} |
|
|
|
public bool IsGradientSupported { get { return true; } } |
|
|
|
public bool IsHessianSupported { get { return true; } } |
|
|
|
|
|
|
|
public bool IsFinished { get; set; } |
|
|
|
|
|
|
|
public IObjectiveFunction ToObjectiveFunction() |
|
|
|
{ |
|
|
|
Tuple<double, Vector<double>, Matrix<double>> function(Vector<double> point) |
|
|
|
{ |
|
|
|
EvaluateFunction(point); |
|
|
|
EvaluateJacobian(point); |
|
|
|
EvaluateAt(point); |
|
|
|
|
|
|
|
return new Tuple<double, Vector<double>, Matrix<double>>(Residue, Gradient, Hessian); |
|
|
|
return new Tuple<double, Vector<double>, Matrix<double>>(Value, Gradient, Hessian); |
|
|
|
} |
|
|
|
|
|
|
|
LowerBound = null; |
|
|
|
UpperBound = null; |
|
|
|
Scales = null; |
|
|
|
IsFixed = null; |
|
|
|
IsBounded = false; |
|
|
|
|
|
|
|
var objective = new GradientHessianObjectiveFunction(function); |
|
|
|
return objective; |
|
|
|
} |
|
|
|
@ -244,50 +269,37 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set observed data to fit.
|
|
|
|
/// </summary>
|
|
|
|
public void SetObserved(double[] observedX, double[] observedY, double[] weights = null) |
|
|
|
{ |
|
|
|
if (observedX == null || observedY == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("The data set can't be null."); |
|
|
|
} |
|
|
|
if (observedX.Length != observedY.Length) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The observed x data can't have different from observed y data."); |
|
|
|
} |
|
|
|
|
|
|
|
var wVector = (weights == null) |
|
|
|
? null |
|
|
|
: Vector<double>.Build.DenseOfArray(weights); |
|
|
|
SetObserved(Vector<double>.Build.DenseOfArray(observedX), Vector<double>.Build.DenseOfArray(observedY), wVector); |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set parameters.
|
|
|
|
/// <para/>
|
|
|
|
/// If bounded, the paramneters will be projected to unconstrained range by the mapping rule from the MINPACK.
|
|
|
|
/// If the projection is not needed, set IsBounded = false befre calling the Minimization method.
|
|
|
|
/// Set parameters and bounds.
|
|
|
|
/// </summary>
|
|
|
|
/// <param name="lowerBound">The lower bounds of parameters.</param>
|
|
|
|
/// <param name="upperBound">The upper bounds of parameters.</param>
|
|
|
|
/// /// <param name="scales">The scaling constants of parameters</param>
|
|
|
|
/// <param name="scales">The scaling constants of parameters</param>
|
|
|
|
/// <param name="isFixed">The list to the parameters fix or free.</param>
|
|
|
|
public void SetParameters(Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null) |
|
|
|
public void SetParameters(Vector<double> initialGuess, Vector<double> lowerBound = null, Vector<double> upperBound = null, Vector<double> scales = null, List<bool> isFixed = null) |
|
|
|
{ |
|
|
|
if (initialGuess == null) |
|
|
|
{ |
|
|
|
throw new ArgumentNullException("initialGuess"); |
|
|
|
} |
|
|
|
coefficients = initialGuess; |
|
|
|
|
|
|
|
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 != initialGuess.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); |
|
|
|
} |
|
|
|
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 && lowerBound != null && upperBound.Count != lowerBound.Count) |
|
|
|
if (upperBound != null && upperBound.Count != initialGuess.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the lower bounds."); |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); |
|
|
|
} |
|
|
|
UpperBound = upperBound; |
|
|
|
|
|
|
|
@ -295,13 +307,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
{ |
|
|
|
throw new ArgumentException("The scales must be finite."); |
|
|
|
} |
|
|
|
if (scales != null && lowerBound != null && scales.Count != lowerBound.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the lower bounds."); |
|
|
|
} |
|
|
|
if (scales != null && upperBound != null && scales.Count != upperBound.Count) |
|
|
|
if (scales != null && scales.Count != initialGuess.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the upper bounds."); |
|
|
|
throw new ArgumentException("The upper bounds can't have different elements from the initial guess."); |
|
|
|
} |
|
|
|
if (scales != null && scales.Count(x => x < 0) > 0) |
|
|
|
{ |
|
|
|
@ -309,48 +317,20 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
} |
|
|
|
Scales = scales; |
|
|
|
|
|
|
|
IsBounded = (LowerBound != null || UpperBound != null || Scales != null); |
|
|
|
|
|
|
|
if (isFixed != null && lowerBound != null && isFixed.Count != lowerBound.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The initial guess can't have different elements from the lower bounds."); |
|
|
|
} |
|
|
|
if (isFixed != null && upperBound != null && isFixed.Count != upperBound.Count) |
|
|
|
if (isFixed != null && isFixed.Count != initialGuess.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The initial guess can't have different elements from the upper bounds."); |
|
|
|
} |
|
|
|
if (isFixed != null && scales != null && isFixed.Count != scales.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("The initial guess can't have different elements from the scales."); |
|
|
|
throw new ArgumentException("The isFixed can't have different elements from the initial guess."); |
|
|
|
} |
|
|
|
if (isFixed != null && isFixed.Count(p => p == true) == isFixed.Count) |
|
|
|
{ |
|
|
|
throw new ArgumentException("All the parameters can't be fixed."); |
|
|
|
} |
|
|
|
IsFixed = isFixed; |
|
|
|
} |
|
|
|
|
|
|
|
/// <summary>
|
|
|
|
/// Set parameters.
|
|
|
|
/// <para/>
|
|
|
|
/// If bounded, the paramneters will be projected to unconstrained range by the mapping rule.
|
|
|
|
/// If the projection is not needed, set IsBounded = false befre calling the Minimization method.
|
|
|
|
/// </summary>
|
|
|
|
/// <param name="lowerBound">The lower bounds of parameters.</param>
|
|
|
|
/// <param name="upperBound">The upper bounds of parameters.</param>
|
|
|
|
/// <param name="scales">The scaling constants of parameters</param>
|
|
|
|
/// <param name="isFixed">The list to the parameters fix or free.</param>
|
|
|
|
public void SetParameters(double[] lowerBound = null, double[] upperBound = null, double[] scales = null, bool[] isFixed = null) |
|
|
|
{ |
|
|
|
var lb = (lowerBound == null) ? null : Vector<double>.Build.DenseOfArray(lowerBound); |
|
|
|
var ub = (upperBound == null) ? null : Vector<double>.Build.DenseOfArray(upperBound); |
|
|
|
var sc = (scales == null) ? null : Vector<double>.Build.DenseOfArray(scales); |
|
|
|
var fp = (isFixed == null) ? null : isFixed.ToList(); |
|
|
|
|
|
|
|
SetParameters(lb, ub, sc, fp); |
|
|
|
isBounded = LowerBound != null || UpperBound != null || Scales != null; |
|
|
|
} |
|
|
|
|
|
|
|
public void EvaluateFunction(Vector<double> parameters) |
|
|
|
public void EvaluateAt(Vector<double> parameters) |
|
|
|
{ |
|
|
|
ValidateParameters(parameters); |
|
|
|
|
|
|
|
@ -386,69 +366,84 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
//
|
|
|
|
// 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.
|
|
|
|
var Pext = (FunctionEvaluations > 0 && this.IsBounded) |
|
|
|
? ProjectParametersToExternal(parameters) |
|
|
|
: parameters.Clone(); |
|
|
|
|
|
|
|
// Project parameters, now this.Parameters are the internal parameters.
|
|
|
|
Parameters = (this.IsBounded) |
|
|
|
Pext = (FunctionEvaluations > 0 && isBounded) |
|
|
|
? ProjectParametersToExternal(parameters) |
|
|
|
: parameters.Clone(); |
|
|
|
Pint = (isBounded) |
|
|
|
? ProjectParametersToInternal(Pext) |
|
|
|
: Pext; |
|
|
|
|
|
|
|
this.coefficients = Pint; |
|
|
|
|
|
|
|
if (IsFinished) |
|
|
|
{ |
|
|
|
this.coefficients = Pext; |
|
|
|
} |
|
|
|
|
|
|
|
hasFunctionValue = false; |
|
|
|
hasJacobianValue = false; |
|
|
|
|
|
|
|
// don't keep references unnecessarily
|
|
|
|
jacobianValue = null; |
|
|
|
gradientValue = null; |
|
|
|
hessianValue = null; |
|
|
|
} |
|
|
|
|
|
|
|
#region Private Methods
|
|
|
|
|
|
|
|
private void EvaluateFunction() |
|
|
|
{ |
|
|
|
// Calculates the residuals, (y[i] - f(x[i]; p)) * L[i]
|
|
|
|
if (Values == null) |
|
|
|
if (ModelValues == null) |
|
|
|
{ |
|
|
|
Values = Vector<double>.Build.Dense(NumberOfObservations); |
|
|
|
ModelValues = Vector<double>.Build.Dense(NumberOfObservations); |
|
|
|
} |
|
|
|
for (int i = 0; i < NumberOfObservations; i++) |
|
|
|
{ |
|
|
|
Values[i] = userFunction(Pext, ObservedX[i]); |
|
|
|
ModelValues[i] = userFunction(Pext, ObservedX[i]); |
|
|
|
} |
|
|
|
FunctionEvaluations++; |
|
|
|
|
|
|
|
// calculate the weighted residuals
|
|
|
|
Residuals = (Weights == null) |
|
|
|
? ObservedY - Values |
|
|
|
: (ObservedY - Values).PointwiseMultiply(L); |
|
|
|
residuals = (Weights == null) |
|
|
|
? ObservedY - ModelValues |
|
|
|
: (ObservedY - ModelValues).PointwiseMultiply(L); |
|
|
|
|
|
|
|
// Calculate the residual sum of squares
|
|
|
|
Residue = Residuals.DotProduct(Residuals); |
|
|
|
functionValue = residuals.DotProduct(residuals); |
|
|
|
|
|
|
|
return; |
|
|
|
} |
|
|
|
|
|
|
|
public void EvaluateJacobian(Vector<double> parameters) |
|
|
|
private void EvaluateJacobian() |
|
|
|
{ |
|
|
|
var Pext = (IsBounded) |
|
|
|
? ProjectParametersToExternal(parameters) |
|
|
|
: parameters.Clone(); |
|
|
|
|
|
|
|
// Calculates the jacobian of x and p.
|
|
|
|
if (userDerivatives != null) |
|
|
|
{ |
|
|
|
// analytical jacobian
|
|
|
|
if (Jacobian == null) |
|
|
|
if (jacobianValue == null) |
|
|
|
{ |
|
|
|
Jacobian = Matrix<double>.Build.Dense(NumberOfObservations, NumberOfParameters); |
|
|
|
jacobianValue = Matrix<double>.Build.Dense(NumberOfObservations, NumberOfParameters); |
|
|
|
} |
|
|
|
for (int i = 0; i < NumberOfObservations; i++) |
|
|
|
{ |
|
|
|
Jacobian.SetRow(i, userDerivatives(Pext, ObservedX[i])); |
|
|
|
jacobianValue.SetRow(i, userDerivatives(Pext, ObservedX[i])); |
|
|
|
} |
|
|
|
JacobianEvaluations++; |
|
|
|
} |
|
|
|
else |
|
|
|
{ |
|
|
|
// numerical jacobian
|
|
|
|
Jacobian = NumericalJacobian(Pext, Values, AccuracyOrder); |
|
|
|
FunctionEvaluations += AccuracyOrder; |
|
|
|
jacobianValue = NumericalJacobian(Pext, ModelValues, accuracyOrder); |
|
|
|
FunctionEvaluations += accuracyOrder; |
|
|
|
} |
|
|
|
|
|
|
|
var scaleFactors = (this.IsBounded) |
|
|
|
? ScaleFactorsOfJacobian(Parameters) |
|
|
|
: Vector<double>.Build.Dense(Parameters.Count, 1.0); |
|
|
|
var scaleFactors = (isBounded && !IsFinished) |
|
|
|
? ScaleFactorsOfJacobian(Pint) |
|
|
|
: Vector<double>.Build.Dense(Pint.Count, 1.0); |
|
|
|
|
|
|
|
// Jint(x; Pint) = Jext(x; Pext) * scale where scale = dPext/dPint
|
|
|
|
// project jacobian: Jint(x; Pint) = Jext(x; Pext) * scale where scale = dPext/dPint
|
|
|
|
for (int i = 0; i < NumberOfObservations; i++) |
|
|
|
{ |
|
|
|
for (int j = 0; j < NumberOfParameters; j++) |
|
|
|
@ -456,58 +451,24 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
if (IsFixed != null && IsFixed[j]) |
|
|
|
{ |
|
|
|
// if j-th parameter is fixed, set J[i, j] = 0
|
|
|
|
Jacobian[i, j] = 0.0; |
|
|
|
jacobianValue[i, j] = 0.0; |
|
|
|
} |
|
|
|
else |
|
|
|
{ |
|
|
|
Jacobian[i, j] = Jacobian[i, j] * scaleFactors[j]; |
|
|
|
jacobianValue[i, j] = jacobianValue[i, j] * scaleFactors[j]; |
|
|
|
} |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
// Gradient, g = -J'W(y − f(x; p)) = -J'L(L'E) = -J'LR
|
|
|
|
Gradient = (Weights == null) |
|
|
|
? -Jacobian.Transpose() * (ObservedY - Values) |
|
|
|
: -Jacobian.Transpose() * Weights * (ObservedY - Values); |
|
|
|
gradientValue = (Weights == null) |
|
|
|
? -jacobianValue.Transpose() * (ObservedY - ModelValues) |
|
|
|
: -jacobianValue.Transpose() * Weights * (ObservedY - ModelValues); |
|
|
|
|
|
|
|
// approximated Hessian, H = J'WJ + ∑LRiHi ~ J'WJ near the minimum
|
|
|
|
Hessian = (Weights == null) |
|
|
|
? Jacobian.Transpose() * Jacobian |
|
|
|
: Jacobian.Transpose() * Weights * Jacobian; |
|
|
|
} |
|
|
|
|
|
|
|
public void EvaluateCovariance(Vector<double> parameters) |
|
|
|
{ |
|
|
|
// convert to bounded(external) parameters
|
|
|
|
var Pext = (IsBounded) |
|
|
|
? ProjectParametersToExternal(parameters) |
|
|
|
: parameters.Clone(); |
|
|
|
|
|
|
|
// set IsBounded = false to get external Parameters and covariance matrix
|
|
|
|
this.IsBounded = false; |
|
|
|
|
|
|
|
EvaluateFunction(Pext); |
|
|
|
EvaluateJacobian(Pext); |
|
|
|
|
|
|
|
// restore isBounded
|
|
|
|
this.IsBounded = (LowerBound != null || UpperBound != null); |
|
|
|
|
|
|
|
if (Hessian == null || Residuals == null || DegreeOfFreedom < 1) |
|
|
|
{ |
|
|
|
Covariance = null; |
|
|
|
Correlation = null; |
|
|
|
return; |
|
|
|
} |
|
|
|
|
|
|
|
var covariance = Hessian.PseudoInverse() * Residuals.DotProduct(Residuals) / DegreeOfFreedom; |
|
|
|
Covariance = covariance; |
|
|
|
|
|
|
|
var correlation = covariance.Clone(); |
|
|
|
var d = correlation.Diagonal().PointwiseSqrt(); |
|
|
|
var dd = d.OuterProduct(d); |
|
|
|
Correlation = correlation.PointwiseDivide(dd); |
|
|
|
|
|
|
|
return; |
|
|
|
hessianValue = (Weights == null) |
|
|
|
? jacobianValue.Transpose() * jacobianValue |
|
|
|
: jacobianValue.Transpose() * Weights * jacobianValue; |
|
|
|
} |
|
|
|
|
|
|
|
private void ValidateParameters(Vector<double> parameters) |
|
|
|
@ -538,16 +499,13 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
} |
|
|
|
} |
|
|
|
|
|
|
|
#region Numerical Derivatives
|
|
|
|
|
|
|
|
// Numerical derivatives by using the central or forward finite difference
|
|
|
|
private Matrix<double> NumericalJacobian(Vector<double> parameters, Vector<double> currentValues, int accuracyOrder = 2) |
|
|
|
private Matrix<double> NumericalJacobian(Vector<double> Pext, Vector<double> currentValues, int accuracyOrder = 2) |
|
|
|
{ |
|
|
|
const double sqrtEpsilon = 1.4901161193847656250E-8; // sqrt(machineEpsilon)
|
|
|
|
|
|
|
|
Matrix<double> derivertives = Matrix<double>.Build.Dense(NumberOfObservations, NumberOfParameters); |
|
|
|
|
|
|
|
var d = 0.000003 * parameters.PointwiseAbs().PointwiseMaximum(sqrtEpsilon); |
|
|
|
var d = 0.000003 * Pext.PointwiseAbs().PointwiseMaximum(sqrtEpsilon); |
|
|
|
|
|
|
|
var h = Vector<double>.Build.Dense(NumberOfParameters); |
|
|
|
for (int i = 0; i < NumberOfObservations; i++) |
|
|
|
@ -560,12 +518,12 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
if (accuracyOrder >= 6) |
|
|
|
{ |
|
|
|
// f'(x) = {- f(x - 3h) + 9f(x - 2h) - 45f(x - h) + 45f(x + h) - 9f(x + 2h) + f(x + 3h)} / 60h + O(h^6)
|
|
|
|
var f1 = userFunction(parameters - 3 * h, x); |
|
|
|
var f2 = userFunction(parameters - 2 * h, x); |
|
|
|
var f3 = userFunction(parameters - h, x); |
|
|
|
var f4 = userFunction(parameters + h, x); |
|
|
|
var f5 = userFunction(parameters + 2 * h, x); |
|
|
|
var f6 = userFunction(parameters + 3 * h, x); |
|
|
|
var f1 = userFunction(Pext - 3 * h, x); |
|
|
|
var f2 = userFunction(Pext - 2 * h, x); |
|
|
|
var f3 = userFunction(Pext - h, x); |
|
|
|
var f4 = userFunction(Pext + h, x); |
|
|
|
var f5 = userFunction(Pext + 2 * h, x); |
|
|
|
var f6 = userFunction(Pext + 3 * h, x); |
|
|
|
|
|
|
|
var prime = (-f1 + 9 * f2 - 45 * f3 + 45 * f4 - 9 * f5 + f6) / (60 * h[j]); |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -574,11 +532,11 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
{ |
|
|
|
// f'(x) = {-137f(x) + 300f(x + h) - 300f(x + 2h) + 200f(x + 3h) - 75f(x + 4h) + 12f(x + 5h)} / 60h + O(h^5)
|
|
|
|
var f1 = currentValues[i]; |
|
|
|
var f2 = userFunction(parameters + h, x); |
|
|
|
var f3 = userFunction(parameters + 2 * h, x); |
|
|
|
var f4 = userFunction(parameters + 3 * h, x); |
|
|
|
var f5 = userFunction(parameters + 4 * h, x); |
|
|
|
var f6 = userFunction(parameters + 5 * h, x); |
|
|
|
var f2 = userFunction(Pext + h, x); |
|
|
|
var f3 = userFunction(Pext + 2 * h, x); |
|
|
|
var f4 = userFunction(Pext + 3 * h, x); |
|
|
|
var f5 = userFunction(Pext + 4 * h, x); |
|
|
|
var f6 = userFunction(Pext + 5 * h, x); |
|
|
|
|
|
|
|
var prime = (-137 * f1 + 300 * f2 - 300 * f3 + 200 * f4 - 75 * f5 + 12 * f6) / (60 * h[j]); |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -586,10 +544,10 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
else if (accuracyOrder == 4) |
|
|
|
{ |
|
|
|
// f'(x) = {f(x - 2h) - 8f(x - h) + 8f(x + h) - f(x + 2h)} / 12h + O(h^4)
|
|
|
|
var f1 = userFunction(parameters - 2 * h, x); |
|
|
|
var f2 = userFunction(parameters - h, x); |
|
|
|
var f3 = userFunction(parameters + h, x); |
|
|
|
var f4 = userFunction(parameters + 2 * h, x); |
|
|
|
var f1 = userFunction(Pext - 2 * h, x); |
|
|
|
var f2 = userFunction(Pext - h, x); |
|
|
|
var f3 = userFunction(Pext + h, x); |
|
|
|
var f4 = userFunction(Pext + 2 * h, x); |
|
|
|
|
|
|
|
var prime = (f1 - 8 * f2 + 8 * f3 - f4) / (12 * h[j]); |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -598,9 +556,9 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
{ |
|
|
|
// f'(x) = {-11f(x) + 18f(x + h) - 9f(x + 2h) + 2f(x + 3h)} / 6h + O(h^3)
|
|
|
|
var f1 = currentValues[i]; |
|
|
|
var f2 = userFunction(parameters + h, x); |
|
|
|
var f3 = userFunction(parameters + 2 * h, x); |
|
|
|
var f4 = userFunction(parameters + 3 * h, x); |
|
|
|
var f2 = userFunction(Pext + h, x); |
|
|
|
var f3 = userFunction(Pext + 2 * h, x); |
|
|
|
var f4 = userFunction(Pext + 3 * h, x); |
|
|
|
|
|
|
|
var prime = (-11 * f1 + 18 * f2 - 9 * f3 + 2 * f4) / (6 * h[j]); |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -608,8 +566,8 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
else if (accuracyOrder == 2) |
|
|
|
{ |
|
|
|
// f'(x) = {f(x + h) - f(x - h)} / 2h + O(h^2)
|
|
|
|
var f1 = userFunction(parameters + h, x); |
|
|
|
var f2 = userFunction(parameters - h, x); |
|
|
|
var f1 = userFunction(Pext + h, x); |
|
|
|
var f2 = userFunction(Pext - h, x); |
|
|
|
|
|
|
|
var prime = (f1 - f2) / (2 * h[j]); |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -618,7 +576,7 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
{ |
|
|
|
// f'(x) = {- f(x) + f(x + h)} / h + O(h)
|
|
|
|
var f1 = currentValues[i]; |
|
|
|
var f2 = userFunction(parameters + h, x); |
|
|
|
var f2 = userFunction(Pext + h, x); |
|
|
|
|
|
|
|
var prime = (-f1 + f2) / h[j]; |
|
|
|
derivertives[i, j] = prime; |
|
|
|
@ -631,10 +589,6 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
return derivertives; |
|
|
|
} |
|
|
|
|
|
|
|
#endregion Numerical Derivatives
|
|
|
|
|
|
|
|
#region Projection
|
|
|
|
|
|
|
|
private Vector<double> ProjectParametersToInternal(Vector<double> Pext) |
|
|
|
{ |
|
|
|
var Pint = Pext.Clone(); |
|
|
|
@ -771,6 +725,6 @@ namespace MathNet.Numerics.Optimization.ObjectiveModels |
|
|
|
return scale; |
|
|
|
} |
|
|
|
|
|
|
|
#endregion Projection
|
|
|
|
#endregion Private Methods
|
|
|
|
} |
|
|
|
} |
|
|
|
|
diluculo
8 years ago