Search in sources :

Example 6 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method calculateSensitivity.

private void calculateSensitivity() {
    if (!settings.calculateSensitivity) {
        return;
    }
    if (!filterAnalysisResult.bestFilter.isEmpty()) {
        final Consumer<String> output = createSensitivityWindow();
        final Ticker ticker = ImageJUtils.createTicker(filterAnalysisResult.bestFilter.size(), 0, "Calculating sensitivity ...");
        for (final String type : filterAnalysisResult.bestFilter.keySet()) {
            final DirectFilter filter = filterAnalysisResult.bestFilter.get(type).getFilter();
            FractionClassificationResult score = scoreFilter(filter, defaultMinimalFilter, fitResultData.resultsList, coordinateStore);
            score = getOriginalScore(score);
            final String message = type + "\t\t\t" + MathUtils.rounded(score.getJaccard(), 4) + "\t\t" + MathUtils.rounded(score.getPrecision(), 4) + "\t\t" + MathUtils.rounded(score.getRecall(), 4);
            output.accept(message);
            // List all the parameters that can be adjusted.
            final int parameters = filter.getNumberOfParameters();
            for (int index = 0; index < parameters; index++) {
                // For each parameter compute as upward + downward delta and get the average gradient
                final DirectFilter higher = (DirectFilter) filter.adjustParameter(index, settings.delta);
                final DirectFilter lower = (DirectFilter) filter.adjustParameter(index, -settings.delta);
                FractionClassificationResult scoreHigher = scoreFilter(higher, defaultMinimalFilter, fitResultData.resultsList, coordinateStore);
                scoreHigher = getOriginalScore(scoreHigher);
                FractionClassificationResult scoreLower = scoreFilter(lower, defaultMinimalFilter, fitResultData.resultsList, coordinateStore);
                scoreLower = getOriginalScore(scoreLower);
                final StringBuilder sb = new StringBuilder();
                sb.append('\t').append(filter.getParameterName(index)).append('\t');
                sb.append(MathUtils.rounded(filter.getParameterValue(index), 4)).append('\t');
                final double dx1 = higher.getParameterValue(index) - filter.getParameterValue(index);
                final double dx2 = filter.getParameterValue(index) - lower.getParameterValue(index);
                addSensitivityScore(sb, score.getJaccard(), scoreHigher.getJaccard(), scoreLower.getJaccard(), dx1, dx2);
                addSensitivityScore(sb, score.getPrecision(), scoreHigher.getPrecision(), scoreLower.getPrecision(), dx1, dx2);
                addSensitivityScore(sb, score.getRecall(), scoreHigher.getRecall(), scoreLower.getRecall(), dx1, dx2);
                output.accept(sb.toString());
            }
            ticker.tick();
        }
        final String message = "-=-=-=-";
        output.accept(message);
        ImageJUtils.finished();
    }
}
Also used : Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) IDirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.IDirectFilter) DirectFilter(uk.ac.sussex.gdsc.smlm.results.filter.DirectFilter) FractionClassificationResult(uk.ac.sussex.gdsc.core.match.FractionClassificationResult)

Example 7 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method scoreFilters.

/**
 * Score filters.
 *
 * @param points the points (must be sorted by duplicate distance)
 * @param createTextResult set to true to create the text result
 * @return the score results
 */
@Nullable
private ParameterScoreResult[] scoreFilters(double[][] points, boolean createTextResult) {
    if (points == null || points.length == 0) {
        return null;
    }
    gaResultsListToScore = gaResultsList;
    gaSubset = false;
    ParameterScoreResult[] scoreResults = new ParameterScoreResult[points.length];
    if (scoreResults.length == 1) {
        // No need to multi-thread this
        final int failCount = (int) Math.round(points[0][0]);
        final double residualsThreshold = points[0][1];
        final double duplicateDistance = points[0][2];
        scoreResults[0] = scoreFilter(searchScoreFilter, defaultMinimalFilter, failCount, residualsThreshold, duplicateDistance, createCoordinateStore(duplicateDistance), createTextResult);
    } else {
        // Multi-thread score all the result
        final int nThreads = getThreads(scoreResults.length);
        final BlockingQueue<ParameterScoreJob> jobs = new ArrayBlockingQueue<>(nThreads * 2);
        final List<Thread> threads = new LinkedList<>();
        final Ticker ticker = ImageJUtils.createTicker(scoreResults.length, nThreads, "Scoring Filters");
        for (int i = 0; i < nThreads; i++) {
            final ParameterScoreWorker worker = new ParameterScoreWorker(jobs, scoreResults, createTextResult, ticker);
            final Thread t = new Thread(worker);
            threads.add(t);
            t.start();
        }
        ticker.start();
        for (int i = 0; i < points.length; i++) {
            if (IJ.escapePressed()) {
                break;
            }
            put(jobs, new ParameterScoreJob(points[i], i));
        }
        // Finish all the worker threads by passing in a null job
        for (int i = 0; i < threads.size(); i++) {
            put(jobs, new ParameterScoreJob(null, -1));
        }
        // Wait for all to finish
        for (int i = 0; i < threads.size(); i++) {
            try {
                threads.get(i).join();
            } catch (final InterruptedException ex) {
                Logger.getLogger(BenchmarkFilterAnalysis.class.getName()).log(Level.WARNING, "Interrupted!", ex);
                Thread.currentThread().interrupt();
                throw new ConcurrentRuntimeException("Unexpected interruption", ex);
            }
        }
        threads.clear();
        ImageJUtils.finished();
        // In case the threads were interrupted
        if (ImageJUtils.isInterrupted()) {
            scoreResults = null;
        }
    }
    finishScoring();
    return scoreResults;
}
Also used : Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) LinkedList(java.util.LinkedList) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue) Nullable(uk.ac.sussex.gdsc.core.annotation.Nullable)

Example 8 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker in project GDSC-SMLM by aherbert.

the class BenchmarkFilterAnalysis method readResults.

private MultiPathFitResults[] readResults() {
    // Extract all the results in memory into a list per frame. This can be cached
    boolean update = false;
    Pair<Integer, TIntObjectHashMap<UniqueIdPeakResult[]>> coords = coordinateCache.get();
    if (coords.getKey() != simulationParameters.id) {
        coords = Pair.of(simulationParameters.id, getCoordinates(results));
        coordinateCache.set(coords);
        update = true;
    }
    actualCoordinates = coords.getValue();
    spotFitResults = BenchmarkSpotFit.getBenchmarkSpotFitResults();
    FitResultData localFitResultData = fitResultDataCache.get();
    final SettingsList scoreSettings = new SettingsList(settings.partialMatchDistance, settings.upperMatchDistance, settings.partialSignalFactor, settings.upperSignalFactor);
    final boolean equalScoreSettings = scoreSettings.equals(localFitResultData.scoreSettings);
    if (update || localFitResultData.fittingId != spotFitResults.id || !equalScoreSettings || localFitResultData.differentSettings(settings)) {
        IJ.showStatus("Reading results ...");
        if (localFitResultData.fittingId < 0) {
            // Copy the settings from the fitter if this is the first run.
            // This just starts the plugin with sensible settings.
            // Q. Should this be per new simulation or fitting result instead?
            final FitEngineConfiguration config = BenchmarkSpotFit.getFitEngineConfiguration();
            settings.failCount = config.getFailuresLimit();
            settings.duplicateDistance = config.getDuplicateDistance();
            settings.duplicateDistanceAbsolute = config.getDuplicateDistanceAbsolute();
            settings.residualsThreshold = (BenchmarkSpotFit.getComputeDoublets()) ? BenchmarkSpotFit.getMultiFilter().residualsThreshold : 1;
        }
        // This functionality is for choosing the optimum filter for the given scoring metric.
        if (!equalScoreSettings) {
            filterAnalysisResult.scores.clear();
        }
        localFitResultData = new FitResultData(spotFitResults.id, scoreSettings, settings);
        // @formatter:off
        // -=-=-=-
        // The scoring is designed to find the best fitter+filter combination for the given spot
        // candidates. The ideal combination would correctly fit+pick all the candidate positions
        // that are close to a localisation.
        // 
        // Use the following scoring scheme for all candidates:
        // 
        // Candidates
        // +----------------------------------------+
        // |   Actual matches                       |
        // |  +-----------+                TN       |
        // |  |  FN       |                         |
        // |  |      +----------                    |
        // |  |      | TP |    | Fitted             |
        // |  +-----------+    | spots              |
        // |         |     FP  |                    |
        // |         +---------+                    |
        // +----------------------------------------+
        // 
        // Candidates     = All the spot candidates
        // Actual matches = Any spot candidate or fitted spot candidate that matches a localisation
        // Fitted spots   = Any spot candidate that was successfully fitted
        // 
        // TP = A spot candidate that was fitted and matches a localisation and is accepted
        // FP = A spot candidate that was fitted but does not match a localisation and is accepted
        // FN = A spot candidate that failed to be fitted but matches a localisation
        // = A spot candidate that was fitted and matches a localisation and is rejected
        // TN = A spot candidate that failed to be fitted and does not match a localisation
        // = A spot candidate that was fitted and does not match a localisation and is rejected
        // 
        // When fitting only produces one result it is possible to compute the TN score.
        // Since unfitted candidates can only be TN or FN we could accumulate these scores and cache
        // them. This was the old method of benchmarking single spot fitting and allowed more scores
        // to be computed.
        // 
        // When fitting produces multiple results then we have to score each fit result against all
        // possible actual results and keep a record of the scores. These can then be assessed when
        // the specific results have been chosen by result filtering.
        // 
        // Using a distance ramped scoring function the degree of match can be varied from 0 to 1.
        // Using a signal-factor ramped scoring function the degree of fitted can be varied from 0
        // to 1. When using ramped scoring functions the fractional allocation of scores using the
        // above scheme is performed, i.e. candidates are treated as if they both match and unmatch.
        // This results in an equivalent to multiple analysis using different thresholds and averaging
        // of the scores.
        // 
        // The totals TP+FP+TN+FN must equal the number of spot candidates. This allows different
        // fitting methods to be compared since the total number of candidates is the same.
        // 
        // Precision = TP / (TP+FP)    : This is always valid as a minimum criteria score
        // Recall    = TP / (TP+FN)    : This is valid between different fitting methods since a
        // method that fits more spots will have a potentially lower FN
        // Jaccard   = TP / (TP+FN+FP) : This is valid between fitting methods
        // 
        // -=-=-=-
        // As an alternative scoring system, different fitting methods can be compared using the same
        // TP value but calculating FN = localisations - TP and FP as Positives - TP. This creates a
        // score against the original number of simulated molecules using everything that was passed
        // through the filter (Positives). This score is comparable when a different spot candidate
        // filter has been used and the total number of candidates is different, e.g. Mean filtering
        // vs. Gaussian filtering
        // -=-=-=-
        // @formatter:on
        final RampedScore distanceScore = RampedScore.of(spotFitResults.distanceInPixels * settings.upperMatchDistance / 100.0, spotFitResults.distanceInPixels * settings.partialMatchDistance / 100.0, false);
        localFitResultData.lowerDistanceInPixels = distanceScore.edge1;
        localFitResultData.distanceInPixels = distanceScore.edge0;
        final double matchDistance = MathUtils.pow2(localFitResultData.distanceInPixels);
        localFitResultData.resultsPrefix3 = "\t" + MathUtils.rounded(distanceScore.edge1 * simulationParameters.pixelPitch) + "\t" + MathUtils.rounded(distanceScore.edge0 * simulationParameters.pixelPitch);
        localFitResultData.limitRange = ", d=" + MathUtils.rounded(distanceScore.edge1 * simulationParameters.pixelPitch) + "-" + MathUtils.rounded(distanceScore.edge0 * simulationParameters.pixelPitch);
        // Signal factor must be greater than 1
        final RampedScore signalScore;
        final double spotSignalFactor = BenchmarkSpotFit.getSignalFactor();
        if (spotSignalFactor > 0 && settings.upperSignalFactor > 0) {
            signalScore = RampedScore.of(spotSignalFactor * settings.upperSignalFactor / 100.0, spotSignalFactor * settings.partialSignalFactor / 100.0, false);
            localFitResultData.lowerSignalFactor = signalScore.edge1;
            localFitResultData.signalFactor = signalScore.edge0;
            localFitResultData.resultsPrefix3 += "\t" + MathUtils.rounded(signalScore.edge1) + "\t" + MathUtils.rounded(signalScore.edge0);
            localFitResultData.limitRange += ", s=" + MathUtils.rounded(signalScore.edge1) + "-" + MathUtils.rounded(signalScore.edge0);
        } else {
            signalScore = null;
            localFitResultData.resultsPrefix3 += "\t0\t0";
            localFitResultData.lowerSignalFactor = localFitResultData.signalFactor = 0;
        }
        // Store all the results
        final ArrayList<MultiPathFitResults> multiPathFitResults = new ArrayList<>(spotFitResults.fitResults.size());
        final List<MultiPathFitResults> syncResults = Collections.synchronizedList(multiPathFitResults);
        // This could be multi-threaded ...
        final int nThreads = getThreads(spotFitResults.fitResults.size());
        final BlockingQueue<Job> jobs = new ArrayBlockingQueue<>(nThreads * 2);
        final List<FitResultsWorker> workers = new LinkedList<>();
        final List<Thread> threads = new LinkedList<>();
        final AtomicInteger uniqueId = new AtomicInteger();
        final CoordinateStore localCoordinateStore = createCoordinateStore();
        final Ticker ticker = ImageJUtils.createTicker(spotFitResults.fitResults.size(), nThreads, null);
        for (int i = 0; i < nThreads; i++) {
            final FitResultsWorker worker = new FitResultsWorker(jobs, syncResults, matchDistance, distanceScore, signalScore, uniqueId, localCoordinateStore.newInstance(), ticker, actualCoordinates);
            final Thread t = new Thread(worker);
            workers.add(worker);
            threads.add(t);
            t.start();
        }
        spotFitResults.fitResults.forEachEntry((frame, candidates) -> {
            put(jobs, new Job(frame, candidates));
            return true;
        });
        // Finish all the worker threads by passing in a null job
        for (int i = 0; i < threads.size(); i++) {
            put(jobs, new Job(0, null));
        }
        // Wait for all to finish
        for (int i = 0; i < threads.size(); i++) {
            try {
                threads.get(i).join();
                final FitResultsWorker worker = workers.get(i);
                localFitResultData.matches += worker.matches;
                localFitResultData.fittedResults += worker.included;
                localFitResultData.totalResults += worker.total;
                localFitResultData.notDuplicateCount += worker.notDuplicateCount;
                localFitResultData.newResultCount += worker.newResultCount;
                localFitResultData.countActual += worker.includedActual;
                if (i == 0) {
                    localFitResultData.depthStats = worker.depthStats;
                    localFitResultData.depthFitStats = worker.depthFitStats;
                    localFitResultData.signalFactorStats = worker.signalFactorStats;
                    localFitResultData.distanceStats = worker.distanceStats;
                } else {
                    localFitResultData.depthStats.add(worker.depthStats);
                    localFitResultData.depthFitStats.add(worker.depthFitStats);
                    localFitResultData.signalFactorStats.add(worker.signalFactorStats);
                    localFitResultData.distanceStats.add(worker.distanceStats);
                }
            } catch (final InterruptedException ex) {
                Thread.currentThread().interrupt();
                throw new ConcurrentRuntimeException("Unexpected interrupt", ex);
            }
        }
        threads.clear();
        ImageJUtils.finished();
        localFitResultData.maxUniqueId = uniqueId.get();
        localFitResultData.resultsList = multiPathFitResults.toArray(new MultiPathFitResults[0]);
        Arrays.sort(localFitResultData.resultsList, (o1, o2) -> Integer.compare(o1.getFrame(), o2.getFrame()));
        MultiPathFilter.resetValidationFlag(localFitResultData.resultsList);
        fitResultDataCache.set(localFitResultData);
    }
    fitResultData = localFitResultData;
    return localFitResultData.resultsList;
}
Also used : ArrayList(java.util.ArrayList) GridCoordinateStore(uk.ac.sussex.gdsc.smlm.results.filter.GridCoordinateStore) CoordinateStore(uk.ac.sussex.gdsc.smlm.results.filter.CoordinateStore) ArrayBlockingQueue(java.util.concurrent.ArrayBlockingQueue) SettingsList(uk.ac.sussex.gdsc.core.utils.SettingsList) FitEngineConfiguration(uk.ac.sussex.gdsc.smlm.engine.FitEngineConfiguration) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) LinkedList(java.util.LinkedList) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) ConcurrentRuntimeException(org.apache.commons.lang3.concurrent.ConcurrentRuntimeException) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) TIntObjectHashMap(gnu.trove.map.hash.TIntObjectHashMap) RampedScore(uk.ac.sussex.gdsc.core.utils.RampedScore) MultiPathFitResults(uk.ac.sussex.gdsc.smlm.results.filter.MultiPathFitResults)

Example 9 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker 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);
}
Also used : UnitSphereSampler(org.apache.commons.rng.sampling.UnitSphereSampler) Color(java.awt.Color) Arrays(java.util.Arrays) ByteProcessor(ij.process.ByteProcessor) Calibration(uk.ac.sussex.gdsc.smlm.data.config.CalibrationProtos.Calibration) IntUnaryOperator(java.util.function.IntUnaryOperator) HistogramPlotBuilder(uk.ac.sussex.gdsc.core.ij.HistogramPlot.HistogramPlotBuilder) IdFramePeakResultComparator(uk.ac.sussex.gdsc.smlm.results.sort.IdFramePeakResultComparator) UnaryOperator(java.util.function.UnaryOperator) RealVector(org.apache.commons.math3.linear.RealVector) Evaluation(org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem.Evaluation) MultivariateJacobianFunction(org.apache.commons.math3.fitting.leastsquares.MultivariateJacobianFunction) VisibleForTesting(uk.ac.sussex.gdsc.core.data.VisibleForTesting) MemoryPeakResults(uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults) NonPositiveDefiniteMatrixException(org.apache.commons.math3.linear.NonPositiveDefiniteMatrixException) LeastSquaresFactory(org.apache.commons.math3.fitting.leastsquares.LeastSquaresFactory) RowSorter(javax.swing.RowSorter) JFrame(javax.swing.JFrame) LutHelper(uk.ac.sussex.gdsc.core.ij.process.LutHelper) KeyStroke(javax.swing.KeyStroke) ConvergenceException(org.apache.commons.math3.exception.ConvergenceException) DistanceUnit(uk.ac.sussex.gdsc.smlm.data.config.UnitProtos.DistanceUnit) KeyEvent(java.awt.event.KeyEvent) WindowAdapter(java.awt.event.WindowAdapter) TextUtils(uk.ac.sussex.gdsc.core.utils.TextUtils) Plot(ij.gui.Plot) TIntHashSet(gnu.trove.set.hash.TIntHashSet) ImagePlus(ij.ImagePlus) DefaultTableCellRenderer(javax.swing.table.DefaultTableCellRenderer) TDoubleArrayList(gnu.trove.list.array.TDoubleArrayList) SumOfSquaredDeviations(uk.ac.sussex.gdsc.core.math.SumOfSquaredDeviations) BasicStroke(java.awt.BasicStroke) RealMatrix(org.apache.commons.math3.linear.RealMatrix) MixtureMultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution) FDistribution(org.apache.commons.math3.distribution.FDistribution) PlugIn(ij.plugin.PlugIn) ActionListener(java.awt.event.ActionListener) MultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution.MultivariateGaussianDistribution) PolygonRoi(ij.gui.PolygonRoi) StoredData(uk.ac.sussex.gdsc.core.utils.StoredData) WindowManager(ij.WindowManager) PeakResult(uk.ac.sussex.gdsc.smlm.results.PeakResult) Supplier(java.util.function.Supplier) PointRoi(ij.gui.PointRoi) Trace(uk.ac.sussex.gdsc.smlm.results.Trace) MultiDialog(uk.ac.sussex.gdsc.core.ij.gui.MultiDialog) UnitSphereSampler(org.apache.commons.rng.sampling.UnitSphereSampler) GenericDialog(ij.gui.GenericDialog) AbstractTableModel(javax.swing.table.AbstractTableModel) SortUtils(uk.ac.sussex.gdsc.core.utils.SortUtils) Overlay(ij.gui.Overlay) IntFunction(java.util.function.IntFunction) Array2DRowRealMatrix(org.apache.commons.math3.linear.Array2DRowRealMatrix) Pair(org.apache.commons.math3.util.Pair) Mean(uk.ac.sussex.gdsc.core.math.Mean) Window(java.awt.Window) LeastSquaresProblem(org.apache.commons.math3.fitting.leastsquares.LeastSquaresProblem) AttributePeakResult(uk.ac.sussex.gdsc.smlm.results.AttributePeakResult) JScrollPane(javax.swing.JScrollPane) ConvergenceChecker(org.apache.commons.math3.optim.ConvergenceChecker) ListSelectionListener(javax.swing.event.ListSelectionListener) PeakResultStoreList(uk.ac.sussex.gdsc.smlm.results.PeakResultStoreList) Optimum(org.apache.commons.math3.fitting.leastsquares.LeastSquaresOptimizer.Optimum) DoubleEquality(uk.ac.sussex.gdsc.core.utils.DoubleEquality) TIntObjectHashMap(gnu.trove.map.hash.TIntObjectHashMap) TIntArrayList(gnu.trove.list.array.TIntArrayList) Mixers(uk.ac.sussex.gdsc.core.utils.rng.Mixers) TextWindow(ij.text.TextWindow) IntConsumer(java.util.function.IntConsumer) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) DataException(uk.ac.sussex.gdsc.core.data.DataException) NonBlockingExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.NonBlockingExtendedGenericDialog) ScreenDimensionHelper(uk.ac.sussex.gdsc.core.ij.gui.ScreenDimensionHelper) MathUtils(uk.ac.sussex.gdsc.core.utils.MathUtils) CalibrationWriter(uk.ac.sussex.gdsc.smlm.data.config.CalibrationWriter) ListSelectionEvent(javax.swing.event.ListSelectionEvent) ArrayRealVector(org.apache.commons.math3.linear.ArrayRealVector) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate) JMenuBar(javax.swing.JMenuBar) BufferedTextWindow(uk.ac.sussex.gdsc.core.ij.BufferedTextWindow) ExtendedGenericDialog(uk.ac.sussex.gdsc.core.ij.gui.ExtendedGenericDialog) JMenu(javax.swing.JMenu) TIntIntHashMap(gnu.trove.map.hash.TIntIntHashMap) MultivariateGaussianMixtureExpectationMaximization(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization) WindowEvent(java.awt.event.WindowEvent) List(java.util.List) SimpleArrayUtils(uk.ac.sussex.gdsc.core.utils.SimpleArrayUtils) JTable(javax.swing.JTable) LUT(ij.process.LUT) TypeConverter(uk.ac.sussex.gdsc.core.data.utils.TypeConverter) Roi(ij.gui.Roi) IntStream(java.util.stream.IntStream) PrecisionResultProcedure(uk.ac.sussex.gdsc.smlm.results.procedures.PrecisionResultProcedure) ParameterValidator(org.apache.commons.math3.fitting.leastsquares.ParameterValidator) TDoubleList(gnu.trove.list.TDoubleList) ValidationUtils(uk.ac.sussex.gdsc.core.utils.ValidationUtils) CompletableFuture(java.util.concurrent.CompletableFuture) LevenbergMarquardtOptimizer(org.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizer) WindowOrganiser(uk.ac.sussex.gdsc.core.ij.plugin.WindowOrganiser) SingularMatrixException(org.apache.commons.math3.linear.SingularMatrixException) AtomicReference(java.util.concurrent.atomic.AtomicReference) SwingConstants(javax.swing.SwingConstants) DoubleUnaryOperator(java.util.function.DoubleUnaryOperator) JMenuItem(javax.swing.JMenuItem) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) DoubleData(uk.ac.sussex.gdsc.core.utils.DoubleData) TFloatArrayList(gnu.trove.list.array.TFloatArrayList) UniformRandomProvider(org.apache.commons.rng.UniformRandomProvider) ChiSquaredDistributionTable(uk.ac.sussex.gdsc.smlm.function.ChiSquaredDistributionTable) LutColour(uk.ac.sussex.gdsc.core.ij.process.LutHelper.LutColour) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) ActionEvent(java.awt.event.ActionEvent) TooManyIterationsException(org.apache.commons.math3.exception.TooManyIterationsException) CalibrationReader(uk.ac.sussex.gdsc.smlm.data.config.CalibrationReader) Consumer(java.util.function.Consumer) ImageWindow(ij.gui.ImageWindow) SimpleRegression(org.apache.commons.math3.stat.regression.SimpleRegression) BinMethod(uk.ac.sussex.gdsc.core.ij.HistogramPlot.BinMethod) HistogramPlot(uk.ac.sussex.gdsc.core.ij.HistogramPlot) ImageJUtils(uk.ac.sussex.gdsc.core.ij.ImageJUtils) TableColumnAdjuster(uk.ac.sussex.gdsc.smlm.ij.gui.TableColumnAdjuster) IJ(ij.IJ) BitSet(java.util.BitSet) DiagonalMatrix(org.apache.commons.math3.linear.DiagonalMatrix) Collections(java.util.Collections) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) UniformRandomProviders(uk.ac.sussex.gdsc.core.utils.rng.UniformRandomProviders) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) CompletableFuture(java.util.concurrent.CompletableFuture) DoubleDoubleBiPredicate(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.DoubleDoubleBiPredicate) MixtureMultivariateGaussianDistribution(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization.MixtureMultivariateGaussianDistribution) AtomicInteger(java.util.concurrent.atomic.AtomicInteger) Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) MultivariateGaussianMixtureExpectationMaximization(uk.ac.sussex.gdsc.smlm.math3.distribution.fitting.MultivariateGaussianMixtureExpectationMaximization)

Example 10 with Ticker

use of uk.ac.sussex.gdsc.core.logging.Ticker 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);
}
Also used : Ticker(uk.ac.sussex.gdsc.core.logging.Ticker) Statistics(uk.ac.sussex.gdsc.core.utils.Statistics) TIntHashSet(gnu.trove.set.hash.TIntHashSet) TIntArrayList(gnu.trove.list.array.TIntArrayList) PeakResult(uk.ac.sussex.gdsc.smlm.results.PeakResult) AttributePeakResult(uk.ac.sussex.gdsc.smlm.results.AttributePeakResult) Trace(uk.ac.sussex.gdsc.smlm.results.Trace) LocalList(uk.ac.sussex.gdsc.core.utils.LocalList) ConvergenceException(org.apache.commons.math3.exception.ConvergenceException) TooManyIterationsException(org.apache.commons.math3.exception.TooManyIterationsException)

Aggregations

Ticker (uk.ac.sussex.gdsc.core.logging.Ticker)33 LinkedList (java.util.LinkedList)16 Future (java.util.concurrent.Future)12 ImageStack (ij.ImageStack)11 ArrayBlockingQueue (java.util.concurrent.ArrayBlockingQueue)10 LocalList (uk.ac.sussex.gdsc.core.utils.LocalList)10 ConcurrentRuntimeException (org.apache.commons.lang3.concurrent.ConcurrentRuntimeException)9 ExecutorService (java.util.concurrent.ExecutorService)8 AtomicInteger (java.util.concurrent.atomic.AtomicInteger)8 ImagePlus (ij.ImagePlus)7 ArrayList (java.util.ArrayList)6 List (java.util.List)6 Statistics (uk.ac.sussex.gdsc.core.utils.Statistics)6 MemoryPeakResults (uk.ac.sussex.gdsc.smlm.results.MemoryPeakResults)6 TIntObjectHashMap (gnu.trove.map.hash.TIntObjectHashMap)5 ImageProcessor (ij.process.ImageProcessor)5 TextWindow (ij.text.TextWindow)5 Point (java.awt.Point)5 Rectangle (java.awt.Rectangle)5 UniformRandomProvider (org.apache.commons.rng.UniformRandomProvider)5