Examples with Gaussian org.apache.commons.math3.analysis.function.Gaussian used on opensource projects

Example 36 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project vcell by virtualcell.

the class FitBleachSpotOp method fitToGaussian.

static GaussianFitResults fitToGaussian(double init_center_i, double init_center_j, double init_radius2, FloatImage image) {
    // 
    // do some optimization on the image (fitting to a Gaussian)
    // set initial guesses from ROI operation.
    // 
    ISize imageSize = image.getISize();
    final int num_i = imageSize.getX();
    final int num_j = imageSize.getY();
    final float[] floatPixels = image.getFloatPixels();
    // 
    // initial guess based on previous fit of ROI
    // do gaussian fit in index space for center and standard deviation (later to translate it back to world coordinates)
    // 
    final int window_size = (int) Math.sqrt(init_radius2) * 4;
    // final int window_min_i = 0;       // (int) Math.max(0, Math.floor(init_center_i - window_size/2));
    // final int window_max_i = num_i-1; // (int) Math.min(num_i-1, Math.ceil(init_center_i + window_size/2));
    // final int window_min_j = 0;       // (int) Math.max(0, Math.floor(init_center_j - window_size/2));
    // final int window_max_j = num_j-1; // (int) Math.min(num_j-1, Math.ceil(init_center_j + window_size/2));
    final int window_min_i = (int) Math.max(0, Math.floor(init_center_i - window_size / 2));
    final int window_max_i = (int) Math.min(num_i - 1, Math.ceil(init_center_i + window_size / 2));
    final int window_min_j = (int) Math.max(0, Math.floor(init_center_j - window_size / 2));
    final int window_max_j = (int) Math.min(num_j - 1, Math.ceil(init_center_j + window_size / 2));
    final int PARAM_INDEX_CENTER_I = 0;
    final int PARAM_INDEX_CENTER_J = 1;
    final int PARAM_INDEX_K = 2;
    final int PARAM_INDEX_HIGH = 3;
    final int PARAM_INDEX_RADIUS_SQUARED = 4;
    final int NUM_PARAMETERS = 5;
    double[] initParameters = new double[NUM_PARAMETERS];
    initParameters[PARAM_INDEX_CENTER_I] = init_center_i;
    initParameters[PARAM_INDEX_CENTER_J] = init_center_j;
    initParameters[PARAM_INDEX_HIGH] = 1.0;
    initParameters[PARAM_INDEX_K] = 10;
    initParameters[PARAM_INDEX_RADIUS_SQUARED] = init_radius2;
    PowellOptimizer optimizer = new PowellOptimizer(1e-4, 1e-1);
    MultivariateFunction func = new MultivariateFunction() {

        @Override
        public double value(double[] point) {
            double center_i = point[PARAM_INDEX_CENTER_I];
            double center_j = point[PARAM_INDEX_CENTER_J];
            double high = point[PARAM_INDEX_HIGH];
            double K = point[PARAM_INDEX_K];
            double radius2 = point[PARAM_INDEX_RADIUS_SQUARED];
            double error2 = 0;
            for (int j = window_min_j; j <= window_max_j; j++) {
                // double y = j - center_j;
                double y = j;
                for (int i = window_min_i; i <= window_max_i; i++) {
                    // double x = i - center_i;
                    double x = i;
                    double modelValue = high - FastMath.exp(-K * FastMath.exp(-2 * (x * x + y * y) / radius2));
                    double imageValue = floatPixels[j * num_i + i];
                    double error = modelValue - imageValue;
                    error2 += error * error;
                }
            }
            System.out.println(new GaussianFitResults(center_i, center_j, radius2, K, high, error2));
            return error2;
        }
    };
    PointValuePair pvp = optimizer.optimize(new ObjectiveFunction(func), new InitialGuess(initParameters), new MaxEval(100000), GoalType.MINIMIZE);
    double[] fittedParamValues = pvp.getPoint();
    double fitted_center_i = fittedParamValues[PARAM_INDEX_CENTER_I];
    double fitted_center_j = fittedParamValues[PARAM_INDEX_CENTER_J];
    double fitted_radius2 = fittedParamValues[PARAM_INDEX_RADIUS_SQUARED];
    double fitted_K = fittedParamValues[PARAM_INDEX_K];
    double fitted_high = fittedParamValues[PARAM_INDEX_HIGH];
    double objectiveFunctionValue = pvp.getValue();
    return new GaussianFitResults(fitted_center_i, fitted_center_j, fitted_radius2, fitted_K, fitted_high, objectiveFunctionValue);
}

Example 37 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project narchy by automenta.

the class MyCMAESOptimizer method updateCovariance.

/**
 * Update of the covariance matrix C.
 *
 * @param hsig Flag indicating a small correction.
 * @param bestArx Fitness-sorted matrix of the argument vectors producing the
 * current offspring.
 * @param arz Unsorted matrix containing the gaussian random values of the
 * current offspring.
 * @param arindex Indices indicating the fitness-order of the current offspring.
 * @param xold xmean matrix of the previous generation.
 */
private void updateCovariance(boolean hsig, final RealMatrix bestArx, final RealMatrix arz, final int[] arindex, final RealMatrix xold) {
    double negccov = 0;
    if (ccov1 + ccovmu > 0) {
        final RealMatrix arpos = bestArx.subtract(repmat(xold, 1, mu)).scalarMultiply(// mu difference vectors
        1 / sigma);
        final RealMatrix roneu = pc.multiply(pc.transpose()).scalarMultiply(// rank one update
        ccov1);
        // minor correction if hsig==false
        double oldFac = hsig ? 0 : ccov1 * cc * (2 - cc);
        oldFac += 1 - ccov1 - ccovmu;
        if (isActiveCMA) {
            // Adapt covariance matrix C active CMA
            negccov = (1 - ccovmu) * 0.25 * mueff / (Math.pow(dimension + 2, 1.5) + 2 * mueff);
            // keep at least 0.66 in all directions, small popsize are most
            // critical
            final double negminresidualvariance = 0.66;
            // where to make up for the variance loss
            final double negalphaold = 0.5;
            // prepare vectors, compute negative updating matrix Cneg
            final int[] arReverseIndex = reverse(arindex);
            RealMatrix arzneg = selectColumns(arz, MathArrays.copyOf(arReverseIndex, mu));
            RealMatrix arnorms = sqrt(sumRows(square(arzneg)));
            final int[] idxnorms = sortedIndices(arnorms.getRow(0));
            final RealMatrix arnormsSorted = selectColumns(arnorms, idxnorms);
            final int[] idxReverse = reverse(idxnorms);
            final RealMatrix arnormsReverse = selectColumns(arnorms, idxReverse);
            arnorms = divide(arnormsReverse, arnormsSorted);
            final int[] idxInv = inverse(idxnorms);
            final RealMatrix arnormsInv = selectColumns(arnorms, idxInv);
            // check and set learning rate negccov
            final double negcovMax = (1 - negminresidualvariance) / square(arnormsInv).multiply(weights).getEntry(0, 0);
            if (negccov > negcovMax) {
                negccov = negcovMax;
            }
            arzneg = times(arzneg, repmat(arnormsInv, dimension, 1));
            final RealMatrix artmp = BD.multiply(arzneg);
            final RealMatrix Cneg = artmp.multiply(diag(weights)).multiply(artmp.transpose());
            oldFac += negalphaold * negccov;
            C = C.scalarMultiply(oldFac).add(// regard old matrix
            roneu).add(// plus rank one update
            arpos.scalarMultiply(// plus rank mu update
            ccovmu + (1 - negalphaold) * negccov).multiply(times(repmat(weights, 1, dimension), arpos.transpose()))).subtract(Cneg.scalarMultiply(negccov));
        } else {
            // Adapt covariance matrix C - nonactive
            C = // regard old matrix
            C.scalarMultiply(oldFac).add(// plus rank one update
            roneu).add(// plus rank mu update
            arpos.scalarMultiply(ccovmu).multiply(times(repmat(weights, 1, dimension), arpos.transpose())));
        }
    }
    updateBD(negccov);
}