Examples of org.apache.commons.math3.fitting.HarmonicFitter$ParameterGuesser

Package org.apache.commons.math3.fitting

Examples of org.apache.commons.math3.fitting.HarmonicFitter$ParameterGuesser

org.apache.commons.math3.exception.DimensionMismatchException
This class guesses harmonic coefficients from a sample.
The algorithm used to guess the coefficients is as follows:

We know f (t) at some sampling points t_i and want to find a, ω and φ such that f (t) = a cos (ω t + φ).

From the analytical expression, we can compute two primitives :
```
 If2  (t) = ∫ f²  = a² × [t + S (t)] / 2 If'2 (t) = ∫ f'² = a² ω² × [t - S (t)] / 2 where S (t) = sin (2 (ω t + φ)) / (2 ω) 
```
We can remove S between these expressions :
```
 If'2 (t) = a² ω² t - ω² If2 (t) 
```
The preceding expression shows that If'2 (t) is a linear combination of both t and If2 (t): If'2 (t) = A × t + B × If2 (t)

From the primitive, we can deduce the same form for definite integrals between t₁ and t_i for each t_i :
```
 If2 (t_i) - If2 (t₁) = A × (t_i - t₁) + B × (If2 (t_i) - If2 (t₁)) 
```
We can find the coefficients A and B that best fit the sample to this linear expression by computing the definite integrals for each sample points.

For a bilinear expression z (x_i, y_i) = A × x_i + B × y_i, the coefficients A and B that minimize a least square criterion ∑ (z_i - z (x_i, y_i))² are given by these expressions:
```
 ∑y_iy_i ∑x_iz_i - ∑x_iy_i ∑y_iz_i A = ------------------------ ∑x_ix_i ∑y_iy_i - ∑x_iy_i ∑x_iy_i ∑x_ix_i ∑y_iz_i - ∑x_iy_i ∑x_iz_i B = ------------------------ ∑x_ix_i ∑y_iy_i - ∑x_iy_i ∑x_iy_i 
```
In fact, we can assume both a and ω are positive and compute them directly, knowing that A = a² ω² and that B = - ω². The complete algorithm is therefore:
```
 for each t_i from t₁ to t_n-1, compute: f  (t_i) f' (t_i) = (f (t_i+1) - f(t_i-1)) / (t_i+1 - t_i-1) x_i = t_i - t₁ y_i = ∫ f² from t₁ to t_i z_i = ∫ f'² from t₁ to t_i update the sums ∑x_ix_i, ∑y_iy_i, ∑x_iy_i, ∑x_iz_i and ∑y_iz_i end for |-------------------------- \  | ∑y_iy_i ∑x_iz_i - ∑x_iy_i ∑y_iz_i a     =  \ | ------------------------ \| ∑x_iy_i ∑x_iz_i - ∑x_ix_i ∑y_iz_i |-------------------------- \  | ∑x_iy_i ∑x_iz_i - ∑x_ix_i ∑y_iz_i ω     =  \ | ------------------------ \| ∑x_ix_i ∑y_iy_i - ∑x_iy_i ∑x_iy_i 
```
Once we know ω, we can compute:
```
 fc = ω f (t) cos (ω t) - f' (t) sin (ω t) fs = ω f (t) sin (ω t) + f' (t) cos (ω t) 
```
It appears that fc = a ω cos (φ) and fs = -a ω sin (φ), so we can use these expressions to compute φ. The best estimate over the sample is given by averaging these expressions.

Since integrals and means are involved in the preceding estimations, these operations run in O(n) time, where n is the number of measurements.

                return false;
            }
            final int    n = FastMath.max(1, (int) FastMath.ceil(FastMath.abs(dt) / maxCheckInterval));
            final double h = dt / n;


            final UnivariateFunction f = new UnivariateFunction() {
                public double value(final double t) throws LocalMaxCountExceededException {
                    try {
                        interpolator.setInterpolatedTime(t);
                        return handler.g(t, getCompleteState(interpolator));
                    } catch (MaxCountExceededException mcee) {
                        throw new LocalMaxCountExceededException(mcee);
                    }
                }
            };


            double ta = t0;
            double ga = g0;
            for (int i = 0; i < n; ++i) {


                // evaluate handler value at the end of the substep
                final double tb = t0 + (i + 1) * h;
                interpolator.setInterpolatedTime(tb);
                final double gb = handler.g(tb, getCompleteState(interpolator));


                // check events occurrence
                if (g0Positive ^ (gb >= 0)) {
                    // there is a sign change: an event is expected during this step


                    // variation direction, with respect to the integration direction
                    increasing = gb >= ga;


                    // find the event time making sure we select a solution just at or past the exact root
                    final double root;
                    if (solver instanceof BracketedUnivariateSolver<?>) {
                        @SuppressWarnings("unchecked")
                        BracketedUnivariateSolver<UnivariateFunction> bracketing =
                                (BracketedUnivariateSolver<UnivariateFunction>) solver;
                        root = forward ?
                               bracketing.solve(maxIterationCount, f, ta, tb, AllowedSolution.RIGHT_SIDE) :
                               bracketing.solve(maxIterationCount, f, tb, ta, AllowedSolution.LEFT_SIDE);
                    } else {
                        final double baseRoot = forward ?
                                                solver.solve(maxIterationCount, f, ta, tb) :
                                                solver.solve(maxIterationCount, f, tb, ta);
                        final int remainingEval = maxIterationCount - solver.getEvaluations();
                        BracketedUnivariateSolver<UnivariateFunction> bracketing =
                                new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
                        root = forward ?
                               UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
                                                                   baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
                               UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
                                                                   baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);
                    }


                    if ((!Double.isNaN(previousEventTime)) &&
                        (FastMath.abs(root - ta) <= convergence) &&
                        (FastMath.abs(root - previousEventTime) <= convergence)) {
                        // we have either found nothing or found (again ?) a past event,
                        // retry the substep excluding this value, and taking care to have the
                        // required sign in case the g function is noisy around its zero and
                        // crosses the axis several times
                        do {
                            ta = forward ? ta + convergence : ta - convergence;
                            ga = f.value(ta);
                        } while ((g0Positive ^ (ga >= 0)) && (forward ^ (ta >= tb)));
                        --i;
                    } else if (Double.isNaN(previousEventTime) ||
                               (FastMath.abs(previousEventTime - root) > convergence)) {
                        pendingEventTime = root;

View Full Code Here

    protected PointValuePair optimizeInternal(final int maxEval,
                                              final MultivariateDifferentiableFunction f,
                                              final GoalType goalType,
                                              final OptimizationData... optData) {
        // Store optimization problem characteristics.
        gradient = new GradientFunction(f);


        // Perform optimization.
        return super.optimizeInternal(maxEval, f, goalType, optData);
    }

View Full Code Here

                        final double baseRoot = forward ?
                                                solver.solve(maxIterationCount, f, ta, tb) :
                                                solver.solve(maxIterationCount, f, tb, ta);
                        final int remainingEval = maxIterationCount - solver.getEvaluations();
                        BracketedUnivariateSolver<UnivariateFunction> bracketing =
                                new PegasusSolver(solver.getRelativeAccuracy(), solver.getAbsoluteAccuracy());
                        root = forward ?
                               UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
                                                                   baseRoot, ta, tb, AllowedSolution.RIGHT_SIDE) :
                               UnivariateSolverUtils.forceSide(remainingEval, f, bracketing,
                                                                   baseRoot, tb, ta, AllowedSolution.LEFT_SIDE);

View Full Code Here


                // tests for termination and stringent tolerances
                if (FastMath.abs(actRed) <= TWO_EPS &&
                    preRed <= TWO_EPS &&
                    ratio <= 2.0) {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE,
                                                   costRelativeTolerance);
                } else if (delta <= TWO_EPS * xNorm) {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE,
                                                   parRelativeTolerance);
                } else if (maxCosine <= TWO_EPS) {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE,
                                                   orthoTolerance);
                }
            }
        }
    }

View Full Code Here

                for (int j = k; j < nR; ++j) {
                    double aki = weightedJacobian[j][permutation[i]];
                    norm2 += aki * aki;
                }
                if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
                    throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
                                                   nR, nC);
                }
                if (norm2 > ak2) {
                    nextColumn = i;
                    ak2        = norm2;

View Full Code Here

                }


                // tests for termination and stringent tolerances
                // (2.2204e-16 is the machine epsilon for IEEE754)
                if ((FastMath.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE,
                                                   costRelativeTolerance);
                } else if (delta <= 2.2204e-16 * xNorm) {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE,
                                                   parRelativeTolerance);
                } else if (maxCosine <= 2.2204e-16)  {
                    throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE,
                                                   orthoTolerance);
                }
            }
        }
    }

View Full Code Here

                for (int j = k; j < nR; ++j) {
                    double aki = weightedJacobian[j][permutation[i]];
                    norm2 += aki * aki;
                }
                if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
                    throw new ConvergenceException(LocalizedFormats.UNABLE_TO_PERFORM_QR_DECOMPOSITION_ON_JACOBIAN,
                                                   nR, nC);
                }
                if (norm2 > ak2) {
                    nextColumn = i;
                    ak2        = norm2;

View Full Code Here

     * length.
     */
    protected double[] computeResiduals(double[] objectiveValue) {
        final double[] target = getTarget();
        if (objectiveValue.length != target.length) {
            throw new DimensionMismatchException(target.length,
                                                 objectiveValue.length);
        }


        final double[] residuals = new double[target.length];
        for (int i = 0; i < target.length; i++) {

View Full Code Here


        /** {@inheritDoc} */
        public RealVector solve(final RealVector b) {
            final int m = lTData.length;
            if (b.getDimension() != m) {
                throw new DimensionMismatchException(b.getDimension(), m);
            }


            final double[] x = b.toArray();


            // Solve LY = b

View Full Code Here


        /** {@inheritDoc} */
        public RealMatrix solve(RealMatrix b) {
            final int m = lTData.length;
            if (b.getRowDimension() != m) {
                throw new DimensionMismatchException(b.getRowDimension(), m);
            }


            final int nColB = b.getColumnDimension();
            final double[][] x = b.getData();

View Full Code Here

0 1 2 3 4 5 6 7 8 9

TOP

Related Classes of org.apache.commons.math3.fitting.HarmonicFitter$ParameterGuesser

org.apache.commons.math3.analysis.differentiation.DerivativeStructure

org.apache.commons.math3.analysis.differentiation.MultivariateDifferentiableVectorFunction

org.apache.commons.math3.analysis.function.HarmonicOscillator

org.apache.commons.math3.analysis.function.Sin

org.apache.commons.math3.analysis.function.Sinc

org.apache.commons.math3.analysis.MultivariateFunction

org.apache.commons.math3.analysis.QuinticFunction

org.apache.commons.math3.analysis.solvers.BrentSolver

org.apache.commons.math3.analysis.UnivariateFunction

org.apache.commons.math3.complex.Complex

All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.