Examples with PointVectorValuePair org.apache.commons.math3.optimization.PointVectorValuePair used on opensource projects

Example 1 with PointVectorValuePair

use of org.apache.commons.math3.optimization.PointVectorValuePair in project ffx by mjschnie.

the class ParticleMeshEwaldCart method leastSquaresPredictor.

/**
 * The least-squares predictor with induced dipole information from 8-10
 * previous steps reduces the number SCF iterations by ~50%.
 */
private void leastSquaresPredictor() {
    if (predictorCount < 2) {
        return;
    }
    try {
        /**
         * The Jacobian and target do not change during the LS optimization,
         * so it's most efficient to update them once before the
         * Least-Squares optimizer starts.
         */
        leastSquaresPredictor.updateJacobianAndTarget();
        int maxEvals = 100;
        fill(leastSquaresPredictor.initialSolution, 0.0);
        leastSquaresPredictor.initialSolution[0] = 1.0;
        PointVectorValuePair optimum = leastSquaresOptimizer.optimize(maxEvals, leastSquaresPredictor, leastSquaresPredictor.calculateTarget(), leastSquaresPredictor.weights, leastSquaresPredictor.initialSolution);
        double[] optimalValues = optimum.getPoint();
        if (logger.isLoggable(Level.FINEST)) {
            logger.finest(String.format("\n LS RMS:            %10.6f", leastSquaresOptimizer.getRMS()));
            logger.finest(String.format(" LS Iterations:     %10d", leastSquaresOptimizer.getEvaluations()));
            logger.finest(String.format(" Jacobian Evals:    %10d", leastSquaresOptimizer.getJacobianEvaluations()));
            logger.finest(String.format(" Chi Square:        %10.6f", leastSquaresOptimizer.getChiSquare()));
            logger.finest(String.format(" LS Coefficients"));
            for (int i = 0; i < predictorOrder - 1; i++) {
                logger.finest(String.format(" %2d  %10.6f", i + 1, optimalValues[i]));
            }
        }
        int mode;
        switch(lambdaMode) {
            case OFF:
            case CONDENSED:
                mode = 0;
                break;
            case CONDENSED_NO_LIGAND:
                mode = 1;
                break;
            case VAPOR:
                mode = 2;
                break;
            default:
                mode = 0;
        }
        /**
         * Initialize a pointer into predictor induced dipole array.
         */
        int index = predictorStartIndex;
        /**
         * Apply the LS coefficients in order to provide an initial guess at
         * the converged induced dipoles.
         */
        for (int k = 0; k < predictorOrder - 1; k++) {
            /**
             * Set the current coefficient.
             */
            double c = optimalValues[k];
            for (int i = 0; i < nAtoms; i++) {
                for (int j = 0; j < 3; j++) {
                    inducedDipole[0][i][j] += c * predictorInducedDipole[mode][index][i][j];
                    inducedDipoleCR[0][i][j] += c * predictorInducedDipoleCR[mode][index][i][j];
                }
            }
            index++;
            if (index >= predictorOrder) {
                index = 0;
            }
        }
    } catch (Exception e) {
        logger.log(Level.WARNING, " Exception computing predictor coefficients", e);
    }
}

Example 2 with PointVectorValuePair

use of org.apache.commons.math3.optimization.PointVectorValuePair in project dawnsci by DawnScience.

the class EllipseFitter method geometricFit.

@Override
public void geometricFit(IDataset x, IDataset y, double[] init) {
    residual = Double.NaN;
    if (x.getSize() < PARAMETERS || y.getSize() < PARAMETERS) {
        throw new IllegalArgumentException("Need " + PARAMETERS + " or more points");
    }
    if (init == null)
        init = quickfit(x, y);
    else if (init.length < PARAMETERS)
        throw new IllegalArgumentException("Need " + PARAMETERS + " parameters");
    EllipseCoordinatesFunction f = (EllipseCoordinatesFunction) getFitFunction(x, y);
    LevenbergMarquardtOptimizer opt = new LevenbergMarquardtOptimizer();
    try {
        PointVectorValuePair result = opt.optimize(new ModelFunction(f), new ModelFunctionJacobian(f.jacobian()), f.getTarget(), f.getWeight(), f.calcAllInitValues(init), new MaxEval(MAX_EVALUATIONS));
        double[] point = result.getPointRef();
        for (int i = 0; i < PARAMETERS; i++) parameters[i] = point[i];
        residual = opt.getRMS();
        logger.trace("Ellipse fit: rms = {}, x^2 = {}", residual, opt.getChiSquare());
    } catch (DimensionMismatchException e) {
    // cannot happen
    } catch (IllegalArgumentException e) {
    // should not happen!
    } catch (TooManyEvaluationsException e) {
        throw new IllegalArgumentException("Problem with optimizer converging");
    }
}

Example 3 with PointVectorValuePair

use of org.apache.commons.math3.optimization.PointVectorValuePair in project gdsc by aherbert.

the class CellOutliner_PlugIn method fitPolygonalCell.

/**
 * Find an ellipse that optimises the fit to the polygon detected edges.
 *
 * @param roi the roi
 * @param params the params
 * @param weightMap the weight map
 * @return the ellipse parameters
 */
@Nullable
private double[] fitPolygonalCell(PolygonRoi roi, double[] params, FloatProcessor weightMap) {
    // Get an estimate of the starting parameters using the current polygon
    final double[] startPoint = estimateStartPoint(roi, weightMap.getWidth(), weightMap.getHeight());
    final int maxEval = 2000;
    final DifferentiableEllipticalFitFunction func = new DifferentiableEllipticalFitFunction(roi, weightMap);
    final double relativeThreshold = 100 * Precision.EPSILON;
    final double absoluteThreshold = 100 * Precision.SAFE_MIN;
    final ConvergenceChecker<PointVectorValuePair> checker = new SimplePointChecker<>(relativeThreshold, absoluteThreshold);
    final double initialStepBoundFactor = 10;
    final double costRelativeTolerance = 1e-10;
    final double parRelativeTolerance = 1e-10;
    final double orthoTolerance = 1e-10;
    final double threshold = Precision.SAFE_MIN;
    final LevenbergMarquardtOptimizer optimiser = new LevenbergMarquardtOptimizer(initialStepBoundFactor, costRelativeTolerance, parRelativeTolerance, orthoTolerance, threshold);
    try {
        // @formatter:off
        final LeastSquaresProblem problem = new LeastSquaresBuilder().maxEvaluations(Integer.MAX_VALUE).maxIterations(maxEval).start(startPoint).target(func.calculateTarget()).weight(new DiagonalMatrix(func.calculateWeights())).model(func, func::jacobian).checkerPair(checker).build();
        // @formatter:on
        final Optimum solution = optimiser.optimize(problem);
        if (debug) {
            IJ.log(String.format("Eval = %d (Iter = %d), RMS = %f", solution.getEvaluations(), solution.getIterations(), solution.getRMS()));
        }
        return solution.getPoint().toArray();
    } catch (final Exception ex) {
        IJ.log("Failed to find an elliptical solution, defaulting to polygon: " + ex.getMessage());
    }
    return null;
}

Example 4 with PointVectorValuePair

use of org.apache.commons.math3.optimization.PointVectorValuePair in project s1tbx by senbox-org.

the class ArcDataIntegration method integrateArcsL2.

/**
 * Integrates arc data to node data using weighted least squares (L2 norm).
 *     This function integrates arc `data` by finding `point_data` (the returned
 *     value), such that `point_data` minimizes:
 *         ||diag(weights) * (A * point_data - data)||_2
 *     where `A` is an incidence matrix of dimensions 'number of arcs' by
 *     'number of nodes'.
 *     Description of the algorithm:
 *         1. Apply weights
 *             A = diag(weights) * A
 *             b = diag(weights) * data
 *         2. Minimize
 *             ||A*x - b||_2
 *
 * @param arcs every row should contain two indices corresponding to the end nodes of a specific arc.
 * @param data data to be integrated.
 * @param weights quality of arc data.
 * @return
 */
public static double[] integrateArcsL2(int[][] arcs, double[] data, double[] weights) {
    // graph size
    int noOfNodes = Arrays.stream(arcs).flatMapToInt(Arrays::stream).max().getAsInt() + 1;
    int noOfArcs = arcs.length;
    SystemUtils.LOG.fine("Number of nodes: " + noOfNodes);
    SystemUtils.LOG.fine("Number of arcs: " + noOfArcs);
    // A: coefficients matrix
    double[][] incidenceMatrix = new double[noOfArcs][noOfNodes];
    // get column and set to zero reference node occurrences
    int[] sourceNodes = getColumnVector(arcs, 0);
    // get column and set to zero reference node occurrences
    int[] targetNodes = getColumnVector(arcs, 1);
    // source nodes
    setValues(incidenceMatrix, repeat(-1.0, noOfArcs), null, sourceNodes);
    // target nodes
    setValues(incidenceMatrix, repeat(1.0, noOfArcs), null, targetNodes);
    // remove reference node column
    if (!referenceNodeIsLast) {
        // remove first column
        shiftColumnsLeft(incidenceMatrix, 1);
    }
    dropLastColumn(incidenceMatrix);
    // problem
    ModelFunction problem = new ModelFunction(new MultivariateVectorFunction() {

        RealMatrix jacobian = new Array2DRowRealMatrix(incidenceMatrix);

        public double[] value(double[] params) {
            double[] values = jacobian.operate(params);
            return values;
        }
    });
    ModelFunctionJacobian problemJacobian = new ModelFunctionJacobian(new MultivariateMatrixFunction() {

        double[][] jacobian = incidenceMatrix;

        public double[][] value(double[] params) {
            return jacobian;
        }
    });
    // optimization
    ConvergenceChecker<PointVectorValuePair> checker = new SimplePointChecker(0.00001, -1.0);
    LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer(checker);
    int maxIterations = 10000;
    PointVectorValuePair solution = optimizer.optimize(new MaxEval(maxIterations), problem, problemJacobian, new Target(data), new Weight(weights), new InitialGuess(new double[noOfNodes - 1]));
    // return solution
    double[] nodeData;
    if (referenceNodeIsLast) {
        nodeData = Arrays.copyOf(solution.getPoint(), noOfNodes);
        // add reference node back (last)
        nodeData[noOfNodes - 1] = 0;
    } else {
        nodeData = new double[noOfNodes];
        // reference node (first) is 0
        System.arraycopy(solution.getPoint(), 0, nodeData, 1, noOfNodes - 1);
    }
    return nodeData;
}

Example 5 with PointVectorValuePair

use of org.apache.commons.math3.optimization.PointVectorValuePair in project dawnsci by DawnScience.

the class CircleFitter method geometricFit.

@Override
public void geometricFit(IDataset x, IDataset y, double[] init) {
    residual = Double.NaN;
    if (x.getSize() < PARAMETERS || y.getSize() < PARAMETERS) {
        throw new IllegalArgumentException("Need " + PARAMETERS + " or more points");
    }
    if (x.getSize() == PARAMETERS) {
        init = quickfit(x, y);
        for (int i = 0; i < PARAMETERS; i++) parameters[i] = init[i];
        residual = 0;
        return;
    }
    if (init == null)
        init = quickfit(x, y);
    else if (init.length < PARAMETERS)
        throw new IllegalArgumentException("Need " + PARAMETERS + " parameters");
    CircleCoordinatesFunction f = (CircleCoordinatesFunction) getFitFunction(x, y);
    LevenbergMarquardtOptimizer opt = new LevenbergMarquardtOptimizer();
    try {
        PointVectorValuePair result = opt.optimize(new ModelFunction(f), new ModelFunctionJacobian(f.jacobian()), f.getTarget(), f.getWeight(), f.calcAllInitValues(init), new MaxEval(MAX_EVALUATIONS));
        double[] point = result.getPointRef();
        for (int i = 0; i < PARAMETERS; i++) parameters[i] = point[i];
        residual = opt.getRMS();
        logger.trace("Circle fit: rms = {}, x^2 = {}", opt.getRMS(), opt.getChiSquare());
    } catch (DimensionMismatchException e) {
    // cannot happen
    } catch (IllegalArgumentException e) {
    // should not happen!
    } catch (TooManyEvaluationsException e) {
        throw new IllegalArgumentException("Problem with optimizer converging");
    }
}