NumCpp  2.9.0
A Templatized Header Only C++ Implementation of the Python NumPy Library
gradient.hpp
Go to the documentation of this file.
1 #pragma once
29 
30 #include <complex>
31 #include <string>
32 
37 #include "NumCpp/Core/Shape.hpp"
38 #include "NumCpp/Core/Types.hpp"
39 #include "NumCpp/NdArray.hpp"
40 
41 namespace nc
42 {
43  //============================================================================
44  // Method Description:
54  template<typename dtype>
56  {
58 
59  switch (inAxis)
60  {
61  case Axis::ROW:
62  {
63  const auto inShape = inArray.shape();
64  if (inShape.rows < 2)
65  {
66  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 row.");
67  }
68 
69  // first do the first and last rows
70  auto returnArray = NdArray<double>(inShape);
71  for (uint32 col = 0; col < inShape.cols; ++col)
72  {
73  returnArray(0, col) = static_cast<double>(inArray(1, col)) - static_cast<double>(inArray(0, col));
74  returnArray(-1, col) =
75  static_cast<double>(inArray(-1, col)) - static_cast<double>(inArray(-2, col));
76  }
77 
78  // then rip through the rest of the array
79  for (uint32 col = 0; col < inShape.cols; ++col)
80  {
81  for (uint32 row = 1; row < inShape.rows - 1; ++row)
82  {
83  returnArray(row, col) =
84  (static_cast<double>(inArray(row + 1, col)) - static_cast<double>(inArray(row - 1, col))) /
85  2.;
86  }
87  }
88 
89  return returnArray;
90  }
91  case Axis::COL:
92  {
93  const auto inShape = inArray.shape();
94  if (inShape.cols < 2)
95  {
96  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 columns.");
97  }
98 
99  // first do the first and last columns
100  auto returnArray = NdArray<double>(inShape);
101  for (uint32 row = 0; row < inShape.rows; ++row)
102  {
103  returnArray(row, 0) = static_cast<double>(inArray(row, 1)) - static_cast<double>(inArray(row, 0));
104  returnArray(row, -1) =
105  static_cast<double>(inArray(row, -1)) - static_cast<double>(inArray(row, -2));
106  }
107 
108  // then rip through the rest of the array
109  for (uint32 row = 0; row < inShape.rows; ++row)
110  {
111  for (uint32 col = 1; col < inShape.cols - 1; ++col)
112  {
113  returnArray(row, col) =
114  (static_cast<double>(inArray(row, col + 1)) - static_cast<double>(inArray(row, col - 1))) /
115  2.;
116  }
117  }
118 
119  return returnArray;
120  }
121  default:
122  {
123  // will return the gradient of the flattened array
124  if (inArray.size() < 2)
125  {
126  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 element.");
127  }
128 
129  auto returnArray = NdArray<double>(1, inArray.size());
130  returnArray[0] = static_cast<double>(inArray[1]) - static_cast<double>(inArray[0]);
131  returnArray[-1] = static_cast<double>(inArray[-1]) - static_cast<double>(inArray[-2]);
132 
133  stl_algorithms::transform(inArray.cbegin() + 2,
134  inArray.cend(),
135  inArray.cbegin(),
136  returnArray.begin() + 1,
137  [](dtype value1, dtype value2) -> double
138  { return (static_cast<double>(value1) - static_cast<double>(value2)) / 2.; });
139 
140  return returnArray;
141  }
142  }
143  }
144 
145  //============================================================================
146  // Method Description:
156  template<typename dtype>
157  NdArray<std::complex<double>> gradient(const NdArray<std::complex<dtype>>& inArray, Axis inAxis = Axis::ROW)
158  {
160 
161  switch (inAxis)
162  {
163  case Axis::ROW:
164  {
165  const auto inShape = inArray.shape();
166  if (inShape.rows < 2)
167  {
168  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 row.");
169  }
170 
171  // first do the first and last rows
172  auto returnArray = NdArray<std::complex<double>>(inShape);
173  for (uint32 col = 0; col < inShape.cols; ++col)
174  {
175  returnArray(0, col) = complex_cast<double>(inArray(1, col)) - complex_cast<double>(inArray(0, col));
176  returnArray(-1, col) =
177  complex_cast<double>(inArray(-1, col)) - complex_cast<double>(inArray(-2, col));
178  }
179 
180  // then rip through the rest of the array
181  for (uint32 col = 0; col < inShape.cols; ++col)
182  {
183  for (uint32 row = 1; row < inShape.rows - 1; ++row)
184  {
185  returnArray(row, col) = (complex_cast<double>(inArray(row + 1, col)) -
186  complex_cast<double>(inArray(row - 1, col))) /
187  2.;
188  }
189  }
190 
191  return returnArray;
192  }
193  case Axis::COL:
194  {
195  const auto inShape = inArray.shape();
196  if (inShape.cols < 2)
197  {
198  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 columns.");
199  }
200 
201  // first do the first and last columns
202  auto returnArray = NdArray<std::complex<double>>(inShape);
203  for (uint32 row = 0; row < inShape.rows; ++row)
204  {
205  returnArray(row, 0) = complex_cast<double>(inArray(row, 1)) - complex_cast<double>(inArray(row, 0));
206  returnArray(row, -1) =
207  complex_cast<double>(inArray(row, -1)) - complex_cast<double>(inArray(row, -2));
208  }
209 
210  // then rip through the rest of the array
211  for (uint32 row = 0; row < inShape.rows; ++row)
212  {
213  for (uint32 col = 1; col < inShape.cols - 1; ++col)
214  {
215  returnArray(row, col) = (complex_cast<double>(inArray(row, col + 1)) -
216  complex_cast<double>(inArray(row, col - 1))) /
217  2.;
218  }
219  }
220 
221  return returnArray;
222  }
223  default:
224  {
225  // will return the gradient of the flattened array
226  if (inArray.size() < 2)
227  {
228  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 element.");
229  }
230 
231  auto returnArray = NdArray<std::complex<double>>(1, inArray.size());
232  returnArray[0] = complex_cast<double>(inArray[1]) - complex_cast<double>(inArray[0]);
233  returnArray[-1] = complex_cast<double>(inArray[-1]) - complex_cast<double>(inArray[-2]);
234 
236  inArray.cbegin() + 2,
237  inArray.cend(),
238  inArray.cbegin(),
239  returnArray.begin() + 1,
240  [](const std::complex<dtype>& value1, const std::complex<dtype>& value2) -> std::complex<double>
241  { return (complex_cast<double>(value1) - complex_cast<double>(value2)) / 2.; });
242 
243  return returnArray;
244  }
245  }
246  }
247 } // namespace nc
#define THROW_INVALID_ARGUMENT_ERROR(msg)
Definition: Error.hpp:36
#define STATIC_ASSERT_ARITHMETIC(dtype)
Definition: StaticAsserts.hpp:37
size_type size() const noexcept
Definition: NdArrayCore.hpp:4105
const_iterator cbegin() const noexcept
Definition: NdArrayCore.hpp:1221
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4092
const_iterator cend() const noexcept
Definition: NdArrayCore.hpp:1529
iterator begin() noexcept
Definition: NdArrayCore.hpp:1171
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction)
Definition: StlAlgorithms.hpp:784
Definition: Coordinate.hpp:45
Axis
Enum To describe an axis.
Definition: Types.hpp:47
NdArray< double > gradient(const NdArray< dtype > &inArray, Axis inAxis=Axis::ROW)
Definition: gradient.hpp:55
std::uint32_t uint32
Definition: Types.hpp:40