commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/fitting/SimpleCurveFitter.java - commons-math - Git at Google

 /*
  * Licensed to the Apache Software Foundation (ASF) under one or more
  * contributor license agreements.  See the NOTICE file distributed with
  * this work for additional information regarding copyright ownership.
  * The ASF licenses this file to You under the Apache License, Version 2.0
  * (the "License"); you may not use this file except in compliance with
  * the License.  You may obtain a copy of the License at
  *
  *      http://www.apache.org/licenses/LICENSE-2.0
  *
  * Unless required by applicable law or agreed to in writing, software
  * distributed under the License is distributed on an "AS IS" BASIS,
  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  * See the License for the specific language governing permissions and
  * limitations under the License.
  */
 package org.apache.commons.math4.legacy.fitting;

 import java.util.Collections;
 import java.util.Collection;
 import java.util.Comparator;
 import java.util.List;
 import java.util.ArrayList;

 import org.apache.commons.math4.legacy.exception.ZeroException;
 import org.apache.commons.math4.legacy.exception.OutOfRangeException;
 import org.apache.commons.math4.legacy.analysis.ParametricUnivariateFunction;
 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresBuilder;
 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem;
 import org.apache.commons.math4.legacy.linear.DiagonalMatrix;

 /**
  * Fits points to a user-defined {@link ParametricUnivariateFunction function}.
  *
  * @since 3.4
  */
 public class SimpleCurveFitter extends AbstractCurveFitter {
     /** Function to fit. */
     private final ParametricUnivariateFunction function;
     /** Initial guess for the parameters. */
     private final double[] initialGuess;
     /** Parameter guesser. */
     private final ParameterGuesser guesser;
     /** Maximum number of iterations of the optimization algorithm. */
     private final int maxIter;

     /**
      * Constructor used by the factory methods.
      *
      * @param function Function to fit.
      * @param initialGuess Initial guess. Cannot be {@code null}. Its length must
      * be consistent with the number of parameters of the {@code function} to fit.
      * @param guesser Method for providing an initial guess (if {@code initialGuess}
      * is {@code null}).
      * @param maxIter Maximum number of iterations of the optimization algorithm.
      */
     protected SimpleCurveFitter(ParametricUnivariateFunction function,
                                 double[] initialGuess,
                                 ParameterGuesser guesser,
                                 int maxIter) {
         this.function = function;
         this.initialGuess = initialGuess;
         this.guesser = guesser;
         this.maxIter = maxIter;
     }

     /**
      * Creates a curve fitter.
      * The maximum number of iterations of the optimization algorithm is set
      * to {@link Integer#MAX_VALUE}.
      *
      * @param f Function to fit.
      * @param start Initial guess for the parameters.  Cannot be {@code null}.
      * Its length must be consistent with the number of parameters of the
      * function to fit.
      * @return a curve fitter.
      *
      * @see #withStartPoint(double[])
      * @see #withMaxIterations(int)
      */
     public static SimpleCurveFitter create(ParametricUnivariateFunction f,
                                            double[] start) {
         return new SimpleCurveFitter(f, start, null, Integer.MAX_VALUE);
     }

     /**
      * Creates a curve fitter.
      * The maximum number of iterations of the optimization algorithm is set
      * to {@link Integer#MAX_VALUE}.
      *
      * @param f Function to fit.
      * @param guesser Method for providing an initial guess.
      * @return a curve fitter.
      *
      * @see #withStartPoint(double[])
      * @see #withMaxIterations(int)
      */
     public static SimpleCurveFitter create(ParametricUnivariateFunction f,
                                            ParameterGuesser guesser) {
         return new SimpleCurveFitter(f, null, guesser, Integer.MAX_VALUE);
     }

     /**
      * Configure the start point (initial guess).
      * @param newStart new start point (initial guess)
      * @return a new instance.
      */
     public SimpleCurveFitter withStartPoint(double[] newStart) {
         return new SimpleCurveFitter(function,
                                      newStart.clone(),
                                      null,
                                      maxIter);
     }

     /**
      * Configure the maximum number of iterations.
      * @param newMaxIter maximum number of iterations
      * @return a new instance.
      */
     public SimpleCurveFitter withMaxIterations(int newMaxIter) {
         return new SimpleCurveFitter(function,
                                      initialGuess,
                                      guesser,
                                      newMaxIter);
     }

     /** {@inheritDoc} */
     @Override
     protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> observations) {
         // Prepare least-squares problem.
         final int len = observations.size();
         final double[] target  = new double[len];
         final double[] weights = new double[len];

         int count = 0;
         for (WeightedObservedPoint obs : observations) {
             target[count]  = obs.getY();
             weights[count] = obs.getWeight();
             ++count;
         }

         final AbstractCurveFitter.TheoreticalValuesFunction model
             = new AbstractCurveFitter.TheoreticalValuesFunction(function,
                                                                 observations);

         final double[] startPoint = initialGuess != null ?
             initialGuess :
             // Compute estimation.
             guesser.guess(observations);

         // Create an optimizer for fitting the curve to the observed points.
         return new LeastSquaresBuilder().
                 maxEvaluations(Integer.MAX_VALUE).
                 maxIterations(maxIter).
                 start(startPoint).
                 target(target).
                 weight(new DiagonalMatrix(weights)).
                 model(model.getModelFunction(), model.getModelFunctionJacobian()).
                 build();
     }

     /**
      * Guesses the parameters.
      */
     public abstract static class ParameterGuesser {
         /** Comparator. */
         private static final Comparator<WeightedObservedPoint> CMP = new Comparator<WeightedObservedPoint>() {
                 /** {@inheritDoc} */
                 @Override
                 public int compare(WeightedObservedPoint p1,
                                    WeightedObservedPoint p2) {
                     if (p1 == null && p2 == null) {
                         return 0;
                     }
                     if (p1 == null) {
                         return -1;
                     }
                     if (p2 == null) {
                         return 1;
                     }
                     int comp = Double.compare(p1.getX(), p2.getX());
                     if (comp != 0) {
                         return comp;
                     }
                     comp = Double.compare(p1.getY(), p2.getY());
                     if (comp != 0) {
                         return comp;
                     }
                     return Double.compare(p1.getWeight(), p2.getWeight());
                 }
             };

         /**
          * Computes an estimation of the parameters.
          *
          * @param obs Observations.
          * @return the guessed parameters.
          */
         public abstract double[] guess(Collection<WeightedObservedPoint> obs);

         /**
          * Sort the observations.
          *
          * @param unsorted Input observations.
          * @return the input observations, sorted.
          */
         protected List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
             final List<WeightedObservedPoint> observations = new ArrayList<>(unsorted);
             Collections.sort(observations, CMP);
             return observations;
         }

         /**
          * Finds index of point in specified points with the largest Y.
          *
          * @param points Points to search.
          * @return the index in specified points array.
          */
         protected int findMaxY(WeightedObservedPoint[] points) {
             int maxYIdx = 0;
             for (int i = 1; i < points.length; i++) {
                 if (points[i].getY() > points[maxYIdx].getY()) {
                     maxYIdx = i;
                 }
             }
             return maxYIdx;
         }

         /**
          * Interpolates using the specified points to determine X at the
          * specified Y.
          *
          * @param points Points to use for interpolation.
          * @param startIdx Index within points from which to start the search for
          * interpolation bounds points.
          * @param idxStep Index step for searching interpolation bounds points.
          * @param y Y value for which X should be determined.
          * @return the value of X for the specified Y.
          * @throws ZeroException if {@code idxStep} is 0.
          * @throws OutOfRangeException if specified {@code y} is not within the
          * range of the specified {@code points}.
          */
         protected double interpolateXAtY(WeightedObservedPoint[] points,
                                          int startIdx,
                                          int idxStep,
                                          double y) {
             if (idxStep == 0) {
                 throw new ZeroException();
             }
             final WeightedObservedPoint[] twoPoints
                 = getInterpolationPointsForY(points, startIdx, idxStep, y);
             final WeightedObservedPoint p1 = twoPoints[0];
             final WeightedObservedPoint p2 = twoPoints[1];
             if (p1.getY() == y) {
                 return p1.getX();
             }
             if (p2.getY() == y) {
                 return p2.getX();
             }
             return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
                                 (p2.getY() - p1.getY()));
         }

         /**
          * Gets the two bounding interpolation points from the specified points
          * suitable for determining X at the specified Y.
          *
          * @param points Points to use for interpolation.
          * @param startIdx Index within points from which to start search for
          * interpolation bounds points.
          * @param idxStep Index step for search for interpolation bounds points.
          * @param y Y value for which X should be determined.
          * @return the array containing two points suitable for determining X at
          * the specified Y.
          * @throws ZeroException if {@code idxStep} is 0.
          * @throws OutOfRangeException if specified {@code y} is not within the
          * range of the specified {@code points}.
          */
         private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
                                                                    int startIdx,
                                                                    int idxStep,
                                                                    double y) {
             if (idxStep == 0) {
                 throw new ZeroException();
             }
             for (int i = startIdx;
                  idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
                  i += idxStep) {
                 final WeightedObservedPoint p1 = points[i];
                 final WeightedObservedPoint p2 = points[i + idxStep];
                 if (isBetween(y, p1.getY(), p2.getY())) {
                     if (idxStep < 0) {
                         return new WeightedObservedPoint[] { p2, p1 };
                     } else {
                         return new WeightedObservedPoint[] { p1, p2 };
                     }
                 }
             }

             // Boundaries are replaced by dummy values because the raised
             // exception is caught and the message never displayed.
             // TODO: Exceptions should not be used for flow control.
             throw new OutOfRangeException(y,
                                           Double.NEGATIVE_INFINITY,
                                           Double.POSITIVE_INFINITY);
         }

         /**
          * Determines whether a value is between two other values.
          *
          * @param value Value to test whether it is between {@code boundary1}
          * and {@code boundary2}.
          * @param boundary1 One end of the range.
          * @param boundary2 Other end of the range.
          * @return {@code true} if {@code value} is between {@code boundary1} and
          * {@code boundary2} (inclusive), {@code false} otherwise.
          */
         private boolean isBetween(double value,
                                   double boundary1,
                                   double boundary2) {
             return (value >= boundary1 && value <= boundary2) ||
                 (value >= boundary2 && value <= boundary1);
         }
     }
 }
	/*
	* Licensed to the Apache Software Foundation (ASF) under one or more
	* contributor license agreements. See the NOTICE file distributed with
	* this work for additional information regarding copyright ownership.
	* The ASF licenses this file to You under the Apache License, Version 2.0
	* (the "License"); you may not use this file except in compliance with
	* the License. You may obtain a copy of the License at
	*
	* http://www.apache.org/licenses/LICENSE-2.0
	*
	* Unless required by applicable law or agreed to in writing, software
	* distributed under the License is distributed on an "AS IS" BASIS,
	* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
	* See the License for the specific language governing permissions and
	* limitations under the License.
	*/
	package org.apache.commons.math4.legacy.fitting;

	import java.util.Collections;
	import java.util.Collection;
	import java.util.Comparator;
	import java.util.List;
	import java.util.ArrayList;

	import org.apache.commons.math4.legacy.exception.ZeroException;
	import org.apache.commons.math4.legacy.exception.OutOfRangeException;
	import org.apache.commons.math4.legacy.analysis.ParametricUnivariateFunction;
	import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresBuilder;
	import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem;
	import org.apache.commons.math4.legacy.linear.DiagonalMatrix;

	/**
	* Fits points to a user-defined {@link ParametricUnivariateFunction function}.
	*
	* @since 3.4
	*/
	public class SimpleCurveFitter extends AbstractCurveFitter {
	/** Function to fit. */
	private final ParametricUnivariateFunction function;
	/** Initial guess for the parameters. */
	private final double[] initialGuess;
	/** Parameter guesser. */
	private final ParameterGuesser guesser;
	/** Maximum number of iterations of the optimization algorithm. */
	private final int maxIter;

	/**
	* Constructor used by the factory methods.
	*
	* @param function Function to fit.
	* @param initialGuess Initial guess. Cannot be {@code null}. Its length must
	* be consistent with the number of parameters of the {@code function} to fit.
	* @param guesser Method for providing an initial guess (if {@code initialGuess}
	* is {@code null}).
	* @param maxIter Maximum number of iterations of the optimization algorithm.
	*/
	protected SimpleCurveFitter(ParametricUnivariateFunction function,
	double[] initialGuess,
	ParameterGuesser guesser,
	int maxIter) {
	this.function = function;
	this.initialGuess = initialGuess;
	this.guesser = guesser;
	this.maxIter = maxIter;
	}

	/**
	* Creates a curve fitter.
	* The maximum number of iterations of the optimization algorithm is set
	* to {@link Integer#MAX_VALUE}.
	*
	* @param f Function to fit.
	* @param start Initial guess for the parameters. Cannot be {@code null}.
	* Its length must be consistent with the number of parameters of the
	* function to fit.
	* @return a curve fitter.
	*
	* @see #withStartPoint(double[])
	* @see #withMaxIterations(int)
	*/
	public static SimpleCurveFitter create(ParametricUnivariateFunction f,
	double[] start) {
	return new SimpleCurveFitter(f, start, null, Integer.MAX_VALUE);
	}

	/**
	* Creates a curve fitter.
	* The maximum number of iterations of the optimization algorithm is set
	* to {@link Integer#MAX_VALUE}.
	*
	* @param f Function to fit.
	* @param guesser Method for providing an initial guess.
	* @return a curve fitter.
	*
	* @see #withStartPoint(double[])
	* @see #withMaxIterations(int)
	*/
	public static SimpleCurveFitter create(ParametricUnivariateFunction f,
	ParameterGuesser guesser) {
	return new SimpleCurveFitter(f, null, guesser, Integer.MAX_VALUE);
	}

	/**
	* Configure the start point (initial guess).
	* @param newStart new start point (initial guess)
	* @return a new instance.
	*/
	public SimpleCurveFitter withStartPoint(double[] newStart) {
	return new SimpleCurveFitter(function,
	newStart.clone(),
	null,
	maxIter);
	}

	/**
	* Configure the maximum number of iterations.
	* @param newMaxIter maximum number of iterations
	* @return a new instance.
	*/
	public SimpleCurveFitter withMaxIterations(int newMaxIter) {
	return new SimpleCurveFitter(function,
	initialGuess,
	guesser,
	newMaxIter);
	}

	/** {@inheritDoc} */
	@Override
	protected LeastSquaresProblem getProblem(Collection<WeightedObservedPoint> observations) {
	// Prepare least-squares problem.
	final int len = observations.size();
	final double[] target = new double[len];
	final double[] weights = new double[len];

	int count = 0;
	for (WeightedObservedPoint obs : observations) {
	target[count] = obs.getY();
	weights[count] = obs.getWeight();
	++count;
	}

	final AbstractCurveFitter.TheoreticalValuesFunction model
	= new AbstractCurveFitter.TheoreticalValuesFunction(function,
	observations);

	final double[] startPoint = initialGuess != null ?
	initialGuess :
	// Compute estimation.
	guesser.guess(observations);

	// Create an optimizer for fitting the curve to the observed points.
	return new LeastSquaresBuilder().
	maxEvaluations(Integer.MAX_VALUE).
	maxIterations(maxIter).
	start(startPoint).
	target(target).
	weight(new DiagonalMatrix(weights)).
	model(model.getModelFunction(), model.getModelFunctionJacobian()).
	build();
	}

	/**
	* Guesses the parameters.
	*/
	public abstract static class ParameterGuesser {
	/** Comparator. */
	private static final Comparator<WeightedObservedPoint> CMP = new Comparator<WeightedObservedPoint>() {
	/** {@inheritDoc} */
	@Override
	public int compare(WeightedObservedPoint p1,
	WeightedObservedPoint p2) {
	if (p1 == null && p2 == null) {
	return 0;
	}
	if (p1 == null) {
	return -1;
	}
	if (p2 == null) {
	return 1;
	}
	int comp = Double.compare(p1.getX(), p2.getX());
	if (comp != 0) {
	return comp;
	}
	comp = Double.compare(p1.getY(), p2.getY());
	if (comp != 0) {
	return comp;
	}
	return Double.compare(p1.getWeight(), p2.getWeight());
	}
	};

	/**
	* Computes an estimation of the parameters.
	*
	* @param obs Observations.
	* @return the guessed parameters.
	*/
	public abstract double[] guess(Collection<WeightedObservedPoint> obs);

	/**
	* Sort the observations.
	*
	* @param unsorted Input observations.
	* @return the input observations, sorted.
	*/
	protected List<WeightedObservedPoint> sortObservations(Collection<WeightedObservedPoint> unsorted) {
	final List<WeightedObservedPoint> observations = new ArrayList<>(unsorted);
	Collections.sort(observations, CMP);
	return observations;
	}

	/**
	* Finds index of point in specified points with the largest Y.
	*
	* @param points Points to search.
	* @return the index in specified points array.
	*/
	protected int findMaxY(WeightedObservedPoint[] points) {
	int maxYIdx = 0;
	for (int i = 1; i < points.length; i++) {
	if (points[i].getY() > points[maxYIdx].getY()) {
	maxYIdx = i;
	}
	}
	return maxYIdx;
	}

	/**
	* Interpolates using the specified points to determine X at the
	* specified Y.
	*
	* @param points Points to use for interpolation.
	* @param startIdx Index within points from which to start the search for
	* interpolation bounds points.
	* @param idxStep Index step for searching interpolation bounds points.
	* @param y Y value for which X should be determined.
	* @return the value of X for the specified Y.
	* @throws ZeroException if {@code idxStep} is 0.
	* @throws OutOfRangeException if specified {@code y} is not within the
	* range of the specified {@code points}.
	*/
	protected double interpolateXAtY(WeightedObservedPoint[] points,
	int startIdx,
	int idxStep,
	double y) {
	if (idxStep == 0) {
	throw new ZeroException();
	}
	final WeightedObservedPoint[] twoPoints
	= getInterpolationPointsForY(points, startIdx, idxStep, y);
	final WeightedObservedPoint p1 = twoPoints[0];
	final WeightedObservedPoint p2 = twoPoints[1];
	if (p1.getY() == y) {
	return p1.getX();
	}
	if (p2.getY() == y) {
	return p2.getX();
	}
	return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
	(p2.getY() - p1.getY()));
	}

	/**
	* Gets the two bounding interpolation points from the specified points
	* suitable for determining X at the specified Y.
	*
	* @param points Points to use for interpolation.
	* @param startIdx Index within points from which to start search for
	* interpolation bounds points.
	* @param idxStep Index step for search for interpolation bounds points.
	* @param y Y value for which X should be determined.
	* @return the array containing two points suitable for determining X at
	* the specified Y.
	* @throws ZeroException if {@code idxStep} is 0.
	* @throws OutOfRangeException if specified {@code y} is not within the
	* range of the specified {@code points}.
	*/
	private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
	int startIdx,
	int idxStep,
	double y) {
	if (idxStep == 0) {
	throw new ZeroException();
	}
	for (int i = startIdx;
	idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
	i += idxStep) {
	final WeightedObservedPoint p1 = points[i];
	final WeightedObservedPoint p2 = points[i + idxStep];
	if (isBetween(y, p1.getY(), p2.getY())) {
	if (idxStep < 0) {
	return new WeightedObservedPoint[] { p2, p1 };
	} else {
	return new WeightedObservedPoint[] { p1, p2 };
	}
	}
	}

	// Boundaries are replaced by dummy values because the raised
	// exception is caught and the message never displayed.
	// TODO: Exceptions should not be used for flow control.
	throw new OutOfRangeException(y,
	Double.NEGATIVE_INFINITY,
	Double.POSITIVE_INFINITY);
	}

	/**
	* Determines whether a value is between two other values.
	*
	* @param value Value to test whether it is between {@code boundary1}
	* and {@code boundary2}.
	* @param boundary1 One end of the range.
	* @param boundary2 Other end of the range.
	* @return {@code true} if {@code value} is between {@code boundary1} and
	* {@code boundary2} (inclusive), {@code false} otherwise.
	*/
	private boolean isBetween(double value,
	double boundary1,
	double boundary2) {
	return (value >= boundary1 && value <= boundary2) \|\|
	(value >= boundary2 && value <= boundary1);
	}
	}
	}