scran_pca
Principal component analysis for single-cell data
Loading...
Searching...
No Matches
blocked_pca.hpp
Go to the documentation of this file.
1#ifndef SCRAN_PCA_BLOCKED_PCA_HPP
2#define SCRAN_PCA_BLOCKED_PCA_HPP
3
4#include <vector>
5#include <cmath>
6#include <algorithm>
7#include <type_traits>
8#include <cstddef>
9
10#include "tatami/tatami.hpp"
11#include "irlba/irlba.hpp"
12#include "irlba/parallel.hpp"
13#include "Eigen/Dense"
15#include "sanisizer/sanisizer.hpp"
16
17#include "utils.hpp"
18
25namespace scran_pca {
26
35 // Avoid throwing an error if too many PCs are requested.
37 }
47 int number = 25;
48
55 bool scale = false;
56
61 bool transpose = true;
62
72 scran_blocks::WeightPolicy block_weight_policy = scran_blocks::WeightPolicy::VARIABLE;
73
79
86
91 bool realize_matrix = true;
92
97 int num_threads = 1;
98
103};
104
108namespace internal {
109
110/*****************************************************
111 ************* Blocking data structures **************
112 *****************************************************/
113
114template<typename Index_, class EigenVector_>
115struct BlockingDetails {
116 std::vector<Index_> block_size;
117
118 bool weighted = false;
119 typedef typename EigenVector_::Scalar Weight;
120
121 // The below should only be used if weighted = true.
122 std::vector<Weight> per_element_weight;
123 Weight total_block_weight = 0;
124 EigenVector_ expanded_weights;
125};
126
127template<class EigenVector_, typename Index_, typename Block_>
128BlockingDetails<Index_, EigenVector_> compute_blocking_details(
129 const Index_ ncells,
130 const Block_* block,
131 const scran_blocks::WeightPolicy block_weight_policy,
132 const scran_blocks::VariableWeightParameters& variable_block_weight_parameters)
133{
134 BlockingDetails<Index_, EigenVector_> output;
135 output.block_size = tatami_stats::tabulate_groups(block, ncells);
136 if (block_weight_policy == scran_blocks::WeightPolicy::NONE) {
137 return output;
138 }
139
140 const auto& block_size = output.block_size;
141 const auto nblocks = block_size.size();
142 output.weighted = true;
143 auto& total_weight = output.total_block_weight;
144 auto& element_weight = output.per_element_weight;
145 sanisizer::resize(element_weight, nblocks);
146
147 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
148 const auto bsize = block_size[b];
149
150 // Computing effective block weights that also incorporate division by the
151 // block size. This avoids having to do the division by block size in the
152 // 'compute_blockwise_mean_and_variance*()' functions.
153 if (bsize) {
154 typename EigenVector_::Scalar block_weight = 1;
155 if (block_weight_policy == scran_blocks::WeightPolicy::VARIABLE) {
156 block_weight = scran_blocks::compute_variable_weight(bsize, variable_block_weight_parameters);
157 }
158
159 element_weight[b] = block_weight / bsize;
160 total_weight += block_weight;
161 } else {
162 element_weight[b] = 0;
163 }
164 }
165
166 // Setting a placeholder value to avoid problems with division by zero.
167 if (total_weight == 0) {
168 total_weight = 1;
169 }
170
171 // Expanding them for multiplication in the IRLBA wrappers.
172 auto sqrt_weights = element_weight;
173 for (auto& s : sqrt_weights) {
174 s = std::sqrt(s);
175 }
176
177 auto& expanded = output.expanded_weights;
178 sanisizer::resize(expanded, ncells);
179 for (Index_ c = 0; c < ncells; ++c) {
180 expanded.coeffRef(c) = sqrt_weights[block[c]];
181 }
182
183 return output;
184}
185
186/*****************************************************************
187 ************ Computing the blockwise mean and variance **********
188 *****************************************************************/
189
190template<typename Num_, typename Value_, typename Index_, typename Block_, typename EigenVector_, typename Float_>
191void compute_sparse_mean_and_variance_blocked(
192 const Num_ num_nonzero,
193 const Value_* values,
194 const Index_* indices,
195 const Block_* block,
196 const BlockingDetails<Index_, EigenVector_>& block_details,
197 Float_* centers,
198 Float_& variance,
199 std::vector<Index_>& block_copy,
200 const Num_ num_all)
201{
202 const auto& block_size = block_details.block_size;
203 const auto nblocks = block_size.size();
204
205 std::fill_n(centers, nblocks, 0);
206 for (Num_ i = 0; i < num_nonzero; ++i) {
207 centers[block[indices[i]]] += values[i];
208 }
209 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
210 auto bsize = block_size[b];
211 if (bsize) {
212 centers[b] /= bsize;
213 }
214 }
215
216 // Computing the variance from the sum of squared differences.
217 // This is technically not the correct variance estimate if we
218 // were to consider the loss of residual d.f. from estimating
219 // the block means, but it's what the PCA sees, so whatever.
220 variance = 0;
221 std::copy(block_size.begin(), block_size.end(), block_copy.begin());
222
223 if (block_details.weighted) {
224 for (Num_ i = 0; i < num_nonzero; ++i) {
225 const Block_ curb = block[indices[i]];
226 const auto diff = values[i] - centers[curb];
227 variance += diff * diff * block_details.per_element_weight[curb];
228 --block_copy[curb];
229 }
230 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
231 const auto val = centers[b];
232 variance += val * val * block_copy[b] * block_details.per_element_weight[b];
233 }
234 } else {
235 for (Num_ i = 0; i < num_nonzero; ++i) {
236 const Block_ curb = block[indices[i]];
237 const auto diff = values[i] - centers[curb];
238 variance += diff * diff;
239 --block_copy[curb];
240 }
241 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
242 const auto val = centers[b];
243 variance += val * val * block_copy[b];
244 }
245 }
246
247 // COMMENT ON DENOMINATOR:
248 // If we're not dealing with weights, we compute the actual sample
249 // variance for easy interpretation (and to match up with the
250 // per-PC calculations in internal::clean_up).
251 //
252 // If we're dealing with weights, the concept of the sample variance
253 // becomes somewhat weird, but we just use the same denominator for
254 // consistency in clean_up_projected. Magnitude doesn't matter when
255 // scaling for internal::process_scale_vector anyway.
256 variance /= num_all - 1;
257}
258
259template<class IrlbaSparseMatrix_, typename Block_, class Index_, class EigenVector_, class EigenMatrix_>
260void compute_blockwise_mean_and_variance_realized_sparse(
261 const IrlbaSparseMatrix_& emat, // this should be column-major with genes in the columns.
262 const Block_* block,
263 const BlockingDetails<Index_, EigenVector_>& block_details,
264 EigenMatrix_& centers,
265 EigenVector_& variances,
266 const int nthreads)
267{
268 const auto ngenes = emat.cols();
269 tatami::parallelize([&](const int, const decltype(I(ngenes)) start, const decltype(I(ngenes)) length) -> void {
270 const auto ncells = emat.rows();
271 const auto& values = emat.get_values();
272 const auto& indices = emat.get_indices();
273 const auto& pointers = emat.get_pointers();
274
275 const auto nblocks = block_details.block_size.size();
276 static_assert(!EigenMatrix_::IsRowMajor);
277 auto block_copy = sanisizer::create<std::vector<Index_> >(nblocks);
278
279 for (decltype(I(start)) g = start, end = start + length; g < end; ++g) {
280 const auto offset = pointers[g];
281 const auto next_offset = pointers[g + 1]; // increment won't overflow as 'g < end' and 'end' is of the same type.
282 compute_sparse_mean_and_variance_blocked(
283 static_cast<decltype(I(ncells))>(next_offset - offset),
284 values.data() + offset,
285 indices.data() + offset,
286 block,
287 block_details,
288 centers.data() + sanisizer::product_unsafe<std::size_t>(g, nblocks),
289 variances[g],
290 block_copy,
291 ncells
292 );
293 }
294 }, ngenes, nthreads);
295}
296
297template<typename Num_, typename Value_, typename Block_, typename Index_, typename EigenVector_, typename Float_>
298void compute_dense_mean_and_variance_blocked(
299 const Num_ number,
300 const Value_* values,
301 const Block_* block,
302 const BlockingDetails<Index_, EigenVector_>& block_details,
303 Float_* centers,
304 Float_& variance)
305{
306 const auto& block_size = block_details.block_size;
307 const auto nblocks = block_size.size();
308 std::fill_n(centers, nblocks, 0);
309 for (Num_ i = 0; i < number; ++i) {
310 centers[block[i]] += values[i];
311 }
312 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
313 const auto& bsize = block_size[b];
314 if (bsize) {
315 centers[b] /= bsize;
316 }
317 }
318
319 variance = 0;
320
321 if (block_details.weighted) {
322 for (Num_ i = 0; i < number; ++i) {
323 const auto curb = block[i];
324 const auto delta = values[i] - centers[curb];
325 variance += delta * delta * block_details.per_element_weight[curb];
326 }
327 } else {
328 for (Num_ i = 0; i < number; ++i) {
329 const auto curb = block[i];
330 const auto delta = values[i] - centers[curb];
331 variance += delta * delta;
332 }
333 }
334
335 variance /= number - 1; // See COMMENT ON DENOMINATOR above.
336}
337
338template<class EigenMatrix_, typename Block_, class Index_, class EigenVector_>
339void compute_blockwise_mean_and_variance_realized_dense(
340 const EigenMatrix_& emat, // this should be column-major with genes in the columns.
341 const Block_* block,
342 const BlockingDetails<Index_, EigenVector_>& block_details,
343 EigenMatrix_& centers,
344 EigenVector_& variances,
345 const int nthreads)
346{
347 const auto ngenes = emat.cols();
348 tatami::parallelize([&](const int, const decltype(I(ngenes)) start, const decltype(I(ngenes)) length) -> void {
349 const auto ncells = emat.rows();
350 static_assert(!EigenMatrix_::IsRowMajor);
351 const auto nblocks = block_details.block_size.size();
352 for (decltype(I(start)) g = start, end = start + length; g < end; ++g) {
353 compute_dense_mean_and_variance_blocked(
354 ncells,
355 emat.data() + sanisizer::product_unsafe<std::size_t>(g, ncells),
356 block,
357 block_details,
358 centers.data() + sanisizer::product_unsafe<std::size_t>(g, nblocks),
359 variances[g]
360 );
361 }
362 }, ngenes, nthreads);
363}
364
365template<typename Value_, typename Index_, typename Block_, class EigenMatrix_, class EigenVector_>
366void compute_blockwise_mean_and_variance_tatami(
367 const tatami::Matrix<Value_, Index_>& mat, // this should have genes in the rows!
368 const Block_* block,
369 const BlockingDetails<Index_, EigenVector_>& block_details,
370 EigenMatrix_& centers,
371 EigenVector_& variances,
372 const int nthreads)
373{
374 const auto& block_size = block_details.block_size;
375 const auto nblocks = block_size.size();
376 const Index_ ngenes = mat.nrow();
377 const Index_ ncells = mat.ncol();
378
379 if (mat.prefer_rows()) {
380 tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
381 static_assert(!EigenMatrix_::IsRowMajor);
382 auto block_copy = sanisizer::create<std::vector<Index_> >(nblocks);
384
385 if (mat.is_sparse()) {
387 auto ext = tatami::consecutive_extractor<true>(mat, true, start, length);
388 for (Index_ g = start, end = start + length; g < end; ++g) {
389 auto range = ext->fetch(vbuffer.data(), ibuffer.data());
390 compute_sparse_mean_and_variance_blocked(
391 range.number,
392 range.value,
393 range.index,
394 block,
395 block_details,
396 centers.data() + sanisizer::product_unsafe<std::size_t>(g, nblocks),
397 variances[g],
398 block_copy,
399 ncells
400 );
401 }
402 } else {
403 auto ext = tatami::consecutive_extractor<false>(mat, true, start, length);
404 for (Index_ g = start, end = start + length; g < end; ++g) {
405 auto ptr = ext->fetch(vbuffer.data());
406 compute_dense_mean_and_variance_blocked(
407 ncells,
408 ptr,
409 block,
410 block_details,
411 centers.data() + sanisizer::product_unsafe<std::size_t>(g, nblocks),
412 variances[g]
413 );
414 }
415 }
416 }, ngenes, nthreads);
417
418 } else {
419 typedef typename EigenVector_::Scalar Scalar;
420 std::vector<std::pair<decltype(I(nblocks)), Scalar> > block_multipliers;
421 block_multipliers.reserve(nblocks);
422
423 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
424 const auto bsize = block_size[b];
425 if (bsize > 1) { // skipping blocks with NaN variances.
426 Scalar mult = bsize - 1; // need to convert variances back into sum of squared differences.
427 if (block_details.weighted) {
428 mult *= block_details.per_element_weight[b];
429 }
430 block_multipliers.emplace_back(b, mult);
431 }
432 }
433
434 tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
435 std::vector<std::vector<Scalar> > re_centers, re_variances;
436 re_centers.reserve(nblocks);
437 re_variances.reserve(nblocks);
438 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
439 re_centers.emplace_back(length);
440 re_variances.emplace_back(length);
441 }
442
444
445 if (mat.is_sparse()) {
446 std::vector<tatami_stats::variances::RunningSparse<Scalar, Value_, Index_> > running;
447 running.reserve(nblocks);
448 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
449 running.emplace_back(length, re_centers[b].data(), re_variances[b].data(), /* skip_nan = */ false, /* subtract = */ start);
450 }
451
453 auto ext = tatami::consecutive_extractor<true>(mat, false, static_cast<Index_>(0), ncells, start, length);
454 for (Index_ c = 0; c < ncells; ++c) {
455 const auto range = ext->fetch(vbuffer.data(), ibuffer.data());
456 running[block[c]].add(range.value, range.index, range.number);
457 }
458
459 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
460 running[b].finish();
461 }
462
463 } else {
464 std::vector<tatami_stats::variances::RunningDense<Scalar, Value_, Index_> > running;
465 running.reserve(nblocks);
466 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
467 running.emplace_back(length, re_centers[b].data(), re_variances[b].data(), /* skip_nan = */ false);
468 }
469
470 auto ext = tatami::consecutive_extractor<false>(mat, false, static_cast<Index_>(0), ncells, start, length);
471 for (Index_ c = 0; c < ncells; ++c) {
472 auto ptr = ext->fetch(vbuffer.data());
473 running[block[c]].add(ptr);
474 }
475
476 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
477 running[b].finish();
478 }
479 }
480
481 static_assert(!EigenMatrix_::IsRowMajor);
482 for (Index_ i = 0; i < length; ++i) {
483 auto mptr = centers.data() + sanisizer::product_unsafe<std::size_t>(start + i, nblocks);
484 for (decltype(I(nblocks)) b = 0; b < nblocks; ++b) {
485 mptr[b] = re_centers[b][i];
486 }
487
488 auto& my_var = variances[start + i];
489 my_var = 0;
490 for (const auto& bm : block_multipliers) {
491 my_var += re_variances[bm.first][i] * bm.second;
492 }
493 my_var /= ncells - 1; // See COMMENT ON DENOMINATOR above.
494 }
495 }, ngenes, nthreads);
496 }
497}
498
499/******************************************************************
500 ************ Project matrices on their rotation vectors **********
501 ******************************************************************/
502
503template<class EigenMatrix_, class EigenVector_>
504const EigenMatrix_& scale_rotation_matrix(const EigenMatrix_& rotation, bool scale, const EigenVector_& scale_v, EigenMatrix_& tmp) {
505 if (scale) {
506 tmp = (rotation.array().colwise() / scale_v.array()).matrix();
507 return tmp;
508 } else {
509 return rotation;
510 }
511}
512
513template<class IrlbaSparseMatrix_, class EigenMatrix_>
514inline void project_matrix_realized_sparse(
515 const IrlbaSparseMatrix_& emat, // cell in rows, genes in the columns, CSC.
516 EigenMatrix_& components, // dims in rows, cells in columns
517 const EigenMatrix_& scaled_rotation, // genes in rows, dims in columns
518 int nthreads)
519{
520 const auto rank = scaled_rotation.cols();
521 const auto ncells = emat.rows();
522 const auto ngenes = emat.cols();
523
524 // Store as transposed for more cache efficiency.
525 components.resize(
526 sanisizer::cast<decltype(I(components.rows()))>(rank),
527 sanisizer::cast<decltype(I(components.cols()))>(ncells)
528 );
529 components.setZero();
530
531 const auto& values = emat.get_values();
532 const auto& indices = emat.get_indices();
533
534 if (nthreads == 1) {
535 const auto& pointers = emat.get_pointers();
536 auto multipliers = sanisizer::create<Eigen::VectorXd>(rank);
537 for (decltype(I(ngenes)) g = 0; g < ngenes; ++g) {
538 multipliers.noalias() = scaled_rotation.row(g);
539 const auto start = pointers[g], end = pointers[g + 1]; // increment is safe as 'g + 1 <= ngenes'.
540 for (auto i = start; i < end; ++i) {
541 components.col(indices[i]).noalias() += values[i] * multipliers;
542 }
543 }
544
545 } else {
546 const auto& row_nonzero_starts = emat.get_secondary_nonzero_starts();
547 irlba::parallelize(nthreads, [&](const int t) -> void {
548 const auto& starts = row_nonzero_starts[t];
549 const auto& ends = row_nonzero_starts[t + 1]; // increment is safe as 't + 1 <= nthreads'.
550 auto multipliers = sanisizer::create<Eigen::VectorXd>(rank);
551
552 for (decltype(I(ngenes)) g = 0; g < ngenes; ++g) {
553 multipliers.noalias() = scaled_rotation.row(g);
554 const auto start = starts[g], end = ends[g];
555 for (auto i = start; i < end; ++i) {
556 components.col(indices[i]).noalias() += values[i] * multipliers;
557 }
558 }
559 });
560 }
561}
562
563template<typename Value_, typename Index_, class EigenMatrix_>
564void project_matrix_transposed_tatami(
565 const tatami::Matrix<Value_, Index_>& mat, // genes in rows, cells in columns
566 EigenMatrix_& components,
567 const EigenMatrix_& scaled_rotation, // genes in rows, dims in columns
568 const int nthreads)
569{
570 const auto rank = scaled_rotation.cols();
571 const auto ngenes = mat.nrow();
572 const auto ncells = mat.ncol();
573 typedef typename EigenMatrix_::Scalar Scalar;
574
575 // Store as transposed for more cache efficiency.
576 components.resize(
577 sanisizer::cast<decltype(I(components.rows()))>(rank),
578 sanisizer::cast<decltype(I(components.cols()))>(ncells)
579 );
580
581 if (mat.prefer_rows()) {
582 tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
583 static_assert(!EigenMatrix_::IsRowMajor);
584 const auto vptr = scaled_rotation.data();
586
587 std::vector<std::vector<Scalar> > local_buffers; // create separate buffers to avoid false sharing.
588 local_buffers.reserve(rank);
589 for (decltype(I(rank)) r = 0; r < rank; ++r) {
590 local_buffers.emplace_back(tatami::cast_Index_to_container_size<decltype(I(local_buffers.front()))>(length));
591 }
592
593 if (mat.is_sparse()) {
595 auto ext = tatami::consecutive_extractor<true>(mat, true, static_cast<Index_>(0), ngenes, start, length);
596 for (Index_ g = 0; g < ngenes; ++g) {
597 const auto range = ext->fetch(vbuffer.data(), ibuffer.data());
598 for (decltype(I(rank)) r = 0; r < rank; ++r) {
599 const auto mult = vptr[sanisizer::nd_offset<std::size_t>(g, ngenes, r)];
600 auto& local_buffer = local_buffers[r];
601 for (Index_ i = 0; i < range.number; ++i) {
602 local_buffer[range.index[i] - start] += range.value[i] * mult;
603 }
604 }
605 }
606
607 } else {
608 auto ext = tatami::consecutive_extractor<false>(mat, true, static_cast<Index_>(0), ngenes, start, length);
609 for (Index_ g = 0; g < ngenes; ++g) {
610 const auto ptr = ext->fetch(vbuffer.data());
611 for (decltype(I(rank)) r = 0; r < rank; ++r) {
612 const auto mult = vptr[sanisizer::nd_offset<std::size_t>(g, ngenes, r)];
613 auto& local_buffer = local_buffers[r];
614 for (Index_ i = 0; i < length; ++i) {
615 local_buffer[i] += ptr[i] * mult;
616 }
617 }
618 }
619 }
620
621 for (decltype(I(rank)) r = 0; r < rank; ++r) {
622 for (Index_ c = 0; c < length; ++c) {
623 components.coeffRef(r, c + start) = local_buffers[r][c];
624 }
625 }
626
627 }, ncells, nthreads);
628
629 } else {
630 tatami::parallelize([&](const int, const Index_ start, const Index_ length) -> void {
631 static_assert(!EigenMatrix_::IsRowMajor);
633
634 if (mat.is_sparse()) {
635 std::vector<Index_> ibuffer(ngenes);
636 auto ext = tatami::consecutive_extractor<true>(mat, false, start, length);
637
638 for (Index_ c = start, end = start + length; c < end; ++c) {
639 const auto range = ext->fetch(vbuffer.data(), ibuffer.data());
640 static_assert(!EigenMatrix_::IsRowMajor);
641 for (decltype(I(rank)) r = 0; r < rank; ++r) {
642 auto& output = components.coeffRef(r, c);
643 output = 0;
644 const auto rotptr = scaled_rotation.data() + sanisizer::product_unsafe<std::size_t>(r, ngenes);
645 for (Index_ i = 0; i < range.number; ++i) {
646 output += rotptr[range.index[i]] * range.value[i];
647 }
648 }
649 }
650
651 } else {
652 auto ext = tatami::consecutive_extractor<false>(mat, false, start, length);
653 for (Index_ c = start, end = start + length; c < end; ++c) {
654 const auto ptr = ext->fetch(vbuffer.data());
655 static_assert(!EigenMatrix_::IsRowMajor);
656 for (decltype(I(rank)) r = 0; r < rank; ++r) {
657 const auto rotptr = scaled_rotation.data() + sanisizer::product_unsafe<std::size_t>(r, ngenes);
658 components.coeffRef(r, c) = std::inner_product(rotptr, rotptr + ngenes, ptr, static_cast<Scalar>(0));
659 }
660 }
661 }
662 }, ncells, nthreads);
663 }
664}
665
666template<class EigenMatrix_, class EigenVector_>
667void clean_up_projected(EigenMatrix_& projected, EigenVector_& D) {
668 // Empirically centering to give nice centered PCs, because we can't
669 // guarantee that the projection is centered in this manner.
670 for (decltype(I(projected.rows())) i = 0, prows = projected.rows(); i < prows; ++i) {
671 projected.row(i).array() -= projected.row(i).sum() / projected.cols();
672 }
673
674 // Just dividing by the number of observations - 1 regardless of weighting.
675 const typename EigenMatrix_::Scalar denom = projected.cols() - 1;
676 for (auto& d : D) {
677 d = d * d / denom;
678 }
679}
680
681/*******************************
682 ***** Residual wrapper ********
683 *******************************/
684
685// This wrapper class mimics multiplication with the residuals,
686// i.e., after subtracting the per-block mean from each cell.
687template<class Matrix_, typename Block_, class EigenMatrix_, class EigenVector_>
688class ResidualWrapper {
689public:
690 ResidualWrapper(const Matrix_& mat, const Block_* block, const EigenMatrix_& means) : my_mat(mat), my_block(block), my_means(means) {}
691
692public:
693 Eigen::Index rows() const { return my_mat.rows(); }
694 Eigen::Index cols() const { return my_mat.cols(); }
695
696public:
697 struct Workspace {
698 template<typename NumBlocks_>
699 Workspace(NumBlocks_ nblocks, irlba::WrappedWorkspace<Matrix_> c) :
700 sub(sanisizer::cast<decltype(I(sub.size()))>(nblocks)),
701 child(std::move(c))
702 {}
703
704 EigenVector_ sub;
705 EigenVector_ holding;
707 };
708
709 Workspace workspace() const {
710 return Workspace(my_means.rows(), irlba::wrapped_workspace(my_mat));
711 }
712
713 template<class Right_>
714 void multiply(const Right_& rhs, Workspace& work, EigenVector_& output) const {
715 const auto& realized_rhs = [&]() -> const auto& {
716 if constexpr(std::is_same<Right_, EigenVector_>::value) {
717 return rhs;
718 } else {
719 work.holding.noalias() = rhs;
720 return work.holding;
721 }
722 }();
723
724 irlba::wrapped_multiply(my_mat, realized_rhs, work.child, output);
725
726 work.sub.noalias() = my_means * realized_rhs;
727 for (decltype(I(output.size())) i = 0, end = output.size(); i < end; ++i) {
728 auto& val = output.coeffRef(i);
729 val -= work.sub.coeff(my_block[i]);
730 }
731 }
732
733public:
734 struct AdjointWorkspace {
735 template<typename NumBlocks_>
736 AdjointWorkspace(NumBlocks_ nblocks, irlba::WrappedAdjointWorkspace<Matrix_> c) :
737 aggr(sanisizer::cast<decltype(I(aggr.size()))>(nblocks)),
738 child(std::move(c))
739 {}
740
741 EigenVector_ aggr;
742 EigenVector_ holding;
744 };
745
746 AdjointWorkspace adjoint_workspace() const {
747 return AdjointWorkspace(my_means.rows(), irlba::wrapped_adjoint_workspace(my_mat));
748 }
749
750 template<class Right_>
751 void adjoint_multiply(const Right_& rhs, AdjointWorkspace& work, EigenVector_& output) const {
752 const auto& realized_rhs = [&]() {
753 if constexpr(std::is_same<Right_, EigenVector_>::value) {
754 return rhs;
755 } else {
756 work.holding.noalias() = rhs;
757 return work.holding;
758 }
759 }();
760
761 irlba::wrapped_adjoint_multiply(my_mat, realized_rhs, work.child, output);
762
763 work.aggr.setZero();
764 for (decltype(I(realized_rhs.size())) i = 0, end = realized_rhs.size(); i < end; ++i) {
765 work.aggr.coeffRef(my_block[i]) += realized_rhs.coeff(i);
766 }
767
768 output.noalias() -= my_means.adjoint() * work.aggr;
769 }
770
771public:
772 template<class EigenMatrix2_>
773 EigenMatrix2_ realize() const {
774 EigenMatrix2_ output = irlba::wrapped_realize<EigenMatrix2_>(my_mat);
775 for (decltype(I(output.rows())) i = 0, end = output.rows(); i < end; ++i) {
776 output.row(i) -= my_means.row(my_block[i]);
777 }
778 return output;
779 }
780
781private:
782 const Matrix_& my_mat;
783 const Block_* my_block;
784 const EigenMatrix_& my_means;
785};
786
787/**************************
788 ***** Dispatchers ********
789 **************************/
790
791template<bool realize_matrix_, bool sparse_, typename Value_, typename Index_, typename Block_, class EigenMatrix_, class EigenVector_>
792void run_blocked(
794 const Block_* block,
795 const BlockingDetails<Index_, EigenVector_>& block_details,
796 const BlockedPcaOptions& options,
797 EigenMatrix_& components,
798 EigenMatrix_& rotation,
799 EigenVector_& variance_explained,
800 EigenMatrix_& center_m,
801 EigenVector_& scale_v,
802 typename EigenVector_::Scalar& total_var,
803 bool& converged)
804{
805 Index_ ngenes = mat.nrow(), ncells = mat.ncol();
806
807 auto emat = [&]{
808 if constexpr(!realize_matrix_) {
809 return internal::TransposedTatamiWrapper<EigenVector_, Value_, Index_>(mat, options.num_threads);
810
811 } else if constexpr(sparse_) {
812 // 'extracted' contains row-major contents... but we implicitly transpose it to CSC with genes in columns.
814 mat,
815 /* row = */ true,
816 [&]{
818 opt.two_pass = false;
819 opt.num_threads = options.num_threads;
820 return opt;
821 }()
822 );
823 return irlba::ParallelSparseMatrix(ncells, ngenes, std::move(extracted.value), std::move(extracted.index), std::move(extracted.pointers), true, options.num_threads);
824
825 } else {
826 // Perform an implicit transposition by performing a row-major extraction into a column-major transposed matrix.
827 EigenMatrix_ emat(
828 sanisizer::cast<decltype(I(std::declval<EigenMatrix_>().rows()))>(ncells),
829 sanisizer::cast<decltype(I(std::declval<EigenMatrix_>().cols()))>(ngenes)
830 );
831 static_assert(!EigenMatrix_::IsRowMajor);
833 mat,
834 /* row_major = */ true,
835 emat.data(),
836 [&]{
837 tatami::ConvertToDenseOptions opt;
838 opt.num_threads = options.num_threads;
839 return opt;
840 }()
841 );
842 return emat;
843 }
844 }();
845
846 const auto nblocks = block_details.block_size.size();
847 center_m.resize(
848 sanisizer::cast<decltype(I(center_m.rows()))>(nblocks),
849 sanisizer::cast<decltype(I(center_m.cols()))>(ngenes)
850 );
851 sanisizer::resize(scale_v, ngenes);
852
853 if constexpr(!realize_matrix_) {
854 compute_blockwise_mean_and_variance_tatami(mat, block, block_details, center_m, scale_v, options.num_threads);
855 } else if constexpr(sparse_) {
856 compute_blockwise_mean_and_variance_realized_sparse(emat, block, block_details, center_m, scale_v, options.num_threads);
857 } else {
858 compute_blockwise_mean_and_variance_realized_dense(emat, block, block_details, center_m, scale_v, options.num_threads);
859 }
860 total_var = internal::process_scale_vector(options.scale, scale_v);
861
862 ResidualWrapper<decltype(I(emat)), Block_, EigenMatrix_, EigenVector_> centered(emat, block, center_m);
863
864 if (block_details.weighted) {
865 if (options.scale) {
866 irlba::Scaled<true, decltype(I(centered)), EigenVector_> scaled(centered, scale_v, /* divide = */ true);
867 irlba::Scaled<false, decltype(I(scaled)), EigenVector_> weighted(scaled, block_details.expanded_weights, /* divide = */ false);
868 auto out = irlba::compute(weighted, options.number, components, rotation, variance_explained, options.irlba_options);
869 converged = out.first;
870 } else {
871 irlba::Scaled<false, decltype(I(centered)), EigenVector_> weighted(centered, block_details.expanded_weights, /* divide = */ false);
872 auto out = irlba::compute(weighted, options.number, components, rotation, variance_explained, options.irlba_options);
873 converged = out.first;
874 }
875
876 EigenMatrix_ tmp;
877 const auto& scaled_rotation = scale_rotation_matrix(rotation, options.scale, scale_v, tmp);
878
879 // This transposes 'components' to be a NDIM * NCELLS matrix.
880 if constexpr(!realize_matrix_) {
881 project_matrix_transposed_tatami(mat, components, scaled_rotation, options.num_threads);
882 } else if constexpr(sparse_) {
883 project_matrix_realized_sparse(emat, components, scaled_rotation, options.num_threads);
884 } else {
885 components.noalias() = (emat * scaled_rotation).adjoint();
886 }
887
888 // Subtracting each block's mean from the PCs.
889 if (options.components_from_residuals) {
890 EigenMatrix_ centering = (center_m * scaled_rotation).adjoint();
891 for (decltype(I(ncells)) c =0 ; c < ncells; ++c) {
892 components.col(c) -= centering.col(block[c]);
893 }
894 }
895
896 clean_up_projected(components, variance_explained);
897 if (!options.transpose) {
898 components.adjointInPlace();
899 }
900
901 } else {
902 if (options.scale) {
903 irlba::Scaled<true, decltype(I(centered)), EigenVector_> scaled(centered, scale_v, /* divide = */ true);
904 const auto out = irlba::compute(scaled, options.number, components, rotation, variance_explained, options.irlba_options);
905 converged = out.first;
906 } else {
907 const auto out = irlba::compute(centered, options.number, components, rotation, variance_explained, options.irlba_options);
908 converged = out.first;
909 }
910
911 if (options.components_from_residuals) {
912 internal::clean_up(mat.ncol(), components, variance_explained);
913 if (options.transpose) {
914 components.adjointInPlace();
915 }
916 } else {
917 EigenMatrix_ tmp;
918 const auto& scaled_rotation = scale_rotation_matrix(rotation, options.scale, scale_v, tmp);
919
920 // This transposes 'components' to be a NDIM * NCELLS matrix.
921 if constexpr(!realize_matrix_) {
922 project_matrix_transposed_tatami(mat, components, scaled_rotation, options.num_threads);
923 } else if constexpr(sparse_) {
924 project_matrix_realized_sparse(emat, components, scaled_rotation, options.num_threads);
925 } else {
926 components.noalias() = (emat * scaled_rotation).adjoint();
927 }
928
929 clean_up_projected(components, variance_explained);
930 if (!options.transpose) {
931 components.adjointInPlace();
932 }
933 }
934 }
935}
936
937}
948template<typename EigenMatrix_, typename EigenVector_>
956 EigenMatrix_ components;
957
963 EigenVector_ variance_explained;
964
969 typename EigenVector_::Scalar total_variance = 0;
970
976 EigenMatrix_ rotation;
977
983 EigenMatrix_ center;
984
989 EigenVector_ scale;
990
994 bool converged = false;
995};
996
1044template<typename Value_, typename Index_, typename Block_, typename EigenMatrix_, class EigenVector_>
1047 auto bdetails = internal::compute_blocking_details<EigenVector_>(mat.ncol(), block, options.block_weight_policy, options.variable_block_weight_parameters);
1048
1049 EigenMatrix_& components = output.components;
1050 EigenMatrix_& rotation = output.rotation;
1051 EigenVector_& variance_explained = output.variance_explained;
1052 EigenMatrix_& center_m = output.center;
1053 EigenVector_& scale_v = output.scale;
1054 auto& total_var = output.total_variance;
1055 bool& converged = output.converged;
1056
1057 if (mat.sparse()) {
1058 if (options.realize_matrix) {
1059 internal::run_blocked<true, true>(mat, block, bdetails, options, components, rotation, variance_explained, center_m, scale_v, total_var, converged);
1060 } else {
1061 internal::run_blocked<false, true>(mat, block, bdetails, options, components, rotation, variance_explained, center_m, scale_v, total_var, converged);
1062 }
1063 } else {
1064 if (options.realize_matrix) {
1065 internal::run_blocked<true, false>(mat, block, bdetails, options, components, rotation, variance_explained, center_m, scale_v, total_var, converged);
1066 } else {
1067 internal::run_blocked<false, false>(mat, block, bdetails, options, components, rotation, variance_explained, center_m, scale_v, total_var, converged);
1068 }
1069 }
1070
1071 if (!options.scale) {
1072 output.scale = EigenVector_();
1073 }
1074}
1075
1095template<typename EigenMatrix_ = Eigen::MatrixXd, class EigenVector_ = Eigen::VectorXd, typename Value_, typename Index_, typename Block_>
1098 blocked_pca(mat, block, options, output);
1099 return output;
1100}
1101
1102}
1103
1104#endif
virtual Index_ ncol() const=0
virtual Index_ nrow() const=0
virtual bool prefer_rows() const=0
virtual bool is_sparse() const=0
virtual std::unique_ptr< MyopicSparseExtractor< Value_, Index_ > > sparse(bool row, const Options &opt) const=0
void wrapped_adjoint_multiply(const Matrix_ &matrix, const Right_ &rhs, WrappedAdjointWorkspace< Matrix_ > &work, EigenVector_ &out)
WrappedWorkspace< Matrix_ > wrapped_workspace(const Matrix_ &matrix)
std::pair< bool, int > compute(const Matrix_ &matrix, Eigen::Index number, EigenMatrix_ &outU, EigenMatrix_ &outV, EigenVector_ &outD, const Options &options)
void parallelize(Task_ num_tasks, Run_ run_task)
void wrapped_multiply(const Matrix_ &matrix, const Right_ &rhs, WrappedWorkspace< Matrix_ > &work, EigenVector_ &out)
WrappedAdjointWorkspace< Matrix_ > wrapped_adjoint_workspace(const Matrix_ &matrix)
typename get_adjoint_workspace< Matrix_ >::type WrappedAdjointWorkspace
EigenMatrix_ wrapped_realize(const Matrix_ &matrix)
typename get_workspace< Matrix_ >::type WrappedWorkspace
double compute_variable_weight(double s, const VariableWeightParameters &params)
Principal component analysis on single-cell data.
void blocked_pca(const tatami::Matrix< Value_, Index_ > &mat, const Block_ *block, const BlockedPcaOptions &options, BlockedPcaResults< EigenMatrix_, EigenVector_ > &output)
Definition blocked_pca.hpp:1045
CompressedSparseContents< StoredValue_, StoredIndex_, StoredPointer_ > retrieve_compressed_sparse_contents(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const RetrieveCompressedSparseContentsOptions &options)
decltype(std::declval< Container_ >().size()) cast_Index_to_container_size(Index_ x)
void parallelize(Function_ fun, Index_ tasks, int threads)
Container_ create_container_of_Index_size(Index_ x, Args_ &&... args)
void convert_to_dense(const Matrix< InputValue_, InputIndex_ > &matrix, bool row_major, StoredValue_ *store, const ConvertToDenseOptions &options)
auto consecutive_extractor(const Matrix< Value_, Index_ > &matrix, bool row, Index_ iter_start, Index_ iter_length, Args_ &&... args)
Options for blocked_pca().
Definition blocked_pca.hpp:30
scran_blocks::VariableWeightParameters variable_block_weight_parameters
Definition blocked_pca.hpp:78
int num_threads
Definition blocked_pca.hpp:97
bool components_from_residuals
Definition blocked_pca.hpp:85
bool realize_matrix
Definition blocked_pca.hpp:91
irlba::Options irlba_options
Definition blocked_pca.hpp:102
bool transpose
Definition blocked_pca.hpp:61
scran_blocks::WeightPolicy block_weight_policy
Definition blocked_pca.hpp:72
int number
Definition blocked_pca.hpp:47
bool scale
Definition blocked_pca.hpp:55
Results of blocked_pca().
Definition blocked_pca.hpp:949
EigenVector_::Scalar total_variance
Definition blocked_pca.hpp:969
bool converged
Definition blocked_pca.hpp:994
EigenMatrix_ components
Definition blocked_pca.hpp:956
EigenMatrix_ rotation
Definition blocked_pca.hpp:976
EigenVector_ scale
Definition blocked_pca.hpp:989
EigenMatrix_ center
Definition blocked_pca.hpp:983
EigenVector_ variance_explained
Definition blocked_pca.hpp:963