Alexandria  2.22.0
Please provide a description of the project.
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 
23 namespace Euclid {
24 namespace NdArray {
25 
26 template <typename T>
27 T sum(const NdArray<T>& array) {
28  return std::accumulate(array.begin(), array.end(), T{});
29 }
30 
31 template <typename T>
32 NdArray<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 
68 template <typename T, typename Iterator>
69 NdArray<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 
116 template <typename T>
117 std::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 
123 template <typename T>
124 std::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 
130 template <typename T>
131 void 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 == 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  qsort_r(
150  &(*array.begin()), count, row_size,
151  [](const void* a, const void* b, void* context) {
152  const T* arow = static_cast<const T*>(a);
153  const T* brow = static_cast<const T*>(b);
154  const auto* attr_idx_ptr = static_cast<const std::vector<size_t>*>(context);
155  for (auto& idx : *attr_idx_ptr) {
156  if (*(arow + idx) < *(brow + idx)) {
157  return -1;
158  } else if (*(arow + idx) > *(brow + idx)) {
159  return 1;
160  }
161  }
162  return 0;
163  },
164  &attr_idx);
165 }
166 
167 } // namespace NdArray
168 } // namespace Euclid
169 
170 #endif