Alexandria 2.31.0
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#ifdef __APPLE__
150 auto compare = [](void* context, const void* a, const void* b) {
151#else
152 auto compare = [](const void* a, const void* b, void* context) {
153#endif
154 const T* arow = static_cast<const T*>(a);
155 const T* brow = static_cast<const T*>(b);
156 const auto* attr_idx_ptr = static_cast<const std::vector<size_t>*>(context);
157 for (auto& idx : *attr_idx_ptr) {
158 if (*(arow + idx) < *(brow + idx)) {
159 return -1;
160 } else if (*(arow + idx) > *(brow + idx)) {
161 return 1;
162 }
163 }
164 return 0;
165 };
166
167#ifdef __APPLE__
168 qsort_r(&(*array.begin()), count, row_size, &attr_idx, compare);
169#else
170 qsort_r(&(*array.begin()), count, row_size, compare, &attr_idx);
171#endif
172}
173
174} // namespace NdArray
175} // namespace Euclid
176
177#endif