Examples with BasePoissonFisherInformation uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation used on opensource projects

Example 1 with BasePoissonFisherInformation

use of uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation in project GDSC-SMLM by aherbert.

the class CameraModelFisherInformationAnalysis method plotFromCache.

private static void plotFromCache() {
    // Build a list of curve stored in the cache
    final String[] names = new String[cache.size()];
    final PoissonFisherInformationData[] datas = new PoissonFisherInformationData[cache.size()];
    int count = 0;
    for (final Entry<FiKey, PoissonFisherInformationData> e : cache.entrySet()) {
        names[count] = e.getKey().toString();
        datas[count] = e.getValue();
        count++;
    }
    final ExtendedGenericDialog gd = new ExtendedGenericDialog(TITLE);
    gd.addChoice("Fisher_information", names, cachePlot);
    gd.addChoice("Plot_point", POINT_OPTION, pointOption);
    gd.showDialog();
    if (gd.wasCanceled()) {
        return;
    }
    final int index = gd.getNextChoiceIndex();
    pointOption = gd.getNextChoiceIndex();
    cachePlot = names[index];
    final PoissonFisherInformationData data = datas[index];
    final FiKey key = new FiKey(data);
    count = 0;
    final double[] exp = new double[data.getAlphaSampleCount()];
    final double[] alpha = new double[data.getAlphaSampleCount()];
    for (final AlphaSample s : data.getAlphaSampleList()) {
        exp[count] = s.getLog10Mean();
        alpha[count] = s.getAlpha();
        count++;
    }
    // Just in case
    // Sort.sortArrays(alpha, exp, true);
    // Test if we can use ImageJ support for a X log scale
    final boolean logScaleX = ((float) Math.pow(10, exp[0]) != 0);
    double[] x = exp;
    String xTitle = "log10(photons)";
    if (logScaleX) {
        final double[] photons = new double[exp.length];
        for (int i = 0; i < photons.length; i++) {
            photons[i] = Math.pow(10, exp[i]);
        }
        x = photons;
        xTitle = "photons";
    }
    // Get interpolation for alpha. Convert to base e.
    final double[] logU = exp.clone();
    final double scale = Math.log(10);
    for (int i = 0; i < logU.length; i++) {
        logU[i] *= scale;
    }
    final BasePoissonFisherInformation if1 = getInterpolatedPoissonFisherInformation(key.getType(), logU, alpha, null);
    // Interpolate with 5 points per sample for smooth curve
    final int n = 5;
    final TDoubleArrayList iexp = new TDoubleArrayList();
    final TDoubleArrayList iphotons = new TDoubleArrayList();
    for (int i = 1; i < exp.length; i++) {
        final int i_1 = i - 1;
        final double h = (exp[i] - exp[i_1]) / n;
        for (int j = 0; j < n; j++) {
            final double e = exp[i_1] + j * h;
            iexp.add(e);
            iphotons.add(Math.pow(10, e));
        }
    }
    iexp.add(exp[exp.length - 1]);
    iphotons.add(Math.pow(10, exp[exp.length - 1]));
    final double[] photons = iphotons.toArray();
    final double[] ix = (logScaleX) ? photons : iexp.toArray();
    final double[] ialpha1 = getAlpha(if1, photons);
    final int pointShape = getPointShape(pointOption);
    final String title = "Cached Relative Fisher Information";
    final Plot plot = new Plot(title, xTitle, "Noise coefficient (alpha)");
    plot.setLimits(x[0], x[x.length - 1], -0.05, 1.05);
    if (logScaleX) {
        plot.setLogScaleX();
    }
    plot.setColor(Color.blue);
    plot.addPoints(ix, ialpha1, Plot.LINE);
    // Option to show nodes
    if (pointShape != -1) {
        plot.addPoints(x, alpha, pointShape);
    }
    plot.setColor(Color.BLACK);
    plot.addLabel(0, 0, cachePlot);
    ImageJUtils.display(title, plot);
}

Example 2 with BasePoissonFisherInformation

use of uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation in project GDSC-SMLM by aherbert.

the class CameraModelFisherInformationAnalysis method analyse.

private void analyse() {
    final CameraType type1 = CameraType.forNumber(settings.getCamera1Type());
    final CameraType type2 = CameraType.forNumber(settings.getCamera2Type());
    final FiKey key1 = new FiKey(type1, settings.getCamera1Gain(), settings.getCamera1Noise());
    final FiKey key2 = new FiKey(type2, settings.getCamera2Gain(), settings.getCamera2Noise());
    final BasePoissonFisherInformation f1 = createPoissonFisherInformation(key1);
    if (f1 == null) {
        return;
    }
    final BasePoissonFisherInformation f2 = createPoissonFisherInformation(key2);
    if (f2 == null) {
        return;
    }
    final double[] exp = createExponents();
    if (exp == null) {
        return;
    }
    final double[] photons = new double[exp.length];
    for (int i = 0; i < photons.length; i++) {
        photons[i] = Math.pow(10, exp[i]);
    }
    // Get the alpha.
    // This may be from the cache or computed shutdown the executor if it was created.
    double[] alpha1;
    double[] alpha2;
    try {
        alpha1 = getAlpha(photons, exp, f1, key1);
        alpha2 = getAlpha(photons, exp, f2, key2);
    } finally {
        if (es != null) {
            es.shutdownNow();
        }
    }
    // Compute the Poisson Fisher information
    final double[] fi1 = getFisherInformation(alpha1, photons);
    final double[] fi2 = getFisherInformation(alpha2, photons);
    // ==============
    if (debug && f2 instanceof PoissonGammaGaussianFisherInformation) {
        final PoissonGammaGaussianFisherInformation pgg = (PoissonGammaGaussianFisherInformation) f2;
        final double t = 200;
        final double fi = pgg.getFisherInformation(t);
        final double alpha = fi * t;
        final double[][] data1 = pgg.getFisherInformationFunction(false);
        final double[][] data2 = pgg.getFisherInformationFunction(true);
        final double[] fif = data1[1];
        int max = 0;
        for (int j = 1; j < fif.length; j++) {
            if (fif[max] < fif[j]) {
                max = j;
            }
        }
        ImageJUtils.log("PGG(p=%g) max=%g", t, data1[0][max]);
        final String title = TITLE + " photons=" + MathUtils.rounded(t) + " alpha=" + MathUtils.rounded(alpha);
        final Plot plot = new Plot(title, "Count", "FI function");
        double yMax = MathUtils.max(data1[1]);
        yMax = MathUtils.maxDefault(yMax, data2[1]);
        plot.setLimits(data2[0][0], data2[0][data2[0].length - 1], 0, yMax);
        plot.setColor(Color.red);
        plot.addPoints(data1[0], data1[1], Plot.LINE);
        plot.setColor(Color.blue);
        plot.addPoints(data2[0], data2[1], Plot.LINE);
        ImageJUtils.display(title, plot);
    }
    // ==============
    final Color color1 = Color.BLUE;
    final Color color2 = Color.RED;
    final WindowOrganiser wo = new WindowOrganiser();
    // Test if we can use ImageJ support for a X log scale
    final boolean logScaleX = ((float) photons[0] != 0);
    final double[] x = (logScaleX) ? photons : exp;
    final String xTitle = (logScaleX) ? "photons" : "log10(photons)";
    // Get interpolation for alpha. Convert to base e.
    final double[] logU = exp.clone();
    final double scale = Math.log(10);
    for (int i = 0; i < logU.length; i++) {
        logU[i] *= scale;
    }
    final BasePoissonFisherInformation if1 = getInterpolatedPoissonFisherInformation(type1, logU, alpha1, f1);
    final BasePoissonFisherInformation if2 = getInterpolatedPoissonFisherInformation(type2, logU, alpha2, f2);
    // Interpolate with 5 points per sample for smooth curve
    final int n = 5 * exp.length;
    final double[] iexp = new double[n + 1];
    final double[] iphotons = new double[iexp.length];
    final double h = (exp[exp.length - 1] - exp[0]) / n;
    for (int i = 0; i <= n; i++) {
        iexp[i] = exp[0] + i * h;
        iphotons[i] = Math.pow(10, iexp[i]);
    }
    final double[] ix = (logScaleX) ? iphotons : iexp;
    final double[] ialpha1 = getAlpha(if1, iphotons);
    final double[] ialpha2 = getAlpha(if2, iphotons);
    final int pointShape = getPointShape(settings.getPointOption());
    final String name1 = getName(key1);
    final String name2 = getName(key2);
    final String legend = name1 + "\n" + name2;
    String title = "Relative Fisher Information";
    Plot plot = new Plot(title, xTitle, "Noise coefficient (alpha)");
    plot.setLimits(x[0], x[x.length - 1], -0.05, 1.05);
    if (logScaleX) {
        plot.setLogScaleX();
    }
    plot.setColor(color1);
    plot.addPoints(ix, ialpha1, Plot.LINE);
    plot.setColor(color2);
    plot.addPoints(ix, ialpha2, Plot.LINE);
    plot.setColor(Color.BLACK);
    plot.addLegend(legend);
    // Option to show nodes
    if (pointShape != -1) {
        plot.setColor(color1);
        plot.addPoints(x, alpha1, pointShape);
        plot.setColor(color2);
        plot.addPoints(x, alpha2, pointShape);
        plot.setColor(Color.BLACK);
    }
    ImageJUtils.display(title, plot, 0, wo);
    // The approximation should not produce an infinite computation
    double[] limits = new double[] { fi2[fi2.length - 1], fi2[fi2.length - 1] };
    limits = limits(limits, fi1);
    limits = limits(limits, fi2);
    // Check if we can use ImageJ support for a Y log scale
    final boolean logScaleY = ((float) limits[1] <= Float.MAX_VALUE);
    if (!logScaleY) {
        for (int i = 0; i < fi1.length; i++) {
            fi1[i] = Math.log10(fi1[i]);
            fi2[i] = Math.log10(fi2[i]);
        }
        limits[0] = Math.log10(limits[0]);
        limits[1] = Math.log10(limits[1]);
    }
    final String yTitle = (logScaleY) ? "Fisher Information" : "log10(Fisher Information)";
    title = "Fisher Information";
    plot = new Plot(title, xTitle, yTitle);
    plot.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
    if (logScaleX) {
        plot.setLogScaleX();
    }
    if (logScaleY) {
        plot.setLogScaleY();
    }
    plot.setColor(color1);
    plot.addPoints(x, fi1, Plot.LINE);
    plot.setColor(color2);
    plot.addPoints(x, fi2, Plot.LINE);
    plot.setColor(Color.BLACK);
    plot.addLegend(legend);
    // // Option to show nodes
    // This gets messy as the lines are straight
    // if (pointShape != -1)
    // {
    // plot.setColor(color1);
    // plot.addPoints(x, pgFI, pointShape);
    // plot.setColor(color3);
    // plot.addPoints(x, pggFI, pointShape);
    // plot.setColor(Color.BLACK);
    // }
    ImageJUtils.display(title, plot, 0, wo);
    wo.tile();
}

Example 3 with BasePoissonFisherInformation

use of uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation in project GDSC-SMLM by aherbert.

the class CameraModelFisherInformationAnalysis method getAlpha.

private double[] getAlpha(final double[] photons, double[] exp, final BasePoissonFisherInformation fi, FiKey key) {
    final CameraType type = key.getType();
    final double[] alpha = new double[photons.length];
    if (!type.isFast()) {
        final int[] index;
        // Try and load from the cache
        final PoissonFisherInformationData data = load(key);
        if (data != null) {
            // Dump the samples
            final TDoubleArrayList meanList = new TDoubleArrayList(data.getAlphaSampleCount());
            final TDoubleArrayList alphalist = new TDoubleArrayList(data.getAlphaSampleCount());
            for (final AlphaSample sample : data.getAlphaSampleList()) {
                meanList.add(sample.getLog10Mean());
                alphalist.add(sample.getAlpha());
            }
            final double[] exp2 = meanList.toArray();
            final double[] alphas = alphalist.toArray();
            SortUtils.sortData(alphas, exp2, true, false);
            // Find any exponent not in the array
            final TIntArrayList list = new TIntArrayList(exp.length);
            for (int i = 0; i < exp.length; i++) {
                // Assume exp2 is sorted
                final int j = Arrays.binarySearch(exp2, exp[i]);
                if (j < 0) {
                    // Add to indices to compute
                    list.add(i);
                } else {
                    // Get alpha
                    alpha[i] = alphas[j];
                }
            }
            index = list.toArray();
        } else {
            // Compute all
            index = SimpleArrayUtils.natural(alpha.length);
        }
        if (index.length > 0) {
            IJ.showStatus("Computing " + getName(key));
            final int nThreads = Prefs.getThreads();
            if (es == null) {
                es = Executors.newFixedThreadPool(nThreads);
            }
            final Ticker ticker = ImageJUtils.createTicker(index.length, nThreads);
            final int nPerThread = (int) Math.ceil((double) index.length / nThreads);
            final LocalList<Future<?>> futures = new LocalList<>(nThreads);
            for (int i = 0; i < index.length; i += nPerThread) {
                final int start = i;
                final int end = Math.min(index.length, i + nPerThread);
                futures.add(es.submit(() -> {
                    final BasePoissonFisherInformation fi2 = fi.copy();
                    for (int ii = start; ii < end; ii++) {
                        final int j = index[ii];
                        alpha[j] = fi2.getAlpha(photons[j]);
                        ticker.tick();
                    }
                }));
            }
            ConcurrencyUtils.waitForCompletionUnchecked(futures);
            ImageJUtils.finished();
            save(key, exp, alpha);
        }
    } else {
        // Simple single threaded method.
        for (int i = 0; i < alpha.length; i++) {
            alpha[i] = fi.getAlpha(photons[i]);
        }
    }
    return alpha;
}

Example 4 with BasePoissonFisherInformation

use of uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation in project GDSC-SMLM by aherbert.

the class CameraModelFisherInformationAnalysis method loadFunction.

/**
 * Load the Poisson Fisher information for the camera type from the cache.
 *
 * @param type the type
 * @param gain the gain
 * @param noise the noise
 * @param relativeError the relative error (used to compare the gain and noise)
 * @return the poisson fisher information (or null)
 */
public static InterpolatedPoissonFisherInformation loadFunction(CameraType type, double gain, double noise, double relativeError) {
    if (type == null || type.isFast()) {
        return null;
    }
    final FiKey key = new FiKey(type, gain, noise);
    PoissonFisherInformationData data = load(key);
    if (data == null && relativeError > 0 && relativeError < 0.1) {
        // Fuzzy matching
        double error = 1;
        for (final PoissonFisherInformationData d : cache.values()) {
            final double e1 = DoubleEquality.relativeError(gain, d.getGain());
            if (e1 < relativeError) {
                final double e2 = DoubleEquality.relativeError(noise, d.getNoise());
                if (e2 < relativeError) {
                    // Combined error - Euclidean distance
                    final double e = e1 * e1 + e2 * e2;
                    if (error > e) {
                        error = e;
                        data = d;
                    }
                }
            }
        }
    }
    if (data == null) {
        return null;
    }
    // Dump the samples. Convert to base e.
    final double scale = Math.log(10);
    final TDoubleArrayList meanList = new TDoubleArrayList(data.getAlphaSampleCount());
    final TDoubleArrayList alphalist = new TDoubleArrayList(data.getAlphaSampleCount());
    for (final AlphaSample sample : data.getAlphaSampleList()) {
        meanList.add(sample.getLog10Mean() * scale);
        alphalist.add(sample.getAlpha());
    }
    final double[] means = meanList.toArray();
    final double[] alphas = alphalist.toArray();
    SortUtils.sortData(alphas, means, true, false);
    final BasePoissonFisherInformation upperf = (type == CameraType.EM_CCD) ? new HalfPoissonFisherInformation() : createPoissonGaussianApproximationFisherInformation(noise / gain);
    return new InterpolatedPoissonFisherInformation(means, alphas, type.isLowerFixedI(), upperf);
}

Example 5 with BasePoissonFisherInformation

use of uk.ac.sussex.gdsc.smlm.function.BasePoissonFisherInformation in project GDSC-SMLM by aherbert.

the class CreateData method createLikelihoodFunction.

/**
 * Creates the likelihood function. This is used for CRLB computation.
 */
private void createLikelihoodFunction() {
    final CameraType cameraType = settings.getCameraType();
    final boolean isCcd = CalibrationProtosHelper.isCcdCameraType(cameraType);
    fiFunction = new BasePoissonFisherInformation[settings.getSize() * settings.getSize()];
    if (isCcd) {
        BasePoissonFisherInformation fi;
        final CreateDataSettingsHelper helper = new CreateDataSettingsHelper(settings);
        final double readNoise = helper.getReadNoiseInCounts();
        if (cameraType == CameraType.EMCCD) {
            // We only want the amplification (without QE applied)
            final double amp = helper.getAmplification();
            // This should be interpolated from a stored curve
            final InterpolatedPoissonFisherInformation i = CameraModelFisherInformationAnalysis.loadFunction(CameraModelFisherInformationAnalysis.CameraType.EM_CCD, amp, readNoise);
            if (i == null) {
                throw new IllegalStateException("No stored Fisher information for EM-CCD camera with gain " + amp + " and noise " + readNoise + "\n \nPlease generate using the " + CameraModelFisherInformationAnalysis.TITLE);
            }
            fi = i;
        } else {
            // This is fast enough to compute dynamically.
            // Read noise is in electrons so use directly.
            fi = new PoissonGaussianFisherInformation(settings.getReadNoise());
        }
        Arrays.fill(fiFunction, fi);
    } else if (cameraType == CameraType.SCMOS) {
        // Build per-pixel likelihood function.
        // Get the normalised variance per pixel.
        final float[] v = cameraModel.getNormalisedVariance(cameraModel.getBounds());
        // Build the function
        for (int i = 0; i < fiFunction.length; i++) {
            fiFunction[i] = new PoissonGaussianFisherInformation(Math.sqrt(v[i]));
        }
    } else {
        throw new IllegalArgumentException("Unsupported camera type: " + CalibrationProtosHelper.getName(cameraType));
    }
}