1 /*
2 * Copyright (C) 2015 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.fitting;
17
18 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.Matrix;
20 import com.irurueta.algebra.SingularValueDecomposer;
21 import com.irurueta.numerical.EvaluationException;
22 import com.irurueta.numerical.NotReadyException;
23
24 /**
25 * Fits provided data (x,y) to a function made of a linear combination of
26 * functions used as a basis (i.e. f(x) = a * f0(x) + b * f1(x) + ...).
27 * Where f0, f1, ... is the function basis which ideally should be formed by
28 * orthogonal function.
29 * This is an improvement on SimpleLinearFitter for cases where basis functions
30 * cannot be clearly distinguished from provided data, to avoid matrix
31 * singularities and obtain better results.
32 * This class is based on the implementation available at Numerical Recipes
33 * 3rd Ed, page 795.
34 */
35 public class SvdSingleDimensionLinearFitter extends SingleDimensionLinearFitter {
36
37 /**
38 * Default tolerance.
39 */
40 public static final double DEFAULT_TOL = 1e-12;
41
42 /**
43 * Tolerance to define convergence threshold for SVD.
44 */
45 private double tol;
46
47 /**
48 * Constructor.
49 */
50 public SvdSingleDimensionLinearFitter() {
51 super();
52 tol = DEFAULT_TOL;
53 }
54
55 /**
56 * Constructor.
57 *
58 * @param x input points x where a linear single dimensional function f(x) =
59 * a * f0(x) + b * f1(x) + ...
60 * @param y result of evaluation of linear single dimensional function f(x)
61 * at provided x points.
62 * @param sig standard deviations of each pair of points (x, y).
63 * @throws IllegalArgumentException if provided arrays don't have the same
64 * length.
65 */
66 public SvdSingleDimensionLinearFitter(final double[] x, final double[] y, final double[] sig) {
67 super(x, y, sig);
68 tol = DEFAULT_TOL;
69 }
70
71 /**
72 * Constructor.
73 *
74 * @param x input points x where a linear single dimensional function f(x) =
75 * a * f0(x) + b * f1(x) + ...
76 * @param y result of evaluation of linear single dimensional function f(x)
77 * at provided x points.
78 * @param sig standard deviation of all pair of points assuming that
79 * standard deviations are constant.
80 * @throws IllegalArgumentException if provided arrays don't have the same
81 * length.
82 */
83 public SvdSingleDimensionLinearFitter(final double[] x, final double[] y, final double sig) {
84 super(x, y, sig);
85 tol = DEFAULT_TOL;
86 }
87
88 /**
89 * Constructor.
90 *
91 * @param evaluator evaluator to evaluate function at provided point and
92 * obtain the evaluation of function basis at such point.
93 * @throws FittingException if evaluation fails.
94 */
95 public SvdSingleDimensionLinearFitter(final LinearFitterSingleDimensionFunctionEvaluator evaluator)
96 throws FittingException {
97 super(evaluator);
98 tol = DEFAULT_TOL;
99 }
100
101 /**
102 * Constructor.
103 *
104 * @param evaluator evaluator to evaluate function at provided point and
105 * obtain the evaluation of function basis at such point.
106 * @param x input points x where a linear single dimensional function f(x) =
107 * a * f0(x) + b * f1(x) + ...
108 * @param y result of evaluation of linear single dimensional function f(x)
109 * at provided x points.
110 * @param sig standard deviation of all pair of points assuming that
111 * standard deviations are constant.
112 * @throws FittingException if evaluation fails.
113 * @throws IllegalArgumentException if provided arrays don't have the same
114 * length.
115 */
116 public SvdSingleDimensionLinearFitter(
117 final LinearFitterSingleDimensionFunctionEvaluator evaluator, final double[] x, final double[] y,
118 final double[] sig) throws FittingException {
119 super(evaluator, x, y, sig);
120 tol = DEFAULT_TOL;
121 }
122
123 /**
124 * Constructor.
125 *
126 * @param evaluator evaluator to evaluate function at provided point and
127 * obtain the evaluation of function basis at such point.
128 * @param x input points x where a linear single dimensional function f(x) =
129 * a * f0(x) + b * f1(x) + ...
130 * @param y result of evaluation of linear single dimensional function f(x)
131 * at provided x points.
132 * @param sig standard deviation of all pair of points assuming that
133 * standard deviations are constant.
134 * @throws FittingException if evaluation fails.
135 * @throws IllegalArgumentException if provided arrays don't have the same
136 * length.
137 */
138 public SvdSingleDimensionLinearFitter(
139 final LinearFitterSingleDimensionFunctionEvaluator evaluator, final double[] x, final double[] y,
140 final double sig) throws FittingException {
141 super(evaluator, x, y, sig);
142 tol = DEFAULT_TOL;
143 }
144
145 /**
146 * Returns tolerance to define convergence threshold for SVD.
147 *
148 * @return tolerance to define convergence threshold for SVD.
149 */
150 public double getTol() {
151 return tol;
152 }
153
154 /**
155 * Sets tolerance to define convergence threshold for SVD.
156 *
157 * @param tol tolerance to define convergence threshold for SVD.
158 */
159 public void setTol(final double tol) {
160 this.tol = tol;
161 }
162
163 /**
164 * Fits a function to provided data so that parameters associated to that
165 * function can be estimated along with their covariance matrix and chi
166 * square value.
167 *
168 * @throws FittingException if fitting fails.
169 * @throws NotReadyException if enough input data has not yet been provided.
170 */
171 @SuppressWarnings("DuplicatedCode")
172 @Override
173 public void fit() throws FittingException, NotReadyException {
174 if (!isReady()) {
175 throw new NotReadyException();
176 }
177
178 try {
179 resultAvailable = false;
180
181 int i;
182 int j;
183 int k;
184 double tmp;
185 final double thresh;
186 double sum;
187 final var aa = new Matrix(ndat, ma);
188 final var b = new double[ndat];
189 for (i = 0; i < ndat; i++) {
190 evaluator.evaluate(x[i], afunc);
191 tmp = 1.0 / sig[i];
192 for (j = 0; j < ma; j++) {
193 aa.setElementAt(i, j, afunc[j] * tmp);
194 }
195 b[i] = y[i] * tmp;
196 }
197
198 final var svd = new SingularValueDecomposer(aa);
199 svd.decompose();
200 thresh = (tol > 0.0 ? tol * svd.getSingularValues()[0] : -1.0);
201 svd.solve(b, thresh, a);
202 chisq = 0.0;
203 for (i = 0; i < ndat; i++) {
204 sum = 0.0;
205 for (j = 0; j < ma; j++) {
206 sum += aa.getElementAt(i, j) * a[j];
207 }
208 chisq += Math.pow(sum - b[i], 2.0);
209 }
210 for (i = 0; i < ma; i++) {
211 for (j = 0; j < i + 1; j++) {
212 sum = 0.0;
213 final var w = svd.getSingularValues();
214 final var tsh = svd.getNegligibleSingularValueThreshold();
215 final var v = svd.getV();
216 for (k = 0; k < ma; k++) {
217 if (w[k] > tsh) {
218 sum += v.getElementAt(i, k) * v.getElementAt(j, k) / Math.pow(w[k], 2.0);
219 }
220 }
221 covar.setElementAt(j, i, sum);
222 covar.setElementAt(i, j, sum);
223 }
224 }
225
226 resultAvailable = true;
227
228 } catch (final AlgebraException | EvaluationException e) {
229 throw new FittingException(e);
230 }
231 }
232 }