Examples with LocalList uk.ac.sussex.gdsc.core.utils.LocalList used on opensource projects

Example 16 with LocalList

use of uk.ac.sussex.gdsc.core.utils.LocalList in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method getCoordinates.

private static TIntObjectHashMap<UniqueIdPeakResult[]> getCoordinates(MemoryPeakResults results) {
    final TIntObjectHashMap<UniqueIdPeakResult[]> coords = new TIntObjectHashMap<>();
    if (results.size() > 0) {
        // Do not use HashMap directly to build the coords object since there
        // will be many calls to getEntry(). Instead sort the results and use
        // a new list for each time point
        results.sort();
        final Counter uniqueId = new Counter();
        final FrameCounter counter = new FrameCounter();
        final LocalList<PeakResult> tmp = new LocalList<>();
        // Add the results to the lists
        results.forEach((PeakResultProcedure) result -> {
            if (counter.advanceAndReset(result.getFrame()) && !tmp.isEmpty()) {
                coords.put(counter.previousFrame(), tmp.toArray(new UniqueIdPeakResult[0]));
                tmp.clear();
            }
            tmp.add(new UniqueIdPeakResult(tmp.size(), uniqueId.getAndIncrement(), result));
        });
        if (!tmp.isEmpty()) {
            coords.put(counter.currentFrame(), tmp.toArray(new UniqueIdPeakResult[0]));
        }
    }
    return coords;
}

Example 17 with LocalList

use of uk.ac.sussex.gdsc.core.utils.LocalList in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method createMixed.

/**
 * Creates the multivariate gaussian mixture as the best of many repeats of the expectation
 * maximisation algorithm.
 *
 * @param data the data
 * @param dimensions the dimensions
 * @param numComponents the number of components
 * @return the multivariate gaussian mixture expectation maximization
 */
private MultivariateGaussianMixtureExpectationMaximization createMixed(final double[][] data, int dimensions, int numComponents) {
    // Fit a mixed multivariate Gaussian with different repeats.
    final UnitSphereSampler sampler = UnitSphereSampler.of(UniformRandomProviders.create(Mixers.stafford13(settings.seed++)), dimensions);
    final LocalList<CompletableFuture<MultivariateGaussianMixtureExpectationMaximization>> results = new LocalList<>(settings.repeats);
    final DoubleDoubleBiPredicate test = createConvergenceTest(settings.relativeError);
    if (settings.debug) {
        ImageJUtils.log("  Fitting %d components", numComponents);
    }
    final Ticker ticker = ImageJUtils.createTicker(settings.repeats, 2, "Fitting...");
    final AtomicInteger failures = new AtomicInteger();
    for (int i = 0; i < settings.repeats; i++) {
        final double[] vector = sampler.sample();
        results.add(CompletableFuture.supplyAsync(() -> {
            final MultivariateGaussianMixtureExpectationMaximization fitter = new MultivariateGaussianMixtureExpectationMaximization(data);
            try {
                // This may also throw the same exceptions due to inversion of the covariance matrix
                final MixtureMultivariateGaussianDistribution initialMixture = MultivariateGaussianMixtureExpectationMaximization.estimate(data, numComponents, point -> {
                    double dot = 0;
                    for (int j = 0; j < dimensions; j++) {
                        dot += vector[j] * point[j];
                    }
                    return dot;
                });
                final boolean result = fitter.fit(initialMixture, settings.maxIterations, test);
                // Log the result. Note: The ImageJ log is synchronized.
                if (settings.debug) {
                    ImageJUtils.log("  Fit: log-likelihood=%s, iter=%d, converged=%b", fitter.getLogLikelihood(), fitter.getIterations(), result);
                }
                return result ? fitter : null;
            } catch (NonPositiveDefiniteMatrixException | SingularMatrixException ex) {
                failures.getAndIncrement();
                if (settings.debug) {
                    ImageJUtils.log("  Fit failed during iteration %d. No variance in a sub-population " + "component (check alpha is not always 1.0).", fitter.getIterations());
                }
            } finally {
                ticker.tick();
            }
            return null;
        }));
    }
    ImageJUtils.finished();
    if (failures.get() != 0 && settings.debug) {
        ImageJUtils.log("  %d component fit failed %d/%d", numComponents, failures.get(), settings.repeats);
    }
    // Collect results and return the best model.
    return results.stream().map(f -> f.join()).filter(f -> f != null).sorted((f1, f2) -> Double.compare(f2.getLogLikelihood(), f1.getLogLikelihood())).findFirst().orElse(null);
}

Example 18 with LocalList

use of uk.ac.sussex.gdsc.core.utils.LocalList in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method extractTrackData.

/**
 * Extract the track data. This extracts different descriptors of the track using a rolling local
 * window.
 *
 * <p>Distances are converted to {@code unit} using the provided converter and time units are
 * converted from frame to seconds (s). The diffusion coefficients is in unit^2/s.
 *
 * <p>If categories are added they are remapped to be a natural sequence starting from 0.
 *
 * @param tracks the tracks
 * @param distanceConverter the distance converter
 * @param deltaT the time step of each frame in seconds (delta T)
 * @param hasCategory if true add the category from the result to the end of the data
 * @return the track data (track lengths and descriptors)
 */
private Pair<int[], double[][]> extractTrackData(List<Trace> tracks, TypeConverter<DistanceUnit> distanceConverter, double deltaT, boolean hasCategory) {
    final List<double[]> data = new LocalList<>(tracks.size());
    double[] x = new double[0];
    double[] y = new double[0];
    final int wm1 = settings.window - 1;
    final int valuesLength = hasCategory ? 5 : 4;
    final int mid = settings.window / 2;
    final MsdFitter msdFitter = new MsdFitter(settings, deltaT);
    final double significance = settings.significance;
    final int[] fitResult = new int[4];
    // Factor for the diffusion coefficient: 1/N * 1/(2*dimensions*deltaT) = 1 / 4Nt
    // with N the number of points to average.
    final double diffusionCoefficientFactor = 1.0 / (4 * wm1 * deltaT);
    // Used for the standard deviations
    final Statistics statsX = new Statistics();
    final Statistics statsY = new Statistics();
    final Ticker ticker = ImageJUtils.createTicker(tracks.size(), 1, "Computing track features...");
    // Collect the track categories. Later these are renumbered to a natural sequence from 0.
    final TIntHashSet catSet = new TIntHashSet();
    // final StoredDataStatistics statsAlpha = new StoredDataStatistics(tracks.size());
    // Process each track
    final TIntArrayList lengths = new TIntArrayList(tracks.size());
    for (final Trace track : tracks) {
        // Get xy coordinates
        final int size = track.size();
        if (x.length < size) {
            x = new double[size];
            y = new double[size];
        }
        for (int i = 0; i < size; i++) {
            final PeakResult peak = track.get(i);
            x[i] = distanceConverter.convert(peak.getXPosition());
            y[i] = distanceConverter.convert(peak.getYPosition());
        }
        final int smwm1 = size - wm1;
        final int previousSize = data.size();
        for (int k = 0; k < smwm1; k++) {
            final double[] values = new double[valuesLength];
            data.add(values);
            // middle of even sized windows is between two localisations.
            if (hasCategory) {
                final int cat = track.get(k + mid).getCategory();
                values[4] = cat;
                catSet.add(cat);
            }
            // First point in window = k
            // Last point in window = k + w - 1 = k + wm1
            final int end = k + wm1;
            // 1. Anomalous exponent.
            msdFitter.setData(x, y);
            try {
                msdFitter.fit(k, null);
                // statsAlpha.add(msdFitter.alpha);
                int fitIndex = msdFitter.positiveSlope ? 0 : 2;
                // If better then this is anomalous diffusion
                final double alpha = msdFitter.lvmSolution2.getPoint().getEntry(2);
                values[0] = alpha;
                if (msdFitter.pValue > significance) {
                    fitIndex++;
                }
                fitResult[fitIndex]++;
            // values[0] = msdFitter.alpha;
            // // Debug
            // if (
            // // msdFitter.pValue < 0.2
            // msdFitter.pValue < 0.2 && values[0] < 0
            // // !msdFitter.positiveSlope
            // ) {
            // final RealVector p = msdFitter.lvmSolution2.getPoint();
            // final String title = "anomalous exponent";
            // final Plot plot = new Plot(title, "time (s)", "MSD (um^2)");
            // final double[] t = SimpleArrayUtils.newArray(msdFitter.s.length, deltaT, deltaT);
            // plot.addLabel(0, 0, msdFitter.lvmSolution2.getPoint().toString() + " p="
            // + msdFitter.pValue + ". " + msdFitter.lvmSolution1.getPoint().toString());
            // plot.addPoints(t, msdFitter.s, Plot.CROSS);
            // plot.addPoints(t, msdFitter.model2.value(p).getFirst().toArray(), Plot.LINE);
            // plot.setColor(Color.BLUE);
            // plot.addPoints(t,
            // msdFitter.model1.value(msdFitter.lvmSolution1.getPoint()).getFirst().toArray(),
            // Plot.LINE);
            // plot.setColor(Color.RED);
            // final double[] yy = Arrays.stream(t).map(msdFitter.reg::predict).toArray();
            // plot.addPoints(t, yy, Plot.CIRCLE);
            // System.out.printf("%s : %s", msdFitter.lvmSolution2.getPoint(), values[0]);
            // ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT);
            // System.out.println();
            // }
            } catch (TooManyIterationsException | ConvergenceException ex) {
                if (settings.debug) {
                    ImageJUtils.log("Failed to fit anomalous exponent: " + ex.getMessage());
                }
                // Ignore this and leave as Brownian motion
                values[0] = 1.0;
            }
            // Referenced papers:
            // Hozé, N. H., D. (2017) Statistical methods for large ensembles of super-resolution
            // stochastic single particle trajectories in cell biology.
            // Annual Review of Statistics and Its Application 4, 189-223
            // 
            // Amitai, A., Seeber, A., Gasser, S. M. & Holcman, D. (2017) Visualization of Chromatin
            // Decompaction and Break Site Extrusion as Predicted by Statistical Polymer
            // Modeling of Single-Locus Trajectories. Cell reports 18, 1200-1214
            // 2. Effective diffusion coefficient (Hozé, eq 10).
            // This is the average squared jump distance between successive points
            // divided by 1 / (2 * dimensions * deltaT), i.e. 1 / 4t.
            double sum = 0;
            for (int i = k; i < end; i++) {
                sum += MathUtils.distance2(x[i], y[i], x[i + 1], y[i + 1]);
            }
            values[1] = sum * diffusionCoefficientFactor;
            // 3. Length of confinement (Amitai et al, eq 1).
            // Compute the average of the standard deviation of the position in each dimension.
            statsX.reset();
            statsY.reset();
            for (int i = k; i <= end; i++) {
                statsX.add(x[i]);
                statsY.add(y[i]);
            }
            values[2] = (statsX.getStandardDeviation() + statsY.getStandardDeviation()) / 2;
            // 4. Magnitude of drift vector (Hozé, eq 9).
            // Note: The drift field is given as the expected distance between successive points, i.e.
            // the average step. Since all track windows are the same length this is the same
            // as the distance between the first and last point divided by the number of points.
            // The drift field in each dimension is combined to create a drift norm, i.e. Euclidean
            // distance.
            values[3] = MathUtils.distance(x[k], y[k], x[end], y[end]) / wm1;
        }
        lengths.add(data.size() - previousSize);
        ticker.tick();
    }
    ImageJUtils.finished();
    if (settings.debug) {
        ImageJUtils.log("  +Slope, significant:   %d", fitResult[0]);
        ImageJUtils.log("  +Slope, insignificant: %d", fitResult[1]);
        ImageJUtils.log("  -Slope, significant:   %d", fitResult[2]);
        ImageJUtils.log("  -Slope, insignificant: %d", fitResult[3]);
    }
    ImageJUtils.log("Insignificant anomalous exponents: %d / %d", fitResult[1] + fitResult[3], data.size());
    // System.out.println(statsAlpha.getStatistics().toString());
    final double[][] trackData = data.toArray(new double[0][0]);
    if (hasCategory) {
        final int[] categories = catSet.toArray();
        Arrays.sort(categories);
        // Only remap if non-compact (i.e. not 0 to n)
        final int max = categories[categories.length - 1];
        if (categories[0] != 0 || max + 1 != categories.length) {
            final int[] remap = new int[max + 1];
            for (int i = 0; i < categories.length; i++) {
                remap[categories[i]] = i;
            }
            final int end = valuesLength - 1;
            for (final double[] values : trackData) {
                values[end] = remap[(int) values[end]];
            }
        }
    }
    return Pair.create(lengths.toArray(), trackData);
}

Example 19 with LocalList

use of uk.ac.sussex.gdsc.core.utils.LocalList in project GDSC-SMLM by aherbert.

the class TrackPopulationAnalysis method getTracks.

/**
 * Gets the tracks. Each track has contiguous frames and the length is enough to fit
 * {@code minTrackLength} overlapping windows of the specified size:
 *
 * <pre>
 * length >= window + minTrackLength - 1
 * </pre>
 *
 * @param combinedResults the combined results
 * @param window the window size
 * @param minTrackLength the minimum track length (assumed to be {@code >= 1})
 * @return the tracks
 */
private static List<Trace> getTracks(List<MemoryPeakResults> combinedResults, int window, int minTrackLength) {
    final LocalList<Trace> tracks = new LocalList<>();
    final Statistics stats = new Statistics();
    final int minSize = window + Math.max(minTrackLength, 1) - 1;
    combinedResults.forEach(results -> {
        final int start = tracks.size();
        // Sort by id then frame
        results = results.copy();
        results.sort(IdFramePeakResultComparator.INSTANCE);
        final int size = results.size();
        // Skip IDs not associated with clustering
        int index = 0;
        while (index < size && results.get(index).getId() < 1) {
            index++;
        }
        // Initialise current id and frame
        int id = results.get(index).getId() - 1;
        int frame = results.get(index).getFrame();
        Trace track = new Trace();
        for (; index < size; index++) {
            final PeakResult result = results.get(index);
            // Same ID and contiguous frames
            if (result.getId() != id || result.getFrame() != frame + 1) {
                addTrack(minSize, tracks, track);
                track = new Trace();
            }
            id = result.getId();
            frame = result.getFrame();
            track.add(result);
        }
        addTrack(minSize, tracks, track);
        stats.reset();
        for (int i = start; i < tracks.size(); i++) {
            stats.add(tracks.unsafeGet(i).size());
        }
        final StringBuilder sb = new StringBuilder(256);
        TextUtils.formatTo(sb, "%s tracks=%d, length=%s +/- %s", results.getName(), stats.getN(), MathUtils.rounded(stats.getMean(), 3), MathUtils.rounded(stats.getStandardDeviation(), 3));
        // Limit of diffusion coefficient from the localisation precision.
        // Just use the entire dataset for simplicity (i.e. not the tracks of min length).
        final PrecisionResultProcedure pp = new PrecisionResultProcedure(results);
        try {
            pp.getPrecision();
            final Mean mean = new Mean();
            for (final double p : pp.precisions) {
                mean.add(p);
            }
            // 2nDt = MSD (n = number of dimensions)
            // D = MSD / 2nt
            final CalibrationReader reader = results.getCalibrationReader();
            final double t = reader.getExposureTime() / 1000.0;
            // Assume computed in nm. Convert to um.
            final double x = mean.getMean() / 1000;
            final double d = x * x / (2 * t);
            TextUtils.formatTo(sb, ", precision=%s nm, D limit=%s um^2/s", MathUtils.rounded(x * 1000, 4), MathUtils.rounded(d, 4));
        } catch (final DataException ex) {
        // No precision
        }
        IJ.log(sb.toString());
    });
    return tracks;
}

Example 20 with LocalList

use of uk.ac.sussex.gdsc.core.utils.LocalList in project GDSC-SMLM by aherbert.

the class BenchmarkFit method getStartPoints.

/**
 * Gets the start points.
 *
 * @param is3D Set to true if 3D
 * @return The starting points for the fitting
 */
private double[][] getStartPoints(boolean is3D) {
    if (startPoints != null) {
        return startPoints;
    }
    final LocalList<double[]> list = new LocalList<>();
    // Set up origin with an offset
    final double[] origin = new double[3];
    switch(originXY) {
        case 2:
            // Offset from the origin in pixels
            origin[0] = offsetX / benchmarkParameters.pixelPitch;
            origin[1] = offsetY / benchmarkParameters.pixelPitch;
            break;
        case 0:
        case 1:
            // No offset
            break;
        default:
            throw new IllegalStateException();
    }
    switch(originZ) {
        case 2:
            // Offset from the origin in pixels
            origin[2] = offsetZ / benchmarkParameters.pixelPitch;
            break;
        case 0:
        case 1:
            // No offset
            break;
        default:
            throw new IllegalStateException();
    }
    if (zeroOffset) {
        list.add(origin.clone());
    }
    if (offsetPoints > 0 && ((offsetRangeX > 0 || offsetRangeY > 0) || is3D && offsetRangeZ > 0)) {
        final double[] min = new double[] { -Math.max(0, offsetRangeX), -Math.max(0, offsetRangeY), -Math.max(0, offsetRangeZ) };
        final double[] range = new double[] { 2 * min[0], 2 * min[1], 2 * min[2] };
        final HaltonSequenceGenerator halton = new HaltonSequenceGenerator((is3D) ? 3 : 2);
        for (int i = 0; i < offsetPoints; i++) {
            final double[] offset = origin.clone();
            final double[] v = halton.nextVector();
            for (int j = 0; j < v.length; j++) {
                offset[j] += v[j] * range[j] + min[j];
            }
            list.add(offset);
        }
    }
    startPoints = list.toArray(new double[0][]);
    return startPoints;
}