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.GaussJordanElimination;
20 import com.irurueta.algebra.Matrix;
21 import com.irurueta.numerical.EvaluationException;
22 import com.irurueta.numerical.NotReadyException;
23
24 import java.util.Arrays;
25
26 /**
27 * Fits provided data (x,y) to a function made of a linear combination of
28 * functions used as a basis (i.e. f(x) = a * f0(x) + b * f1(x) + ...).
29 * Where f0, f1, ... is the function basis which ideally should be formed by
30 * orthogonal function.
31 * This class is based on the implementation available at Numerical Recipes
32 * 3rd Ed, page 791
33 */
34 public class SimpleSingleDimensionLinearFitter extends SingleDimensionLinearFitter {
35
36 /**
37 * Determines which parameters can be modified during estimation (if true)
38 * and which ones are locked (if false)
39 */
40 private boolean[] ia;
41
42
43 /**
44 * Constructor.
45 */
46 public SimpleSingleDimensionLinearFitter() {
47 super();
48 }
49
50 /**
51 * Constructor.
52 *
53 * @param x input points x where a linear single dimensional function f(x) =
54 * a * f0(x) + b * f1(x) + ...
55 * @param y result of evaluation of linear single dimensional function f(x)
56 * at provided x points.
57 * @param sig standard deviations of each pair of points (x, y).
58 * @throws IllegalArgumentException if provided arrays don't have the same
59 * length.
60 */
61 public SimpleSingleDimensionLinearFitter(final double[] x, final double[] y, final double[] sig) {
62 super(x, y, sig);
63 }
64
65 /**
66 * Constructor.
67 *
68 * @param x input points x where a linear single dimensional function f(x) =
69 * a * f0(x) + b * f1(x) + ...
70 * @param y result of evaluation of linear single dimensional function f(x)
71 * at provided x points.
72 * @param sig standard deviation of all pair of points assuming that
73 * standard deviations are constant.
74 * @throws IllegalArgumentException if provided arrays don't have the same
75 * length.
76 */
77 public SimpleSingleDimensionLinearFitter(final double[] x, final double[] y, final double sig) {
78 super(x, y, sig);
79 }
80
81 /**
82 * Constructor.
83 *
84 * @param evaluator evaluator to evaluate function at provided point and
85 * obtain the evaluation of function basis at such point.
86 * @throws FittingException if evaluation fails.
87 */
88 public SimpleSingleDimensionLinearFitter(
89 final LinearFitterSingleDimensionFunctionEvaluator evaluator) throws FittingException {
90 super();
91 setFunctionEvaluator(evaluator);
92 }
93
94 /**
95 * Constructor.
96 *
97 * @param evaluator evaluator to evaluate function at provided point and
98 * obtain the evaluation of function basis at such point.
99 * @param x input points x where a linear single dimensional function f(x) =
100 * a * f0(x) + b * f1(x) + ...
101 * @param y result of evaluation of linear single dimensional function f(x)
102 * at provided x points.
103 * @param sig standard deviation of all pair of points assuming that
104 * standard deviations are constant.
105 * @throws FittingException if evaluation fails.
106 * @throws IllegalArgumentException if provided arrays don't have the same
107 * length .
108 */
109 public SimpleSingleDimensionLinearFitter(
110 final LinearFitterSingleDimensionFunctionEvaluator evaluator, final double[] x, final double[] y,
111 final double[] sig) throws FittingException {
112 super(x, y, sig);
113 setFunctionEvaluator(evaluator);
114 }
115
116 /**
117 * Constructor.
118 *
119 * @param evaluator evaluator to evaluate function at provided point and
120 * obtain the evaluation of function basis at such point.
121 * @param x input points x where a linear single dimensional function f(x) =
122 * a * f0(x) + b * f1(x) + ...
123 * @param y result of evaluation of linear single dimensional function f(x)
124 * at provided x points.
125 * @param sig standard deviation of all pair of points assuming that
126 * standard deviations are constant.
127 * @throws FittingException if evaluation fails.
128 * @throws IllegalArgumentException if provided arrays don't have the same
129 * length.
130 */
131 public SimpleSingleDimensionLinearFitter(
132 final LinearFitterSingleDimensionFunctionEvaluator evaluator, final double[] x, final double[] y,
133 final double sig) throws FittingException {
134 super(x, y, sig);
135 setFunctionEvaluator(evaluator);
136 }
137
138 /**
139 * Sets function evaluator to evaluate function at a given point and obtain
140 * the evaluation of function basis at such point.
141 *
142 * @param evaluator function evaluator.
143 * @throws FittingException if evaluation fails.
144 */
145 @Override
146 public final void setFunctionEvaluator(
147 final LinearFitterSingleDimensionFunctionEvaluator evaluator) throws FittingException {
148 super.setFunctionEvaluator(evaluator);
149 if (ma > 0) {
150 ia = new boolean[ma];
151 Arrays.fill(ia, true);
152 }
153 }
154
155 /**
156 * Fits a function to provided data so that parameters associated to that
157 * function can be estimated along with their covariance matrix and chi
158 * square value.
159 *
160 * @throws FittingException if fitting fails.
161 * @throws NotReadyException if enough input data has not yet been provided.
162 */
163 @Override
164 @SuppressWarnings("Duplicates")
165 public void fit() throws FittingException, NotReadyException {
166 if (!isReady()) {
167 throw new NotReadyException();
168 }
169
170 try {
171 resultAvailable = false;
172
173 int i;
174 int j;
175 int k;
176 int l;
177 int m;
178 var mfit = 0;
179 double ym;
180 double wt;
181 double sum;
182 double sig2i;
183 for (j = 0; j < ma; j++) {
184 if (ia[j]) {
185 mfit++;
186 }
187 }
188 if (mfit == 0) {
189 throw new FittingException("lfit: no parameters to be fitted");
190 }
191 final var temp = new Matrix(mfit, mfit);
192 final var beta = new Matrix(mfit, 1);
193 for (i = 0; i < ndat; i++) {
194 evaluator.evaluate(x[i], afunc);
195 ym = y[i];
196 if (mfit < ma) {
197 for (j = 0; j < ma; j++) {
198 if (!ia[j]) ym -= a[j] * afunc[j];
199 }
200 }
201 sig2i = 1.0 / Math.pow(sig[i], 2.0);
202 for (j = 0, l = 0; l < ma; l++) {
203 if (ia[l]) {
204 wt = afunc[l] * sig2i;
205 int index;
206 for (k = 0, m = 0; m <= l; m++) {
207 if (ia[m]) {
208 index = temp.getIndex(j, k++);
209 temp.getBuffer()[index] += wt * afunc[m];
210 }
211 }
212 index = beta.getIndex(j++, 0);
213 beta.getBuffer()[index] += ym * wt;
214 }
215 }
216 }
217 for (j = 1; j < mfit; j++) {
218 for (k = 0; k < j; k++) {
219 temp.setElementAt(k, j, temp.getElementAt(j, k));
220 }
221 }
222 GaussJordanElimination.process(temp, beta);
223 for (j = 0, l = 0; l < ma; l++) {
224 if (ia[l]) {
225 a[l] = beta.getElementAt(j++, 0);
226 }
227 }
228 chisq = 0.0;
229 for (i = 0; i < ndat; i++) {
230 evaluator.evaluate(x[i], afunc);
231 sum = 0.0;
232 for (j = 0; j < ma; j++) {
233 sum += a[j] * afunc[j];
234 }
235 chisq += Math.pow((y[i] - sum) / sig[i], 2.0);
236 }
237 for (j = 0; j < mfit; j++) {
238 for (k = 0; k < mfit; k++) {
239 covar.setElementAt(j, k, temp.getElementAt(j, k));
240 }
241 }
242 for (i = mfit; i < ma; i++) {
243 for (j = 0; j < i + 1; j++) {
244 covar.setElementAt(i, j, 0.0);
245 covar.setElementAt(j, i, 0.0);
246 }
247 }
248 k = mfit - 1;
249 for (j = ma - 1; j >= 0; j--) {
250 if (ia[j]) {
251 for (i = 0; i < ma; i++) {
252 swap(covar.getBuffer(), covar.getBuffer(),
253 covar.getIndex(i, k), covar.getIndex(i, j));
254 }
255 for (i = 0; i < ma; i++) {
256 swap(covar.getBuffer(), covar.getBuffer(),
257 covar.getIndex(k, i), covar.getIndex(j, i));
258 }
259 k--;
260 }
261 }
262
263 resultAvailable = true;
264
265 } catch (final AlgebraException | EvaluationException e) {
266 throw new FittingException(e);
267 }
268 }
269
270 /**
271 * Prevents parameter at position i of linear combination of basis functions
272 * to be modified during function fitting.
273 *
274 * @param i position of parameter to be retained.
275 * @param val value to be set for parameter at position i.
276 */
277 public void hold(final int i, final double val) {
278 ia[i] = false;
279 a[i] = val;
280 }
281
282 /**
283 * Releases parameter at position i of linear combination of basis functions,
284 * so it can be modified again if needed.
285 *
286 * @param i position of parameter to be released.
287 */
288 public void free(final int i) {
289 ia[i] = true;
290 }
291
292 /**
293 * Swaps values of arrays at provided positions.
294 *
295 * @param array1 1st array.
296 * @param array2 2nd array.
297 * @param pos1 1st position.
298 * @param pos2 2nd position.
299 */
300 private static void swap(final double[] array1, final double[] array2, final int pos1, final int pos2) {
301 final var value1 = array1[pos1];
302 final var value2 = array2[pos2];
303 array1[pos1] = value2;
304 array2[pos2] = value1;
305 }
306 }