by Trent Guidry
2. October 2009 06:03
In this post I develop the code for the cubic region elimination class.
This algorithm uses a working set of four points. For each iteration of the cubic region elimination, a third order polynomial (a cubic) is fitted to the four points and a new X value is estimated using the algorithm posted previously at Fitting a cubic to four points and finding the local minimum in C#. If the new X value is not already in the set of points, the function value of the new point is computed, added to the set, and the point with the largest Y value is removed.
The source code for the cubic region elimination class is given below. The code for the region elimination base class was posted previously.
public class CubicRegionElimination : RegionEliminationBase
{
public CubicRegionElimination()
:
base()
{
_numStartingPoints = 4;
}
public override RegionEliminationResults RegionEliminate(Func<double, double> function, SearchPointCollection startPoint, double eliminationFraction)
{
Debug.Assert(startPoint.Count >= 2);
bool succeeded = true;
Cancel = false;
SearchPointCollection searchPoints = startPoint.Clone(4);
CalculateMissingValues(searchPoints, function);
while (searchPoints.Count < 4)
{
double midPoint = (searchPoints[0].X + searchPoints[1].X) / 2.0;
double functionValue = function(midPoint);
searchPoints.Add(new SearchPoint(midPoint, functionValue));
}
double stopRegionSize = (searchPoints[3].X - searchPoints[0].X) * eliminationFraction;
while (!Cancel)
{
try
{
double newX = searchPoints.CubicMinYEstimate();
if (!double.IsNaN(newX))
{
double functionValue = function(newX);
SearchPoint newPoint = new SearchPoint(newX, functionValue);
if (searchPoints.CanAdd(newPoint))
{
searchPoints.Add(newPoint);
}
else
{
Cancel = true;
}
}
else
{
Cancel = true;
succeeded = false;
}
if (searchPoints[3].X - searchPoints[0].X < stopRegionSize)
{
Cancel = true;
}
}
catch (SingularMatrixException)
{
Cancel = true;
}
}
SearchPoint minimumPoint = searchPoints.FindMinY();
return new RegionEliminationResults(minimumPoint.X, minimumPoint.Y, succeeded);
}
Subscribe