Closed lavakin closed 1 year ago
Hi, have you looked e.g. at my implementation of the TAI calculation?
It uses the RcppThread
library which comes with its own progressbar and a parallel for loop:
Rcpp::List rcpp_tei_parallel(const arma::sp_mat& expression,
Rcpp::NumericVector ps,
int ncores = 1){
std::vector< std::string > psnames = ps.attr("names");
int n_col = expression.n_cols;
Rcpp::NumericVector sumx(n_col);
Rcpp::NumericVector teisum(n_col);
Rcpp::NumericVector tei(n_col);
RcppThread::ProgressBar bar(n_col, 1);
RcppThread::parallelFor(0, n_col, [&] (int j) {
for (size_t i = 0; i < expression.n_rows; i++) {
sumx[j] += expression(i, j);
teisum[j] += expression(i, j) * ps[i];
tei[j] = teisum[j]/sumx[j];
}
}, ncores);
return Rcpp::List::create(Rcpp::Named("sumx") = sumx,
Rcpp::Named("teisum") = teisum, Rcpp::Named("tei") = tei);
}
Of ocurse I do not know if it works with the RcppEigen
Best regards
Kristian
The RcppThread
also provide some benchmarks, but I have not looked deeply into it:
Dear Kristian,
Thank you for your question and comments! If I understand correctly, your function computes the TAI value solely from the expression dataset and phylostrata, with the implementation of parallelization and a progress bar.
Regarding my TAI implementation, the parallelization is implicitly handled by the Eigen library, and it runs so quickly that adding a progress bar would not be necessary.
However, thanks to your comment, I looked into the possibility of using rcppthread
to compute the permutations and display a progress bar. Unfortunately, the progress update within the for-loop is atomic and occurs after each iteration. Since generating a single permutation is fast, the threads would mostly wait for each other to update the progress counter, resulting in a loss of the benefits of parallelization. In my current, implementation each thread updates the progress bar after 200 iterations, and it might even be more beneficial to increase this number further.
While it would be technically possible to modify my code so that each iteration computes multiple permutations and then use the progress bar you suggested, such a solution would sacrifice elegance.
Once again, thank you for your comments and for bringing the rcppthread
library to my attention.
Best, Nikola
Dear @lavakin
Thank you so much for finding this elegant solution for the Eigen
library! This will be a super useful addition to the package.
@kullrich
Thank you also for the excellent discussion regarding rcppthread
. If you do find an additional way to use this library with Eigen or even speed up some of the Eigen code for faster permutations then your suggestions would be greatly appreciated!
With many thanks and very best wishes, Hajk
Just a small comment: generating permutations is independent from Eigen (except for adding the generated permutations to the eigen matrix). Eigen is a library optimized for linear algebraic operations, not a general numerical library such as numpy.
Nikola
Hi @lavakin,
the intention of the rcpp_tei_parallel
function is to work with sparse data from single-cell expression.
Sometimes this data has million of cells, and the problem with the TAI.R
function is the following line of code:
TAIProfile <- cpp_TAI(as.matrix(ExpressionMatrix), as.vector(Phylostratum))
So the ExpressionMatrix will be converted into a full matrix before the cpp_TAI
function takes over, which normally would crash any personal machine with low memory profiles.
This is so far handled with the TEI.R
function, which split and distribute the splitted sparseMatrix.
Best regards
Kristian
Dear both,
That's an excellent point @kullrich, then let's also adjust cpp_TAI()
accordingly, since at the time when I wrote it, there was no single-cell data.
@kullrich Would you like to give it a go or does anything speak against it @lavakin ?
With many thanks and very best wishes, Hajk
Dear all,
Thank you for the discussions and the work that has been put in to speed up (my favourite) package. I have installed the specified dependencies (brew install llvm libomp
) on my Mac and loaded the new update to myTAI with devtools::load_all()
.
Unfortunately, I was unable to load myTAI and was left with an error.
18 warnings generated.
clang++ -mmacosx-version-min=10.13 -std=gnu++11 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/usr/local/lib -o myTAI.so RcppExports.o code.o cpp11.o rcpp_funcs.o -lomp -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation
ld: library not found for -lomp
clang: error: linker command failed with exit code 1 (use -v to see invocation)
make: *** [myTAI.so] Error 1
ERROR: compilation failed for package βmyTAIβ
β removing β/private/var/folders/cx/w3nwmg5x1yl4l9pr2kv3mrwh0000gn/T/RtmpMtcQ8p/devtools_install_c4768643e2b/myTAIβ
Error in `(function (command = NULL, args = character(), error_on_status = TRUE, β¦`:
! System command 'R' failed
---
Exit status: 1
stdout & stderr: <printed>
---
Type .Last.error to see the more details.
This was not the case in my last pull request.
With the assumption that users (desiring the developer version of myTAI
) will refer solely to the instructions in README.md
, this issue should be resolved as soon as possible.
Do you have any thoughts on how this can be resolved? On my part, I recommend reverting back to the commit before these changes have been implemented, e.g. https://github.com/drostlab/myTAI/commit/f671be8, if the fix will take some time, lest other users be inconvenienced.
Many thanks in advance.
Best, Sodai
Hi, I resolved with this:
brew install lomp
cd /usr/local/lib
ln -s /usr/local/opt/libomp/lib/libomp.dylib ./libomp.dylib
Best regards Kristian
Interesting, I tried it in a separated environment and it compiled without problems.
please let me know if it works with Kristians fix and if yes, Iβll add it to the README.
If it would be an issue in the future, we can also think about building docker, which would ensure same behaviour through distributions.
Thanks a lot, Kristian for fixing the issue.
Done! I added this to the README
. Thank you so much, @kullrich!
@LotharukpongJS Does this work for you now?
I agree that building a docker
would be great! But we need to also add this to the conda
recepe.
Many thanks, Hajk
Dear both,
That's an excellent point @kullrich, then let's also adjust
cpp_TAI()
accordingly, since at the time when I wrote it, there was no single-cell data.@kullrich Would you like to give it a go or does anything speak against it @lavakin ?
With many thanks and very best wishes, Hajk
Eigen has its own SparseMatrix class, so we can also have a look if there is an easy way to switch to that. (The documentation webpage doesnt work rn, but will have a look once its working again)
Dear @lavakin
I also get:
Plot signature: ' TAI ' and test statistic: ' FlatLineTest ' running 20000 permutations.
Computing permutations
[=========================================] 100%
Computing variances
Number of Eigen threads: 8
Computing permutations
[=========================================] 100%
Computing variances
Number of Eigen threads: 8
Time: 12.3361010551453
Significance status of signature: significant.
Now run 'FlatLineTest(..., permutations = 20000 , plotHistogram = TRUE)' to analyse the permutation test performance
Warning message:
In ks.test.default(filtered_vars, "pgamma", shape = shape, rate = rate) :
ties should not be present for the Kolmogorov-Smirnov test
Would it be possible to find a quick solution to this warning?
Warning message:
In ks.test.default(filtered_vars, "pgamma", shape = shape, rate = rate) :
ties should not be present for the Kolmogorov-Smirnov test
With many thanks, Hajk
I looked into it already at some point and as far as I understand it, it's a property of myTAI. Since we have a lot of repeated vaules in the phylostrata (because there are not many possible values), with a higher number of permutations, we will unavoidably get the same variance multiple times, because the same phylostrata will be assigned to the same genes (eventhough the permutations are different).
This might be also one of the reasons for a bad fit. There is a way to prevent this, but then the permutations would not be randomly drawn anymore, so I don't think that would be a good idea.
But if you have a suggestion, what to do about it, I will be more than happy to fix it :).
Hi @lavakin,
regarding the sparseMatrix
from the RcppEigen
package.
It would be great if one could just use it directly.
I created a small benchmark with the Matrix
package and some matrix crossprod
calculation, since given the PhyloExpressionData
the dimensions do not change.
Here, if the RcppEigen
does not directyl work with the sparseMatrix
one could also write a R
wrapper function to decide wether to use cpp_TAI
or e.g. r_TAI
.
Matrix
packagelibrary(myTAI)
library(Matrix)
r_TAI <- function(eM, Phylostratum){
Divisor <- colSums(eM)
dM <- Matrix::Matrix(0,
ncol=length(Divisor),
nrow=length(Divisor),
sparse=TRUE)
diag(dM) <- 1/Divisor
psM <- Matrix::Matrix(0,
ncol=length(Phylostratum),
nrow=length(Phylostratum),
sparse=TRUE)
diag(psM) <- Phylostratum
# crossprod magic
total <- colSums(t(t(eM) %*% psM) %*% dM)
return(total)
}
And here an example creating sparse data:
# load example data
data(PhyloExpressionSetExample)
ExpressionMatrix <- PhyloExpressionSetExample[,-c(1,2)]
Phylostratum <- PhyloExpressionSetExample[,1]
# create sparse example
ExpressionMatrix[ExpressionMatrix>800]<-0
# create more expression data
test <- cbind(PhyloExpressionSetExample[,c(1,2)],
do.call(cbind, rep(ExpressionMatrix, 2000)))
colnames(test)[-c(1,2)] <- seq(from=1, to=dim(test)[2]-2)
dim(test)
ExpressionMatrix <- as.matrix(test[,-c(1,2)])
ExpressionMatrix_sparse <- Matrix::Matrix(ExpressionMatrix, sparse=TRUE)
str(ExpressionMatrix_sparse)
and the benchmark:
# benchmark
system.time(r_TAI_results <- r_TAI(ExpressionMatrix, Phylostratum))
User System verstrichen
26.203 22.533 60.688
system.time(r_TAI_results_sparse <- r_TAI(ExpressionMatrix_sparse, Phylostratum))
User System verstrichen
5.945 3.670 11.641
system.time(cpp_TAI_results <- myTAI:::cpp_TAI(ExpressionMatrix, Phylostratum))
User System verstrichen
3.011 11.372 19.785
Note: All three different function have slightly different results due to rounding in the used packages and how R handles matrix
versus sparseMatrix
colSums
I guess.
as.character(r_TAI_results[1])
[1] "3.77721412365141"
as.character(r_TAI_results_sparse[1])
[1] "3.77721412365141"
as.character(cpp_TAI_results[1])
[1] "3.77721412365142"
Best regards
Kristian
Hi all,
@kullrich thank you for the solution.
@HajkD regarding the docker, I am not sure if it would make it complicated for some users (biologists) who just want to run myTAI to analyse their new RNA-seq data? One idea is to proceed with the updated multithreading only if OpenMP
is detected in the system. Let me know what you think.
I notice some weird behaviours too. I see %
signs in places that shouldn't exist.
> data("PhyloExpressionSetExample")
> ExprExample <- tf(PhyloExpressionSetExample, log2)
> PlotSignature(ExprExample, permutations = 20000)
Plot signature: ' TAI ' and test statistic: ' FlatLineTest ' running 20000 permutations.
Computing permutations
[=========================================] 100%
Computing variances
Number of Eigen threads: 8
Computing permutations
[=========================================] 100% 60% % % % ][=======================> ]
Computing variances
Number of Eigen threads: 8
Time: 1.83384306828181
Significance status of signature: significant.
Now run 'FlatLineTest(..., permutations = 20000 , plotHistogram = TRUE)' to analyse the permutation test performanceWarning message:
In ks.test.default(filtered_vars, "pgamma", shape = shape, rate = rate) :
ties should not be present for the Kolmogorov-Smirnov test
Furthermore, myTAI
doesn't pass devtools::test()
and devtools::check()
. @lavakin , do you mind having a look?
devtools::test()
[ FAIL 0 | WARN 3480 | SKIP 3 | PASS 200 ]
devtools::check()
ββ R CMD check results βββββββββββββββββββββ myTAI 1.0.1.9000 ββββ
Duration: 4m 59.6s
β― checking compilation flags in Makevars ... WARNING
Non-portable flags in variable 'PKG_CXXFLAGS':
-fopenmp
β― checking for GNU extensions in Makefiles ... WARNING
Found the following file(s) containing GNU extensions:
src/Makevars
Portable Makefiles do not use GNU extensions such as +=, :=, $(shell),
$(wildcard), ifeq ... endif, .NOTPARALLEL See section βWriting portable
packagesβ in the βWriting R Extensionsβ manual.
β― checking package dependencies ... NOTE
Packages suggested but not available for checking: 'taxize', 'Seurat'
β― checking installed package size ... NOTE
installed size is 5.6Mb
sub-directories of 1Mb or more:
data 2.0Mb
doc 1.2Mb
help 1.6Mb
β― checking top-level files ... NOTE
Non-standard file/directory found at top level:
βpkgdownβ
β― checking compiled code ... NOTE
File βmyTAI/libs/myTAI.soβ:
Found β__ZNSt3__14coutEβ, possibly from βstd::coutβ (C++)
Object: βrcpp_funcs.oβ
Compiled code should not call entry points which might terminate R nor
write to stdout/stderr instead of to the console, nor use Fortran I/O
nor system RNGs.
See βWriting portable packagesβ in the βWriting R Extensionsβ manual.
0 errors β | 2 warnings β | 4 notes β
Error: R CMD check found WARNINGs
Execution halted
Exited with status 1.
Hi @LotharukpongJS,
I very much like the idea that the openmpi code is only used upon internal reference check whether or not a user has OpenMP
installed or not. This way, less computer oriented users do not have to bother installing the OpenMP libraries.
@kullrich Thank you very much for this great suggestion! Let's also build this. I like it a lot!
With many thanks and very best wishes, Hajk
Configured cpp to use
OpenMP
, so that Eigen can run in parallel and the permutation function can be parallelized usingomp
library. This might require additional software. Instructions have been added to README.Created a progress bar to show the progress during computing the permutations. Seems like I managed to handle all the pitfalls of multiprocessing and it didn't slow it down π₯³. (Not the nicest implementation but still good enough I hope.)
Changed the default number of permutations to 10000 so that people will (be forced to) admire my progress bar for longer. (And the fit is more accurate ofc)