-
Notifications
You must be signed in to change notification settings - Fork 23
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
fix massive bottleneck #296
Conversation
Codecov ReportAll modified and coverable lines are covered by tests ✅
Additional details and impacted files@@ Coverage Diff @@
## main #296 +/- ##
==========================================
- Coverage 88.71% 88.71% -0.01%
==========================================
Files 112 112
Lines 9122 9115 -7
==========================================
- Hits 8093 8086 -7
Misses 1029 1029
☔ View full report in Codecov by Sentry. |
Any idea why |
The R CMD Check (macOS-latest (release)) job is still "expected", but... hanging? |
src/Distributions.cpp
Outdated
@@ -415,7 +415,7 @@ NumericMatrix C_Vec_WeightedDiscretePdf(NumericVector x, NumericMatrix data, | |||
for (int i = 0; i < nc; i++) { | |||
for (int k = 0; k < n; k++) { | |||
for (int j = 0; j < nr; j++) { | |||
if (data(j, i) == x[k]) { | |||
if (data(j) == x[k]) { |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
C++ not my strong point, but since data
is a NumericVector
as well, shouln't it be data[j]
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Checking with https://teuder.github.io/rcpp4everyone_en/080_vector.html#accessing-vector-elements that should be correct. I wonder what actually happens when using data(j, i)
in this case 🤔
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
still it would be better to have it consistent, since both x
and data
are of the same type
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Verification question: if we never have the statement data(j) == x[k]
TRUE, then the mat(k, i)
(pdf at a new time point not included in the original x
) will be zero and that's what we want, right?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yup
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
- Would it make sense to speed this up by removing the third loop? You pretty much search for the index where the time/data points match every time you scan the vector
data
(again C++ not my strongest point, maybe with thebreak
and scanning withfor
you are already doing pretty good) - There are other places in the code where we can substitute () with [] to make it more consistent - pushed a commit for this
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
How would you remove the third loop? Can you demonstrate with R code?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Was writing the previous response late and didn't think it through, so what I initially thought is not going to work. But I asked ChatGTP and got something similar: an optimization that reduces the time complexity from O(nr * nc * n)
to O(nr + nc * n)
by precomputing the indices of elements in the 'data' array using a hashmap (!):
#include <unordered_map>
// Assuming mat, pdf, data, and x are appropriately defined
// Create a hashmap to store the indices of elements in the 'data' array
std::unordered_map<int, int> dataIndices;
for (int j = 0; j < nr; ++j) {
dataIndices[data[j]] = j;
}
// Use the hashmap to populate the 'mat' matrix efficiently
for (int i = 0; i < nc; ++i) {
for (int k = 0; k < n; ++k) {
auto dataIndexIterator = dataIndices.find(x[k]);
if (dataIndexIterator != dataIndices.end()) {
int j = dataIndexIterator->second;
mat(k, i) = pdf(j, i);
}
}
}
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Heh, looks like your ChatGPT-foo is better than mine, I didn't get too far 😬
I thought that using an ordered set should be the way to go based on how indices are, well, ordered, so surely using an ordered data structure has some sort of benefit?
But I also didn't think this through entirely because I think I still haven't wrapped my head around this fully.
@@ -310,30 +310,28 @@ Arrdist <- R6Class("Arrdist", | |||
.pdf = function(x, log = FALSE) { | |||
"pdf, data, wc" %=% gprm(self, c("pdf", "x", "which.curve")) | |||
mat <- .extCurve(pdf, wc) | |||
out <- t(C_Vec_WeightedDiscretePdf( | |||
x, matrix(data, ncol(mat), private$.ndists), t(mat))) | |||
out <- t(C_Vec_WeightedDiscretePdf(x, data, t(mat))) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@RaphaelS1 so I think I understand what the C_Vec_WeightedDiscretePdf
does => gets a pdf
matrix and subsets it to the x
times (points) instead of the data
it has. It's just the names of x
and data
are not so informative in this context. Some points for discussion:
- Should we consider some name changing? I think
data
is like axnew
? - These matrix tranposes maybe slow down things? Wouldn't be more sense if the code worked as
C_Vec_WeightedDiscretePdf(x, data, mat))
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes maybe but not in this PR
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think (1) would make sense to do now, (2) might be trickier
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
for (1) => #297
Because I have outdated branch protection settings |
Fixes a huge bottleneck in cpp functions for vectorised WeightDisc (includes Matdist and Arrdist)