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.math4.legacy.fitting.leastsquares;
18
19 import org.apache.commons.math4.legacy.analysis.MultivariateMatrixFunction;
20 import org.apache.commons.math4.legacy.analysis.MultivariateVectorFunction;
21 import org.apache.commons.math4.legacy.exception.MathIllegalStateException;
22 import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
23 import org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem.Evaluation;
24 import org.apache.commons.math4.legacy.linear.Array2DRowRealMatrix;
25 import org.apache.commons.math4.legacy.linear.ArrayRealVector;
26 import org.apache.commons.math4.legacy.linear.DiagonalMatrix;
27 import org.apache.commons.math4.legacy.linear.EigenDecomposition;
28 import org.apache.commons.math4.legacy.linear.RealMatrix;
29 import org.apache.commons.math4.legacy.linear.RealVector;
30 import org.apache.commons.math4.legacy.optim.AbstractOptimizationProblem;
31 import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
32 import org.apache.commons.math4.legacy.optim.PointVectorValuePair;
33 import org.apache.commons.math4.core.jdkmath.JdkMath;
34 import org.apache.commons.math4.legacy.core.IntegerSequence;
35 import org.apache.commons.math4.legacy.core.Pair;
36
37 /**
38 * A Factory for creating {@link LeastSquaresProblem}s.
39 *
40 * @since 3.3
41 */
42 public final class LeastSquaresFactory {
43
44 /** Prevent instantiation. */
45 private LeastSquaresFactory() {}
46
47 /**
48 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem}
49 * from the given elements.
50 *
51 * @param model the model function. Produces the computed values.
52 * @param observed the observed (target) values
53 * @param start the initial guess.
54 * @param weight the weight matrix
55 * @param checker convergence checker
56 * @param maxEvaluations the maximum number of times to evaluate the model
57 * @param maxIterations the maximum number to times to iterate in the algorithm
58 * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
59 * will defer the evaluation until access to the value is requested.
60 * @param paramValidator Model parameters validator.
61 * @return the specified General Least Squares problem.
62 *
63 * @since 3.4
64 */
65 public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
66 final RealVector observed,
67 final RealVector start,
68 final RealMatrix weight,
69 final ConvergenceChecker<Evaluation> checker,
70 final int maxEvaluations,
71 final int maxIterations,
72 final boolean lazyEvaluation,
73 final ParameterValidator paramValidator) {
74 final LeastSquaresProblem p = new LocalLeastSquaresProblem(model,
75 observed,
76 start,
77 checker,
78 maxEvaluations,
79 maxIterations,
80 lazyEvaluation,
81 paramValidator);
82 if (weight != null) {
83 return weightMatrix(p, weight);
84 } else {
85 return p;
86 }
87 }
88
89 /**
90 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem}
91 * from the given elements. There will be no weights applied (unit weights).
92 *
93 * @param model the model function. Produces the computed values.
94 * @param observed the observed (target) values
95 * @param start the initial guess.
96 * @param checker convergence checker
97 * @param maxEvaluations the maximum number of times to evaluate the model
98 * @param maxIterations the maximum number to times to iterate in the algorithm
99 * @return the specified General Least Squares problem.
100 */
101 public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
102 final RealVector observed,
103 final RealVector start,
104 final ConvergenceChecker<Evaluation> checker,
105 final int maxEvaluations,
106 final int maxIterations) {
107 return create(model,
108 observed,
109 start,
110 null,
111 checker,
112 maxEvaluations,
113 maxIterations,
114 false,
115 null);
116 }
117
118 /**
119 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem}
120 * from the given elements.
121 *
122 * @param model the model function. Produces the computed values.
123 * @param observed the observed (target) values
124 * @param start the initial guess.
125 * @param weight the weight matrix
126 * @param checker convergence checker
127 * @param maxEvaluations the maximum number of times to evaluate the model
128 * @param maxIterations the maximum number to times to iterate in the algorithm
129 * @return the specified General Least Squares problem.
130 */
131 public static LeastSquaresProblem create(final MultivariateJacobianFunction model,
132 final RealVector observed,
133 final RealVector start,
134 final RealMatrix weight,
135 final ConvergenceChecker<Evaluation> checker,
136 final int maxEvaluations,
137 final int maxIterations) {
138 return weightMatrix(create(model,
139 observed,
140 start,
141 checker,
142 maxEvaluations,
143 maxIterations),
144 weight);
145 }
146
147 /**
148 * Create a {@link org.apache.commons.math4.legacy.fitting.leastsquares.LeastSquaresProblem}
149 * from the given elements.
150 * <p>
151 * This factory method is provided for continuity with previous interfaces. Newer
152 * applications should use {@link #create(MultivariateJacobianFunction, RealVector,
153 * RealVector, ConvergenceChecker, int, int)}, or {@link #create(MultivariateJacobianFunction,
154 * RealVector, RealVector, RealMatrix, ConvergenceChecker, int, int)}.
155 *
156 * @param model the model function. Produces the computed values.
157 * @param jacobian the jacobian of the model with respect to the parameters
158 * @param observed the observed (target) values
159 * @param start the initial guess.
160 * @param weight the weight matrix
161 * @param checker convergence checker
162 * @param maxEvaluations the maximum number of times to evaluate the model
163 * @param maxIterations the maximum number to times to iterate in the algorithm
164 * @return the specified General Least Squares problem.
165 */
166 public static LeastSquaresProblem create(final MultivariateVectorFunction model,
167 final MultivariateMatrixFunction jacobian,
168 final double[] observed,
169 final double[] start,
170 final RealMatrix weight,
171 final ConvergenceChecker<Evaluation> checker,
172 final int maxEvaluations,
173 final int maxIterations) {
174 return create(model(model, jacobian),
175 new ArrayRealVector(observed, false),
176 new ArrayRealVector(start, false),
177 weight,
178 checker,
179 maxEvaluations,
180 maxIterations);
181 }
182
183 /**
184 * Apply a dense weight matrix to the {@link LeastSquaresProblem}.
185 *
186 * @param problem the unweighted problem
187 * @param weights the matrix of weights
188 * @return a new {@link LeastSquaresProblem} with the weights applied. The original
189 * {@code problem} is not modified.
190 */
191 public static LeastSquaresProblem weightMatrix(final LeastSquaresProblem problem,
192 final RealMatrix weights) {
193 final RealMatrix weightSquareRoot = squareRoot(weights);
194 return new LeastSquaresAdapter(problem) {
195 /** {@inheritDoc} */
196 @Override
197 public Evaluation evaluate(final RealVector point) {
198 return new DenseWeightedEvaluation(super.evaluate(point), weightSquareRoot);
199 }
200 };
201 }
202
203 /**
204 * Apply a diagonal weight matrix to the {@link LeastSquaresProblem}.
205 *
206 * @param problem the unweighted problem
207 * @param weights the diagonal of the weight matrix
208 * @return a new {@link LeastSquaresProblem} with the weights applied. The original
209 * {@code problem} is not modified.
210 */
211 public static LeastSquaresProblem weightDiagonal(final LeastSquaresProblem problem,
212 final RealVector weights) {
213 // TODO more efficient implementation
214 return weightMatrix(problem, new DiagonalMatrix(weights.toArray()));
215 }
216
217 /**
218 * Count the evaluations of a particular problem. The {@code counter} will be
219 * incremented every time {@link LeastSquaresProblem#evaluate(RealVector)} is called on
220 * the <em>returned</em> problem.
221 *
222 * @param problem the problem to track.
223 * @param counter the counter to increment.
224 * @return a least squares problem that tracks evaluations
225 */
226 public static LeastSquaresProblem countEvaluations(final LeastSquaresProblem problem,
227 final IntegerSequence.Incrementor counter) {
228 return new LeastSquaresAdapter(problem) {
229
230 /** {@inheritDoc} */
231 @Override
232 public Evaluation evaluate(final RealVector point) {
233 counter.increment();
234 return super.evaluate(point);
235 }
236
237 // Delegate the rest.
238 };
239 }
240
241 /**
242 * View a convergence checker specified for a {@link PointVectorValuePair} as one
243 * specified for an {@link Evaluation}.
244 *
245 * @param checker the convergence checker to adapt.
246 * @return a convergence checker that delegates to {@code checker}.
247 */
248 public static ConvergenceChecker<Evaluation> evaluationChecker(final ConvergenceChecker<PointVectorValuePair> checker) {
249 return new ConvergenceChecker<Evaluation>() {
250 /** {@inheritDoc} */
251 @Override
252 public boolean converged(final int iteration,
253 final Evaluation previous,
254 final Evaluation current) {
255 return checker.converged(
256 iteration,
257 new PointVectorValuePair(
258 previous.getPoint().toArray(),
259 previous.getResiduals().toArray(),
260 false),
261 new PointVectorValuePair(
262 current.getPoint().toArray(),
263 current.getResiduals().toArray(),
264 false)
265 );
266 }
267 };
268 }
269
270 /**
271 * Computes the square-root of the weight matrix.
272 *
273 * @param m Symmetric, positive-definite (weight) matrix.
274 * @return the square-root of the weight matrix.
275 */
276 private static RealMatrix squareRoot(final RealMatrix m) {
277 if (m instanceof DiagonalMatrix) {
278 final int dim = m.getRowDimension();
279 final RealMatrix sqrtM = new DiagonalMatrix(dim);
280 for (int i = 0; i < dim; i++) {
281 sqrtM.setEntry(i, i, JdkMath.sqrt(m.getEntry(i, i)));
282 }
283 return sqrtM;
284 } else {
285 final EigenDecomposition dec = new EigenDecomposition(m);
286 return dec.getSquareRoot();
287 }
288 }
289
290 /**
291 * Combine a {@link MultivariateVectorFunction} with a {@link
292 * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
293 *
294 * @param value the vector value function
295 * @param jacobian the Jacobian function
296 * @return a function that computes both at the same time
297 */
298 public static MultivariateJacobianFunction model(final MultivariateVectorFunction value,
299 final MultivariateMatrixFunction jacobian) {
300 return new LocalValueAndJacobianFunction(value, jacobian);
301 }
302
303 /**
304 * Combine a {@link MultivariateVectorFunction} with a {@link
305 * MultivariateMatrixFunction} to produce a {@link MultivariateJacobianFunction}.
306 */
307 private static final class LocalValueAndJacobianFunction
308 implements ValueAndJacobianFunction {
309 /** Model. */
310 private final MultivariateVectorFunction value;
311 /** Model's Jacobian. */
312 private final MultivariateMatrixFunction jacobian;
313
314 /**
315 * @param value Model function.
316 * @param jacobian Model's Jacobian function.
317 */
318 LocalValueAndJacobianFunction(final MultivariateVectorFunction value,
319 final MultivariateMatrixFunction jacobian) {
320 this.value = value;
321 this.jacobian = jacobian;
322 }
323
324 /** {@inheritDoc} */
325 @Override
326 public Pair<RealVector, RealMatrix> value(final RealVector point) {
327 //TODO get array from RealVector without copying?
328 final double[] p = point.toArray();
329
330 // Evaluate.
331 return new Pair<>(computeValue(p), computeJacobian(p));
332 }
333
334 /** {@inheritDoc} */
335 @Override
336 public RealVector computeValue(final double[] params) {
337 return new ArrayRealVector(value.value(params), false);
338 }
339
340 /** {@inheritDoc} */
341 @Override
342 public RealMatrix computeJacobian(final double[] params) {
343 return new Array2DRowRealMatrix(jacobian.value(params), false);
344 }
345 }
346
347
348 /**
349 * A private, "field" immutable (not "real" immutable) implementation of {@link
350 * LeastSquaresProblem}.
351 * @since 3.3
352 */
353 private static final class LocalLeastSquaresProblem
354 extends AbstractOptimizationProblem<Evaluation>
355 implements LeastSquaresProblem {
356
357 /** Target values for the model function at optimum. */
358 private final RealVector target;
359 /** Model function. */
360 private final MultivariateJacobianFunction model;
361 /** Initial guess. */
362 private final RealVector start;
363 /** Whether to use lazy evaluation. */
364 private final boolean lazyEvaluation;
365 /** Model parameters validator. */
366 private final ParameterValidator paramValidator;
367
368 /**
369 * Create a {@link LeastSquaresProblem} from the given data.
370 *
371 * @param model the model function
372 * @param target the observed data
373 * @param start the initial guess
374 * @param checker the convergence checker
375 * @param maxEvaluations the allowed evaluations
376 * @param maxIterations the allowed iterations
377 * @param lazyEvaluation Whether the call to {@link Evaluation#evaluate(RealVector)}
378 * will defer the evaluation until access to the value is requested.
379 * @param paramValidator Model parameters validator.
380 */
381 LocalLeastSquaresProblem(final MultivariateJacobianFunction model,
382 final RealVector target,
383 final RealVector start,
384 final ConvergenceChecker<Evaluation> checker,
385 final int maxEvaluations,
386 final int maxIterations,
387 final boolean lazyEvaluation,
388 final ParameterValidator paramValidator) {
389 super(maxEvaluations, maxIterations, checker);
390 this.target = target;
391 this.model = model;
392 this.start = start;
393 this.lazyEvaluation = lazyEvaluation;
394 this.paramValidator = paramValidator;
395
396 if (lazyEvaluation &&
397 !(model instanceof ValueAndJacobianFunction)) {
398 // Lazy evaluation requires that value and Jacobian
399 // can be computed separately.
400 throw new MathIllegalStateException(LocalizedFormats.INVALID_IMPLEMENTATION,
401 model.getClass().getName());
402 }
403 }
404
405 /** {@inheritDoc} */
406 @Override
407 public int getObservationSize() {
408 return target.getDimension();
409 }
410
411 /** {@inheritDoc} */
412 @Override
413 public int getParameterSize() {
414 return start.getDimension();
415 }
416
417 /** {@inheritDoc} */
418 @Override
419 public RealVector getStart() {
420 return start == null ? null : start.copy();
421 }
422
423 /** {@inheritDoc} */
424 @Override
425 public Evaluation evaluate(final RealVector point) {
426 // Copy so optimizer can change point without changing our instance.
427 final RealVector p = paramValidator == null ?
428 point.copy() :
429 paramValidator.validate(point.copy());
430
431 if (lazyEvaluation) {
432 return new LazyUnweightedEvaluation((ValueAndJacobianFunction) model,
433 target,
434 p);
435 } else {
436 // Evaluate value and jacobian in one function call.
437 final Pair<RealVector, RealMatrix> value = model.value(p);
438 return new UnweightedEvaluation(value.getFirst(),
439 value.getSecond(),
440 target,
441 p);
442 }
443 }
444
445 /**
446 * Container with the model evaluation at a particular point.
447 */
448 private static final class UnweightedEvaluation extends AbstractEvaluation {
449 /** Point of evaluation. */
450 private final RealVector point;
451 /** Derivative at point. */
452 private final RealMatrix jacobian;
453 /** Computed residuals. */
454 private final RealVector residuals;
455
456 /**
457 * Create an {@link Evaluation} with no weights.
458 *
459 * @param values the computed function values
460 * @param jacobian the computed function Jacobian
461 * @param target the observed values
462 * @param point the abscissa
463 */
464 private UnweightedEvaluation(final RealVector values,
465 final RealMatrix jacobian,
466 final RealVector target,
467 final RealVector point) {
468 super(target.getDimension());
469 this.jacobian = jacobian;
470 this.point = point;
471 this.residuals = target.subtract(values);
472 }
473
474 /** {@inheritDoc} */
475 @Override
476 public RealMatrix getJacobian() {
477 return jacobian;
478 }
479
480 /** {@inheritDoc} */
481 @Override
482 public RealVector getPoint() {
483 return point;
484 }
485
486 /** {@inheritDoc} */
487 @Override
488 public RealVector getResiduals() {
489 return residuals;
490 }
491 }
492
493 /**
494 * Container with the model <em>lazy</em> evaluation at a particular point.
495 */
496 private static final class LazyUnweightedEvaluation extends AbstractEvaluation {
497 /** Point of evaluation. */
498 private final RealVector point;
499 /** Model and Jacobian functions. */
500 private final ValueAndJacobianFunction model;
501 /** Target values for the model function at optimum. */
502 private final RealVector target;
503
504 /**
505 * Create an {@link Evaluation} with no weights.
506 *
507 * @param model the model function
508 * @param target the observed values
509 * @param point the abscissa
510 */
511 private LazyUnweightedEvaluation(final ValueAndJacobianFunction model,
512 final RealVector target,
513 final RealVector point) {
514 super(target.getDimension());
515 // Safe to cast as long as we control usage of this class.
516 this.model = model;
517 this.point = point;
518 this.target = target;
519 }
520
521 /** {@inheritDoc} */
522 @Override
523 public RealMatrix getJacobian() {
524 return model.computeJacobian(point.toArray());
525 }
526
527 /** {@inheritDoc} */
528 @Override
529 public RealVector getPoint() {
530 return point;
531 }
532
533 /** {@inheritDoc} */
534 @Override
535 public RealVector getResiduals() {
536 return target.subtract(model.computeValue(point.toArray()));
537 }
538 }
539 }
540 }
541