-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathstack.cpp
More file actions
150 lines (136 loc) · 6.26 KB
/
stack.cpp
File metadata and controls
150 lines (136 loc) · 6.26 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include <unordered_map>
#include <string>
#include <Rcpp.h>
using namespace Rcpp;
using namespace std;
template<int RTYPE> Vector<RTYPE> cpp_stack_impl(List array_list, int along, Vector<RTYPE> fill, bool ovr) {
auto dimnames = vector<vector<String>>(along); // dim: names along
auto axmap = vector<unordered_map<string, int>>(along); // dim: element name->index
auto ax_unnamed = vector<int>(along); // counter for unnamed dimension elements
auto a2r = vector<vector<vector<int>>>(array_list.size()); // index array>dim>element
// create lookup tables for all present dimension names
for (int ai=0; ai<Rf_xlength(array_list); ai++) { // array index
auto a = as<Vector<RTYPE>>(array_list[ai]);
auto dim = as<vector<int>>(a.attr("dim"));
auto dn = as<List>(a.attr("dimnames"));
if (along == dim.size()+1) { // along introduces new dimension
if (array_list.attr("names") == R_NilValue)
dn.push_back(CharacterVector::create(NA_STRING));
else
dn.push_back(as<vector<string>>(array_list.attr("names"))[ai]);
dim.push_back(1);
}
a2r[ai] = vector<vector<int>>(dim.size());
if (dimnames.size() < dim.size()) {
dimnames.resize(dim.size());
axmap.resize(dim.size());
ax_unnamed.resize(dim.size());
}
for (int d=0; d<dim.size(); d++) { // dimension in array
if (dn.size() == 0 || dn[d] == R_NilValue) {
for (int e=0; e<dim[d]; e++) {
// Rprintf("array %i dim %i: %i -> %i\n", ai, d, e, axmap[d].size() + ax_unnamed[d]);
a2r[ai][d].push_back(axmap[d].size() + ax_unnamed[d]++);
dimnames[d].push_back(NA_STRING);
}
} else {
auto dni = as<vector<string>>(dn[d]);
for (int e=0; e<dim[d]; e++) { // element in dimension
auto it = axmap[d].find(dni[e]);
if (it == axmap[d].end()) {
it = axmap[d].emplace(dni[e], axmap[d].size() + ax_unnamed[d]).first;
dimnames[d].push_back(dni[e]);
}
a2r[ai][d].push_back(it->second);
// Rprintf("array %i dim %i: %s -> %i\n", ai, d, dni[e].c_str(), it->second);
}
}
}
}
/* for (int ai=0; ai<a2r.size(); ai++)
for (int di=0; di<a2r[ai].size(); di++) {
cout << "*** array " << ai << " dim " << di << ": ";
copy(a2r[ai][di].begin(), a2r[ai][di].end(), ostream_iterator<int>(cout, " "));
cout << "\n";
}
*/
// create result array with attributes
auto rdim = IntegerVector(dimnames.size());
for (int ai=0; ai<Rf_xlength(array_list); ai++)
if (a2r[ai].size() != rdim.size())
stop("Names are required for all dimensions except the one stacked along.\n%s",
" Use bind() if you want to just bind together arrays without names.");
auto rdnames = List(dimnames.size());
for (int i=0; i<dimnames.size(); i++) {
rdim[i] = dimnames[i].size();
auto rdni = CharacterVector(dimnames[i].size());
for (int j=0; j<rdni.size(); j++)
rdni[j] = dimnames[i][j];
if (all(is_na(rdni)))
rdnames[i] = R_NilValue;
else
rdnames[i] = rdni;
}
auto n = accumulate(rdim.begin(), rdim.end(), 1, multiplies<int>());
auto result = Vector<RTYPE>(n, fill[0]);
result.attr("dim") = rdim;
result.attr("dimnames") = rdnames;
// fill the result array
int maxdim = rdim.size() - 1;
for (int ai=0; ai<Rf_xlength(array_list); ai++) { // each original array
auto a = as<Vector<RTYPE>>(array_list[ai]);
auto it = vector<vector<int>::iterator>(a2r[ai].size()); // each result dim
for (int d=0; d<it.size(); d++)
it[d] = a2r[ai][d].begin();
int dim_offset;
bool new_offset = true;
for (int aidx=0; aidx<a.size(); aidx++) { // element in original array
if (new_offset) { // dimension jump
dim_offset = 0;
int dim_mult = 1;
for (int d=0; d<maxdim; d++) {
dim_mult *= rdim[d];
dim_offset += dim_mult * *it[d+1];
}
new_offset = false;
}
// Rprintf("result[%i] = a[%i][%i]\n", *it[0] + dim_offset, ai, aidx);
int ri = *it[0] + dim_offset;
auto newval = a[aidx];
// same-value comparisons below: catch double NaN equality
if (ovr || (result[ri] == fill[0]) || (!(result[ri] == result[ri])))
result[ri] = newval;
else
if (!((newval == fill[0]) || (result[ri] == newval) || (!(newval == newval))))
stop("Different values on same position and allow_overwrite=FALSE");
it[0]++;
for (int d=0; d<maxdim; d++) { // check if we're jumping dimensions
if (it[d] != a2r[ai][d].end()) // dim+1 jump only if dim jump
break;
new_offset = true;
it[d] = a2r[ai][d].begin();
it[d+1]++;
}
};
}
return result;
}
// [[Rcpp::export]]
SEXP cpp_stack(List array_list, int along, SEXP fill, bool ovr) {
auto max_type = NILSXP;
for (int ai=0; ai<array_list.size(); ai++) {
int cur_type = TYPEOF(array_list[ai]);
if (cur_type < LGLSXP || cur_type > STRSXP)
stop("Invalid type: %d %s\n", cur_type, type2name(array_list[ai]));
if (cur_type > max_type)
max_type = cur_type;
}
switch(max_type) {
case LGLSXP: return cpp_stack_impl<LGLSXP>(array_list, along, as<LogicalVector>(fill), ovr);
case INTSXP: return cpp_stack_impl<INTSXP>(array_list, along, as<IntegerVector>(fill), ovr);
case REALSXP: return cpp_stack_impl<REALSXP>(array_list, along, as<NumericVector>(fill), ovr);
case CPLXSXP: return cpp_stack_impl<CPLXSXP>(array_list, along, as<ComplexVector>(fill), ovr);
case STRSXP: return cpp_stack_impl<STRSXP>(array_list, along, as<CharacterVector>(fill), ovr);
default: return R_NilValue; // this should not happen
}
}