tatami
C++ API for different matrix representations
All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends Pages
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 output_p.resize(sanisizer::sum<decltype(output_p.size())>(primary, 1));
389
390 if (!options.two_pass) {
391 // Doing a single fragmented run and then concatenating everything together.
393 matrix,
394 row,
395 [&]{
397 roptions.num_threads = options.num_threads;
398 return roptions;
399 }()
400 );
401 const auto& store_v = frag.value;
402 const auto& store_i = frag.index;
403
404 for (InputIndex_ p = 0; p < primary; ++p) {
405 output_p[p + 1] = output_p[p] + store_v[p].size();
406 }
407
408 output_v.reserve(output_p.back());
409 output_i.reserve(output_p.back());
410 for (InputIndex_ p = 0; p < primary; ++p) {
411 output_v.insert(output_v.end(), store_v[p].begin(), store_v[p].end());
412 output_i.insert(output_i.end(), store_i[p].begin(), store_i[p].end());
413 }
414
415 } else if (row == matrix.prefer_rows()) {
416 // First pass to figure out how many non-zeros there are.
417 output_p.resize(sanisizer::sum<decltype(output_p.size())>(primary, 1));
418 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_consistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
419 for (InputIndex_ i = 1; i <= primary; ++i) {
420 output_p[i] += output_p[i - 1];
421 }
422
423 // Second pass to actually fill our vectors.
424 output_v.resize(output_p.back());
425 output_i.resize(output_p.back());
426 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_consistent(
427 matrix,
428 primary,
429 secondary,
430 row,
431 output_p.data(),
432 output_v.data(),
433 output_i.data(),
434 options.num_threads
435 );
436
437 } else {
438 // First pass to figure out how many non-zeros there are.
439 output_p.resize(sanisizer::sum<decltype(output_p.size())>(primary, 1));
440 convert_to_compressed_sparse_internal::count_compressed_sparse_non_zeros_inconsistent(matrix, primary, secondary, row, output_p.data() + 1, options.num_threads);
441 for (InputIndex_ i = 1; i <= primary; ++i) {
442 output_p[i] += output_p[i - 1];
443 }
444
445 // Second pass to actually fill our vectors.
446 output_v.resize(output_p.back());
447 output_i.resize(output_p.back());
448 convert_to_compressed_sparse_internal::fill_compressed_sparse_matrix_inconsistent(
449 matrix,
450 primary,
451 secondary,
452 row,
453 output_p.data(),
454 output_v.data(),
455 output_i.data(),
456 options.num_threads
457 );
458 }
459
460 return output;
461}
462
471 bool two_pass = false;
472
476 int num_threads = 1;
477};
478
496template<
497 typename Value_,
498 typename Index_,
499 typename StoredValue_ = Value_,
500 typename StoredIndex_ = Index_,
501 typename StoredPointer_ = std::size_t,
502 typename InputValue_,
503 typename InputIndex_
504>
505std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>& matrix, bool row, const ConvertToCompressedSparseOptions& options) {
507 matrix,
508 row,
509 [&]{
511 ropt.two_pass = options.two_pass;
512 ropt.num_threads = options.num_threads;
513 return ropt;
514 }()
515 );
516 return std::shared_ptr<Matrix<Value_, Index_> >(
518 Value_,
519 Index_,
520 std::vector<StoredValue_>,
521 std::vector<StoredIndex_>,
522 std::vector<StoredPointer_>
523 >(
524 matrix.nrow(),
525 matrix.ncol(),
526 std::move(comp.value),
527 std::move(comp.index),
528 std::move(comp.pointers),
529 row,
530 []{
531 CompressedSparseMatrixOptions copt;
532 copt.check = false; // no need for checks, as we guarantee correctness.
533 return copt;
534 }()
535 )
536 );
537}
538
542// Backwards compatbility.
543template<typename Value_, typename Index_, typename Count_>
544void count_compressed_sparse_non_zeros(const tatami::Matrix<Value_, Index_>* matrix, bool row, Count_* output, int threads) {
546 *matrix,
547 row,
548 output,
549 [&]{
550 CountCompressedSparseNonZerosOptions copt;
551 copt.num_threads = threads;
552 return copt;
553 }()
554 );
555}
556
557template<typename InputValue_, typename InputIndex_, typename Pointer_, typename StoredValue_, typename StoredIndex_>
559 bool row,
560 const Pointer_* pointers,
561 StoredValue_* output_value,
562 StoredIndex_* output_index,
563 int threads)
564{
566 *matrix,
567 row,
568 pointers,
569 output_value,
570 output_index,
571 [&]{
572 FillCompressedSparseContentsOptions fopt;
573 fopt.num_threads = threads;
574 return fopt;
575 }()
576 );
577}
578
579template<typename StoredValue_, typename StoredIndex_, typename StoredPointer_ = std::size_t, typename InputValue_, typename InputIndex_>
580CompressedSparseContents<StoredValue_, StoredIndex_, StoredPointer_> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass, int threads = 1) {
582 *matrix,
583 row,
584 [&]{
585 RetrieveCompressedSparseContentsOptions opt;
586 opt.two_pass = two_pass;
587 opt.num_threads = threads;
588 return opt;
589 }()
590 );
591}
592
593template<typename Value_ = double, typename Index_ = int, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
594std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool row, bool two_pass = false, int threads = 1) {
596 *matrix,
597 row,
598 [&]{
599 ConvertToCompressedSparseOptions opt;
600 opt.two_pass = two_pass;
601 opt.num_threads = threads;
602 return opt;
603 }()
604 );
605}
606
607template <bool row_, typename Value_, typename Index_, typename InputValue_, typename InputIndex_>
608CompressedSparseContents<Value_, Index_, std::size_t> retrieve_compressed_sparse_contents(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass, int threads = 1) {
609 return retrieve_compressed_sparse_contents<Value_, Index_>(matrix, row_, two_pass, threads);
610}
611
612template <bool row_, typename Value_, typename Index_, typename StoredValue_ = Value_, typename StoredIndex_ = Index_, typename InputValue_, typename InputIndex_>
613std::shared_ptr<Matrix<Value_, Index_> > convert_to_compressed_sparse(const Matrix<InputValue_, InputIndex_>* matrix, bool two_pass = false, int threads = 1) {
615}
620}
621
622#endif
Compressed sparse matrix representation.
Compressed sparse matrix representation.
Definition CompressedSparseMatrix.hpp:505
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:505
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:466
bool two_pass
Definition convert_to_compressed_sparse.hpp:471
int num_threads
Definition convert_to_compressed_sparse.hpp:476
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