From e577669a3220b27b45ede4ffd9dbe6db14bc94fb Mon Sep 17 00:00:00 2001 From: Nathaniel Sherry Date: Sat, 23 Mar 2024 23:23:45 -0400 Subject: [PATCH] Optimize the OptimizingFittingSolver with caching --- .../cyclops/spectrum/ArraySpectrum.java | 36 +++-- .../framework/cyclops/spectrum/Spectrum.java | 20 ++- .../MultisamplingOptimizingFittingSolver.java | 12 +- .../solver/OptimizingFittingSolver.java | 140 ++++++++++++++++-- 4 files changed, 175 insertions(+), 33 deletions(-) diff --git a/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/ArraySpectrum.java b/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/ArraySpectrum.java index b2a837803..1822eee37 100644 --- a/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/ArraySpectrum.java +++ b/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/ArraySpectrum.java @@ -138,37 +138,42 @@ public ArraySpectrum(SpectrumView copy) } - /** - * Copies the values from the given spectrum into this one. - * Values copied will be in the range of 0 .. min(size(), s.size()) exclusive - * @param s - */ + @Override public void copy(SpectrumView s) { - copy(s.backingArrayCopy()); + copy(((Spectrum) s)); } + + @Override + public void copy(Spectrum s) + { + copy(s.backingArray(), 0, Math.min(this.size, s.size()) - 1); + } + /** * Copies the values from the given spectrum into this one. - * Values copied will be in the range of 0 .. min(size(), s.size()) exclusive * @param s + * @param first the first index to copy + * @param last the last index to copy */ @Override - public void copy(Spectrum s) + public void copy(Spectrum s, int first, int last) { - copy(s.backingArray()); + copy(s.backingArray(), first, last); } /** * Copies the given array into the start of this Spectrum's data array - * up to min(this.size(), array.length) exclusive; * @param array the array to copy from + * @param first the first index to copy + * @param last the last index to copy */ - private void copy(float[] array) { + private void copy(float[] array, int first, int last) { int length = Math.min(this.data.length, array.length); - System.arraycopy(array, 0, this.data, 0, length); - maxIndex = Math.max(maxIndex, array.length-1); + System.arraycopy(array, first, this.data, first, last - first + 1); + maxIndex = Math.max(maxIndex, last); } /** @@ -435,6 +440,11 @@ public void zero() { Arrays.fill(data, 0f); } + @Override + public void zero(int first, int last) { + Arrays.fill(data, first, last+1, 0f); + } + public float sum() { float sum = 0; for (int i = 0; i < size; i++) { diff --git a/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/Spectrum.java b/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/Spectrum.java index 71d52865c..dea9b44b9 100644 --- a/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/Spectrum.java +++ b/Framework/Cyclops/Cyclops/src/main/java/org/peakaboo/framework/cyclops/spectrum/Spectrum.java @@ -18,7 +18,18 @@ public interface Spectrum extends SpectrumView { * Values copied will be in the range of 0 .. min(size(), s.size()) exclusive * @param s */ - void copy(Spectrum s); + default void copy(Spectrum s) { + this.copy(s, 0, Math.min(this.size(), s.size()) - 1); + } + + /** + * Copies the values from the given spectrum into this one. + * Values copied will be in the range of 0 .. min(size(), s.size()) exclusive + * @param s + * @param first the first index to copy + * @param last the last index to copy + */ + void copy(Spectrum s, int first, int last); /** * Adds a value to the Spectrum. When a new spectrum is created @@ -58,6 +69,13 @@ public interface Spectrum extends SpectrumView { */ void zero(); + /** + * Populate all entries in a spectrum with 0s + * This will not impact the add cursor's location + */ + void zero(int first, int last); + + /** * Converts a list of integers (points) to a spectrum with matching indices set to 1 diff --git a/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/MultisamplingOptimizingFittingSolver.java b/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/MultisamplingOptimizingFittingSolver.java index 62c218a2a..e4359f21c 100644 --- a/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/MultisamplingOptimizingFittingSolver.java +++ b/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/MultisamplingOptimizingFittingSolver.java @@ -7,7 +7,6 @@ import org.apache.commons.math3.analysis.MultivariateFunction; import org.apache.commons.math3.optim.PointValuePair; import org.peakaboo.curvefit.curve.fitting.CurveView; -import org.peakaboo.curvefit.curve.fitting.FittingResultSetView; public class MultisamplingOptimizingFittingSolver extends OptimizingFittingSolver { @@ -38,17 +37,15 @@ public String pluginUUID() { return "87eeb1e0-6c4e-4f80-9cf7-8c19a07423b5"; } + @Override - public FittingResultSetView solve(FittingSolverContext inputCtx) { + public double[] calculateWeights(OptimizingFittingSolver.Context inputCtx) { int size = inputCtx.fittings.getVisibleCurves().size(); - if (size == 0) { - return getEmptyResult(inputCtx); - } // Create a shallow copy of the input context and then make a deep copy of the // curve list so that we can permute it without impacting the original - FittingSolverContext permCtx = new FittingSolverContext(inputCtx); + OptimizingFittingSolver.Context permCtx = new OptimizingFittingSolver.Context(inputCtx); permCtx.curves = new ArrayList<>(permCtx.curves); int counter = 0; @@ -82,8 +79,7 @@ public FittingResultSetView solve(FittingSolverContext inputCtx) { scalings[i] /= counter; } - return evaluate(scalings, inputCtx); - + return scalings; } diff --git a/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/OptimizingFittingSolver.java b/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/OptimizingFittingSolver.java index b634343ac..bf287fea0 100644 --- a/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/OptimizingFittingSolver.java +++ b/LibPeakaboo/src/main/java/org/peakaboo/curvefit/curve/fitting/solver/OptimizingFittingSolver.java @@ -18,10 +18,13 @@ import org.peakaboo.curvefit.curve.fitting.FittingResultSet; import org.peakaboo.curvefit.curve.fitting.FittingResultSetView; import org.peakaboo.curvefit.curve.fitting.FittingResultView; +import org.peakaboo.curvefit.curve.fitting.FittingSetView; +import org.peakaboo.curvefit.curve.fitting.fitter.CurveFitter; import org.peakaboo.curvefit.curve.fitting.fitter.CurveFitter.CurveFitterContext; import org.peakaboo.framework.cyclops.spectrum.ArraySpectrum; import org.peakaboo.framework.cyclops.spectrum.Spectrum; import org.peakaboo.framework.cyclops.spectrum.SpectrumCalculations; +import org.peakaboo.framework.cyclops.spectrum.SpectrumView; public class OptimizingFittingSolver implements FittingSolver { @@ -60,16 +63,17 @@ public FittingResultSetView solve(FittingSolverContext ctx) { return getEmptyResult(ctx); } - double[] weights = calculateWeights(ctx); + Context octx = new Context(ctx); + double[] weights = calculateWeights(octx); return evaluate(weights, ctx); } - public double[] calculateWeights(FittingSolverContext ctx) { + public double[] calculateWeights(Context ctx) { EvaluationSpace eval = new EvaluationSpace(ctx.data.size()); MultivariateFunction cost = getCostFunction(ctx, eval); double[] guess = getInitialGuess(ctx); - + PointValuePair result = optimizeCostFunction(cost, guess, costFnPrecision); return result.getPoint(); @@ -98,8 +102,18 @@ protected PointValuePair optimizeCostFunction(MultivariateFunction cost, double[ GoalType.MINIMIZE); - //265 reps on test session but occasionally dies? - //optimizer = new BOBYQAOptimizer(Math.max(size+2, size*2)); + //? reps on test session but occasionally dies? +// int size = guess.length; +// var optimizer = new BOBYQAOptimizer(Math.max(size+2, size*2)); +// return optimizer.optimize( +// new ObjectiveFunction(cost), +// new InitialGuess(guess), +// new MaxIter(1000000), +// new MaxEval(1000000), +// new NonNegativeConstraint(true), +// SimpleBounds.unbounded(size), +// GoalType.MINIMIZE +// ); //-1 for rel means don't use rel, just use abs difference // PointValuePair result = new SimplexOptimizer(-1d, 1d).optimize( @@ -137,7 +151,7 @@ public static double[] getInitialGuess(FittingSolverContext ctx) { - protected MultivariateFunction getCostFunction(FittingSolverContext ctx, EvaluationSpace eval) { + protected MultivariateFunction getCostFunction(Context ctx, EvaluationSpace eval) { return new MultivariateFunction() { @Override @@ -151,7 +165,7 @@ public double value(double[] point) { } } - test(point, ctx, eval); + testPowell(point, ctx, eval); float score = score(point, ctx, eval.residual); if (containsNegatives > 0) { return score * (1f+containsNegatives); @@ -166,9 +180,10 @@ public double value(double[] point) { /** - * Calculate the residual from data (signal) and total (fittings). Store the result in residual + * Calculate the residual from data (signal) and total (fittings). Store the result in residual. + * This is a generic version of the test method, untuned for any specific optimizing method */ - private void test(double[] weights, FittingSolverContext ctx, EvaluationSpace eval) { + private void testGeneric(double[] weights, FittingSolverContext ctx, EvaluationSpace eval) { Spectrum total = eval.total; total.zero(); @@ -188,6 +203,87 @@ private void test(double[] weights, FittingSolverContext ctx, EvaluationSpace ev SpectrumCalculations.subtractLists_target(ctx.data, eval.total, eval.residual, first, last); } + + /** + * Calculate the residual from data (signal) and total (fittings). Store the result in residual. + * Note that this version of the test function has been tuned for the PowellOptimizer. + */ + private void testPowell(double[] weights, Context ctx, EvaluationSpace eval) { + + //When there are no intense channels to consider, the residual will be equal to the data + if (ctx.channels.length == 0) { + SpectrumCalculations.subtractFromList_target(ctx.data, eval.residual, 0f); + return; + } + + int curvesLength = weights.length; + double[] lastWeights = ctx.lastWeights; + + // First we check if the weights up to partialIndex are a match + int lastMatchingIndex = -1; + for (int i = 0; i < curvesLength; i++) { + if (weights[i] == lastWeights[i]) { + lastMatchingIndex = i; + } else { + break; + } + } + + int startingCurveIndex = 0; + Spectrum total = eval.total; + int first = ctx.channels[0]; + int last = ctx.channels[ctx.channels.length-1]; + if (lastMatchingIndex >= ctx.partialIndex && ctx.partialIndex > -1) { + // If the weights are a match up to and including the index at which we have + // cached a partial sum, then we can skip some of the work and jump part way + // into the calculations + total.copy(ctx.partial, first, last); + startingCurveIndex = ctx.partialIndex+1; + } else { + // If the weights aren't a match up to and including the index at which we have + // a cached partial value, then we must discard the cached value and start over + ctx.partialIndex = -1; + total.zero(first, last); + } + + // Do all the weights up to the current index match the last run through + boolean matching = true; + List curves = ctx.curves; + for (int i = startingCurveIndex; i < curvesLength; i++) { + double weight = weights[i]; + + // Does this weight match the corresponding weight from the last run through? + boolean match = weight == lastWeights[i]; + + if (matching && !match && i > startingCurveIndex) { + // This is the first weight which does not match the last run through, so we + // save the totals up to but not including this point before adding this curve. + // NB: We skip doing this if this is the first curve in this (partial?) loop, as + // we'll have nothing to add to what's already cached + ctx.partial.copy(total, first, last); + ctx.partialIndex = i-1; + } + + // Update matching with this rounds weight match check + matching &= match; + + // Scale the curve by the weight and add it to the total + curves.get(i).scaleOnto((float)weight, total, first, last); + + } + + // NB: We explicitly don't cache a full match here. The Optimizer will sometimes + // repeat a set of weights. This is not often enough to merit caching every full + // result but often enough that we don't want to invalidate our partial cached + // results + + // Update the weights to the ones from this pass so that the next loop through + // will be able to compare its values to these this runs values + System.arraycopy(weights, 0, ctx.lastWeights, 0, weights.length); + + SpectrumCalculations.subtractLists_target(ctx.data, eval.total, eval.residual, first, last); + + } /** @@ -236,14 +332,36 @@ public static FittingResultSetView evaluate(double[] point, FittingSolverContext } public static class EvaluationSpace { - public Spectrum scratch; public Spectrum total; public Spectrum residual; public EvaluationSpace(int size) { - this.scratch = new ArraySpectrum(size); this.total = new ArraySpectrum(size); this.residual = new ArraySpectrum(size); } } + + public static class Context extends FittingSolverContext { + + public Spectrum partial; + double[] lastWeights; + int partialIndex = -1; + + public Context(SpectrumView data, FittingSetView fittings, CurveFitter fitter) { + super(data, fittings, fitter); + init(); + } + + public Context(FittingSolverContext copy) { + super(copy); + init(); + } + + private void init() { + this.lastWeights = new double[curves.size()]; + this.partial = new ArraySpectrum(data.size()); + } + + } + }