NumCpp  2.8.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.0;
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.0;
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.0;
139  });
140 
141  return returnArray;
142  }
143  }
144  }
145 
146  //============================================================================
147  // Method Description:
157  template<typename dtype>
158  NdArray<std::complex<double>> gradient(const NdArray<std::complex<dtype>>& inArray, Axis inAxis = Axis::ROW)
159  {
161 
162  switch (inAxis)
163  {
164  case Axis::ROW:
165  {
166  const auto inShape = inArray.shape();
167  if (inShape.rows < 2)
168  {
169  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 row.");
170  }
171 
172  // first do the first and last rows
173  auto returnArray = NdArray<std::complex<double>>(inShape);
174  for (uint32 col = 0; col < inShape.cols; ++col)
175  {
176  returnArray(0, col) = complex_cast<double>(inArray(1, col)) - complex_cast<double>(inArray(0, col));
177  returnArray(-1, col) =
178  complex_cast<double>(inArray(-1, col)) - complex_cast<double>(inArray(-2, col));
179  }
180 
181  // then rip through the rest of the array
182  for (uint32 col = 0; col < inShape.cols; ++col)
183  {
184  for (uint32 row = 1; row < inShape.rows - 1; ++row)
185  {
186  returnArray(row, col) = (complex_cast<double>(inArray(row + 1, col)) -
187  complex_cast<double>(inArray(row - 1, col))) /
188  2.0;
189  }
190  }
191 
192  return returnArray;
193  }
194  case Axis::COL:
195  {
196  const auto inShape = inArray.shape();
197  if (inShape.cols < 2)
198  {
199  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 columns.");
200  }
201 
202  // first do the first and last columns
203  auto returnArray = NdArray<std::complex<double>>(inShape);
204  for (uint32 row = 0; row < inShape.rows; ++row)
205  {
206  returnArray(row, 0) = complex_cast<double>(inArray(row, 1)) - complex_cast<double>(inArray(row, 0));
207  returnArray(row, -1) =
208  complex_cast<double>(inArray(row, -1)) - complex_cast<double>(inArray(row, -2));
209  }
210 
211  // then rip through the rest of the array
212  for (uint32 row = 0; row < inShape.rows; ++row)
213  {
214  for (uint32 col = 1; col < inShape.cols - 1; ++col)
215  {
216  returnArray(row, col) = (complex_cast<double>(inArray(row, col + 1)) -
217  complex_cast<double>(inArray(row, col - 1))) /
218  2.0;
219  }
220  }
221 
222  return returnArray;
223  }
224  default:
225  {
226  // will return the gradient of the flattened array
227  if (inArray.size() < 2)
228  {
229  THROW_INVALID_ARGUMENT_ERROR("input array must have more than 1 element.");
230  }
231 
232  auto returnArray = NdArray<std::complex<double>>(1, inArray.size());
233  returnArray[0] = complex_cast<double>(inArray[1]) - complex_cast<double>(inArray[0]);
234  returnArray[-1] = complex_cast<double>(inArray[-1]) - complex_cast<double>(inArray[-2]);
235 
237  inArray.cbegin() + 2,
238  inArray.cend(),
239  inArray.cbegin(),
240  returnArray.begin() + 1,
241  [](const std::complex<dtype>& value1, const std::complex<dtype>& value2) -> std::complex<double>
242  { return (complex_cast<double>(value1) - complex_cast<double>(value2)) / 2.0; });
243 
244  return returnArray;
245  }
246  }
247  }
248 } // 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:4289
const_iterator cbegin() const noexcept
Definition: NdArrayCore.hpp:1221
Shape shape() const noexcept
Definition: NdArrayCore.hpp:4276
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