Alexandria 2.31.4
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
Operations.icpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2021 Euclid Science Ground Segment
3 *
4 * This library is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the Free
6 * Software Foundation; either version 3.0 of the License, or (at your option)
7 * any later version.
8 *
9 * This library is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU Lesser General Public License
15 * along with this library; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
19#ifdef NDARRAY_OPS_IMPL
20
21#include <cstdlib>
22
23namespace Euclid {
24namespace NdArray {
25
26template <typename T>
27T sum(const NdArray<T>& array) {
28 return std::accumulate(array.begin(), array.end(), T{});
29}
30
31template <typename T>
32NdArray<T> sum(const NdArray<T>& array, int axis) {
33 if (axis < 0) {
34 axis = array.shape().size() + axis;
35 }
36 assert(axis >= 0);
37 if (static_cast<std::size_t>(axis) >= array.shape().size()) {
38 throw std::out_of_range("Invalid axis");
39 }
40 // Single axis returns a scalar
41 if (array.shape().size() == 1) {
42 NdArray<T> output{{1}};
43 output.at(0) = sum(array);
44 return output;
45 }
46 // Allocate output
47 auto output_shape = array.shape();
48 auto axis_size = output_shape[axis];
49
50 output_shape.erase(output_shape.begin() + axis);
51 NdArray<T> output(output_shape);
52 // Sum along axis
53 for (std::size_t out_i = 0; out_i < output.size(); ++out_i) {
54 auto out_coordinates = unravel_index(out_i, output_shape);
55 std::vector<std::size_t> in_coordinates(out_coordinates);
56 in_coordinates.insert(in_coordinates.begin() + axis, 0ul);
57
58 auto& output_val = output.at(out_coordinates);
59 for (std::size_t axis_i = 0; axis_i < axis_size; ++axis_i) {
60 in_coordinates[axis] = axis_i;
61 output_val += array.at(in_coordinates);
62 }
63 }
64 // Return
65 return output;
66}
67
68template <typename T, typename Iterator>
69NdArray<T> trapz(const NdArray<T>& array, const Iterator kbegin, const Iterator kend, int axis) {
70 if (axis < 0) {
71 axis = array.shape().size() + axis;
72 }
73 assert(axis >= 0);
74 if (static_cast<std::size_t>(axis) >= array.shape().size()) {
75 throw std::out_of_range("Invalid axis");
76 }
77 // Allocate output
78 auto output_shape = array.shape();
79 auto axis_size = output_shape[axis];
80 // Verify knots size
81 if (static_cast<std::size_t>(kend - kbegin) != axis_size) {
82 throw std::length_error("Integration axis value does not match the size of the array axis");
83 }
84 // Single axis returns a scalar
85 bool output_scalar = (array.shape().size() == 1);
86 if (output_scalar) {
87 output_shape.clear();
88 output_shape.emplace_back(1);
89 } else {
90 output_shape.erase(output_shape.begin() + axis);
91 }
92 NdArray<T> output(output_shape);
93 // Integrate along axis
94 for (std::size_t out_i = 0; out_i < output.size(); ++out_i) {
95 auto out_coordinates = unravel_index(out_i, output_shape);
96 std::vector<std::size_t> in_coordinates(out_coordinates);
97 if (!output_scalar) {
98 in_coordinates.insert(in_coordinates.begin() + axis, 0ul);
99 }
100
101 auto& output_val = output.at(out_coordinates);
102 auto ki = kbegin + 1;
103 for (std::size_t axis_i = 1; axis_i < axis_size; ++axis_i, ++ki) {
104 auto dx = (*ki) - *(ki - 1);
105 in_coordinates[axis] = axis_i;
106 auto b = array.at(in_coordinates);
107 --in_coordinates[axis];
108 auto a = array.at(in_coordinates);
109 output_val += ((a + b) * dx) / 2;
110 }
111 }
112 // Return
113 return output;
114}
115
116template <typename T>
117std::vector<std::size_t> argmax(const NdArray<T>& array) {
118 auto max_iter = std::max_element(array.begin(), array.end());
119 auto max_pos = max_iter - array.begin();
120 return unravel_index(max_pos, array.shape());
121}
122
123template <typename T>
124std::vector<std::size_t> argmin(const NdArray<T>& array) {
125 auto max_iter = std::min_element(array.begin(), array.end());
126 auto max_pos = max_iter - array.begin();
127 return unravel_index(max_pos, array.shape());
128}
129
130template <typename T>
131void sort(NdArray<T>& array, const std::vector<std::string>& attrs) {
132 if (array.shape().size() != 2) {
133 throw std::invalid_argument("Only 2D arrays can be sorted by attribute names");
134 }
135
136 std::vector<size_t> attr_idx(attrs.size());
137 const auto& array_attrs = array.attributes();
138 for (size_t i = 0; i < attrs.size(); ++i) {
139 auto pos = std::find(array_attrs.begin(), array_attrs.end(), attrs[i]);
140 if (pos == array_attrs.end()) {
141 throw std::out_of_range(attrs[i]);
142 }
143 attr_idx[i] = pos - array_attrs.begin();
144 }
145
146 size_t count = array.shape()[0];
147 size_t row_size = array.shape()[1] * sizeof(T);
148
149
150#ifdef __APPLE__
151 auto compare = [](void* context, const void* a, const void* b) {
152#else
153 auto compare = [](const void* a, const void* b, void* context) {
154#endif
155 const T* arow = static_cast<const T*>(a);
156 const T* brow = static_cast<const T*>(b);
157 const auto* attr_idx_ptr = static_cast<const std::vector<size_t>*>(context);
158 for (auto& idx : *attr_idx_ptr) {
159 if (*(arow + idx) < *(brow + idx)) {
160 return -1;
161 } else if (*(arow + idx) > *(brow + idx)) {
162 return 1;
163 }
164 }
165 return 0;
166 };
167
168#ifdef __APPLE__
169 qsort_r(&(*array.begin()), count, row_size, &attr_idx, compare);
170#else
171 qsort_r(&(*array.begin()), count, row_size, compare, &attr_idx);
172#endif
173}
174
175template <typename T1, typename T2>
176bool sameShape(const NdArray<T1>& array_1, const NdArray<T2>& array_2) {
177 auto shape_1 = array_1.shape();
178 auto shape_2 = array_2.shape();
179 if (shape_1.size()!=shape_2.size()) {
180 return false;
181 }
182
183 for (size_t index = 0; index < shape_1.size(); ++index) {
184 if (shape_1[index] != shape_2[index]) {
185 return false;
186 }
187 }
188
189 return true;
190}
191
192
193template <typename T1, typename T2>
194NdArray<T1> multiplyElements(const NdArray<T1>& array_1, const NdArray<T2>& array_2) {
195 static_assert(((std::is_integral<T1>::value && std::is_integral<T2>::value && sizeof(T1) >= sizeof(T2)) ||
196 (std::is_floating_point<T1>::value && std::is_floating_point<T2>::value && sizeof(T1) >= sizeof(T2)) ||
197 (std::is_floating_point<T1>::value && std::is_integral<T2>::value)),
198 "Provided NdArray base Types cannot be multiplyed or stored in the first NdArray base type"
199 );
200 if (!sameShape(array_1, array_2)) {
201 throw std::invalid_argument("Provided NdArrays do not have the same shape.");
202 }
203
204 auto result = array_1.copy();
205 auto iter_2 = array_2.begin();
206 for (auto iter_res = result.begin(); iter_res != result.end(); ++iter_res) {
207 (*iter_res)*=(*iter_2);
208 ++iter_2;
209 }
210
211 return result;
212}
213
214template <typename T1, typename T2>
215NdArray<T1> addElements(const NdArray<T1>& array_1, const NdArray<T2>& array_2) {
216 static_assert(((std::is_integral<T1>::value && std::is_integral<T2>::value && sizeof(T1) >= sizeof(T2)) ||
217 (std::is_floating_point<T1>::value && std::is_floating_point<T2>::value && sizeof(T1) >= sizeof(T2)) ||
218 (std::is_floating_point<T1>::value && std::is_integral<T2>::value)),
219 "Provided NdArray base Types cannot be multiplyed or stored in the first NdArray base type"
220 );
221 if (!sameShape(array_1, array_2)) {
222 throw std::invalid_argument("Provided NdArrays do not have the same shape.");
223 }
224
225 auto result = array_1.copy();
226 auto iter_2 = array_2.begin();
227 for (auto iter_res = result.begin(); iter_res != result.end(); ++iter_res) {
228 (*iter_res)+=(*iter_2);
229 ++iter_2;
230 }
231
232 return result;
233}
234
235} // namespace NdArray
236} // namespace Euclid
237
238#endif
T accumulate(T... args)
T count(T... args)
T find(T... args)
T max_element(T... args)
T min_element(T... args)
std::vector< std::size_t > argmax(const NdArray< T > &array)
NdArray< T1 > multiplyElements(const NdArray< T1 > &array_1, const NdArray< T2 > &array_2)
bool sameShape(const NdArray< T1 > &array_1, const NdArray< T2 > &array_2)
NdArray< T1 > addElements(const NdArray< T1 > &array_1, const NdArray< T2 > &array_2)
void sort(NdArray< T > &array, const std::vector< std::string > &attrs)
std::vector< std::size_t > argmin(const NdArray< T > &array)
T sum(const NdArray< T > &array)
std::vector< std::size_t > unravel_index(std::size_t index, const std::vector< std::size_t > &shape)
NdArray< T > trapz(const NdArray< T > &array, const Iterator kbegin, const Iterator kend, int axis)
T size(T... args)