//
// Math.NET Numerics, part of the Math.NET Project
// http://numerics.mathdotnet.com
// http://github.com/mathnet/mathnet-numerics
//
// Copyright (c) 2009-2014 Math.NET
//
// Permission is hereby granted, free of charge, to any person
// obtaining a copy of this software and associated documentation
// files (the "Software"), to deal in the Software without
// restriction, including without limitation the rights to use,
// copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following
// conditions:
//
// The above copyright notice and this permission notice shall be
// included in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
// OTHER DEALINGS IN THE SOFTWARE.
//
using System;
using System.Diagnostics;
namespace IStation.Numerics.LinearAlgebra.Solvers
{
///
/// Defines an that monitors residuals as stop criterion.
///
public sealed class ResidualStopCriterion : IIterationStopCriterion where T : struct, IEquatable, IFormattable
{
///
/// The maximum value for the residual below which the calculation is considered converged.
///
double _maximum;
///
/// The minimum number of iterations for which the residual has to be below the maximum before
/// the calculation is considered converged.
///
int _minimumIterationsBelowMaximum;
///
/// The status of the calculation
///
IterationStatus _status = IterationStatus.Continue;
///
/// The number of iterations since the residuals got below the maximum.
///
int _iterationCount;
///
/// The iteration number of the last iteration.
///
int _lastIteration = -1;
///
/// Initializes a new instance of the class with the specified
/// maximum residual and minimum number of iterations.
///
///
/// The maximum value for the residual below which the calculation is considered converged.
///
///
/// The minimum number of iterations for which the residual has to be below the maximum before
/// the calculation is considered converged.
///
public ResidualStopCriterion(double maximum, int minimumIterationsBelowMaximum = 0)
{
if (maximum < 0)
{
throw new ArgumentOutOfRangeException(nameof(maximum));
}
if (minimumIterationsBelowMaximum < 0)
{
throw new ArgumentOutOfRangeException(nameof(minimumIterationsBelowMaximum));
}
_maximum = maximum;
_minimumIterationsBelowMaximum = minimumIterationsBelowMaximum;
}
///
/// Gets or sets the maximum value for the residual below which the calculation is considered
/// converged.
///
/// Thrown if the Maximum is set to a negative value.
public double Maximum
{
[DebuggerStepThrough]
get => _maximum;
[DebuggerStepThrough]
set
{
if (value < 0)
{
throw new ArgumentOutOfRangeException(nameof(value));
}
_maximum = value;
}
}
///
/// Gets or sets the minimum number of iterations for which the residual has to be
/// below the maximum before the calculation is considered converged.
///
/// Thrown if the BelowMaximumFor is set to a value less than 1.
public int MinimumIterationsBelowMaximum
{
[DebuggerStepThrough]
get => _minimumIterationsBelowMaximum;
[DebuggerStepThrough]
set
{
if (value < 0)
{
throw new ArgumentOutOfRangeException(nameof(value));
}
_minimumIterationsBelowMaximum = value;
}
}
///
/// Determines the status of the iterative calculation based on the stop criteria stored
/// by the current . Result is set into Status field.
///
/// The number of iterations that have passed so far.
/// The vector containing the current solution values.
/// The right hand side vector.
/// The vector containing the current residual vectors.
///
/// The individual stop criteria may internally track the progress of the calculation based
/// on the invocation of this method. Therefore this method should only be called if the
/// calculation has moved forwards at least one step.
///
public IterationStatus DetermineStatus(int iterationNumber, Vector solutionVector, Vector sourceVector, Vector residualVector)
{
if (iterationNumber < 0)
{
throw new ArgumentOutOfRangeException(nameof(iterationNumber));
}
if (solutionVector.Count != sourceVector.Count)
{
throw new ArgumentException("All vectors must have the same dimensionality.", nameof(sourceVector));
}
if (solutionVector.Count != residualVector.Count)
{
throw new ArgumentException("All vectors must have the same dimensionality.", nameof(residualVector));
}
// Store the infinity norms of both the solution and residual vectors
// These values will be used to calculate the relative drop in residuals
// later on.
var residualNorm = residualVector.InfinityNorm();
// This is criterion 1 from Templates for the solution of linear systems.
// The problem with this criterion is that it's not limiting enough. For now
// we won't use it. Later on we might get back to it.
// return mMaximumResidual * (System.Math.Abs(mMatrixNorm) * System.Math.Abs(solutionNorm) + System.Math.Abs(mVectorNorm));
// For now use criterion 2 from Templates for the solution of linear systems. See page 60.
// Check the residuals by calculating:
// ||r_i|| <= stop_tol * ||b||
var stopCriterion = _maximum*sourceVector.InfinityNorm();
// First check that we have real numbers not NaN's.
// NaN's can occur when the iterative process diverges so we
// stop if that is the case.
if (double.IsNaN(stopCriterion) || double.IsNaN(residualNorm))
{
_iterationCount = 0;
_status = IterationStatus.Diverged;
return _status;
}
// ||r_i|| <= stop_tol * ||b||
// Stop the calculation if it's clearly smaller than the tolerance
if (residualNorm < stopCriterion)
{
if (_lastIteration <= iterationNumber)
{
_iterationCount = iterationNumber - _lastIteration;
_status = _iterationCount >= _minimumIterationsBelowMaximum ? IterationStatus.Converged : IterationStatus.Continue;
}
}
else
{
_iterationCount = 0;
_status = IterationStatus.Continue;
}
_lastIteration = iterationNumber;
return _status;
}
///
/// Gets the current calculation status.
///
public IterationStatus Status
{
[DebuggerStepThrough]
get => _status;
}
///
/// Resets the to the pre-calculation state.
///
public void Reset()
{
_status = IterationStatus.Continue;
_iterationCount = 0;
_lastIteration = -1;
}
///
/// Clones the current and its settings.
///
/// A new instance of the class.
public IIterationStopCriterion Clone()
{
return new ResidualStopCriterion(_maximum, _minimumIterationsBelowMaximum);
}
}
}