From e01e237ecd687d190cb3e3ffa34a0921d8536088 Mon Sep 17 00:00:00 2001 From: Alex Harker Date: Tue, 1 Mar 2022 13:49:19 +0000 Subject: [PATCH 01/12] Correct indexing bugs due to incorrect refactor in robustIteration() --- include/algorithms/util/ARModel.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 76296b350..973f715cd 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -204,8 +204,8 @@ class ARModel // Iterate to find new filtered input for (index i = 0; i < size; i++) { - const double prediction = fowardPrediction(estimates); - estimates[0] = + const double prediction = fowardPrediction(estimates + i); + estimates[i] = prediction + robustResidual(input[i], prediction, robustFactor * sqrt(mVariance)); } From 81f4538a9a9c4e96d15b1a8f8118aa8a6bc73647 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 2 Mar 2022 09:25:05 +0000 Subject: [PATCH 02/12] ARModel: fix model stability by reverting mistaken array indexing change --- include/algorithms/util/ARModel.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 973f715cd..255578760 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -163,7 +163,7 @@ class ARModel directEstimate(input, size, true); // Initialise Estimates - for (index i = mParameters.size(); i < mParameters.size() + size; i++) + for (index i = 0; i < mParameters.size() + size; i++) estimates[asUnsigned(i)] = input[i - mParameters.size()]; // Variance From 0612acd4f369a3b6fddeaf1cb1a6098874b5bca3 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 2 Mar 2022 17:09:20 +0000 Subject: [PATCH 03/12] ARModel: make interface use FluidTensorView --- .../algorithms/public/TransientExtraction.hpp | 46 +++-- include/algorithms/util/ARModel.hpp | 158 +++++++++++------- 2 files changed, 127 insertions(+), 77 deletions(-) diff --git a/include/algorithms/public/TransientExtraction.hpp b/include/algorithms/public/TransientExtraction.hpp index 8bb39fb02..824a233ec 100644 --- a/include/algorithms/public/TransientExtraction.hpp +++ b/include/algorithms/public/TransientExtraction.hpp @@ -130,7 +130,8 @@ class TransientExtraction void analyze() { mModel.setMinVariance(0.0000001); - mModel.estimate(mInput.data() + modelOrder(), analysisSize()); + mModel.estimate(FluidTensorView( + mInput.data(), modelOrder(), analysisSize())); } void detection() @@ -139,12 +140,22 @@ 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); + + assert(mInput.size() > modelOrder() + padSize() + + blockSize() + mDetectHalfWindow + 1 + modelOrder()); + + auto inputView = FluidTensorView( + mInput.data(), modelOrder() + padSize(), + blockSize() + mDetectHalfWindow + 1 + modelOrder()); + auto fwdError = FluidTensorView( + mForwardError.data(), 0, blockSize() + mDetectHalfWindow + 1); + auto backError = FluidTensorView( + mBackwardError.data(), 0, blockSize() + mDetectHalfWindow + 1); + + errorCalculation<&ARModel::forwardErrorArray>(inputView, fwdError, + normFactor); + errorCalculation<&ARModel::backwardErrorArray>(inputView, backError, + normFactor); // Window error functions (brute force convolution) windowError(mForwardWindowedError.data(), @@ -284,7 +295,7 @@ class TransientExtraction residual[i] = input[i + order]; } - if (mRefine) refine(residual, size, Au, u); + if (mRefine) refine(FluidTensorView(residual,0, size), Au, u); for (index i = 0; i < (size - order); i++) transients[i] = input[i + order] - residual[i]; @@ -294,17 +305,19 @@ class TransientExtraction mInput.data() + padSize() + order + order); } - void refine(double* io, index size, Eigen::MatrixXd& Au, Eigen::MatrixXd& ls) + void refine(FluidTensorView 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,14 +354,15 @@ class TransientExtraction return sum; } - template - void errorCalculation(double* error, const double* input, index size, - double normFactor) + template , + FluidTensorView)> + void errorCalculation(FluidTensorView input, + FluidTensorView 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; } diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 255578760..784cb445e 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -14,6 +14,8 @@ under the European Union’s Horizon 2020 research and innovation programme #include "Toeplitz.hpp" #include "../public/WindowFuncs.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidTensor.hpp" +#include "FluidEigenMappings.hpp" #include #include #include @@ -36,80 +38,111 @@ class ARModel double variance() const { return mVariance; } index order() const { return mParameters.size(); } - void estimate(const double* input, index size, index nIterations = 3, + void estimate(FluidTensorView input, index nIterations = 3, double robustFactor = 3.0) { if (nIterations > 0) - robustEstimate(input, size, nIterations, robustFactor); + robustEstimate(input, nIterations, robustFactor); else - directEstimate(input, size, true); + directEstimate(input, true); } - - double fowardPrediction(const double* input) + + + double fowardPrediction(FluidTensorView input) { - return modelPredict>(input); + double prediction; + modelPredict(input, FluidTensorView(&prediction, 0, 1), + std::negate{}, Predict{}); + return prediction; } - double backwardPrediction(const double* input) + double backwardPrediction(FluidTensorView input) { - struct Identity - { - index operator()(index a) { return a; } - }; - return modelPredict(input); + double prediction; + modelPredict(input, FluidTensorView(&prediction, 0, 1), + Identity{}, Predict{}); + return prediction; } - double forwardError(const double* input) + double forwardError(FluidTensorView input) { - return modelError<&ARModel::fowardPrediction>(input); + double error; + forwardErrorArray(input,FluidTensorView(&error,0,1)); + return error; } - double backwardError(const double* input) + double backwardError(FluidTensorView input) { - return modelError<&ARModel::backwardPrediction>(input); + double error; + backwardErrorArray(input,FluidTensorView(&error,0,1)); + return error; } - void forwardErrorArray(double* errors, const double* input, index size) + void forwardErrorArray(FluidTensorView input,FluidTensorView errors) { - modelErrorArray<&ARModel::forwardError>(errors, input, size); + modelPredict(input, errors, std::negate{},Error{}); } - void backwardErrorArray(double* errors, const double* input, index size) + void backwardErrorArray(FluidTensorView input, FluidTensorView errors) { - modelErrorArray<&ARModel::backwardError>(errors, input, size); + modelPredict(input, errors, Identity{}, Error{}); } void setMinVariance(double variance) { mMinVariance = variance; } private: - template - double modelPredict(const double* input) + + struct Identity { - double estimate = 0.0; - - for (index i = 0; i < mParameters.size(); i++) - estimate += mParameters(i) * input[Op()(i + 1)]; - - return estimate; - } - - template - double modelError(const double* input) + index operator()(index a) { return a; } + }; + + struct Predict { - return input[0] - (this->*Method)(input); - } - - template - void modelErrorArray(double* errors, const double* input, index size) + double operator()(double, double estimate) {return estimate;} + }; + + struct Error{ + double operator()(double input, double estimate) {return input - estimate;} + }; + + + /// \pre Op(numPredictons + mParameters.size()) < input.size() && Op(mParameters.size()) >= -input.descriptor().start + template + void modelPredict(FluidTensorView input, + FluidTensorView output, Indexer f_idx, + OutputFn f_out) { - for (index i = 0; i < size; i++) errors[i] = (this->*Method)(input + i); + + index numPredictions = output.size(); + +// std::cout << ((numPredictions - 1) + f_idx(mParameters.size())) << '\t' << input.size() << '\n'; +// std::cout << f_idx(mParameters.size()) << '\t' << -input.descriptor().start << '\n'; + + assert(((numPredictions - 1) + f_idx(mParameters.size())) <= input.size() && + "array bounds error in AR model prediction: input too short"); + assert(f_idx(mParameters.size()) >= -input.descriptor().start && + "array bounds error in AR model prediction: input offset too small"); + + const double* input_ptr = input.data(); + + for(index p = 0; p < numPredictions; p++) + { + double estimate = 0; + for (index i = 0; i < mParameters.size(); i++) + estimate += mParameters(i) * (input_ptr + p)[f_idx(i + 1)]; + output[p] = f_out(input_ptr[p], estimate); + } } - void directEstimate(const double* input, index size, bool updateVariance) + void directEstimate(FluidTensorView input, bool updateVariance) { + + index size = input.size(); + // copy input to a 32 byte aligned block (otherwise risk segfaults on Linux) - VectorXd frame = Eigen::Map(input, size); - + VectorXd frame = _impl::asEigen(input); + if (mUseWindow) { if (mWindow.size() != size) @@ -154,25 +187,26 @@ class ARModel } } - void robustEstimate(const double* input, index size, index nIterations, + void robustEstimate(FluidTensorView input, index nIterations, double robustFactor) { - std::vector estimates(asUnsigned(size + mParameters.size())); - + FluidTensor estimates(input.size() + mParameters.size()); + // Calculate an initial estimate of parameters - directEstimate(input, size, true); + directEstimate(input, true); + assert(input.descriptor().start >= mParameters.size()&&"too little offset into input data"); + // Initialise Estimates - for (index i = 0; i < mParameters.size() + size; i++) - estimates[asUnsigned(i)] = input[i - mParameters.size()]; + for (index i = 0; i < input.size() + mParameters.size(); i++) + estimates[i] = input.data()[i - mParameters.size()]; // Variance - robustVariance(estimates.data() + mParameters.size(), input, size, - robustFactor); + robustVariance(estimates(Slice(mParameters.size())), input, robustFactor); // Iterate for (index iterations = nIterations; iterations--;) - robustIteration(estimates.data() + mParameters.size(), input, size, + robustIteration(estimates(Slice(mParameters.size())), input, robustFactor); } @@ -181,38 +215,40 @@ class ARModel return cs * psiFunction((input - prediction) / cs); } - void robustVariance(double* estimates, const double* input, index size, - double robustFactor) + void robustVariance(FluidTensorView estimates, + FluidTensorView input, + double robustFactor) { double residualSqSum = 0.0; // Iterate to find new filtered input - for (index i = 0; i < size; i++) + for (index i = 0; i < input.size(); i++) { const double residual = - robustResidual(input[i], fowardPrediction(estimates + i), + robustResidual(input[i], fowardPrediction(estimates(Slice(i))), robustFactor * sqrt(mVariance)); residualSqSum += residual * residual; } - setVariance(residualSqSum / size); + setVariance(residualSqSum / input.size()); } - void robustIteration(double* estimates, const double* input, index size, - double robustFactor) + void robustIteration(FluidTensorView estimates, + FluidTensorView input, + double robustFactor) { // Iterate to find new filtered input - for (index i = 0; i < size; i++) + for (index i = 0; i < input.size(); i++) { - const double prediction = fowardPrediction(estimates + i); + const double prediction = fowardPrediction(estimates(Slice(i))); estimates[i] = prediction + robustResidual(input[i], prediction, robustFactor * sqrt(mVariance)); } // New parameters - directEstimate(estimates, size, false); - robustVariance(estimates, input, size, robustFactor); + directEstimate(estimates(Slice(0,input.size())), false); + robustVariance(estimates(Slice(0,input.size())), input, robustFactor); } void setVariance(double variance) From a621a4519afd1f93c543ffca4b4772d94de39005 Mon Sep 17 00:00:00 2001 From: Alex Harker Date: Mon, 14 Mar 2022 15:24:34 +0000 Subject: [PATCH 04/12] Formatting --- include/algorithms/util/ARModel.hpp | 36 ++++++++++++----------------- 1 file changed, 15 insertions(+), 21 deletions(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 784cb445e..bd238505c 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -47,7 +47,6 @@ class ARModel directEstimate(input, true); } - double fowardPrediction(FluidTensorView input) { double prediction; @@ -66,26 +65,26 @@ class ARModel double forwardError(FluidTensorView input) { - double error; - forwardErrorArray(input,FluidTensorView(&error,0,1)); - return error; + double error; + forwardErrorArray(input,FluidTensorView(&error,0,1)); + return error; } double backwardError(FluidTensorView input) { - double error; - backwardErrorArray(input,FluidTensorView(&error,0,1)); - return error; + double error; + backwardErrorArray(input,FluidTensorView(&error,0,1)); + return error; } void forwardErrorArray(FluidTensorView input,FluidTensorView errors) { - modelPredict(input, errors, std::negate{},Error{}); + modelPredict(input, errors, std::negate{},Error{}); } void backwardErrorArray(FluidTensorView input, FluidTensorView errors) { - modelPredict(input, errors, Identity{}, Error{}); + modelPredict(input, errors, Identity{}, Error{}); } void setMinVariance(double variance) { mMinVariance = variance; } @@ -106,14 +105,11 @@ class ARModel double operator()(double input, double estimate) {return input - estimate;} }; - - /// \pre Op(numPredictons + mParameters.size()) < input.size() && Op(mParameters.size()) >= -input.descriptor().start template void modelPredict(FluidTensorView input, - FluidTensorView output, Indexer f_idx, - OutputFn f_out) + FluidTensorView output, Indexer fIdx, + OutputFn fOut) { - index numPredictions = output.size(); // std::cout << ((numPredictions - 1) + f_idx(mParameters.size())) << '\t' << input.size() << '\n'; @@ -126,18 +122,17 @@ class ARModel const double* input_ptr = input.data(); - for(index p = 0; p < numPredictions; p++) + for (index p = 0; p < numPredictions; p++) { double estimate = 0; for (index i = 0; i < mParameters.size(); i++) - estimate += mParameters(i) * (input_ptr + p)[f_idx(i + 1)]; - output[p] = f_out(input_ptr[p], estimate); + estimate += mParameters(i) * (input_ptr + p)[fIdx(i + 1)]; + output[p] = fOut(input_ptr[p], estimate); } } void directEstimate(FluidTensorView input, bool updateVariance) { - index size = input.size(); // copy input to a 32 byte aligned block (otherwise risk segfaults on Linux) @@ -154,7 +149,6 @@ class ARModel frame.array() *= mWindow; } - VectorXd autocorrelation(size); algorithm::autocorrelateReal(autocorrelation.data(), frame.data(), asUnsigned(size)); @@ -247,8 +241,8 @@ class ARModel } // New parameters - directEstimate(estimates(Slice(0,input.size())), false); - robustVariance(estimates(Slice(0,input.size())), input, robustFactor); + directEstimate(estimates(Slice(0, input.size())), false); + robustVariance(estimates(Slice(0, input.size())), input, robustFactor); } void setVariance(double variance) From 19a6d03624ef7d14a3010f8e189e9a781cedf786 Mon Sep 17 00:00:00 2001 From: Alex Harker Date: Mon, 14 Mar 2022 15:39:36 +0000 Subject: [PATCH 05/12] Test ranges --- include/algorithms/util/ARModel.hpp | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index bd238505c..b665b7613 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -112,14 +112,23 @@ class ARModel { index numPredictions = output.size(); -// std::cout << ((numPredictions - 1) + f_idx(mParameters.size())) << '\t' << input.size() << '\n'; -// std::cout << f_idx(mParameters.size()) << '\t' << -input.descriptor().start << '\n'; - - assert(((numPredictions - 1) + f_idx(mParameters.size())) <= input.size() && - "array bounds error in AR model prediction: input too short"); - assert(f_idx(mParameters.size()) >= -input.descriptor().start && - "array bounds error in AR model prediction: input offset too small"); - + assert(fIdx(1) >= -input.descriptor().start && + "array bounds error in AR model prediction: input too short"); + assert(fIdx(1) < input.size() && + "array bounds error in AR model prediction: input too short"); + assert(fIdx(mParameters.size()) >= -input.descriptor().start && + "array bounds error in AR model prediction: input too short"); + assert(fIdx(mParameters.size()) < input.size() && + "array bounds error in AR model prediction: input too short"); + assert(((numPredictions - 1) + fIdx(1)) >= -input.descriptor().start && + "array bounds error in AR model prediction: input too short"); + assert(((numPredictions - 1) + fIdx(1)) < input.size() && + "array bounds error in AR model prediction: input too short"); + assert(((numPredictions - 1) + fIdx(mParameters.size())) >= -input.descriptor().start && + "array bounds error in AR model prediction: input too short"); + assert(((numPredictions - 1) + fIdx(mParameters.size())) < input.size() && + "array bounds error in AR model prediction: input too short"); + const double* input_ptr = input.data(); for (index p = 0; p < numPredictions; p++) From 05a2fb834261c76afbce8a30695c1f8ae18343c9 Mon Sep 17 00:00:00 2001 From: Alex Harker Date: Tue, 1 Mar 2022 13:49:19 +0000 Subject: [PATCH 06/12] Correct indexing bugs due to incorrect refactor in robustIteration() --- include/algorithms/util/ARModel.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 76296b350..973f715cd 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -204,8 +204,8 @@ class ARModel // Iterate to find new filtered input for (index i = 0; i < size; i++) { - const double prediction = fowardPrediction(estimates); - estimates[0] = + const double prediction = fowardPrediction(estimates + i); + estimates[i] = prediction + robustResidual(input[i], prediction, robustFactor * sqrt(mVariance)); } From 100c747d2871f644d9a941aeae95811082376ac0 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 2 Mar 2022 09:25:05 +0000 Subject: [PATCH 07/12] ARModel: fix model stability by reverting mistaken array indexing change --- include/algorithms/util/ARModel.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 973f715cd..255578760 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -163,7 +163,7 @@ class ARModel directEstimate(input, size, true); // Initialise Estimates - for (index i = mParameters.size(); i < mParameters.size() + size; i++) + for (index i = 0; i < mParameters.size() + size; i++) estimates[asUnsigned(i)] = input[i - mParameters.size()]; // Variance From e05933bef61d4ccda42a67d9a1357e64b0d35f7e Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 16 Mar 2022 15:36:37 +0000 Subject: [PATCH 08/12] TransientExtraction: safer indexing on windowing, with added asserts. More `FliudTensorView`, fewer pointers --- .../algorithms/public/TransientExtraction.hpp | 50 +++++++++++-------- 1 file changed, 28 insertions(+), 22 deletions(-) diff --git a/include/algorithms/public/TransientExtraction.hpp b/include/algorithms/public/TransientExtraction.hpp index 824a233ec..efe13fecd 100644 --- a/include/algorithms/public/TransientExtraction.hpp +++ b/include/algorithms/public/TransientExtraction.hpp @@ -82,6 +82,8 @@ 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 +93,8 @@ 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 +105,8 @@ 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; @@ -141,16 +146,12 @@ class TransientExtraction // Forward and backward error const double normFactor = 1.0 / sqrt(mModel.variance()); - assert(mInput.size() > modelOrder() + padSize() + - blockSize() + mDetectHalfWindow + 1 + modelOrder()); - auto inputView = FluidTensorView( - mInput.data(), modelOrder() + padSize(), - blockSize() + mDetectHalfWindow + 1 + modelOrder()); + mInput.data(), modelOrder() + padSize(), blockSize() + modelOrder()); auto fwdError = FluidTensorView( - mForwardError.data(), 0, blockSize() + mDetectHalfWindow + 1); + mForwardError.data(), 0, blockSize() + mDetectHalfWindow); auto backError = FluidTensorView( - mBackwardError.data(), 0, blockSize() + mDetectHalfWindow + 1); + mBackwardError.data(), 0, blockSize() + mDetectHalfWindow); errorCalculation<&ARModel::forwardErrorArray>(inputView, fwdError, normFactor); @@ -158,12 +159,15 @@ class TransientExtraction normFactor); // Window error functions (brute force convolution) - windowError(mForwardWindowedError.data(), - mForwardError.data() + modelOrder(), hopSize()); - windowError(mBackwardWindowedError.data(), - mBackwardError.data() + modelOrder(), hopSize()); - - // Detection + auto fwdWindowedError = FluidTensorView( + mForwardWindowedError.data(), 0, hopSize()); + auto backWindowedError = FluidTensorView( + mBackwardWindowedError.data(), 0, hopSize()); + + windowError(fwdError(Slice(modelOrder())),fwdWindowedError); + windowError(backError(Slice(modelOrder())),backWindowedError); + + // Detection index count = 0; const double hiThresh = mDetectThreshHi; const double loThresh = mDetectThreshLo; @@ -369,30 +373,32 @@ class TransientExtraction // 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 error, FluidTensorView 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); From b832f8a0a62ce436edbd483aa5235bb721839bf0 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 16 Mar 2022 15:37:57 +0000 Subject: [PATCH 09/12] TransientExtraction: formatting --- .../algorithms/public/TransientExtraction.hpp | 56 +++++++++++-------- 1 file changed, 32 insertions(+), 24 deletions(-) diff --git a/include/algorithms/public/TransientExtraction.hpp b/include/algorithms/public/TransientExtraction.hpp index efe13fecd..7773db1ad 100644 --- a/include/algorithms/public/TransientExtraction.hpp +++ b/include/algorithms/public/TransientExtraction.hpp @@ -83,7 +83,8 @@ 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"); + assert(modelOrder() >= mDetectHalfWindow && + "model order needs to be >= half filter size"); frame(input, inSize); analyze(); detection(); @@ -94,7 +95,8 @@ class TransientExtraction RealVectorView residual) { assert(mInitialized && "Call init() before processing"); - assert(modelOrder() >= mDetectHalfWindow && "model order needs to be >= half filter size"); + assert(modelOrder() >= mDetectHalfWindow && + "model order needs to be >= half filter size"); index inSize = input.extent(0); frame(input.data(), inSize); analyze(); @@ -106,7 +108,8 @@ class TransientExtraction RealVectorView transients, RealVectorView residual) { assert(mInitialized && "Call init() before processing"); - assert(modelOrder() >= mDetectHalfWindow && "model order needs to be >= half filter size"); + 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; @@ -148,8 +151,8 @@ class TransientExtraction auto inputView = FluidTensorView( mInput.data(), modelOrder() + padSize(), blockSize() + modelOrder()); - auto fwdError = FluidTensorView( - mForwardError.data(), 0, blockSize() + mDetectHalfWindow); + auto fwdError = FluidTensorView(mForwardError.data(), 0, + blockSize() + mDetectHalfWindow); auto backError = FluidTensorView( mBackwardError.data(), 0, blockSize() + mDetectHalfWindow); @@ -159,15 +162,15 @@ class TransientExtraction normFactor); // Window error functions (brute force convolution) - auto fwdWindowedError = FluidTensorView( - mForwardWindowedError.data(), 0, hopSize()); - auto backWindowedError = FluidTensorView( - mBackwardWindowedError.data(), 0, hopSize()); - - windowError(fwdError(Slice(modelOrder())),fwdWindowedError); - windowError(backError(Slice(modelOrder())),backWindowedError); - - // Detection + auto fwdWindowedError = + FluidTensorView(mForwardWindowedError.data(), 0, hopSize()); + auto backWindowedError = + FluidTensorView(mBackwardWindowedError.data(), 0, hopSize()); + + windowError(fwdError(Slice(modelOrder())), fwdWindowedError); + windowError(backError(Slice(modelOrder())), backWindowedError); + + // Detection index count = 0; const double hiThresh = mDetectThreshHi; const double loThresh = mDetectThreshLo; @@ -179,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; @@ -299,7 +304,7 @@ class TransientExtraction residual[i] = input[i + order]; } - if (mRefine) refine(FluidTensorView(residual,0, size), Au, u); + if (mRefine) refine(FluidTensorView(residual, 0, size), Au, u); for (index i = 0; i < (size - order); i++) transients[i] = input[i + order] - residual[i]; @@ -309,14 +314,14 @@ class TransientExtraction mInput.data() + padSize() + order + order); } - void refine(FluidTensorView io, Eigen::MatrixXd& Au, Eigen::MatrixXd& ls) + void refine(FluidTensorView 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) @@ -363,7 +368,7 @@ class TransientExtraction void errorCalculation(FluidTensorView input, FluidTensorView error, double normFactor) { - (mModel.*Method)(input,error); + (mModel.*Method)(input, error); // Take absolutes and normalise for (index i = 0; i < error.size(); i++) @@ -373,11 +378,14 @@ class TransientExtraction // Triangle window double calcWindow(double norm) { return std::min(norm, 1.0 - norm); } - void windowError(FluidTensorView error, FluidTensorView errorWindowed) + void windowError(FluidTensorView error, + FluidTensorView errorWindowed) { - assert(error.descriptor().start >= mDetectHalfWindow && "insufficient offset for filter size"); - assert(error.size() > errorWindowed.size() - 1 + (mDetectHalfWindow) && "insufficient input for filter size"); - + 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; From b2d4d9f725f5cc6283fa4baba281cd79248a14bc Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 16 Mar 2022 15:41:25 +0000 Subject: [PATCH 10/12] ARModel: Formatting --- include/algorithms/util/ARModel.hpp | 58 ++++++++++++++++------------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index b665b7613..6b8607782 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -11,11 +11,11 @@ under the European Union’s Horizon 2020 research and innovation programme #pragma once #include "ConvolutionTools.hpp" +#include "FluidEigenMappings.hpp" #include "Toeplitz.hpp" #include "../public/WindowFuncs.hpp" #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" -#include "FluidEigenMappings.hpp" #include #include #include @@ -46,7 +46,7 @@ class ARModel else directEstimate(input, true); } - + double fowardPrediction(FluidTensorView input) { double prediction; @@ -66,23 +66,25 @@ class ARModel double forwardError(FluidTensorView input) { double error; - forwardErrorArray(input,FluidTensorView(&error,0,1)); + forwardErrorArray(input, FluidTensorView(&error, 0, 1)); return error; } double backwardError(FluidTensorView input) { double error; - backwardErrorArray(input,FluidTensorView(&error,0,1)); + backwardErrorArray(input, FluidTensorView(&error, 0, 1)); return error; } - void forwardErrorArray(FluidTensorView input,FluidTensorView errors) + void forwardErrorArray(FluidTensorView input, + FluidTensorView errors) { - modelPredict(input, errors, std::negate{},Error{}); + modelPredict(input, errors, std::negate{}, Error{}); } - void backwardErrorArray(FluidTensorView input, FluidTensorView errors) + void backwardErrorArray(FluidTensorView input, + FluidTensorView errors) { modelPredict(input, errors, Identity{}, Error{}); } @@ -90,28 +92,31 @@ class ARModel void setMinVariance(double variance) { mMinVariance = variance; } private: - struct Identity { index operator()(index a) { return a; } }; - + struct Predict { - double operator()(double, double estimate) {return estimate;} + double operator()(double, double estimate) { return estimate; } }; - - struct Error{ - double operator()(double input, double estimate) {return input - estimate;} + + struct Error + { + double operator()(double input, double estimate) + { + return input - estimate; + } }; - + template void modelPredict(FluidTensorView input, FluidTensorView output, Indexer fIdx, OutputFn fOut) { index numPredictions = output.size(); - + assert(fIdx(1) >= -input.descriptor().start && "array bounds error in AR model prediction: input too short"); assert(fIdx(1) < input.size() && @@ -124,29 +129,31 @@ class ARModel "array bounds error in AR model prediction: input too short"); assert(((numPredictions - 1) + fIdx(1)) < input.size() && "array bounds error in AR model prediction: input too short"); - assert(((numPredictions - 1) + fIdx(mParameters.size())) >= -input.descriptor().start && + assert(((numPredictions - 1) + fIdx(mParameters.size())) >= + -input.descriptor().start && "array bounds error in AR model prediction: input too short"); assert(((numPredictions - 1) + fIdx(mParameters.size())) < input.size() && "array bounds error in AR model prediction: input too short"); - + const double* input_ptr = input.data(); - + for (index p = 0; p < numPredictions; p++) { double estimate = 0; for (index i = 0; i < mParameters.size(); i++) - estimate += mParameters(i) * (input_ptr + p)[fIdx(i + 1)]; + estimate += mParameters(i) * (input_ptr + p)[fIdx(i + 1)]; output[p] = fOut(input_ptr[p], estimate); } } - void directEstimate(FluidTensorView input, bool updateVariance) + void directEstimate(FluidTensorView input, + bool updateVariance) { index size = input.size(); - + // copy input to a 32 byte aligned block (otherwise risk segfaults on Linux) VectorXd frame = _impl::asEigen(input); - + if (mUseWindow) { if (mWindow.size() != size) @@ -194,12 +201,13 @@ class ARModel double robustFactor) { FluidTensor estimates(input.size() + mParameters.size()); - + // Calculate an initial estimate of parameters directEstimate(input, true); - assert(input.descriptor().start >= mParameters.size()&&"too little offset into input data"); - + assert(input.descriptor().start >= mParameters.size() && + "too little offset into input data"); + // Initialise Estimates for (index i = 0; i < input.size() + mParameters.size(); i++) estimates[i] = input.data()[i - mParameters.size()]; From c00accb914935b42fc285c3754785ecfeba288d2 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 16 Mar 2022 16:06:14 +0000 Subject: [PATCH 11/12] TestTransientSlice: change to approval tests --- .../algorithms/public/TestTransientSlice.cpp | 40 +++++---------- ..._predictable_on_real_material.approved.txt | 50 +++++++++++++++++++ ..._material_with_heavy_settings.approved.txt | 22 ++++++++ ...dictable_on_sharp_sine_bursts.approved.txt | 8 +++ 4 files changed, 93 insertions(+), 27 deletions(-) create mode 100644 tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material.approved.txt create mode 100644 tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material_with_heavy_settings.approved.txt create mode 100644 tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_sharp_sine_bursts.approved.txt diff --git a/tests/algorithms/public/TestTransientSlice.cpp b/tests/algorithms/public/TestTransientSlice.cpp index b50f50104..55f9915e8 100644 --- a/tests/algorithms/public/TestTransientSlice.cpp +++ b/tests/algorithms/public/TestTransientSlice.cpp @@ -1,12 +1,16 @@ -#define CATCH_CONFIG_MAIN +#define APPROVALS_CATCH +#define CATCH_CONFIG_MAIN +// #include +#include #include "SlicerTestHarness.hpp" + #include -#include #include #include #include #include + #include #include #include @@ -24,6 +28,9 @@ using WindowSize = StrongType; using ClumpLength = StrongType; using MinSliceLength = StrongType; +auto directoryDisposer = + ApprovalTests::Approvals::useApprovalsSubdirectory("approval_tests"); + struct TestParams { Order order; @@ -107,20 +114,16 @@ TEST_CASE("TransientSlice is predictable on sharp sine bursts", "[TransientSlice][slicers]") { auto source = testsignals::sharpSines(); - auto expected = std::vector{1000, 22050, 33075}; - + auto params = TestParams{Order(20), BlockSize(256), Padding(128), Skew(0), ThreshFwd(2), ThreshBack(1.1), WindowSize(14), ClumpLength(25), MinSliceLength(1000)}; - auto matcher = Catch::Matchers::Approx(expected); - index margin = 8; - matcher.margin(margin); auto result = runTest(source, params); - CHECK_THAT(result, matcher); + ApprovalTests::Approvals::verifyAll("Detected",result); } @@ -128,48 +131,31 @@ TEST_CASE("TransientSlice is predictable on real material", "[TransientSlice][slicers]") { auto source = testsignals::monoEurorackSynth(); - auto expected = std::vector{ - 144, 19188, 34706, 47223, 49465, 58299, 68185, 86942, 105689, - 106751, 117438, 139521, 152879, 161525, 167573, 179045, 186295, 205049, - 223795, 248985, 250356, 256304, 263609, 280169, 297483, 306502, 310674, - 312505, 319114, 327659, 335217, 346778, 364673, 368356, 384718, 400937, - 431226, 433295, 434501, 435764, 439536, 441625, 444028, 445795, 452031, - 453392, 465467, 481514, 494518, 496119, 505754, 512477, 514270}; auto params = TestParams{Order(20), BlockSize(256), Padding(128), Skew(0), ThreshFwd(2), ThreshBack(1.1), WindowSize(14), ClumpLength(25), MinSliceLength(1000)}; - auto matcher = Catch::Matchers::Approx(expected); - index margin = 1; - matcher.margin(margin); auto result = runTest(source, params); - CHECK_THAT(result, matcher); + ApprovalTests::Approvals::verifyAll("Detected",result); } TEST_CASE("TransientSlice is predictable on real material with heavy settings", "[TransientSlice][slicers]") { auto source = testsignals::monoEurorackSynth(); - auto expected = std::vector{ - 140, 19182, 34704, 47217, 58297, 68182, 86941, 105688, 117356, - 122134, 139498, 150485, 161516, 167571, 179043, 186293, 205047, 220493}; auto params = TestParams{Order(200), BlockSize(2048), Padding(1024), Skew(1), ThreshFwd(3), ThreshBack(1), WindowSize(15), ClumpLength(30), MinSliceLength(4410)}; - auto matcher = Catch::Matchers::Approx(expected); - index margin = 1; - matcher.margin(margin); - auto result = runTest(source(Slice(0, 220500)), params); - CHECK_THAT(result, matcher); + ApprovalTests::Approvals::verifyAll("Detected",result); } diff --git a/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material.approved.txt b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material.approved.txt new file mode 100644 index 000000000..04615e8fb --- /dev/null +++ b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material.approved.txt @@ -0,0 +1,50 @@ +Detected + + +[0] = 146 +[1] = 34710 +[2] = 47223 +[3] = 58304 +[4] = 68190 +[5] = 86944 +[6] = 105695 +[7] = 106756 +[8] = 117800 +[9] = 122140 +[10] = 139529 +[11] = 152884 +[12] = 161527 +[13] = 167575 +[14] = 186297 +[15] = 205055 +[16] = 223797 +[17] = 227228 +[18] = 238074 +[19] = 248989 +[20] = 250361 +[21] = 256309 +[22] = 306846 +[23] = 310675 +[24] = 312508 +[25] = 319117 +[26] = 327660 +[27] = 335219 +[28] = 346780 +[29] = 364675 +[30] = 400940 +[31] = 429155 +[32] = 431273 +[33] = 433296 +[34] = 434505 +[35] = 439538 +[36] = 441627 +[37] = 452033 +[38] = 453397 +[39] = 465470 +[40] = 474529 +[41] = 481516 +[42] = 494519 +[43] = 496122 +[44] = 505753 +[45] = 514268 + diff --git a/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material_with_heavy_settings.approved.txt b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material_with_heavy_settings.approved.txt new file mode 100644 index 000000000..77c003331 --- /dev/null +++ b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_real_material_with_heavy_settings.approved.txt @@ -0,0 +1,22 @@ +Detected + + +[0] = 149 +[1] = 19183 +[2] = 34706 +[3] = 47222 +[4] = 58300 +[5] = 68183 +[6] = 86942 +[7] = 105690 +[8] = 117439 +[9] = 122135 +[10] = 139500 +[11] = 150486 +[12] = 161522 +[13] = 167572 +[14] = 179048 +[15] = 186294 +[16] = 205050 +[17] = 220493 + diff --git a/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_sharp_sine_bursts.approved.txt b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_sharp_sine_bursts.approved.txt new file mode 100644 index 000000000..b6568c5db --- /dev/null +++ b/tests/algorithms/public/approval_tests/TestTransientSlice.TransientSlice_is_predictable_on_sharp_sine_bursts.approved.txt @@ -0,0 +1,8 @@ +Detected + + +[0] = 995 +[1] = 11023 +[2] = 22045 +[3] = 33070 + From 8ed2985ea7c8a85f00a459f66af6af6fae182003 Mon Sep 17 00:00:00 2001 From: weefuzzy Date: Wed, 16 Mar 2022 16:10:43 +0000 Subject: [PATCH 12/12] Hope to see meanginful output on CI test failures --- .github/workflows/flucoma-core-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/flucoma-core-ci.yml b/.github/workflows/flucoma-core-ci.yml index d56dda8ef..8cfd964e8 100644 --- a/.github/workflows/flucoma-core-ci.yml +++ b/.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