Examples with AlphaSample uk.ac.sussex.gdsc.smlm.data.config.FisherProtos.AlphaSample used on opensource projects

Example 1 with AlphaSample

use of uk.ac.sussex.gdsc.smlm.data.config.FisherProtos.AlphaSample 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 AlphaSample

use of uk.ac.sussex.gdsc.smlm.data.config.FisherProtos.AlphaSample 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 3 with AlphaSample

use of uk.ac.sussex.gdsc.smlm.data.config.FisherProtos.AlphaSample in project GDSC-SMLM by aherbert.

the class CameraModelFisherInformationAnalysis method save.

/**
 * Save the data to the cache.
 *
 * @param key the key
 * @param log10photons the log 10 photons
 * @param alpha the alpha
 */
private static void save(FiKey key, double[] log10photons, double[] alpha) {
    final CameraType type = key.getType();
    if (type.isFast()) {
        return;
    }
    PoissonFisherInformationData data = cache.get(key);
    if (data != null) {
        // This should only be called if new values have been computed.
        // so assume we must merge the lists.
        // Note: The lists must be sorted.
        final AlphaSample[] list1 = data.getAlphaSampleList().toArray(new AlphaSample[0]);
        final LocalList<AlphaSample> list = new LocalList<>(list1.length + log10photons.length);
        int i1 = 0;
        int i2 = 0;
        final AlphaSample.Builder a2 = AlphaSample.newBuilder();
        while (i1 < list1.length && i2 < log10photons.length) {
            final AlphaSample a1 = list1[i1];
            final double mean = log10photons[i2];
            if (a1.getLog10Mean() == mean) {
                list.add(a1);
                i1++;
                i2++;
            } else if (a1.getLog10Mean() < mean) {
                list.add(a1);
                i1++;
            } else {
                // (a1.getLog10Mean() > mean)
                a2.setLog10Mean(mean);
                a2.setAlpha(alpha[i2++]);
                list.add(a2.build());
            }
        }
        while (i1 < list1.length) {
            list.add(list1[i1++]);
        }
        while (i2 < log10photons.length) {
            a2.setLog10Mean(log10photons[i2]);
            a2.setAlpha(alpha[i2++]);
            list.add(a2.build());
        }
        final PoissonFisherInformationData.Builder b = data.toBuilder();
        b.clearAlphaSample();
        b.addAllAlphaSample(list);
        data = b.build();
    } else {
        final PoissonFisherInformationData.Builder b = PoissonFisherInformationData.newBuilder();
        final AlphaSample.Builder sample = AlphaSample.newBuilder();
        b.setType(key.type);
        b.setGain(key.gain);
        b.setNoise(key.noise);
        for (int i = 0; i < log10photons.length; i++) {
            sample.setLog10Mean(log10photons[i]);
            sample.setAlpha(alpha[i]);
            b.addAlphaSample(sample);
        }
        data = b.build();
    }
    cache.put(key, data);
    // Also save EM-CCD data to file
    if (type == CameraType.EM_CCD) {
        final int t = type.ordinal();
        // Get the EM-CCD keys
        final LocalList<FiKey> list = new LocalList<>(cache.size());
        for (final FiKey k : cache.keySet()) {
            if (k.type == t) {
                list.add(k);
            }
        }
        if (list.size() > MAX_DATA_TO_FILE) {
            // Sort by age - Youngest first
            list.sort((o1, o2) -> Long.compare(o2.timeStamp, o1.timeStamp));
        }
        final PoissonFisherInformationCache.Builder cacheData = PoissonFisherInformationCache.newBuilder();
        for (int i = Math.min(list.size(), MAX_DATA_TO_FILE); i-- > 0; ) {
            cacheData.addData(cache.get(list.unsafeGet(i)));
        }
        SettingsManager.writeSettings(cacheData);
    }
}

Example 4 with AlphaSample

use of uk.ac.sussex.gdsc.smlm.data.config.FisherProtos.AlphaSample 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);
}