Fix/armodel transients

.github/workflows/flucoma-core-ci.yml

@@ -35,5 +35,5 @@ jobs:
 
     - name: Test
       working-directory: ${{github.workspace}}/build
-      run: ctest -C ${{env.BUILD_TYPE}} -j3
+      run: ctest -C ${{env.BUILD_TYPE}} -j3 --output-on-failure
 

include/algorithms/public/TransientExtraction.hpp

@@ -82,6 +82,9 @@ class TransientExtraction
 
   index detect(const double* input, index inSize)
   {
+    assert(mInitialized && "Call init() before processing");
+    assert(modelOrder() >= mDetectHalfWindow &&
+           "model order needs to be >= half filter size");
     frame(input, inSize);
     analyze();
     detection();
@@ -91,7 +94,9 @@ class TransientExtraction
   void process(const RealVectorView input, RealVectorView transients,
                RealVectorView residual)
   {
-    assert(mInitialized);
+    assert(mInitialized && "Call init() before processing");
+    assert(modelOrder() >= mDetectHalfWindow &&
+           "model order needs to be >= half filter size");
     index inSize = input.extent(0);
     frame(input.data(), inSize);
     analyze();
@@ -102,6 +107,9 @@ class TransientExtraction
   void process(const RealVectorView input, const RealVectorView unknowns,
                RealVectorView transients, RealVectorView residual)
   {
+    assert(mInitialized && "Call init() before processing");
+    assert(modelOrder() >= mDetectHalfWindow &&
+           "model order needs to be >= half filter size");
     index inSize = input.extent(0);
     std::copy(unknowns.data(), unknowns.data() + hopSize(), mDetect.data());
     mCount = 0;
@@ -130,7 +138,8 @@ class TransientExtraction
   void analyze()
   {
     mModel.setMinVariance(0.0000001);
-    mModel.estimate(mInput.data() + modelOrder(), analysisSize());
+    mModel.estimate(FluidTensorView<const double, 1>(
+        mInput.data(), modelOrder(), analysisSize()));
   }
 
   void detection()
@@ -139,18 +148,27 @@ class TransientExtraction
 
     // Forward and backward error
     const double normFactor = 1.0 / sqrt(mModel.variance());
-    errorCalculation<&ARModel::forwardErrorArray>(
-        mForwardError.data(), input, blockSize() + mDetectHalfWindow + 1,
-        normFactor);
-    errorCalculation<&ARModel::backwardErrorArray>(
-        mBackwardError.data(), input, blockSize() + mDetectHalfWindow + 1,
-        normFactor);
+
+    auto inputView = FluidTensorView<const double, 1>(
+        mInput.data(), modelOrder() + padSize(), blockSize() + modelOrder());
+    auto fwdError = FluidTensorView<double, 1>(mForwardError.data(), 0,
+                                               blockSize() + mDetectHalfWindow);
+    auto backError = FluidTensorView<double, 1>(
+        mBackwardError.data(), 0, blockSize() + mDetectHalfWindow);
+
+    errorCalculation<&ARModel::forwardErrorArray>(inputView, fwdError,
+                                                  normFactor);
+    errorCalculation<&ARModel::backwardErrorArray>(inputView, backError,
+                                                   normFactor);
 
     // Window error functions (brute force convolution)
-    windowError(mForwardWindowedError.data(),
-                mForwardError.data() + modelOrder(), hopSize());
-    windowError(mBackwardWindowedError.data(),
-                mBackwardError.data() + modelOrder(), hopSize());
+    auto fwdWindowedError =
+        FluidTensorView<double, 1>(mForwardWindowedError.data(), 0, hopSize());
+    auto backWindowedError =
+        FluidTensorView<double, 1>(mBackwardWindowedError.data(), 0, hopSize());
+
+    windowError(fwdError(Slice(modelOrder())), fwdWindowedError);
+    windowError(backError(Slice(modelOrder())), backWindowedError);
 
     // Detection
     index        count = 0;
@@ -164,7 +182,9 @@ class TransientExtraction
     {
       if (!click && (mBackwardWindowedError[asUnsigned(i)] > loThresh) &&
           (mForwardWindowedError[asUnsigned(i)] > hiThresh))
-      { click = true; }
+      {
+        click = true;
+      }
       else if (click && (mBackwardWindowedError[asUnsigned(i)] < loThresh))
       {
         click = false;
@@ -284,7 +304,7 @@ class TransientExtraction
         residual[i] = input[i + order];
     }
 
-    if (mRefine) refine(residual, size, Au, u);
+    if (mRefine) refine(FluidTensorView<double, 1>(residual, 0, size), Au, u);
 
     for (index i = 0; i < (size - order); i++)
       transients[i] = input[i + order] - residual[i];
@@ -294,17 +314,19 @@ class TransientExtraction
               mInput.data() + padSize() + order + order);
   }
 
-  void refine(double* io, index size, Eigen::MatrixXd& Au, Eigen::MatrixXd& ls)
+  void refine(FluidTensorView<double, 1> io, Eigen::MatrixXd& Au,
+              Eigen::MatrixXd& ls)
   {
     const double energy = mModel.variance() * mCount;
     double       energyLS = 0.0;
     index        order = modelOrder();
+    index        size = io.size();
 
     for (index i = 0; i < (size - order); i++)
     {
       if (mDetect[asUnsigned(i)] != 0)
       {
-        const double error = mModel.forwardError(io + i);
+        const double error = mModel.forwardError(io(Slice(i)));
         energyLS += error * error;
       }
     }
@@ -341,44 +363,50 @@ class TransientExtraction
     return sum;
   }
 
-  template <void (ARModel::*Method)(double*, const double*, index)>
-  void errorCalculation(double* error, const double* input, index size,
-                        double normFactor)
+  template <void (ARModel::*Method)(FluidTensorView<const double, 1>,
+                                    FluidTensorView<double, 1>)>
+  void errorCalculation(FluidTensorView<const double, 1> input,
+                        FluidTensorView<double, 1> error, double normFactor)
   {
-    (mModel.*Method)(error, input, size);
+    (mModel.*Method)(input, error);
 
     // Take absolutes and normalise
-    for (index i = 0; i < size; i++)
+    for (index i = 0; i < error.size(); i++)
       error[i] = std::fabs(error[i]) * normFactor;
   }
 
   // Triangle window
   double calcWindow(double norm) { return std::min(norm, 1.0 - norm); }
 
-  void windowError(double* errorWindowed, const double* error, index size)
+  void windowError(FluidTensorView<const double, 1> error,
+                   FluidTensorView<double, 1>       errorWindowed)
   {
-    const index  windowSize = mDetectHalfWindow * 2 + 1;
-    const index  windowOffset = mDetectHalfWindow;
+    assert(error.descriptor().start >= mDetectHalfWindow &&
+           "insufficient offset for filter size");
+    assert(error.size() > errorWindowed.size() - 1 + (mDetectHalfWindow) &&
+           "insufficient input for filter size");
+
+    const index  windowSize = mDetectHalfWindow * 2 + 3;
+    const index  windowOffset = mDetectHalfWindow + 1;
     const double powFactor = mDetectPowerFactor;
 
     // Calculate window normalisation factor
     double windowNormFactor = 0.0;
 
     for (index j = 0; j < windowSize; j++)
-      windowNormFactor += calcWindow((double) j / windowSize);
+      windowNormFactor += calcWindow((double) j / (windowSize - 1));
 
     windowNormFactor = 1.0 / windowNormFactor;
 
     // Do window processing
-    for (index i = 0; i < size; i++)
+    for (index i = 0; i < errorWindowed.size(); i++)
     {
       double windowed = 0.0;
 
-      for (index j = 1; j < windowSize; j++)
+      for (index j = 1; j < windowSize - 1; j++)
       {
         const double value = pow(fabs(error[i - windowOffset + j]), powFactor);
-        windowed += value * calcWindow((double) j / windowSize);
-        ;
+        windowed += value * calcWindow((double) j / (windowSize - 1));
       }
 
       errorWindowed[i] = pow((windowed * windowNormFactor), 1.0 / powFactor);