Examples with RandomGenerator org.apache.commons.math3.random.RandomGenerator used on opensource projects

Example 46 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project GDSC-SMLM by aherbert.

the class PulseActivationAnalysis method simulateMolecules.

private float[][] simulateMolecules(RandomDataGenerator rdg, int c) {
    int n = sim_nMolecules[c];
    float[][] molecules = new float[n][];
    if (n == 0)
        return molecules;
    // Draw the shapes
    Shape[] shapes = createShapes(rdg, c);
    // Sample positions from within the shapes
    boolean canSample = shapes[0].canSample();
    RandomGenerator rand = rdg.getRandomGenerator();
    while (n-- > 0) {
        float[] coords;
        if (canSample) {
            int next = rand.nextInt(shapes.length);
            coords = shapes[next].sample(rand);
        } else {
            coords = shapes[n % shapes.length].getPosition();
        }
        // Avoid out-of-bounds positions
        if (outOfBounds(coords[0]) || outOfBounds(coords[1]))
            n++;
        else
            molecules[n] = coords;
    }
    return molecules;
}

Example 47 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project GDSC-SMLM by aherbert.

the class PCPALMFitting method runBoundedOptimiser.

private PointValuePair runBoundedOptimiser(double[][] gr, double[] initialSolution, double[] lB, double[] uB, SumOfSquaresModelFunction function) {
    // Create the functions to optimise
    ObjectiveFunction objective = new ObjectiveFunction(new SumOfSquaresMultivariateFunction(function));
    ObjectiveFunctionGradient gradient = new ObjectiveFunctionGradient(new SumOfSquaresMultivariateVectorFunction(function));
    final boolean debug = false;
    // Try a BFGS optimiser since this will produce a deterministic solution and can respect bounds.
    PointValuePair optimum = null;
    boundedEvaluations = 0;
    final MaxEval maxEvaluations = new MaxEval(2000);
    MultivariateOptimizer opt = null;
    for (int iteration = 0; iteration <= fitRestarts; iteration++) {
        try {
            opt = new BFGSOptimizer();
            final double relativeThreshold = 1e-6;
            // Configure maximum step length for each dimension using the bounds
            double[] stepLength = new double[lB.length];
            for (int i = 0; i < stepLength.length; i++) stepLength[i] = (uB[i] - lB[i]) * 0.3333333;
            // The GoalType is always minimise so no need to pass this in
            optimum = opt.optimize(maxEvaluations, gradient, objective, new InitialGuess((optimum == null) ? initialSolution : optimum.getPointRef()), new SimpleBounds(lB, uB), new BFGSOptimizer.GradientTolerance(relativeThreshold), new BFGSOptimizer.StepLength(stepLength));
            if (debug)
                System.out.printf("BFGS Iter %d = %g (%d)\n", iteration, optimum.getValue(), opt.getEvaluations());
        } catch (TooManyEvaluationsException e) {
            // No need to restart
            break;
        } catch (RuntimeException e) {
            // No need to restart
            break;
        } finally {
            boundedEvaluations += opt.getEvaluations();
        }
    }
    // Try a CMAES optimiser which is non-deterministic. To overcome this we perform restarts.
    // CMAESOptimiser based on Matlab code:
    // https://www.lri.fr/~hansen/cmaes.m
    // Take the defaults from the Matlab documentation
    //Double.NEGATIVE_INFINITY;
    double stopFitness = 0;
    boolean isActiveCMA = true;
    int diagonalOnly = 0;
    int checkFeasableCount = 1;
    //Well19937c();
    RandomGenerator random = new Well44497b();
    boolean generateStatistics = false;
    ConvergenceChecker<PointValuePair> checker = new SimpleValueChecker(1e-6, 1e-10);
    // The sigma determines the search range for the variables. It should be 1/3 of the initial search region.
    double[] range = new double[lB.length];
    for (int i = 0; i < lB.length; i++) range[i] = (uB[i] - lB[i]) / 3;
    OptimizationData sigma = new CMAESOptimizer.Sigma(range);
    OptimizationData popSize = new CMAESOptimizer.PopulationSize((int) (4 + Math.floor(3 * Math.log(initialSolution.length))));
    SimpleBounds bounds = new SimpleBounds(lB, uB);
    opt = new CMAESOptimizer(maxEvaluations.getMaxEval(), stopFitness, isActiveCMA, diagonalOnly, checkFeasableCount, random, generateStatistics, checker);
    // Restart the optimiser several times and take the best answer.
    for (int iteration = 0; iteration <= fitRestarts; iteration++) {
        try {
            // Start from the initial solution
            PointValuePair constrainedSolution = opt.optimize(new InitialGuess(initialSolution), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations);
            if (debug)
                System.out.printf("CMAES Iter %d initial = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations());
            boundedEvaluations += opt.getEvaluations();
            if (optimum == null || constrainedSolution.getValue() < optimum.getValue()) {
                optimum = constrainedSolution;
            }
        } catch (TooManyEvaluationsException e) {
        } catch (TooManyIterationsException e) {
        } finally {
            boundedEvaluations += maxEvaluations.getMaxEval();
        }
        if (optimum == null)
            continue;
        try {
            // Also restart from the current optimum
            PointValuePair constrainedSolution = opt.optimize(new InitialGuess(optimum.getPointRef()), objective, GoalType.MINIMIZE, bounds, sigma, popSize, maxEvaluations);
            if (debug)
                System.out.printf("CMAES Iter %d restart = %g (%d)\n", iteration, constrainedSolution.getValue(), opt.getEvaluations());
            if (constrainedSolution.getValue() < optimum.getValue()) {
                optimum = constrainedSolution;
            }
        } catch (TooManyEvaluationsException e) {
        } catch (TooManyIterationsException e) {
        } finally {
            boundedEvaluations += maxEvaluations.getMaxEval();
        }
    }
    return optimum;
}

Example 48 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project GDSC-SMLM by aherbert.

the class GaussianPSFModel method sample.

private double[] sample(final int n, final double mu, final double sigma) {
    final double[] x = new double[n];
    final RandomGenerator random = rand.getRandomGenerator();
    for (int i = 0; i < n; i++) {
        x[i] = sigma * random.nextGaussian() + mu;
    }
    return x;
}

Example 49 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project GDSC-SMLM by aherbert.

the class ImageModel method setRandomGenerator.

/**
	 * Set the random generator for creating the image
	 * 
	 * @param random
	 */
public void setRandomGenerator(RandomGenerator random) {
    if (random == null)
        throw new NullPointerException("Random generator must not be null");
    this.random = random;
    this.randomGenerator = new RandomDataGenerator(random);
}

Example 50 with RandomGenerator

use of org.apache.commons.math3.random.RandomGenerator in project GDSC-SMLM by aherbert.

the class ImagePSFModel method sample.

private double[][] sample(final int n, double x0, double x1, double x2) {
    final int slice = getSlice(x2);
    if (slice < 0 || slice >= sumImage.length)
        return new double[][] { null, null };
    final double[] sumPsf = cumulativeImage[slice];
    final RandomGenerator randomX, randomY;
    // Use the input generator
    randomX = rand.getRandomGenerator();
    randomY = rand.getRandomGenerator();
    //// Debugging - Use a uniform distribution to sample x
    //randomX = new AbstractRandomGenerator()
    //{
    //	int pos = 0;
    //
    //	@Override
    //	public double nextDouble()
    //	{
    //		double p = (double) pos / n;
    //		if (pos++ >= n)
    //			pos = 0;
    //		return p;
    //	}
    //
    //	@Override
    //	public void setSeed(long seed)
    //	{
    //		pos = Math.abs((int) seed) % n;
    //	}
    //};
    //// Debugging - Use a fixed distribution to sample y
    //randomY = new AbstractRandomGenerator()
    //{
    //	public double nextDouble()
    //	{
    //		return 0.5;
    //	}
    //
    //	@Override
    //	public void setSeed(long seed)
    //	{
    //	}
    //};
    // Ensure the generated index is adjusted to the correct position
    // The index will be generated at 0,0 of a pixel in the PSF image.
    // We must subtract the PSF centre so that the middle coords are zero.
    x0 -= xyCentre[slice][0] * unitsPerPixel;
    x1 -= xyCentre[slice][1] * unitsPerPixel;
    //x0 -= 0.5 * psfWidth * unitsPerPixel;
    //x1 -= 0.5 * psfWidth * unitsPerPixel;
    final double max = sumPsf[sumPsf.length - 1];
    double[] x = new double[n];
    double[] y = new double[n];
    int count = 0;
    double sx = 0, sy = 0, s = 0;
    for (int i = 0; i < n; i++) {
        final double p = randomX.nextDouble();
        // If outside the observed PSF then skip 
        if (p > max)
            continue;
        final int index = findIndex(sumPsf, p);
        // Interpolate xi using the fraction of the pixel
        double xi = index % psfWidth;
        xi += (p - sumPsf[index]) / (sumPsf[index + 1] - sumPsf[index]);
        // Add random dither within pixel for y
        final double yi = randomY.nextDouble() + (index / psfWidth);
        if (COM_CHECK) {
            final double v = 1;
            sx += xi * v;
            sy += yi * v;
            s += v;
        }
        x[count] = x0 + (xi * this.unitsPerPixel);
        y[count] = x1 + (yi * this.unitsPerPixel);
        count++;
    }
    if (COM_CHECK) {
        sx = sx / s - xyCentre[slice][0];
        sy = sy / s - xyCentre[slice][1];
        System.out.printf("%dx%d sample centre [ %f %f ] ( %f %f )\n", psfWidth, psfWidth, sx, sy, sx / unitsPerPixel, sy / unitsPerPixel);
    }
    x = Arrays.copyOf(x, count);
    y = Arrays.copyOf(y, count);
    return new double[][] { x, y };
}