Examples with Max org.apache.commons.math3.stat.descriptive.rank.Max used on opensource projects

Example 41 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project cassandra by apache.

the class LongSharedExecutorPoolTest method testPromptnessOfExecution.

private void testPromptnessOfExecution(long intervalNanos, float loadIncrement) throws InterruptedException, ExecutionException {
    final int executorCount = 4;
    int threadCount = 8;
    int maxQueued = 1024;
    final WeibullDistribution workTime = new WeibullDistribution(3, 200000);
    final long minWorkTime = TimeUnit.MICROSECONDS.toNanos(1);
    final long maxWorkTime = TimeUnit.MILLISECONDS.toNanos(1);
    final int[] threadCounts = new int[executorCount];
    final WeibullDistribution[] workCount = new WeibullDistribution[executorCount];
    final ExecutorService[] executors = new ExecutorService[executorCount];
    for (int i = 0; i < executors.length; i++) {
        executors[i] = SharedExecutorPool.SHARED.newExecutor(threadCount, maxQueued, "test" + i, "test" + i);
        threadCounts[i] = threadCount;
        workCount[i] = new WeibullDistribution(2, maxQueued);
        threadCount *= 2;
        maxQueued *= 2;
    }
    long runs = 0;
    long events = 0;
    final TreeSet<Batch> pending = new TreeSet<>();
    final BitSet executorsWithWork = new BitSet(executorCount);
    long until = 0;
    // (beyond max queued size), and longer operations
    for (float multiplier = 0f; multiplier < 2.01f; ) {
        if (System.nanoTime() > until) {
            System.out.println(String.format("Completed %.0fK batches with %.1fM events", runs * 0.001f, events * 0.000001f));
            events = 0;
            until = System.nanoTime() + intervalNanos;
            multiplier += loadIncrement;
            System.out.println(String.format("Running for %ds with load multiplier %.1f", TimeUnit.NANOSECONDS.toSeconds(intervalNanos), multiplier));
        }
        // wait a random amount of time so we submit new tasks in various stages of
        long timeout;
        if (pending.isEmpty())
            timeout = 0;
        else if (Math.random() > 0.98)
            timeout = Long.MAX_VALUE;
        else if (pending.size() == executorCount)
            timeout = pending.first().timeout;
        else
            timeout = (long) (Math.random() * pending.last().timeout);
        while (!pending.isEmpty() && timeout > System.nanoTime()) {
            Batch first = pending.first();
            boolean complete = false;
            try {
                for (Result result : first.results.descendingSet()) result.future.get(timeout - System.nanoTime(), TimeUnit.NANOSECONDS);
                complete = true;
            } catch (TimeoutException e) {
            }
            if (!complete && System.nanoTime() > first.timeout) {
                for (Result result : first.results) if (!result.future.isDone())
                    throw new AssertionError();
                complete = true;
            }
            if (complete) {
                pending.pollFirst();
                executorsWithWork.clear(first.executorIndex);
            }
        }
        // if we've emptied the executors, give all our threads an opportunity to spin down
        if (timeout == Long.MAX_VALUE)
            Uninterruptibles.sleepUninterruptibly(10, TimeUnit.MILLISECONDS);
        // submit a random batch to the first free executor service
        int executorIndex = executorsWithWork.nextClearBit(0);
        if (executorIndex >= executorCount)
            continue;
        executorsWithWork.set(executorIndex);
        ExecutorService executor = executors[executorIndex];
        TreeSet<Result> results = new TreeSet<>();
        int count = (int) (workCount[executorIndex].sample() * multiplier);
        long targetTotalElapsed = 0;
        long start = System.nanoTime();
        long baseTime;
        if (Math.random() > 0.5)
            baseTime = 2 * (long) (workTime.sample() * multiplier);
        else
            baseTime = 0;
        for (int j = 0; j < count; j++) {
            long time;
            if (baseTime == 0)
                time = (long) (workTime.sample() * multiplier);
            else
                time = (long) (baseTime * Math.random());
            if (time < minWorkTime)
                time = minWorkTime;
            if (time > maxWorkTime)
                time = maxWorkTime;
            targetTotalElapsed += time;
            Future<?> future = executor.submit(new WaitTask(time));
            results.add(new Result(future, System.nanoTime() + time));
        }
        long end = start + (long) Math.ceil(targetTotalElapsed / (double) threadCounts[executorIndex]) + TimeUnit.MILLISECONDS.toNanos(100L);
        long now = System.nanoTime();
        if (runs++ > executorCount && now > end)
            throw new AssertionError();
        events += results.size();
        pending.add(new Batch(results, end, executorIndex));
    //            System.out.println(String.format("Submitted batch to executor %d with %d items and %d permitted millis", executorIndex, count, TimeUnit.NANOSECONDS.toMillis(end - start)));
    }
}

Example 42 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project GDSC-SMLM by aherbert.

the class CreateData method createUniformDistribution.

/**
	 * Create distribution within an XY border
	 * 
	 * @param border
	 * @return
	 */
private UniformDistribution createUniformDistribution(double border) {
    double depth = (settings.fixedDepth) ? settings.depth / settings.pixelPitch : settings.depth / (2 * settings.pixelPitch);
    // Ensure the focal plane is in the middle of the zDepth
    double[] max = new double[] { settings.size / 2 - border, settings.size / 2 - border, depth };
    double[] min = new double[3];
    for (int i = 0; i < 3; i++) min[i] = -max[i];
    if (settings.fixedDepth)
        min[2] = max[2];
    // Try using different distributions:
    final RandomGenerator rand1 = createRandomGenerator();
    if (settings.distribution.equals(DISTRIBUTION[UNIFORM_HALTON])) {
        return new UniformDistribution(min, max, rand1.nextInt());
    }
    if (settings.distribution.equals(DISTRIBUTION[UNIFORM_SOBOL])) {
        SobolSequenceGenerator rvg = new SobolSequenceGenerator(3);
        rvg.skipTo(rand1.nextInt());
        return new UniformDistribution(min, max, rvg);
    }
    // Create a distribution using random generators for each dimension 
    UniformDistribution distribution = new UniformDistribution(min, max, this);
    return distribution;
}

Example 43 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project GDSC-SMLM by aherbert.

the class PCPALMMolecules method calculateAveragePrecision.

/**
	 * Calculate the average precision by fitting a skewed Gaussian to the histogram of the precision distribution.
	 * <p>
	 * A simple mean and SD of the histogram is computed. If the mean of the Skewed Gaussian does not fit within 3 SDs
	 * of the simple mean then the simple mean is returned.
	 * 
	 * @param molecules
	 * @param title
	 *            the plot title (null if no plot should be displayed)
	 * @param histogramBins
	 * @param logFitParameters
	 *            Record the fit parameters to the ImageJ log
	 * @param removeOutliers
	 *            The distribution is created using all values within 1.5x the inter-quartile range (IQR) of the data
	 * @return The average precision
	 */
public double calculateAveragePrecision(ArrayList<Molecule> molecules, String title, int histogramBins, boolean logFitParameters, boolean removeOutliers) {
    // Plot histogram of the precision
    float[] data = new float[molecules.size()];
    DescriptiveStatistics stats = new DescriptiveStatistics();
    double yMin = Double.NEGATIVE_INFINITY, yMax = 0;
    for (int i = 0; i < data.length; i++) {
        data[i] = (float) molecules.get(i).precision;
        stats.addValue(data[i]);
    }
    // Set the min and max y-values using 1.5 x IQR 
    if (removeOutliers) {
        double lower = stats.getPercentile(25);
        double upper = stats.getPercentile(75);
        if (Double.isNaN(lower) || Double.isNaN(upper)) {
            if (logFitParameters)
                Utils.log("Error computing IQR: %f - %f", lower, upper);
        } else {
            double iqr = upper - lower;
            yMin = FastMath.max(lower - iqr, stats.getMin());
            yMax = FastMath.min(upper + iqr, stats.getMax());
            if (logFitParameters)
                Utils.log("  Data range: %f - %f. Plotting 1.5x IQR: %f - %f", stats.getMin(), stats.getMax(), yMin, yMax);
        }
    }
    if (yMin == Double.NEGATIVE_INFINITY) {
        yMin = stats.getMin();
        yMax = stats.getMax();
        if (logFitParameters)
            Utils.log("  Data range: %f - %f", yMin, yMax);
    }
    float[][] hist = Utils.calcHistogram(data, yMin, yMax, histogramBins);
    Plot2 plot = null;
    if (title != null) {
        plot = new Plot2(title, "Precision", "Frequency");
        float[] xValues = hist[0];
        float[] yValues = hist[1];
        if (xValues.length > 0) {
            double xPadding = 0.05 * (xValues[xValues.length - 1] - xValues[0]);
            plot.setLimits(xValues[0] - xPadding, xValues[xValues.length - 1] + xPadding, 0, Maths.max(yValues) * 1.05);
        }
        plot.addPoints(xValues, yValues, Plot2.BAR);
        Utils.display(title, plot);
    }
    // Extract non-zero data
    float[] x = Arrays.copyOf(hist[0], hist[0].length);
    float[] y = hist[1];
    int count = 0;
    float dx = (x[1] - x[0]) * 0.5f;
    for (int i = 0; i < y.length; i++) if (y[i] > 0) {
        x[count] = x[i] + dx;
        y[count] = y[i];
        count++;
    }
    x = Arrays.copyOf(x, count);
    y = Arrays.copyOf(y, count);
    // Sense check to fitted data. Get mean and SD of histogram
    double[] stats2 = Utils.getHistogramStatistics(x, y);
    double mean = stats2[0];
    if (logFitParameters)
        log("  Initial Statistics: %f +/- %f", stats2[0], stats2[1]);
    // Standard Gaussian fit
    double[] parameters = fitGaussian(x, y);
    if (parameters == null) {
        log("  Failed to fit initial Gaussian");
        return mean;
    }
    double newMean = parameters[1];
    double error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    mean = newMean;
    if (logFitParameters)
        log("  Initial Gaussian: %f @ %f +/- %f", parameters[0], parameters[1], parameters[2]);
    double[] initialSolution = new double[] { parameters[0], parameters[1], parameters[2], -1 };
    // Fit to a skewed Gaussian (or appropriate function)
    double[] skewParameters = fitSkewGaussian(x, y, initialSolution);
    if (skewParameters == null) {
        log("  Failed to fit Skewed Gaussian");
        return mean;
    }
    SkewNormalFunction sn = new SkewNormalFunction(skewParameters);
    if (logFitParameters)
        log("  Skewed Gaussian: %f @ %f +/- %f (a = %f) => %f +/- %f", skewParameters[0], skewParameters[1], skewParameters[2], skewParameters[3], sn.getMean(), Math.sqrt(sn.getVariance()));
    newMean = sn.getMean();
    error = Math.abs(stats2[0] - newMean) / stats2[1];
    if (error > 3) {
        log("  Failed to fit Skewed Gaussian: %f standard deviations from histogram mean", error);
        return mean;
    }
    if (newMean < yMin || newMean > yMax) {
        log("  Failed to fit Skewed Gaussian: %f outside data range %f - %f", newMean, yMin, yMax);
        return mean;
    }
    // Use original histogram x-axis to maintain all the bins
    if (plot != null) {
        x = hist[0];
        for (int i = 0; i < y.length; i++) x[i] += dx;
        plot.setColor(Color.red);
        addToPlot(plot, x, skewParameters, Plot2.LINE);
        plot.setColor(Color.black);
        Utils.display(title, plot);
    }
    // Return the average precision from the fitted curve
    return newMean;
}

Example 44 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max 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 };
}

Example 45 with Max

use of org.apache.commons.math3.stat.descriptive.rank.Max in project GDSC-SMLM by aherbert.

the class PoissonFunctionTest method cumulativeProbabilityIsOneWithRealAbove4.

private void cumulativeProbabilityIsOneWithRealAbove4(final double gain, final double mu, int min, int max) {
    final double o = mu * gain;
    final PoissonFunction f = new PoissonFunction(1.0 / gain, true);
    double p = 0;
    SimpsonIntegrator in = new SimpsonIntegrator(3, 30);
    try {
        p = in.integrate(Integer.MAX_VALUE, new UnivariateFunction() {

            public double value(double x) {
                return f.likelihood(x, o);
            }
        }, min, max);
        System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertEquals(String.format("g=%f, mu=%f", gain, mu), 1, p, 0.02);
    } catch (TooManyEvaluationsException e) {
        //double inc = max / 20000.0;
        //for (double x = 0; x <= max; x += inc)
        //{
        //	final double pp = f.likelihood(x, o);
        //	//System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, pp);
        //	p += pp;
        //}
        //System.out.printf("g=%f, mu=%f, o=%f, p=%f\n", gain, mu, o, p);
        Assert.assertFalse(e.getMessage(), true);
    }
}