List of usage examples for org.apache.commons.math3.exception ConvergenceException ConvergenceException
public ConvergenceException()
From source file:edu.uchc.octane.GaussianFitAstigmatism.java
/** * calculate the Z coordinate based on astigmatism *//* w w w . j av a2 s .c om*/ void calculateZ() { if (calibration_ == null) { z_ = 0; return; } UnivariateFunction func = new UnivariateFunction() { @Override public double value(double z) { double sigmax = FastMath.sqrt(pvp_.getPoint()[3] / 2); double sigmay = FastMath.sqrt(pvp_.getPoint()[4] / 2); double vx = calibration_[0] + (z - calibration_[1]) * (z - calibration_[1]) * calibration_[2] - sigmax; double vy = calibration_[3] + (z - calibration_[4]) * (z - calibration_[4]) * calibration_[5] - sigmay; return vx * vx + vy * vy; } }; BrentOptimizer optimizer = new BrentOptimizer(1e-4, 1e-4); UnivariatePointValuePair upvp = optimizer.optimize(new UnivariateObjectiveFunction(func), GoalType.MINIMIZE, MaxEval.unlimited(), new SearchInterval(2 * z0min_ - z0max_, 2 * z0max_ - z0min_)); if (upvp.getValue() > errTol_) { throw (new ConvergenceException()); } z_ = upvp.getPoint(); }
From source file:com.itemanalysis.psychometrics.factoranalysis.GPArotation.java
/** * Conducts orthogonal rotation of factor loadings. * * @param A matrix of orthogonal factor loadings * @return a matrix of rotated factor loadings. * @throws ConvergenceException//from w ww . j ava 2s. c o m */ private RotationResults GPForth(RealMatrix A, boolean normalize, int maxIter, double eps) throws ConvergenceException { int ncol = A.getColumnDimension(); if (normalize) { //elementwise division by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value / W.getEntry(row, column); } }); } RealMatrix Tmat = new IdentityMatrix(ncol); double alpha = 1; RealMatrix L = A.multiply(Tmat); gpFunction.computeValues(L); double f = gpFunction.getValue(); RealMatrix VgQ = gpFunction.getGradient(); RealMatrix G = A.transpose().multiply(VgQ); double VgQtF = gpFunction.getValue(); RealMatrix VgQt = gpFunction.getGradient(); RealMatrix Tmatt = null; int iter = 0; double s = eps + 0.5; double s2 = 0; int innnerIter = 11; while (iter < maxIter) { RealMatrix M = Tmat.transpose().multiply(G); RealMatrix S = (M.add(M.transpose())); S = S.scalarMultiply(0.5); RealMatrix Gp = G.subtract(Tmat.multiply(S)); s = Math.sqrt((Gp.transpose().multiply(Gp)).getTrace()); s2 = Math.pow(s, 2); if (s < eps) break; alpha *= 2.0; for (int j = 0; j < innnerIter; j++) { Gp = Gp.scalarMultiply(alpha); RealMatrix X = (Tmat.subtract(Gp)); SingularValueDecomposition SVD = new SingularValueDecomposition(X); Tmatt = SVD.getU().multiply(SVD.getV().transpose()); L = A.multiply(Tmatt); gpFunction.computeValues(L); VgQt = gpFunction.getGradient(); VgQtF = gpFunction.getValue(); if (VgQtF < f - 0.5 * s2 * alpha) { break; } alpha /= 2.0; } Tmat = Tmatt; f = VgQtF; G = A.transpose().multiply(VgQt); iter++; } boolean convergence = s < eps; if (!convergence) { throw new ConvergenceException(); } if (normalize) { //elementwise multiplication by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value * W.getEntry(row, column); } }); } RealMatrix Phi = Tmat.transpose().multiply(Tmat); RotationResults result = new RotationResults(gpFunction.getValue(), L, Phi, Tmat, rotationMethod); return result; }
From source file:gedi.util.math.stat.distributions.NormalMixtureDistribution.java
public static NormalMixtureDistribution fit(NormalMixtureDistribution initialMixture, double[] data, final int maxIterations, final double threshold) { if (maxIterations < 1) { throw new NotStrictlyPositiveException(maxIterations); }/* w w w . j a v a 2s . c om*/ if (threshold < Double.MIN_VALUE) { throw new NotStrictlyPositiveException(threshold); } final int n = data.length; final int k = initialMixture.getNumComponents(); if (k == 1) return new NormalMixtureDistribution(new NormalDistribution[] { new NormalDistribution(new Mean().evaluate(data), new StandardDeviation().evaluate(data)) }, new double[] { 1 }); int numIterations = 0; double previousLogLikelihood = 0d; double logLikelihood = Double.NEGATIVE_INFINITY; // Initialize model to fit to initial mixture. NormalMixtureDistribution fittedModel = new NormalMixtureDistribution(initialMixture.components, initialMixture.mixing); while (numIterations++ <= maxIterations && FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { previousLogLikelihood = logLikelihood; logLikelihood = 0d; // E-step: compute the data dependent parameters of the expectation // function. // The percentage of row's total density between a row and a // component final double[][] gamma = new double[n][k]; // Sum of gamma for each component final double[] gammaSums = new double[k]; for (int i = 0; i < n; i++) { final double rowDensity = fittedModel.density(data[i]); logLikelihood += FastMath.log(rowDensity); for (int j = 0; j < k; j++) { gamma[i][j] = fittedModel.mixing[j] * fittedModel.components[j].density(data[i]) / rowDensity; gammaSums[j] += gamma[i][j]; } } logLikelihood /= n; // System.out.println(logLikelihood); // M-step: compute the new parameters based on the expectation // function. final double[] newWeights = gammaSums.clone(); ArrayUtils.mult(newWeights, 1.0 / n); NormalDistribution[] comp = new NormalDistribution[k]; for (int j = 0; j < k; j++) { double m = 0; for (int i = 0; i < n; i++) { m += gamma[i][j] * data[i]; } m /= gammaSums[j]; double var = 0; for (int i = 0; i < n; i++) { double d = m - data[i]; var += gamma[i][j] * d * d; } var /= gammaSums[j]; comp[j] = new NormalDistribution(m, Math.sqrt(var)); } // Update current model fittedModel = new NormalMixtureDistribution(comp, newWeights); } if (FastMath.abs(previousLogLikelihood - logLikelihood) > threshold) { // Did not converge before the maximum number of iterations throw new ConvergenceException(); } return fittedModel; }
From source file:de.tuberlin.uebb.jdae.llmsl.Block.java
public void exec() { views[0].loadD(point, variables);/*w w w . j a v a 2 s . co m*/ forceCompute(); int steps = options.maxIterations; while (!writeNegResidual(residual.data)) { if (steps-- <= 0) { views[0].loadD(point, variables); forceCompute(); writeJacobian(); System.out.println("Convergence Error:"); System.out.println(jacobianMatrix.toString()); throw new ConvergenceException(); } writeJacobian(); if (!solver.setA(jacobianMatrix)) throw new ConvergenceException(); solver.solve(residual, x); for (int i = 0; i < point.length; ++i) point[i] += x.data[i]; forceCompute(); } views[0].set(0, variables, point); }
From source file:com.itemanalysis.psychometrics.factoranalysis.GPArotation.java
private RotationResults GPFoblq(RealMatrix A, boolean normalize, int maxIter, double eps) throws ConvergenceException { int ncol = A.getColumnDimension(); RealMatrix Tinner = null;/*from ww w. j a v a 2s. co m*/ RealMatrix TinnerInv = null; RealMatrix Tmat = new IdentityMatrix(ncol); if (normalize) { //elementwise division by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value / W.getEntry(row, column); } }); } RealMatrix TmatInv = new LUDecomposition(Tmat).getSolver().getInverse(); RealMatrix L = A.multiply(TmatInv.transpose()); //compute gradientAt and function value gpFunction.computeValues(L); RealMatrix VgQ = gpFunction.getGradient(); RealMatrix VgQt = VgQ; double f = gpFunction.getValue(); double ft = f; RealMatrix G = ((L.transpose().multiply(VgQ).multiply(TmatInv)).transpose()).scalarMultiply(-1.0); int iter = 0; double alpha = 1.0; double s = eps + 0.5; double s2 = Math.pow(s, 2); int innerMaxIter = 10; int innerCount = 0; IdentityMatrix I = new IdentityMatrix(G.getRowDimension()); RealMatrix V1 = MatrixUtils.getVector(ncol, 1.0); while (iter < maxIter) { RealMatrix M = MatrixUtils.multiplyElements(Tmat, G); RealMatrix diagP = new DiagonalMatrix(V1.multiply(M).getRow(0)); RealMatrix Gp = G.subtract(Tmat.multiply(diagP)); s = Math.sqrt(Gp.transpose().multiply(Gp).getTrace()); s2 = Math.pow(s, 2); if (s < eps) { break; } alpha = 2.0 * alpha; innerCount = 0; for (int i = 0; i < innerMaxIter; i++) { RealMatrix X = Tmat.subtract(Gp.scalarMultiply(alpha)); RealMatrix X2 = MatrixUtils.multiplyElements(X, X); RealMatrix V = V1.multiply(X2); V.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return 1.0 / Math.sqrt(value); } }); //compute new value of T, its inverse, and the rotated loadings RealMatrix diagV = new DiagonalMatrix(V.getRow(0)); Tinner = X.multiply(diagV); TinnerInv = new LUDecomposition(Tinner).getSolver().getInverse(); L = A.multiply(TinnerInv.transpose()); //compute new values of the gradientAt and the rotation criteria gpFunction.computeValues(L); VgQt = gpFunction.getGradient(); ft = gpFunction.getValue(); innerCount++; if (ft < f - 0.5 * s2 * alpha) { break; } alpha = alpha / 2.0; } // System.out.println(iter + " " + f + " " + s + " " + Math.log10(s) + " " + alpha + " " + innerCount); Tmat = Tinner; f = ft; G = (L.transpose().multiply(VgQt).multiply(TinnerInv)).transpose().scalarMultiply(-1.0); iter++; } boolean convergence = s < eps; if (!convergence) { throw new ConvergenceException(); } if (normalize) { //elementwise multiplication by normalizing weights final RealMatrix W = getNormalizingWeights(A, true); A.walkInRowOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return value * W.getEntry(row, column); } }); } RealMatrix Phi = Tmat.transpose().multiply(Tmat); RotationResults result = new RotationResults(gpFunction.getValue(), L, Phi, Tmat, rotationMethod); return result; }
From source file:org.orekit.orbits.CartesianOrbit.java
/** Computes the eccentric in-plane argument from the mean in-plane argument. * @param thetaM = mean in-plane argument (rad) * @param ex first component of eccentricity vector * @param ey second component of eccentricity vector * @return the eccentric in-plane argument. *///from www . j av a 2 s. c o m private double meanToEccentric(final double thetaM, final double ex, final double ey) { // Generalization of Kepler equation to in-plane parameters // with thetaE = eta + E and // thetaM = eta + M = thetaE - ex.sin(thetaE) + ey.cos(thetaE) // and eta being counted from an arbitrary reference in the orbital plane double thetaE = thetaM; double thetaEMthetaM = 0.0; int iter = 0; do { final double cosThetaE = FastMath.cos(thetaE); final double sinThetaE = FastMath.sin(thetaE); final double f2 = ex * sinThetaE - ey * cosThetaE; final double f1 = 1.0 - ex * cosThetaE - ey * sinThetaE; final double f0 = thetaEMthetaM - f2; final double f12 = 2.0 * f1; final double shift = f0 * f12 / (f1 * f12 - f0 * f2); thetaEMthetaM -= shift; thetaE = thetaM + thetaEMthetaM; if (FastMath.abs(shift) <= 1.0e-12) { return thetaE; } } while (++iter < 50); throw new ConvergenceException(); }