tatami
C++ API for different matrix representations
Loading...
Searching...
No Matches
convert_to_compressed_sparse.hpp
Go to the documentation of this file.
1#ifndef TATAMI_CONVERT_TO_COMPRESSED_SPARSE_H
2#define TATAMI_CONVERT_TO_COMPRESSED_SPARSE_H
3
8#include "../utils/Index_to_container.hpp"
9
10#include <memory>
11#include <vector>
12#include <cstddef>
13
20namespace tatami {
21
25namespace convert_to_compressed_sparse_internal {
26
27template<typename Value_, typename Index_, typename Count_>
28void count_compressed_sparse_non_zeros_consistent(const tatami::Matrix<Value_, Index_>& matrix, Index_ primary, Index_ secondary, bool row, Count_* output, int threads) {
29 if (matrix.is_sparse()) {
30 Options opt;
31 opt.sparse_extract_value = false;
32 opt.sparse_extract_index = false;
33 opt.sparse_ordered_index = false;
34
35 parallelize([&](int, Index_ start, Index_ length) -> void {
36 auto wrk = consecutive_extractor<true>(matrix, row, start, length, opt);
37 for (Index_ x = 0; x < length; ++x) {
38 auto range = wrk->fetch(NULL, NULL);
39 output[start + x] = range.number;
40 }
41 }, primary, threads);
42
43 } else {
44 parallelize([&](int, Index_ start, Index_ length) -> void {
45 std::vector<Value_> buffer_v(secondary);
46 auto wrk = consecutive_extractor<false>(matrix, row, start, length);
47 for (Index_ p = start, pe = start + length; p < pe; ++p) {
48 auto ptr = wrk->fetch(buffer_v.data());
49 Count_ count = 0;
50 for (Index_ s = 0; s < secondary; ++s) {
51 count += (ptr[s] != 0);
52 }
53 output[p] = count;
54 }
55 }, primary, threads);
56 }
57}
58
59template<typename Value_, typename Index_, typename Count_>
60void count_compressed_sparse_non_zeros_inconsistent(const tatami::Matrix<Value_, Index_>& matrix, Index_ primary, Index_ secondary, bool row, Count_* output, int threads) {
61 auto nz_counts = sanisizer::create<std::vector<std::vector<Count_> > >(threads - 1);
62 for (auto& x : nz_counts) {
63 x.resize(primary);
64 }
65
66 if (matrix.is_sparse()) {
67 Options opt;
68 opt.sparse_extract_value = false;
69 opt.sparse_ordered_index = false;
70
71 parallelize([&](int t, Index_ start, Index_ length) -> void {
72 auto wrk = consecutive_extractor<true>(matrix, !row, start, length, opt);
73 auto buffer_i = create_container_of_Index_size<std::vector<Index_> >(primary);
74 auto my_counts = (t > 0 ? nz_counts[t - 1].data() : output);
75
76 for (Index_ x = 0; x < length; ++x) {
77 auto range = wrk->fetch(NULL, buffer_i.data());
78 for (Index_ i = 0; i < range.number; ++i) {
79 ++my_counts[range.index[i]];
80 }
81 }
82 }, secondary, threads);
83
84 } else {
85 parallelize([&](int t, Index_ start, Index_ length) -> void {
86 auto wrk = consecutive_extractor<false>(matrix, !row, start, length);
87 auto buffer_v = create_container_of_Index_size<std::vector<Value_> >(primary);
88 auto my_counts = (t > 0 ? nz_counts[t - 1].data() : output);
89
90 for (Index_ x = 0; x < length; ++x) {
91 auto ptr = wrk->fetch(buffer_v.data());
92 for (Index_ p = 0; p < primary; ++p) {
93 my_counts[p] += (ptr[p] != 0);
94 }
95 }
96 }, secondary, threads);
97 }
98
99 for (auto& y : nz_counts) {
100 for (Index_ p = 0; p < primary; ++p) {
101 output[p] += y[p];
102 }
103 }
104}
105
106template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
107void fill_compressed_sparse_matrix_consistent(
109 InputIndex_ primary,
110 InputIndex_ secondary,
111 bool row,
112 const Pointer_* pointers,
113 StoredValue_* output_value,
114 StoredIndex_* output_index,
115 int threads)
116{
117 if (matrix.is_sparse()) {
118 Options opt;
119 opt.sparse_ordered_index = false;
120
121 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
122 auto wrk = consecutive_extractor<true>(matrix, row, start, length, opt);
123 auto buffer_v = create_container_of_Index_size<std::vector<InputValue_> >(secondary);
124 auto buffer_i = create_container_of_Index_size<std::vector<InputIndex_> >(secondary);
125
126 for (InputIndex_ p = start, pe = start + length; p < pe; ++p) {
127 // Resist the urge to `fetch()` straight into 'output_v'
128 // and 'output_i', as implementations may assume that they
129 // have the entire 'length' length to play with, and the
130 // output vectors only have whatever is allocated from the
131 // first pass (which might be nothing for an all-zero matrix).
132 auto range = wrk->fetch(buffer_v.data(), buffer_i.data());
133 auto offset = pointers[p];
134 std::copy_n(range.value, range.number, output_value + offset);
135 std::copy_n(range.index, range.number, output_index + offset);
136 }
137 }, primary, threads);
138
139 } else {
140 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
141 std::vector<InputValue_> buffer_v(secondary);
142 auto wrk = consecutive_extractor<false>(matrix, row, start, length);
143
144 for (InputIndex_ p = start, pe = start + length; p < pe; ++p) {
145 auto ptr = wrk->fetch(buffer_v.data());
146 auto offset = pointers[p];
147 for (InputIndex_ s = 0; s < secondary; ++s) {
148 auto val = ptr[s];
149 if (val != 0) {
150 output_value[offset] = val;
151 output_index[offset] = s;
152 ++offset;
153 }
154 }
155 }
156 }, primary, threads);
157 }
158}
159
160template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
161void fill_compressed_sparse_matrix_inconsistent(
163 InputIndex_ primary,
164 InputIndex_ secondary,
165 bool row,
166 const Pointer_* pointers,
167 StoredValue_* output_value,
168 StoredIndex_* output_index,
169 int threads)
170{
171 if (matrix.is_sparse()) {
172 Options opt;
173 opt.sparse_ordered_index = false;
174
175 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
176 auto wrk = consecutive_extractor<true>(matrix, !row, static_cast<InputIndex_>(0), secondary, start, length, opt);
177 auto buffer_v = create_container_of_Index_size<std::vector<InputValue_> >(length);
178 auto buffer_i = create_container_of_Index_size<std::vector<InputIndex_> >(length);
179 std::vector<Pointer_> offset_copy(pointers + start, pointers + start + length);
180
181 for (InputIndex_ x = 0; x < secondary; ++x) {
182 auto range = wrk->fetch(buffer_v.data(), buffer_i.data());
183 for (InputIndex_ i = 0; i < range.number; ++i) {
184 auto& pos = offset_copy[range.index[i] - start];
185 output_value[pos] = range.value[i];
186 output_index[pos] = x;
187 ++pos;
188 }
189 }
190 }, primary, threads);
191
192 } else {
193 parallelize([&](int, InputIndex_ start, InputIndex_ length) -> void {
194 auto wrk = consecutive_extractor<false>(matrix, !row, static_cast<InputIndex_>(0), secondary, start, length);
195 auto buffer_v = create_container_of_Index_size<std::vector<InputValue_> >(length);
196 std::vector<Pointer_> offset_copy(pointers + start, pointers + start + length);
197
198 for (InputIndex_ x = 0; x < secondary; ++x) {
199 auto ptr = wrk->fetch(buffer_v.data());
200 for (InputIndex_ p = 0; p < length; ++p) {
201 auto val = ptr[p];
202 if (val != 0) {
203 auto& pos = offset_copy[p];
204 output_value[pos] = val;
205 output_index[pos] = x;
206 ++pos;
207 }
208 }
209 }
210 }, primary, threads);
211 }
212}
213
214}
228
244template<typename Value_, typename Index_, typename Count_>
246 Index_ NR = matrix.nrow();
247 Index_ NC = matrix.ncol();
248 Index_ primary = (row ? NR : NC);
249 Index_ secondary = (row ? NC : NR);
250 std::fill_n(output, primary, 0);
251
252 if (row == matrix.prefer_rows()) {
253 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_consistent(matrix, primary, secondary, row, output, options.num_threads);
254 } else {
255 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_inconsistent(matrix, primary, secondary, row, output, options.num_threads);
256 }
257}
258
268
288template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
291 bool row,
292 const Pointer_* pointers,
293 StoredValue_* output_value,
294 StoredIndex_* output_index,
296{
297 InputIndex_ NR = matrix.nrow();
298 InputIndex_ NC = matrix.ncol();
299 InputIndex_ primary = (row ? NR : NC);
300 InputIndex_ secondary = (row ? NC : NR);
301
302 if (row == matrix.prefer_rows()) {
303 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_consistent(matrix, primary, secondary, row, pointers, output_value, output_index, options.num_threads);
304 } else {
305 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_inconsistent(matrix, primary, secondary, row, pointers, output_value, output_index, options.num_threads);
306 }
307}
308
319template<typename Value_, typename Index_, typename Pointer_>
324 std::vector<Value_> value;
325
329 std::vector<Index_> index;
330
334 std::vector<Pointer_> pointers;
335};
336
352
370template<typename StoredValue_, typename StoredIndex_, typename StoredPointer_ = std::size_t, typename InputValue_, typename InputIndex_>
373 bool row,
375{
376 // We use size_t as the default pointer type here, as our output consists of vectors
377 // with the default allocator, for which the size_type is unlikely to be bigger than size_t.
379 auto& output_v = output.value;
380 auto& output_i = output.index;
381 auto& output_p = output.pointers;
382
383 InputIndex_ NR = matrix.nrow();
384 InputIndex_ NC = matrix.ncol();
385 InputIndex_ primary = (row ? NR : NC);
386 InputIndex_ secondary = (row ? NC : NR);
387
388 if (!options.two_pass) {
389 // Doing a single fragmented run and then concatenating everything together.
391 matrix,
392 row,
393 [&]{
395 roptions.num_threads = options.num_threads;
396 return roptions;
397 }()
398 );
399 const auto& store_v = frag.value;
400 const auto& store_i = frag.index;
401
402 output_p.resize(static_cast<std::size_t>(primary) + 1);
403 for (InputIndex_ p = 0; p < primary; ++p) {
404 output_p[p + 1] = output_p[p] + store_v[p].size();
405 }
406
407 output_v.reserve(output_p.back());
408 output_i.reserve(output_p.back());
409 for (InputIndex_ p = 0; p < primary; ++p) {
410 output_v.insert(output_v.end(), store_v[p].begin(), store_v[p].end());
411 output_i.insert(output_i.end(), store_i[p].begin(), store_i[p].end());
412 }
413
414 } else if (row == matrix.prefer_rows()) {
415 // First pass to figure out how many non-zeros there are.
416 output_p.resize(static_cast<std::size_t>(primary) + 1);
417 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_consistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
418 for (InputIndex_ i = 1; i <= primary; ++i) {
419 output_p[i] += output_p[i - 1];
420 }
421
422 // Second pass to actually fill our vectors.
423 output_v.resize(output_p.back());
424 output_i.resize(output_p.back());
425 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_consistent(
426 matrix,
427 primary,
428 secondary,
429 row,
430 output_p.data(),
431 output_v.data(),
432 output_i.data(),
433 options.num_threads
434 );
435
436 } else {
437 // First pass to figure out how many non-zeros there are.
438 output_p.resize(static_cast<std::size_t>(primary) + 1);
439 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_inconsistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
440 for (InputIndex_ i = 1; i <= primary; ++i) {
441 output_p[i] += output_p[i - 1];
442 }
443
444 // Second pass to actually fill our vectors.
445 output_v.resize(output_p.back());
446 output_i.resize(output_p.back());
447 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_inconsistent(
448 matrix,
449 primary,
450 secondary,
451 row,
452 output_p.data(),
453 output_v.data(),
454 output_i.data(),
455 options.num_threads
456 );
457 }
458
459 return output;
460}
461
470 bool two_pass = false;
471
475 int num_threads = 1;
476};
477
495template<
496 typename Value_,
497 typename Index_,
498 typename StoredValue_ = Value_,
499 typename StoredIndex_ = Index_,
500 typename StoredPointer_ = std::size_t,
501 typename InputValue_,
502 typename InputIndex_
503>
504std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>& matrix, bool row, const ConvertToCompressedSparseOptions& options) {
506 matrix,
507 row,
508 [&]{
510 ropt.two_pass = options.two_pass;
511 ropt.num_threads = options.num_threads;
512 return ropt;
513 }()
514 );
515 return std::shared_ptr<Matrix<Value_, Index_> >(
517 Value_,
518 Index_,
519 std::vector<StoredValue_>,
520 std::vector<StoredIndex_>,
521 std::vector<StoredPointer_>
522 >(
523 matrix.nrow(),
524 matrix.ncol(),
525 std::move(comp.value),
526 std::move(comp.index),
527 std::move(comp.pointers),
528 row,
529 []{
530 CompressedSparseMatrixOptions copt;
531 copt.check = false; // no need for checks, as we guarantee correctness.
532 return copt;
533 }()
534 )
535 );
536}
537
541// Backwards compatbility.
542template<typename Value_, typename Index_, typename Count_>
543void count_compressed_sparse_non_zeros(const tatami::Matrix<Value_, Index_>* matrix, bool row, Count_* output, int threads) {
545 *matrix,
546 row,
547 output,
548 [&]{
549 CountCompressedSparseNonZerosOptions copt;
550 copt.num_threads = threads;
551 return copt;
552 }()
553 );
554}
555
556template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
558 bool row,
559 const Pointer_* pointers,
560 StoredValue_* output_value,
561 StoredIndex_* output_index,
562 int threads)
563{
565 *matrix,
566 row,
567 pointers,
568 output_value,
569 output_index,
570 [&]{
571 FillCompressedSparseContentsOptions fopt;
572 fopt.num_threads = threads;
573 return fopt;
574 }()
575 );
576}
577
578template<typename StoredValue_, typename StoredIndex_, typename StoredPointer_ = std::size_t, typename InputValue_, typename InputIndex_>
579CompressedSparseContents<StoredValue_, StoredIndex_, StoredPointer_> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass, int threads = 1) {
581 *matrix,
582 row,
583 [&]{
584 RetrieveCompressedSparseContentsOptions opt;
585 opt.two_pass = two_pass;
586 opt.num_threads = threads;
587 return opt;
588 }()
589 );
590}
591
592template<typename Value_ = double, typename Index_ = int, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
593std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass = false, int threads = 1) {
595 *matrix,
596 row,
597 [&]{
598 ConvertToCompressedSparseOptions opt;
599 opt.two_pass = two_pass;
600 opt.num_threads = threads;
601 return opt;
602 }()
603 );
604}
605
606template <bool row_, typename Value_, typename Index_, typename InputValue_, typename InputIndex_>
607CompressedSparseContents<Value_, Index_, std::size_t> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass, int threads = 1) {
608 return retrieve_compressed_sparse_contents<Value_, Index_>(matrix, row_, two_pass, threads);
609}
610
611template <bool row_, typename Value_, typename Index_, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
612std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass = false, int threads = 1) {
614}
619}
620
621#endif
Compressed sparse matrix representation.
Compressed sparse matrix representation.
Definition CompressedSparseMatrix.hpp:503
Virtual class for a matrix.
Definition Matrix.hpp:59
virtual Index_ ncol() const =0
virtual Index_ nrow() const =0
virtual bool prefer_rows() const =0
virtual bool is_sparse() const =0
Templated construction of a new consecutive extractor.
Convert a matrix into a fragmented sparse format.
Flexible representations for matrix data.
Definition Extractor.hpp:15
CompressedSparseContents< StoredValue_, StoredIndex_, StoredPointer_ > retrieve_compressed_sparse_contents(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const RetrieveCompressedSparseContentsOptions &options)
Definition convert_to_compressed_sparse.hpp:371
std::shared_ptr< Matrix< Value_, Index_ > > convert_to_compressed_sparse(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const ConvertToCompressedSparseOptions &options)
Definition convert_to_compressed_sparse.hpp:504
void parallelize(Function_ fun, Index_ tasks, int threads)
Definition parallelize.hpp:42
FragmentedSparseContents< StoredValue_, StoredIndex_ > retrieve_fragmented_sparse_contents(const Matrix< InputValue_, InputIndex_ > &matrix, bool row, const RetrieveFragmentedSparseContentsOptions &options)
Definition convert_to_fragmented_sparse.hpp:81
void fill_compressed_sparse_contents(const tatami::Matrix< InputValue_, InputIndex_ > &matrix, bool row, const Pointer_ *pointers, StoredValue_ *output_value, StoredIndex_ *output_index, const FillCompressedSparseContentsOptions &options)
Definition convert_to_compressed_sparse.hpp:289
void count_compressed_sparse_non_zeros(const tatami::Matrix< Value_, Index_ > &matrix, bool row, Count_ *output, const CountCompressedSparseNonZerosOptions &options)
Definition convert_to_compressed_sparse.hpp:245
Parallelized iteration over a tatami::Matrix.
Compressed sparse contents.
Definition convert_to_compressed_sparse.hpp:320
std::vector< Index_ > index
Definition convert_to_compressed_sparse.hpp:329
std::vector< Value_ > value
Definition convert_to_compressed_sparse.hpp:324
std::vector< Pointer_ > pointers
Definition convert_to_compressed_sparse.hpp:334
Options for convert_to_compressed_sparse().
Definition convert_to_compressed_sparse.hpp:465
bool two_pass
Definition convert_to_compressed_sparse.hpp:470
int num_threads
Definition convert_to_compressed_sparse.hpp:475
Options for count_compressed_sparse_non_zeros().
Definition convert_to_compressed_sparse.hpp:222
int num_threads
Definition convert_to_compressed_sparse.hpp:226
Options for fill_compressed_sparse_contents().
Definition convert_to_compressed_sparse.hpp:262
int num_threads
Definition convert_to_compressed_sparse.hpp:266
Options for retrieve_compressed_sparse_contents().
Definition convert_to_compressed_sparse.hpp:340
int num_threads
Definition convert_to_compressed_sparse.hpp:350
bool two_pass
Definition convert_to_compressed_sparse.hpp:345
Options for retrieve_fragmented_sparse_contents().
Definition convert_to_fragmented_sparse.hpp:61
int num_threads
Definition convert_to_fragmented_sparse.hpp:65