Optimize the OptimizingFittingSolver with caching

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++) {

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

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;
 
 	}
 

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<CurveView> 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());
+		}
+
+	}
+
 }