1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.analysis.interpolation; 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.io.Serializable; 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Arrays; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.analysis.polynomials.PolynomialSplineFunction; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.Localizable; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Implements the <a href="http://en.wikipedia.org/wiki/Local_regression"> 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Local Regression Algorithm</a> (also Loess, Lowess) for interpolation of 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * real univariate functions. 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p/> 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * For reference, see 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <a href="http://www.math.tau.ac.il/~yekutiel/MA seminar/Cleveland 1979.pdf"> 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * William S. Cleveland - Robust Locally Weighted Regression and Smoothing 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Scatterplots</a> 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p/> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This class implements both the loess method and serves as an interpolation 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * adapter to it, allowing to build a spline on the obtained loess fit. 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.0 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class LoessInterpolator 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond implements UnivariateRealInterpolator, Serializable { 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Default value of the bandwidth parameter. */ 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public static final double DEFAULT_BANDWIDTH = 0.3; 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Default value of the number of robustness iterations. */ 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public static final int DEFAULT_ROBUSTNESS_ITERS = 2; 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Default value for accuracy. 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public static final double DEFAULT_ACCURACY = 1e-12; 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** serializable version identifier. */ 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static final long serialVersionUID = 5204927143605193821L; 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The bandwidth parameter: when computing the loess fit at 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a particular point, this fraction of source points closest 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to the current point is taken into account for computing 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a least-squares regression. 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p/> 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0.25 to 0.5. 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double bandwidth; 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The number of robustness iterations parameter: this many 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * robustness iterations are done. 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p/> 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0 (just the initial fit without any 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * robustness iterations) to 4. 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final int robustnessIters; 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If the median residual at a certain robustness iteration 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * is less than this amount, no more iterations are done. 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double accuracy; 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Constructs a new {@link LoessInterpolator} 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * with a bandwidth of {@link #DEFAULT_BANDWIDTH}, 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #DEFAULT_ROBUSTNESS_ITERS} robustness iterations 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * and an accuracy of {#link #DEFAULT_ACCURACY}. 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See {@link #LoessInterpolator(double, int, double)} for an explanation of 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the parameters. 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LoessInterpolator() { 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.bandwidth = DEFAULT_BANDWIDTH; 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.robustnessIters = DEFAULT_ROBUSTNESS_ITERS; 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.accuracy = DEFAULT_ACCURACY; 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Constructs a new {@link LoessInterpolator} 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * with given bandwidth and number of robustness iterations. 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p> 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calling this constructor is equivalent to calling {link {@link 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * #LoessInterpolator(double, int, double) LoessInterpolator(bandwidth, 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * robustnessIters, LoessInterpolator.DEFAULT_ACCURACY)} 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </p> 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param bandwidth when computing the loess fit at 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a particular point, this fraction of source points closest 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to the current point is taken into account for computing 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a least-squares regression.</br> 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0.25 to 0.5, the default value is 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #DEFAULT_BANDWIDTH}. 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param robustnessIters This many robustness iterations are done.</br> 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0 (just the initial fit without any 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * robustness iterations) to 4, the default value is 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #DEFAULT_ROBUSTNESS_ITERS}. 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if bandwidth does not lie in the interval [0,1] 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * or if robustnessIters is negative. 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see #LoessInterpolator(double, int, double) 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LoessInterpolator(double bandwidth, int robustnessIters) throws MathException { 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this(bandwidth, robustnessIters, DEFAULT_ACCURACY); 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Constructs a new {@link LoessInterpolator} 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * with given bandwidth, number of robustness iterations and accuracy. 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param bandwidth when computing the loess fit at 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a particular point, this fraction of source points closest 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * to the current point is taken into account for computing 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a least-squares regression.</br> 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0.25 to 0.5, the default value is 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #DEFAULT_BANDWIDTH}. 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param robustnessIters This many robustness iterations are done.</br> 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A sensible value is usually 0 (just the initial fit without any 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * robustness iterations) to 4, the default value is 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link #DEFAULT_ROBUSTNESS_ITERS}. 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param accuracy If the median residual at a certain robustness iteration 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * is less than this amount, no more iterations are done. 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if bandwidth does not lie in the interval [0,1] 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * or if robustnessIters is negative. 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see #LoessInterpolator(double, int) 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public LoessInterpolator(double bandwidth, int robustnessIters, double accuracy) throws MathException { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (bandwidth < 0 || bandwidth > 1) { 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.BANDWIDTH_OUT_OF_INTERVAL, 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond bandwidth); 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.bandwidth = bandwidth; 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (robustnessIters < 0) { 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.NEGATIVE_ROBUSTNESS_ITERATIONS, robustnessIters); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.robustnessIters = robustnessIters; 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.accuracy = accuracy; 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute an interpolating function by performing a loess fit 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * on the data at the original abscissae and then building a cubic spline 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * with a 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * {@link org.apache.commons.math.analysis.interpolation.SplineInterpolator} 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * on the resulting fit. 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param xval the arguments for the interpolation points 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yval the values for the interpolation points 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return A cubic spline built upon a loess fit to the data at the original abscissae 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if some of the following conditions are false: 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <ul> 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> Arguments and values are of the same size that is greater than zero</li> 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> The arguments are in a strictly increasing order</li> 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> All arguments and values are finite real numbers</li> 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </ul> 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final PolynomialSplineFunction interpolate( 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] xval, final double[] yval) throws MathException { 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new SplineInterpolator().interpolate(xval, smooth(xval, yval)); 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute a weighted loess fit on the data at the original abscissae. 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param xval the arguments for the interpolation points 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yval the values for the interpolation points 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param weights point weights: coefficients by which the robustness weight of a point is multiplied 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return values of the loess fit at corresponding original abscissae 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if some of the following conditions are false: 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <ul> 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> Arguments and values are of the same size that is greater than zero</li> 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> The arguments are in a strictly increasing order</li> 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> All arguments and values are finite real numbers</li> 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </ul> 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 2.1 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final double[] smooth(final double[] xval, final double[] yval, final double[] weights) 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MathException { 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xval.length != yval.length) { 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.MISMATCHED_LOESS_ABSCISSA_ORDINATE_ARRAYS, 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond xval.length, yval.length); 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = xval.length; 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (n == 0) { 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.LOESS_EXPECTS_AT_LEAST_ONE_POINT); 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond checkAllFiniteReal(xval, LocalizedFormats.NON_REAL_FINITE_ABSCISSA); 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond checkAllFiniteReal(yval, LocalizedFormats.NON_REAL_FINITE_ORDINATE); 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond checkAllFiniteReal(weights, LocalizedFormats.NON_REAL_FINITE_WEIGHT); 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond checkStrictlyIncreasing(xval); 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (n == 1) { 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new double[]{yval[0]}; 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (n == 2) { 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new double[]{yval[0], yval[1]}; 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int bandwidthInPoints = (int) (bandwidth * n); 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (bandwidthInPoints < 2) { 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.TOO_SMALL_BANDWIDTH, 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond n, 2.0 / n, bandwidth); 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] res = new double[n]; 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] residuals = new double[n]; 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] sortedResiduals = new double[n]; 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] robustnessWeights = new double[n]; 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Do an initial fit and 'robustnessIters' robustness iterations. 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // This is equivalent to doing 'robustnessIters+1' robustness iterations 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // starting with all robustness weights set to 1. 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Arrays.fill(robustnessWeights, 1); 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int iter = 0; iter <= robustnessIters; ++iter) { 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int[] bandwidthInterval = {0, bandwidthInPoints - 1}; 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // At each x, compute a local weighted linear regression 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x = xval[i]; 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Find out the interval of source points on which 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // a regression is to be made. 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (i > 0) { 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond updateBandwidthInterval(xval, weights, i, bandwidthInterval); 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int ileft = bandwidthInterval[0]; 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int iright = bandwidthInterval[1]; 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Compute the point of the bandwidth interval that is 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // farthest from x 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int edge; 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xval[i] - xval[ileft] > xval[iright] - xval[i]) { 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond edge = ileft; 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond edge = iright; 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Compute a least-squares linear fit weighted by 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // the product of robustness weights and the tricube 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // weight function. 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // See http://en.wikipedia.org/wiki/Linear_regression 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // (section "Univariate linear case") 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // and http://en.wikipedia.org/wiki/Weighted_least_squares 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // (section "Weighted least squares") 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumWeights = 0; 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumX = 0; 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumXSquared = 0; 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumY = 0; 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double sumXY = 0; 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double denom = FastMath.abs(1.0 / (xval[edge] - x)); 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = ileft; k <= iright; ++k) { 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xk = xval[k]; 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yk = yval[k]; 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double dist = (k < i) ? x - xk : xk - x; 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double w = tricube(dist * denom) * robustnessWeights[k] * weights[k]; 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double xkw = xk * w; 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumWeights += w; 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumX += xkw; 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumXSquared += xk * xkw; 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumY += yk * w; 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond sumXY += yk * xkw; 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double meanX = sumX / sumWeights; 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double meanY = sumY / sumWeights; 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double meanXY = sumXY / sumWeights; 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double meanXSquared = sumXSquared / sumWeights; 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double beta; 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.sqrt(FastMath.abs(meanXSquared - meanX * meanX)) < accuracy) { 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond beta = 0; 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond beta = (meanXY - meanX * meanY) / (meanXSquared - meanX * meanX); 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double alpha = meanY - beta * meanX; 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond res[i] = beta * x + alpha; 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond residuals[i] = FastMath.abs(yval[i] - res[i]); 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // No need to recompute the robustness weights at the last 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // iteration, they won't be needed anymore 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (iter == robustnessIters) { 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond break; 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Recompute the robustness weights. 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Find the median residual. 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // An arraycopy and a sort are completely tractable here, 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // because the preceding loop is a lot more expensive 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond System.arraycopy(residuals, 0, sortedResiduals, 0, n); 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Arrays.sort(sortedResiduals); 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double medianResidual = sortedResiduals[n / 2]; 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (FastMath.abs(medianResidual) < accuracy) { 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond break; 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < n; ++i) { 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double arg = residuals[i] / (6 * medianResidual); 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (arg >= 1) { 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond robustnessWeights[i] = 0; 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } else { 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double w = 1 - arg * arg; 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond robustnessWeights[i] = w * w; 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return res; 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute a loess fit on the data at the original abscissae. 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param xval the arguments for the interpolation points 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param yval the values for the interpolation points 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return values of the loess fit at corresponding original abscissae 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if some of the following conditions are false: 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <ul> 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> Arguments and values are of the same size that is greater than zero</li> 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> The arguments are in a strictly increasing order</li> 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <li> All arguments and values are finite real numbers</li> 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * </ul> 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public final double[] smooth(final double[] xval, final double[] yval) 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MathException { 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (xval.length != yval.length) { 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.MISMATCHED_LOESS_ABSCISSA_ORDINATE_ARRAYS, 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond xval.length, yval.length); 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] unitWeights = new double[xval.length]; 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Arrays.fill(unitWeights, 1.0); 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return smooth(xval, yval, unitWeights); 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Given an index interval into xval that embraces a certain number of 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * points closest to xval[i-1], update the interval so that it embraces 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the same number of points closest to xval[i], ignoring zero weights. 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param xval arguments array 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param weights weights array 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param i the index around which the new interval should be computed 380dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param bandwidthInterval a two-element array {left, right} such that: <p/> 381dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <tt>(left==0 or xval[i] - xval[left-1] > xval[right] - xval[i])</tt> 382dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p/> and also <p/> 383dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <tt>(right==xval.length-1 or xval[right+1] - xval[i] > xval[i] - xval[left])</tt>. 384dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The array will be updated. 385dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 386dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static void updateBandwidthInterval(final double[] xval, final double[] weights, 387dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int i, 388dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int[] bandwidthInterval) { 389dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int left = bandwidthInterval[0]; 390dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int right = bandwidthInterval[1]; 391dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 392dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // The right edge should be adjusted if the next point to the right 393dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // is closer to xval[i] than the leftmost point of the current interval 394dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int nextRight = nextNonzero(weights, right); 395dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (nextRight < xval.length && xval[nextRight] - xval[i] < xval[i] - xval[left]) { 396dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int nextLeft = nextNonzero(weights, bandwidthInterval[0]); 397dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond bandwidthInterval[0] = nextLeft; 398dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond bandwidthInterval[1] = nextRight; 399dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 400dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 401dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 402dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 403dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Returns the smallest index j such that j > i && (j==weights.length || weights[j] != 0) 404dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param weights weights array 405dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param i the index from which to start search; must be < weights.length 406dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return the smallest index j such that j > i && (j==weights.length || weights[j] != 0) 407dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 408dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static int nextNonzero(final double[] weights, final int i) { 409dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int j = i + 1; 410dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond while(j < weights.length && weights[j] == 0) { 411dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond j++; 412dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 413dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return j; 414dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 415dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 416dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 417dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Compute the 418dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <a href="http://en.wikipedia.org/wiki/Local_regression#Weight_function">tricube</a> 419dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * weight function 420dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 421dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param x the argument 422dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return (1-|x|^3)^3 423dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 424dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static double tricube(final double x) { 425dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double tmp = 1 - x * x * x; 426dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return tmp * tmp * tmp; 427dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 428dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 429dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 430dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Check that all elements of an array are finite real numbers. 431dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 432dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param values the values array 433dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param pattern pattern of the error message 434dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if one of the values is not a finite real number 435dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 436dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static void checkAllFiniteReal(final double[] values, final Localizable pattern) 437dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MathException { 438dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < values.length; i++) { 439dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double x = values[i]; 440dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (Double.isInfinite(x) || Double.isNaN(x)) { 441dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(pattern, i, x); 442dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 443dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 444dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 445dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 446dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 447dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Check that elements of the abscissae array are in a strictly 448dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * increasing order. 449dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 450dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param xval the abscissae array 451dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws MathException if the abscissae array 452dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * is not in a strictly increasing order 453dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 454dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static void checkStrictlyIncreasing(final double[] xval) 455dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws MathException { 456dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < xval.length; ++i) { 457dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (i >= 1 && xval[i - 1] >= xval[i]) { 458dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new MathException(LocalizedFormats.OUT_OF_ORDER_ABSCISSA_ARRAY, 459dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond i - 1, xval[i - 1], i, xval[i]); 460dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 461dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 462dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 463dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 464