GaussianFitter.java

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

  19. import java.util.Arrays;
  20. import java.util.Comparator;
  21. import org.apache.commons.math3.analysis.function.Gaussian;
  22. import org.apache.commons.math3.exception.NullArgumentException;
  23. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  24. import org.apache.commons.math3.exception.OutOfRangeException;
  25. import org.apache.commons.math3.exception.ZeroException;
  26. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  27. import org.apache.commons.math3.exception.util.LocalizedFormats;
  28. import org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer;
  29. import org.apache.commons.math3.util.FastMath;

  31. /**
  32.  * Fits points to a {@link
  33.  * org.apache.commons.math3.analysis.function.Gaussian.Parametric Gaussian} function.
  34.  * <p>
  35.  * Usage example:
  36.  * <pre>
  37.  *   GaussianFitter fitter = new GaussianFitter(
  38.  *     new LevenbergMarquardtOptimizer());
  39.  *   fitter.addObservedPoint(4.0254623,  531026.0);
  40.  *   fitter.addObservedPoint(4.03128248, 984167.0);
  41.  *   fitter.addObservedPoint(4.03839603, 1887233.0);
  42.  *   fitter.addObservedPoint(4.04421621, 2687152.0);
  43.  *   fitter.addObservedPoint(4.05132976, 3461228.0);
  44.  *   fitter.addObservedPoint(4.05326982, 3580526.0);
  45.  *   fitter.addObservedPoint(4.05779662, 3439750.0);
  46.  *   fitter.addObservedPoint(4.0636168,  2877648.0);
  47.  *   fitter.addObservedPoint(4.06943698, 2175960.0);
  48.  *   fitter.addObservedPoint(4.07525716, 1447024.0);
  49.  *   fitter.addObservedPoint(4.08237071, 717104.0);
  50.  *   fitter.addObservedPoint(4.08366408, 620014.0);
  51.  *   double[] parameters = fitter.fit();
  52.  * </pre>
  53.  *
  54.  * @since 2.2
  55.  * @deprecated As of 3.3. Please use {@link GaussianCurveFitter} and
  56.  * {@link WeightedObservedPoints} instead.
  57.  */
  58. @Deprecated
  59. public class GaussianFitter extends CurveFitter<Gaussian.Parametric> {
  60.     /**
  61.      * Constructs an instance using the specified optimizer.
  62.      *
  63.      * @param optimizer Optimizer to use for the fitting.
  64.      */
  65.     public GaussianFitter(MultivariateVectorOptimizer optimizer) {
  66.         super(optimizer);
  67.     }

  69.     /**
  70.      * Fits a Gaussian function to the observed points.
  71.      *
  72.      * @param initialGuess First guess values in the following order:
  73.      * <ul>
  74.      *  <li>Norm</li>
  75.      *  <li>Mean</li>
  76.      *  <li>Sigma</li>
  77.      * </ul>
  78.      * @return the parameters of the Gaussian function that best fits the
  79.      * observed points (in the same order as above).
  80.      * @since 3.0
  81.      */
  82.     public double[] fit(double[] initialGuess) {
  83.         final Gaussian.Parametric f = new Gaussian.Parametric() {
  84.                 @Override
  85.                 public double value(double x, double ... p) {
  86.                     double v = Double.POSITIVE_INFINITY;
  87.                     try {
  88.                         v = super.value(x, p);
  89.                     } catch (NotStrictlyPositiveException e) { // NOPMD
  90.                         // Do nothing.
  91.                     }
  92.                     return v;
  93.                 }

  95.                 @Override
  96.                 public double[] gradient(double x, double ... p) {
  97.                     double[] v = { Double.POSITIVE_INFINITY,
  98.                                    Double.POSITIVE_INFINITY,
  99.                                    Double.POSITIVE_INFINITY };
  100.                     try {
  101.                         v = super.gradient(x, p);
  102.                     } catch (NotStrictlyPositiveException e) { // NOPMD
  103.                         // Do nothing.
  104.                     }
  105.                     return v;
  106.                 }
  107.             };

  109.         return fit(f, initialGuess);
  110.     }

  112.     /**
  113.      * Fits a Gaussian function to the observed points.
  114.      *
  115.      * @return the parameters of the Gaussian function that best fits the
  116.      * observed points (in the same order as above).
  117.      */
  118.     public double[] fit() {
  119.         final double[] guess = (new ParameterGuesser(getObservations())).guess();
  120.         return fit(guess);
  121.     }

  123.     /**
  124.      * Guesses the parameters {@code norm}, {@code mean}, and {@code sigma}
  125.      * of a {@link org.apache.commons.math3.analysis.function.Gaussian.Parametric}
  126.      * based on the specified observed points.
  127.      */
  128.     public static class ParameterGuesser {
  129.         /** Normalization factor. */
  130.         private final double norm;
  131.         /** Mean. */
  132.         private final double mean;
  133.         /** Standard deviation. */
  134.         private final double sigma;

  136.         /**
  137.          * Constructs instance with the specified observed points.
  138.          *
  139.          * @param observations Observed points from which to guess the
  140.          * parameters of the Gaussian.
  141.          * @throws NullArgumentException if {@code observations} is
  142.          * {@code null}.
  143.          * @throws NumberIsTooSmallException if there are less than 3
  144.          * observations.
  145.          */
  146.         public ParameterGuesser(WeightedObservedPoint[] observations) {
  147.             if (observations == null) {
  148.                 throw new NullArgumentException(LocalizedFormats.INPUT_ARRAY);
  149.             }
  150.             if (observations.length < 3) {
  151.                 throw new NumberIsTooSmallException(observations.length, 3, true);
  152.             }

  154.             final WeightedObservedPoint[] sorted = sortObservations(observations);
  155.             final double[] params = basicGuess(sorted);

  157.             norm = params[0];
  158.             mean = params[1];
  159.             sigma = params[2];
  160.         }

  162.         /**
  163.          * Gets an estimation of the parameters.
  164.          *
  165.          * @return the guessed parameters, in the following order:
  166.          * <ul>
  167.          *  <li>Normalization factor</li>
  168.          *  <li>Mean</li>
  169.          *  <li>Standard deviation</li>
  170.          * </ul>
  171.          */
  172.         public double[] guess() {
  173.             return new double[] { norm, mean, sigma };
  174.         }

  176.         /**
  177.          * Sort the observations.
  178.          *
  179.          * @param unsorted Input observations.
  180.          * @return the input observations, sorted.
  181.          */
  182.         private WeightedObservedPoint[] sortObservations(WeightedObservedPoint[] unsorted) {
  183.             final WeightedObservedPoint[] observations = unsorted.clone();
  184.             final Comparator<WeightedObservedPoint> cmp
  185.                 = new Comparator<WeightedObservedPoint>() {
  186.                 public int compare(WeightedObservedPoint p1,
  187.                                    WeightedObservedPoint p2) {
  188.                     if (p1 == null && p2 == null) {
  189.                         return 0;
  190.                     }
  191.                     if (p1 == null) {
  192.                         return -1;
  193.                     }
  194.                     if (p2 == null) {
  195.                         return 1;
  196.                     }
  197.                     if (p1.getX() < p2.getX()) {
  198.                         return -1;
  199.                     }
  200.                     if (p1.getX() > p2.getX()) {
  201.                         return 1;
  202.                     }
  203.                     if (p1.getY() < p2.getY()) {
  204.                         return -1;
  205.                     }
  206.                     if (p1.getY() > p2.getY()) {
  207.                         return 1;
  208.                     }
  209.                     if (p1.getWeight() < p2.getWeight()) {
  210.                         return -1;
  211.                     }
  212.                     if (p1.getWeight() > p2.getWeight()) {
  213.                         return 1;
  214.                     }
  215.                     return 0;
  216.                 }
  217.             };

  219.             Arrays.sort(observations, cmp);
  220.             return observations;
  221.         }

  223.         /**
  224.          * Guesses the parameters based on the specified observed points.
  225.          *
  226.          * @param points Observed points, sorted.
  227.          * @return the guessed parameters (normalization factor, mean and
  228.          * sigma).
  229.          */
  230.         private double[] basicGuess(WeightedObservedPoint[] points) {
  231.             final int maxYIdx = findMaxY(points);
  232.             final double n = points[maxYIdx].getY();
  233.             final double m = points[maxYIdx].getX();

  235.             double fwhmApprox;
  236.             try {
  237.                 final double halfY = n + ((m - n) / 2);
  238.                 final double fwhmX1 = interpolateXAtY(points, maxYIdx, -1, halfY);
  239.                 final double fwhmX2 = interpolateXAtY(points, maxYIdx, 1, halfY);
  240.                 fwhmApprox = fwhmX2 - fwhmX1;
  241.             } catch (OutOfRangeException e) {
  242.                 // TODO: Exceptions should not be used for flow control.
  243.                 fwhmApprox = points[points.length - 1].getX() - points[0].getX();
  244.             }
  245.             final double s = fwhmApprox / (2 * FastMath.sqrt(2 * FastMath.log(2)));

  247.             return new double[] { n, m, s };
  248.         }

  250.         /**
  251.          * Finds index of point in specified points with the largest Y.
  252.          *
  253.          * @param points Points to search.
  254.          * @return the index in specified points array.
  255.          */
  256.         private int findMaxY(WeightedObservedPoint[] points) {
  257.             int maxYIdx = 0;
  258.             for (int i = 1; i < points.length; i++) {
  259.                 if (points[i].getY() > points[maxYIdx].getY()) {
  260.                     maxYIdx = i;
  261.                 }
  262.             }
  263.             return maxYIdx;
  264.         }

  266.         /**
  267.          * Interpolates using the specified points to determine X at the
  268.          * specified Y.
  269.          *
  270.          * @param points Points to use for interpolation.
  271.          * @param startIdx Index within points from which to start the search for
  272.          * interpolation bounds points.
  273.          * @param idxStep Index step for searching interpolation bounds points.
  274.          * @param y Y value for which X should be determined.
  275.          * @return the value of X for the specified Y.
  276.          * @throws ZeroException if {@code idxStep} is 0.
  277.          * @throws OutOfRangeException if specified {@code y} is not within the
  278.          * range of the specified {@code points}.
  279.          */
  280.         private double interpolateXAtY(WeightedObservedPoint[] points,
  281.                                        int startIdx,
  282.                                        int idxStep,
  283.                                        double y)
  284.             throws OutOfRangeException {
  285.             if (idxStep == 0) {
  286.                 throw new ZeroException();
  287.             }
  288.             final WeightedObservedPoint[] twoPoints
  289.                 = getInterpolationPointsForY(points, startIdx, idxStep, y);
  290.             final WeightedObservedPoint p1 = twoPoints[0];
  291.             final WeightedObservedPoint p2 = twoPoints[1];
  292.             if (p1.getY() == y) {
  293.                 return p1.getX();
  294.             }
  295.             if (p2.getY() == y) {
  296.                 return p2.getX();
  297.             }
  298.             return p1.getX() + (((y - p1.getY()) * (p2.getX() - p1.getX())) /
  299.                                 (p2.getY() - p1.getY()));
  300.         }

  302.         /**
  303.          * Gets the two bounding interpolation points from the specified points
  304.          * suitable for determining X at the specified Y.
  305.          *
  306.          * @param points Points to use for interpolation.
  307.          * @param startIdx Index within points from which to start search for
  308.          * interpolation bounds points.
  309.          * @param idxStep Index step for search for interpolation bounds points.
  310.          * @param y Y value for which X should be determined.
  311.          * @return the array containing two points suitable for determining X at
  312.          * the specified Y.
  313.          * @throws ZeroException if {@code idxStep} is 0.
  314.          * @throws OutOfRangeException if specified {@code y} is not within the
  315.          * range of the specified {@code points}.
  316.          */
  317.         private WeightedObservedPoint[] getInterpolationPointsForY(WeightedObservedPoint[] points,
  318.                                                                    int startIdx,
  319.                                                                    int idxStep,
  320.                                                                    double y)
  321.             throws OutOfRangeException {
  322.             if (idxStep == 0) {
  323.                 throw new ZeroException();
  324.             }
  325.             for (int i = startIdx;
  326.                  idxStep < 0 ? i + idxStep >= 0 : i + idxStep < points.length;
  327.                  i += idxStep) {
  328.                 final WeightedObservedPoint p1 = points[i];
  329.                 final WeightedObservedPoint p2 = points[i + idxStep];
  330.                 if (isBetween(y, p1.getY(), p2.getY())) {
  331.                     if (idxStep < 0) {
  332.                         return new WeightedObservedPoint[] { p2, p1 };
  333.                     } else {
  334.                         return new WeightedObservedPoint[] { p1, p2 };
  335.                     }
  336.                 }
  337.             }

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

  347.         /**
  348.          * Determines whether a value is between two other values.
  349.          *
  350.          * @param value Value to test whether it is between {@code boundary1}
  351.          * and {@code boundary2}.
  352.          * @param boundary1 One end of the range.
  353.          * @param boundary2 Other end of the range.
  354.          * @return {@code true} if {@code value} is between {@code boundary1} and
  355.          * {@code boundary2} (inclusive), {@code false} otherwise.
  356.          */
  357.         private boolean isBetween(double value,
  358.                                   double boundary1,
  359.                                   double boundary2) {
  360.             return (value >= boundary1 && value <= boundary2) ||
  361.                 (value >= boundary2 && value <= boundary1);
  362.         }
  363.     }
  364. }
Created with JaCoCo 0.7.2.201409121644