NumCpp  2.12.1
A Templatized Header Only C++ Implementation of the Python NumPy Library
gradient.hpp
Go to the documentation of this file.
1
28#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
41namespace 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:37
#define STATIC_ASSERT_ARITHMETIC(dtype)
Definition: StaticAsserts.hpp:39
size_type size() const noexcept
Definition: NdArrayCore.hpp:4524
const_iterator cbegin() const noexcept
Definition: NdArrayCore.hpp:1365
const Shape & shape() const noexcept
Definition: NdArrayCore.hpp:4511
const_iterator cend() const noexcept
Definition: NdArrayCore.hpp:1673
iterator begin() noexcept
Definition: NdArrayCore.hpp:1315
OutputIt transform(InputIt first, InputIt last, OutputIt destination, UnaryOperation unaryFunction)
Definition: StlAlgorithms.hpp:775
Definition: Cartesian.hpp:40
Axis
Enum To describe an axis.
Definition: Enums.hpp:36
NdArray< double > gradient(const NdArray< dtype > &inArray, Axis inAxis=Axis::ROW)
Definition: gradient.hpp:55
std::uint32_t uint32
Definition: Types.hpp:40