Author: erans Date: Fri Nov 9 14:30:49 2012 New Revision: 1407467 URL: http://svn.apache.org/viewvc?rev=1407467&view=rev Log: MATH-887 In "LevenbergMarquardtOptimizer", removed usage of deprecated fields and methods from its base class.
Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizer.java commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java?rev=1407467&r1=1407466&r2=1407467&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizer.java Fri Nov 9 14:30:49 2012 @@ -255,6 +255,15 @@ public abstract class AbstractLeastSquar } /** + * Gets the square-root of the weight matrix. + * + * @return the square-root of the weight matrix. + */ + public RealMatrix getWeightSquareRoot() { + return weightMatrixSqrt.copy(); + } + + /** * Sets the cost. * * @param cost Cost value. Modified: commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizer.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizer.java?rev=1407467&r1=1407466&r2=1407467&view=diff ============================================================================== --- commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizer.java (original) +++ commons/proper/math/trunk/src/main/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizer.java Fri Nov 9 14:30:49 2012 @@ -22,6 +22,7 @@ import org.apache.commons.math3.exceptio import org.apache.commons.math3.exception.util.LocalizedFormats; import org.apache.commons.math3.optimization.PointVectorValuePair; import org.apache.commons.math3.optimization.ConvergenceChecker; +import org.apache.commons.math3.linear.RealMatrix; import org.apache.commons.math3.util.Precision; import org.apache.commons.math3.util.FastMath; @@ -132,6 +133,10 @@ public class LevenbergMarquardtOptimizer private final double orthoTolerance; /** Threshold for QR ranking. */ private final double qrRankingThreshold; + /** Weighted residuals. */ + private double[] weightedResidual; + /** Weighted Jacobian. */ + private double[][] weightedJacobian; /** * Build an optimizer for least squares problems with default values @@ -271,66 +276,75 @@ public class LevenbergMarquardtOptimizer /** {@inheritDoc} */ @Override protected PointVectorValuePair doOptimize() { + final int nR = getTarget().length; // Number of observed data. + final double[] currentPoint = getStartPoint(); + final int nC = currentPoint.length; // Number of parameters. + // arrays shared with the other private methods - solvedCols = FastMath.min(rows, cols); - diagR = new double[cols]; - jacNorm = new double[cols]; - beta = new double[cols]; - permutation = new int[cols]; - lmDir = new double[cols]; + solvedCols = FastMath.min(nR, nC); + diagR = new double[nC]; + jacNorm = new double[nC]; + beta = new double[nC]; + permutation = new int[nC]; + lmDir = new double[nC]; // local point double delta = 0; double xNorm = 0; - double[] diag = new double[cols]; - double[] oldX = new double[cols]; - double[] oldRes = new double[rows]; - double[] oldObj = new double[rows]; - double[] qtf = new double[rows]; - double[] work1 = new double[cols]; - double[] work2 = new double[cols]; - double[] work3 = new double[cols]; - - // evaluate the function at the starting point and calculate its norm - updateResidualsAndCost(); + double[] diag = new double[nC]; + double[] oldX = new double[nC]; + double[] oldRes = new double[nR]; + double[] oldObj = new double[nR]; + double[] qtf = new double[nR]; + double[] work1 = new double[nC]; + double[] work2 = new double[nC]; + double[] work3 = new double[nC]; + + final RealMatrix weightMatrixSqrt = getWeightSquareRoot(); + + // Evaluate the function at the starting point and calculate its norm. + double[] currentObjective = computeObjectiveValue(currentPoint); + double[] currentResiduals = computeResiduals(currentObjective); + PointVectorValuePair current = new PointVectorValuePair(currentPoint, currentObjective); + double currentCost = computeCost(currentResiduals); - // outer loop + // Outer loop. lmPar = 0; boolean firstIteration = true; - PointVectorValuePair current = new PointVectorValuePair(point, objective); int iter = 0; final ConvergenceChecker<PointVectorValuePair> checker = getConvergenceChecker(); while (true) { ++iter; + final PointVectorValuePair previous = current; - for (int i=0;i<rows;i++) { - qtf[i]=weightedResiduals[i]; - } + // QR decomposition of the jacobian matrix + qrDecomposition(computeJacobian(currentPoint)); - // compute the Q.R. decomposition of the jacobian matrix - PointVectorValuePair previous = current; - updateJacobian(); - qrDecomposition(); + weightedResidual = weightMatrixSqrt.operate(currentResiduals); + for (int i = 0; i < nR; i++) { + qtf[i] = weightedResidual[i]; + } // compute Qt.res qTy(qtf); + // now we don't need Q anymore, // so let jacobian contain the R matrix with its diagonal elements for (int k = 0; k < solvedCols; ++k) { int pk = permutation[k]; - weightedResidualJacobian[k][pk] = diagR[pk]; + weightedJacobian[k][pk] = diagR[pk]; } if (firstIteration) { // scale the point according to the norms of the columns // of the initial jacobian xNorm = 0; - for (int k = 0; k < cols; ++k) { + for (int k = 0; k < nC; ++k) { double dk = jacNorm[k]; if (dk == 0) { dk = 1.0; } - double xk = dk * point[k]; + double xk = dk * currentPoint[k]; xNorm += xk * xk; diag[k] = dk; } @@ -342,45 +356,44 @@ public class LevenbergMarquardtOptimizer // check orthogonality between function vector and jacobian columns double maxCosine = 0; - if (cost != 0) { + if (currentCost != 0) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double s = jacNorm[pj]; if (s != 0) { double sum = 0; for (int i = 0; i <= j; ++i) { - sum += weightedResidualJacobian[i][pj] * qtf[i]; + sum += weightedJacobian[i][pj] * qtf[i]; } - maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * cost)); + maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * currentCost)); } } } if (maxCosine <= orthoTolerance) { - // convergence has been reached - updateResidualsAndCost(); - current = new PointVectorValuePair(point, objective); + // Convergence has been reached. + setCost(currentCost); return current; } // rescale if necessary - for (int j = 0; j < cols; ++j) { + for (int j = 0; j < nC; ++j) { diag[j] = FastMath.max(diag[j], jacNorm[j]); } - // inner loop + // Inner loop. for (double ratio = 0; ratio < 1.0e-4;) { // save the state for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; - oldX[pj] = point[pj]; + oldX[pj] = currentPoint[pj]; } - double previousCost = cost; - double[] tmpVec = weightedResiduals; - weightedResiduals = oldRes; + final double previousCost = currentCost; + double[] tmpVec = weightedResidual; + weightedResidual = oldRes; oldRes = tmpVec; - tmpVec = objective; - objective = oldObj; + tmpVec = currentObjective; + currentObjective = oldObj; oldObj = tmpVec; // determine the Levenberg-Marquardt parameter @@ -391,7 +404,7 @@ public class LevenbergMarquardtOptimizer for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; lmDir[pj] = -lmDir[pj]; - point[pj] = oldX[pj] + lmDir[pj]; + currentPoint[pj] = oldX[pj] + lmDir[pj]; double s = diag[pj] * lmDir[pj]; lmNorm += s * s; } @@ -401,13 +414,16 @@ public class LevenbergMarquardtOptimizer delta = FastMath.min(delta, lmNorm); } - // evaluate the function at x + p and calculate its norm - updateResidualsAndCost(); + // Evaluate the function at x + p and calculate its norm. + currentObjective = computeObjectiveValue(currentPoint); + currentResiduals = computeResiduals(currentObjective); + current = new PointVectorValuePair(currentPoint, currentObjective); + currentCost = computeCost(currentResiduals); // compute the scaled actual reduction double actRed = -1.0; - if (0.1 * cost < previousCost) { - double r = cost / previousCost; + if (0.1 * currentCost < previousCost) { + double r = currentCost / previousCost; actRed = 1.0 - r * r; } @@ -418,7 +434,7 @@ public class LevenbergMarquardtOptimizer double dirJ = lmDir[pj]; work1[j] = 0; for (int i = 0; i <= j; ++i) { - work1[i] += weightedResidualJacobian[i][pj] * dirJ; + work1[i] += weightedJacobian[i][pj] * dirJ; } } double coeff1 = 0; @@ -438,7 +454,7 @@ public class LevenbergMarquardtOptimizer if (ratio <= 0.25) { double tmp = (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; - if ((0.1 * cost >= previousCost) || (tmp < 0.1)) { + if ((0.1 * currentCost >= previousCost) || (tmp < 0.1)) { tmp = 0.1; } delta = tmp * FastMath.min(delta, 10.0 * lmNorm); @@ -453,33 +469,35 @@ public class LevenbergMarquardtOptimizer // successful iteration, update the norm firstIteration = false; xNorm = 0; - for (int k = 0; k < cols; ++k) { - double xK = diag[k] * point[k]; + for (int k = 0; k < nC; ++k) { + double xK = diag[k] * currentPoint[k]; xNorm += xK * xK; } xNorm = FastMath.sqrt(xNorm); - current = new PointVectorValuePair(point, objective); // tests for convergence. if (checker != null) { // we use the vectorial convergence checker if (checker.converged(iter, previous, current)) { + setCost(currentCost); return current; } } } else { // failed iteration, reset the previous values - cost = previousCost; + currentCost = previousCost; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; - point[pj] = oldX[pj]; + currentPoint[pj] = oldX[pj]; } - tmpVec = weightedResiduals; - weightedResiduals = oldRes; + tmpVec = weightedResidual; + weightedResidual = oldRes; oldRes = tmpVec; - tmpVec = objective; - objective = oldObj; + tmpVec = currentObjective; + currentObjective = oldObj; oldObj = tmpVec; + // Reset "current" to previous values. + current = new PointVectorValuePair(currentPoint, currentObjective); } // Default convergence criteria. @@ -487,6 +505,7 @@ public class LevenbergMarquardtOptimizer preRed <= costRelativeTolerance && ratio <= 2.0) || delta <= parRelativeTolerance * xNorm) { + setCost(currentCost); return current; } @@ -494,13 +513,13 @@ public class LevenbergMarquardtOptimizer // (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); + costRelativeTolerance); } else if (delta <= 2.2204e-16 * xNorm) { throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE, - parRelativeTolerance); + parRelativeTolerance); } else if (maxCosine <= 2.2204e-16) { throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE, - orthoTolerance); + orthoTolerance); } } } @@ -529,21 +548,22 @@ public class LevenbergMarquardtOptimizer * @param work3 work array */ private void determineLMParameter(double[] qy, double delta, double[] diag, - double[] work1, double[] work2, double[] work3) { + double[] work1, double[] work2, double[] work3) { + final int nC = weightedJacobian[0].length; // compute and store in x the gauss-newton direction, if the // jacobian is rank-deficient, obtain a least squares solution for (int j = 0; j < rank; ++j) { lmDir[permutation[j]] = qy[j]; } - for (int j = rank; j < cols; ++j) { + for (int j = rank; j < nC; ++j) { lmDir[permutation[j]] = 0; } for (int k = rank - 1; k >= 0; --k) { int pk = permutation[k]; double ypk = lmDir[pk] / diagR[pk]; for (int i = 0; i < k; ++i) { - lmDir[permutation[i]] -= ypk * weightedResidualJacobian[i][pk]; + lmDir[permutation[i]] -= ypk * weightedJacobian[i][pk]; } lmDir[pk] = ypk; } @@ -579,7 +599,7 @@ public class LevenbergMarquardtOptimizer int pj = permutation[j]; double sum = 0; for (int i = 0; i < j; ++i) { - sum += weightedResidualJacobian[i][pj] * work1[permutation[i]]; + sum += weightedJacobian[i][pj] * work1[permutation[i]]; } double s = (work1[pj] - sum) / diagR[pj]; work1[pj] = s; @@ -594,7 +614,7 @@ public class LevenbergMarquardtOptimizer int pj = permutation[j]; double sum = 0; for (int i = 0; i <= j; ++i) { - sum += weightedResidualJacobian[i][pj] * qy[i]; + sum += weightedJacobian[i][pj] * qy[i]; } sum /= diag[pj]; sum2 += sum * sum; @@ -654,7 +674,7 @@ public class LevenbergMarquardtOptimizer work1[pj] /= work2[j]; double tmp = work1[pj]; for (int i = j + 1; i < solvedCols; ++i) { - work1[permutation[i]] -= weightedResidualJacobian[i][pj] * tmp; + work1[permutation[i]] -= weightedJacobian[i][pj] * tmp; } } sum2 = 0; @@ -698,14 +718,14 @@ public class LevenbergMarquardtOptimizer * @param work work array */ private void determineLMDirection(double[] qy, double[] diag, - double[] lmDiag, double[] work) { + double[] lmDiag, double[] work) { // copy R and Qty to preserve input and initialize s // in particular, save the diagonal elements of R in lmDir for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; for (int i = j + 1; i < solvedCols; ++i) { - weightedResidualJacobian[i][pj] = weightedResidualJacobian[j][permutation[i]]; + weightedJacobian[i][pj] = weightedJacobian[j][permutation[i]]; } lmDir[j] = diagR[pj]; work[j] = qy[j]; @@ -736,7 +756,7 @@ public class LevenbergMarquardtOptimizer final double sin; final double cos; - double rkk = weightedResidualJacobian[k][pk]; + double rkk = weightedJacobian[k][pk]; if (FastMath.abs(rkk) < FastMath.abs(lmDiag[k])) { final double cotan = rkk / lmDiag[k]; sin = 1.0 / FastMath.sqrt(1.0 + cotan * cotan); @@ -749,25 +769,25 @@ public class LevenbergMarquardtOptimizer // compute the modified diagonal element of R and // the modified element of (Qty,0) - weightedResidualJacobian[k][pk] = cos * rkk + sin * lmDiag[k]; + weightedJacobian[k][pk] = cos * rkk + sin * lmDiag[k]; final double temp = cos * work[k] + sin * qtbpj; qtbpj = -sin * work[k] + cos * qtbpj; work[k] = temp; // accumulate the tranformation in the row of s for (int i = k + 1; i < solvedCols; ++i) { - double rik = weightedResidualJacobian[i][pk]; + double rik = weightedJacobian[i][pk]; final double temp2 = cos * rik + sin * lmDiag[i]; lmDiag[i] = -sin * rik + cos * lmDiag[i]; - weightedResidualJacobian[i][pk] = temp2; + weightedJacobian[i][pk] = temp2; } } } // store the diagonal element of s and restore // the corresponding diagonal element of R - lmDiag[j] = weightedResidualJacobian[j][permutation[j]]; - weightedResidualJacobian[j][permutation[j]] = lmDir[j]; + lmDiag[j] = weightedJacobian[j][permutation[j]]; + weightedJacobian[j][permutation[j]] = lmDir[j]; } // solve the triangular system for z, if the system is @@ -786,7 +806,7 @@ public class LevenbergMarquardtOptimizer int pj = permutation[j]; double sum = 0; for (int i = j + 1; i < nSing; ++i) { - sum += weightedResidualJacobian[i][pj] * work[i]; + sum += weightedJacobian[i][pj] * work[i]; } work[j] = (work[j] - sum) / lmDiag[j]; } @@ -818,36 +838,41 @@ public class LevenbergMarquardtOptimizer * are performed in non-increasing columns norms order thanks to columns * pivoting. The diagonal elements of the R matrix are therefore also in * non-increasing absolute values order.</p> + * + * @param jacobian Weighte Jacobian matrix at the current point. * @exception ConvergenceException if the decomposition cannot be performed */ - private void qrDecomposition() throws ConvergenceException { + private void qrDecomposition(RealMatrix jacobian) throws ConvergenceException { + weightedJacobian = jacobian.getData(); + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; // initializations - for (int k = 0; k < cols; ++k) { + for (int k = 0; k < nC; ++k) { permutation[k] = k; double norm2 = 0; - for (int i = 0; i < weightedResidualJacobian.length; ++i) { - double akk = weightedResidualJacobian[i][k]; + for (int i = 0; i < nR; ++i) { + double akk = weightedJacobian[i][k]; norm2 += akk * akk; } jacNorm[k] = FastMath.sqrt(norm2); } // transform the matrix column after column - for (int k = 0; k < cols; ++k) { + for (int k = 0; k < nC; ++k) { // select the column with the greatest norm on active components int nextColumn = -1; double ak2 = Double.NEGATIVE_INFINITY; - for (int i = k; i < cols; ++i) { + for (int i = k; i < nC; ++i) { double norm2 = 0; - for (int j = k; j < weightedResidualJacobian.length; ++j) { - double aki = weightedResidualJacobian[j][permutation[i]]; + 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, - rows, cols); + nR, nC); } if (norm2 > ak2) { nextColumn = i; @@ -863,24 +888,24 @@ public class LevenbergMarquardtOptimizer permutation[k] = pk; // choose alpha such that Hk.u = alpha ek - double akk = weightedResidualJacobian[k][pk]; + double akk = weightedJacobian[k][pk]; double alpha = (akk > 0) ? -FastMath.sqrt(ak2) : FastMath.sqrt(ak2); double betak = 1.0 / (ak2 - akk * alpha); beta[pk] = betak; // transform the current column diagR[pk] = alpha; - weightedResidualJacobian[k][pk] -= alpha; + weightedJacobian[k][pk] -= alpha; // transform the remaining columns - for (int dk = cols - 1 - k; dk > 0; --dk) { + for (int dk = nC - 1 - k; dk > 0; --dk) { double gamma = 0; - for (int j = k; j < weightedResidualJacobian.length; ++j) { - gamma += weightedResidualJacobian[j][pk] * weightedResidualJacobian[j][permutation[k + dk]]; + for (int j = k; j < nR; ++j) { + gamma += weightedJacobian[j][pk] * weightedJacobian[j][permutation[k + dk]]; } gamma *= betak; - for (int j = k; j < weightedResidualJacobian.length; ++j) { - weightedResidualJacobian[j][permutation[k + dk]] -= gamma * weightedResidualJacobian[j][pk]; + for (int j = k; j < nR; ++j) { + weightedJacobian[j][permutation[k + dk]] -= gamma * weightedJacobian[j][pk]; } } } @@ -893,15 +918,18 @@ public class LevenbergMarquardtOptimizer * @param y vector to multiply (will be overwritten with the result) */ private void qTy(double[] y) { - for (int k = 0; k < cols; ++k) { + final int nR = weightedJacobian.length; + final int nC = weightedJacobian[0].length; + + for (int k = 0; k < nC; ++k) { int pk = permutation[k]; double gamma = 0; - for (int i = k; i < rows; ++i) { - gamma += weightedResidualJacobian[i][pk] * y[i]; + for (int i = k; i < nR; ++i) { + gamma += weightedJacobian[i][pk] * y[i]; } gamma *= beta[pk]; - for (int i = k; i < rows; ++i) { - y[i] -= gamma * weightedResidualJacobian[i][pk]; + for (int i = k; i < nR; ++i) { + y[i] -= gamma * weightedJacobian[i][pk]; } } } Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java?rev=1407467&r1=1407466&r2=1407467&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/AbstractLeastSquaresOptimizerAbstractTest.java Fri Nov 9 14:30:49 2012 @@ -353,9 +353,9 @@ public abstract class AbstractLeastSquar circle.addPoint( 35.0, 15.0); circle.addPoint( 45.0, 97.0); AbstractLeastSquaresOptimizer optimizer = createOptimizer(); - PointVectorValuePair optimum = - optimizer.optimize(100, circle, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 }, - new double[] { 98.680, 47.345 }); + PointVectorValuePair optimum + = optimizer.optimize(100, circle, new double[] { 0, 0, 0, 0, 0 }, new double[] { 1, 1, 1, 1, 1 }, + new double[] { 98.680, 47.345 }); Assert.assertTrue(optimizer.getEvaluations() < 10); Assert.assertTrue(optimizer.getJacobianEvaluations() < 10); double rms = optimizer.getRMS(); @@ -399,8 +399,8 @@ public abstract class AbstractLeastSquar circle.addPoint(points[i][0], points[i][1]); } AbstractLeastSquaresOptimizer optimizer = createOptimizer(); - PointVectorValuePair optimum = - optimizer.optimize(100, circle, target, weights, new double[] { -12, -12 }); + PointVectorValuePair optimum + = optimizer.optimize(100, circle, target, weights, new double[] { -12, -12 }); Vector2D center = new Vector2D(optimum.getPointRef()[0], optimum.getPointRef()[1]); Assert.assertTrue(optimizer.getEvaluations() < 25); Assert.assertTrue(optimizer.getJacobianEvaluations() < 20); Modified: commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java URL: http://svn.apache.org/viewvc/commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java?rev=1407467&r1=1407466&r2=1407467&view=diff ============================================================================== --- commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java (original) +++ commons/proper/math/trunk/src/test/java/org/apache/commons/math3/optimization/general/LevenbergMarquardtOptimizerTest.java Fri Nov 9 14:30:49 2012 @@ -118,10 +118,10 @@ public class LevenbergMarquardtOptimizer }, new double[] { 1, 1, 1 }); AbstractLeastSquaresOptimizer optimizer = createOptimizer(); - optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 }); + PointVectorValuePair optimum = optimizer.optimize(100, problem, problem.target, new double[] { 1, 1, 1 }, new double[] { 0, 0, 0 }); Assert.assertTrue(FastMath.sqrt(problem.target.length) * optimizer.getRMS() > 0.6); - optimizer.getCovariances(1.5e-14); + optimizer.computeCovariances(optimum.getPoint(), 1.5e-14); } @Test @@ -223,14 +223,14 @@ public class LevenbergMarquardtOptimizer final LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(); - final PointVectorValuePair optimum = - optimizer.optimize(100, problem, dataPoints[1], weights, + final PointVectorValuePair optimum + = optimizer.optimize(100, problem, dataPoints[1], weights, new double[] { 10, 900, 80, 27, 225 }); final double[] solution = optimum.getPoint(); final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 }; - final double[][] covarMatrix = optimizer.getCovariances(); + final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14); final double[][] expectedCovarMatrix = { { 3.38, -3.69, 27.98, -2.34, -49.24 }, { -3.69, 2492.26, 81.89, -69.21, -8.9 }, @@ -292,7 +292,7 @@ public class LevenbergMarquardtOptimizer final double[] paramFound = optimum.getPoint(); // Retrieve errors estimation. - final double[][] covMatrix = optimizer.getCovariances(); + final double[][] covMatrix = optimizer.computeCovariances(paramFound, 1e-14); final double[] asymptoticStandardErrorFound = optimizer.guessParametersErrors(); final double[] sigmaFound = new double[covMatrix.length]; for (int i = 0; i < covMatrix.length; i++) {