Skip to content

Commit 5d82bd0

Browse files
Virginia BrancatovbrancatLiang YuLiang Yugmgunter
authored andcommitted
Pybind bindings for Poly2D (#873)
* Create Poly2D .cpp and .h files * Adding Poly2D to core.cpp * Add Poly2d to source CMake * Constructor for Poly2d * add dunder methods and read-only properties * Including eval function * Add read_only functions for range and azimuth orders * Add Poly2d unit test * Add poly2d unit test to CMakefile * Evaluate values for each x, y positions in eval function * Change check in x and y size in eval dunder method * removed extra elements in vector eval and formatted * added unit test for vector input and minor formatting * Remove iterator in eval function * More Pythonic __getitem__ dunder methods * Fix bug in eval function to work with x/y arrays * updates to poly2d vector eval replaced std containers with eigen now returns (m,n) array from evaluating poly2d with (m,1) and (n,1) vectors * updates to poly1d vector eval replaced std containers with eigen now returns (m,n) array from evaluating poly1d from (m,n) input array * Revert "updates to poly1d vector eval" This reverts commit 8b760e388914a2ea304501bc8c2e5849a89783d1. * Correct typo (comma) and remove __str__ method * Correct missing & in __call__ method, remove __len__ method * Move return statement in outer loop. Replace py::list with Eigen:ArrayXXd * Replace np.ones with np.full in poly2d unit test * Use Eigen ArrayXXd for Poly2d coefficients * Remove __call__ dunder method * Add coeffs property, remove __getitem__, modify constructor * correct range/azimuth order in constructor, transpose coeffs in unit test * address PR comments and fixes eval operates over meshes instead of vectors coeff property order fixed * refactored eval vector test to eval mesh test * added grid eval, renamed mesh eval, and added docstrings * added unit test for grid eval and consolidated common code * clean up * Correct typo in Poly2d.h and add doc string to pybind implementation * Swap x/y in eval and replace counter size_t with Eigen::Index * Replace range/azimuth with x/y * Correct missing +1 in return coeffs * Change eval2d function to overrid method eval * Modify constructor and remove __getitem__ function * Change docstring of eval function * Remove x_order and y_order from constructor doc string * updated Poly2d init calls * fix x, y order parameter init * changed range/azimuth prefixes to x/y respectively * changed range/azimuth prefixes to x/y respectively * address PR comments fixed overload eval method made y & x args consistent with C++ changed range/azimuth prefixes to x/y respectively * update inits, test grid/mesh with np.polynomial * fix incorrect order * address PR comments fix numerous typos in comments fix erroneous eval arg order add coeff comparison test clang-formatt'd pybind code * Update Poly2d.eval() bindings for array input Previously, Poly2d.eval() could handle input arrays of arbitrary rank by taking input `py::array_t<double>` arguments. However, due to the limited functionality of py::array_t, the implementation was clunky and could have the unintended side-effect of modifying the input arrays. This update replaces the previous Poly2d.eval() implementation with two overloads using Eigen data structures (one for 1-D arrays, one for 2-D arrays) in order to avoid troublesome implementation issues. Co-authored-by: vbrancat <[email protected]> Co-authored-by: Liang Yu <[email protected]> Co-authored-by: Liang Yu <[email protected]> Co-authored-by: Geoffrey M Gunter <[email protected]>
1 parent 9756ef3 commit 5d82bd0

File tree

13 files changed

+463
-144
lines changed

13 files changed

+463
-144
lines changed

cxx/isce3/core/Poly2d.cpp

Lines changed: 10 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,42 +1,35 @@
1-
//
2-
// Author: Joshua Cohen
3-
// Copyright 2017
4-
//
5-
61
#include <iostream>
72
#include "Constants.h"
83
#include "Poly2d.h"
94

10-
/**
5+
/**
116
* @param[in] azi azimuth or y value
127
* @param[in] rng range or x value*/
138
double isce3::core::Poly2d::
14-
eval(double azi, double rng) const {
9+
eval(double y, double x) const {
1510

16-
double xval = (rng - rangeMean) / rangeNorm;
17-
double yval = (azi - azimuthMean) / azimuthNorm;
11+
double xval = (x - xMean) / xNorm;
12+
double yval = (y - yMean) / yNorm;
1813

1914
double scalex;
2015
double scaley = 1.;
2116
double val = 0.;
22-
for (int i=0; i<=azimuthOrder; i++,scaley*=yval) {
17+
for (int i=0; i<=yOrder; i++,scaley*=yval) {
2318
scalex = 1.;
24-
for (int j=0; j<=rangeOrder; j++,scalex*=xval) {
25-
val += scalex * scaley * coeffs[IDX1D(i,j,rangeOrder+1)];
19+
for (int j=0; j<=xOrder; j++,scalex*=xval) {
20+
val += scalex * scaley * coeffs[IDX1D(i,j,xOrder+1)];
2621
}
2722
}
2823
return val;
2924
}
3025

3126
void isce3::core::Poly2d::
3227
printPoly() const {
33-
std::cout << "Polynomial Order: " << azimuthOrder << " - by - " << rangeOrder << std::endl;
34-
for (int i=0; i<=azimuthOrder; i++) {
35-
for (int j=0; j<=rangeOrder; j++) {
28+
std::cout << "Polynomial Order: " << yOrder << " - by - " << xOrder << std::endl;
29+
for (int i=0; i<=yOrder; i++) {
30+
for (int j=0; j<=xOrder; j++) {
3631
std::cout << getCoeff(i,j) << " ";
3732
}
3833
std::cout << std::endl;
3934
}
4035
}
41-
42-
// end of file

cxx/isce3/core/Poly2d.h

Lines changed: 51 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,3 @@
1-
//
2-
// Author: Joshua Cohen
3-
// Copyright 2017
4-
//
5-
61
#pragma once
72

83
#include "forward.h"
@@ -12,60 +7,60 @@
127
#include <vector>
138
#include "Constants.h"
149

15-
/** Data structure for representing 1D polynomials
10+
/** Data structure for representing 2D polynomials
1611
*
17-
* Poly2D is function of the form
12+
* Poly2D is function of the form
1813
* \f[
19-
* f\left( y, x \right) = \sum_{i=0}^{N_y} \sum_{j=0}^{N_x} a_{ij} \cdot \left( \frac{y-\mu_y}{\sigma_y} \right)^i
20-
\cdot \left( \frac{x-\mu_x}{\sigma_x} \right)^j
14+
* f\left( y, x \right) = \sum_{i=0}^{N_y} \sum_{j=0}^{N_x} a_{ij} \cdot \left( \frac{y-\mu_y}{\sigma_y} \right)^i
15+
\cdot \left( \frac{x-\mu_x}{\sigma_x} \right)^j
2116
* \f]
2217
*
23-
* where \f$a_ij\f$ represents the coefficients, \f$\mu_x\f$ and \f$\mu_y\f$ represent the means and
18+
* where \f$a_ij\f$ represents the coefficients, \f$\mu_x\f$ and \f$\mu_y\f$ represent the means and
2419
* \f$\sigma_x\f$ and \f$\sigma_y\f$ represent the norms*/
2520
class isce3::core::Poly2d {
2621
public:
2722

2823
/** Order of polynomial in range or x*/
29-
int rangeOrder;
24+
int xOrder;
3025
/** Order of polynomial in azimuth or y*/
31-
int azimuthOrder;
26+
int yOrder;
3227
/** Mean in range or x direction*/
33-
double rangeMean;
28+
double xMean;
3429
/** Mean in azimuth or y direction*/
35-
double azimuthMean;
30+
double yMean;
3631
/** Norm in range or x direction*/
37-
double rangeNorm;
32+
double xNorm;
3833
/** Norm in azimuth or y direction*/
39-
double azimuthNorm;
34+
double yNorm;
4035
/**Linearized vector of coefficients in row-major format*/
4136
std::vector<double> coeffs;
4237

4338
/** Simple constructor
4439
*
45-
* @param[in] ro Range Order
46-
* @param[in] ao Azimuth Order
47-
* @param[in] rm Range Mean
48-
* @param[in] am Azimuth Mean
49-
* @param[in] rn Range Norm
50-
* @param[in] an Azimuth Norm*/
51-
Poly2d(int ro, int ao, double rm, double am, double rn, double an) : rangeOrder(ro),
52-
azimuthOrder(ao),
53-
rangeMean(rm),
54-
azimuthMean(am),
55-
rangeNorm(rn),
56-
azimuthNorm(an),
57-
coeffs((ro+1)*(ao+1))
40+
* @param[in] xo x/Range Order
41+
* @param[in] yo y/Azimuth Order
42+
* @param[in] xm x/Range Mean
43+
* @param[in] ym y/Azimuth Mean
44+
* @param[in] xn x/Range Norm
45+
* @param[in] yn y/Azimuth Norm*/
46+
Poly2d(int xo, int yo, double xm, double ym, double xn, double yn) : xOrder(xo),
47+
yOrder(yo),
48+
xMean(xm),
49+
yMean(ym),
50+
xNorm(xn),
51+
yNorm(yn),
52+
coeffs((xo+1)*(yo+1))
5853
{}
5954

6055
/** Empty constructor*/
6156
Poly2d() : Poly2d(-1,-1,0.,0.,1.,1.) {}
6257

63-
/** Copy constructor
58+
/** Copy constructor
6459
*
6560
* @param[in] p Poly2D object*/
66-
Poly2d(const Poly2d &p) : rangeOrder(p.rangeOrder), azimuthOrder(p.azimuthOrder),
67-
rangeMean(p.rangeMean), azimuthMean(p.azimuthMean),
68-
rangeNorm(p.rangeNorm), azimuthNorm(p.azimuthNorm),
61+
Poly2d(const Poly2d &p) : xOrder(p.xOrder), yOrder(p.yOrder),
62+
xMean(p.xMean), yMean(p.yMean),
63+
xNorm(p.xNorm), yNorm(p.yNorm),
6964
coeffs(p.coeffs) {}
7065

7166
/** Assignment operator*/
@@ -77,60 +72,60 @@ class isce3::core::Poly2d {
7772
/**Get coefficient by indices*/
7873
inline double getCoeff(int row, int col) const;
7974

80-
/**Evaluate polynomial at given y,x*/
81-
double eval(double azi, double rng) const;
75+
/**Evaluate polynomial at given y/azimuth/row ,x/range/col*/
76+
double eval(double y, double x) const;
8277

8378
/**Printing for debugging*/
8479
void printPoly() const;
8580
};
8681

8782
isce3::core::Poly2d & isce3::core::Poly2d::
8883
operator=(const Poly2d &rhs) {
89-
rangeOrder = rhs.rangeOrder;
90-
azimuthOrder = rhs.azimuthOrder;
91-
rangeMean = rhs.rangeMean;
92-
azimuthMean = rhs.azimuthMean;
93-
rangeNorm = rhs.rangeNorm;
94-
azimuthNorm = rhs.azimuthNorm;
84+
xOrder = rhs.xOrder;
85+
yOrder = rhs.yOrder;
86+
xMean = rhs.xMean;
87+
yMean = rhs.yMean;
88+
xNorm = rhs.xNorm;
89+
yNorm = rhs.yNorm;
9590
coeffs = rhs.coeffs;
9691
return *this;
9792
}
9893

9994
/**
100-
* @param[in] row azimuth/y index
95+
* @param[in] row azimuth/y index
10196
* @param[in] col range/x index
10297
* @param[in] val Coefficient value*/
10398
void isce3::core::Poly2d::
10499
setCoeff(int row, int col, double val) {
105-
if ((row < 0) || (row > azimuthOrder)) {
106-
std::string errstr = "Poly2d::setCoeff - Trying to set coefficient for row " +
107-
std::to_string(row+1) + " out of " +
108-
std::to_string(azimuthOrder+1);
100+
if ((row < 0) || (row > yOrder)) {
101+
std::string errstr = "Poly2d::setCoeff - Trying to set coefficient for row " +
102+
std::to_string(row+1) + " out of " +
103+
std::to_string(yOrder+1);
109104
throw std::out_of_range(errstr);
110105
}
111-
if ((col < 0) || (col > rangeOrder)) {
106+
if ((col < 0) || (col > xOrder)) {
112107
std::string errstr = "Poly2d::setCoeff - Trying to set coefficient for col " +
113-
std::to_string(col+1) + " out of " + std::to_string(rangeOrder+1);
108+
std::to_string(col+1) + " out of " + std::to_string(xOrder+1);
114109
throw std::out_of_range(errstr);
115110
}
116-
coeffs[IDX1D(row,col,rangeOrder+1)] = val;
111+
coeffs[IDX1D(row,col,xOrder+1)] = val;
117112
}
118113

119114
/**
120115
* @param[in] row azimuth/y index
121116
* @param[in] col range/x index*/
122117
double isce3::core::Poly2d::
123118
getCoeff(int row, int col) const {
124-
if ((row < 0) || (row > azimuthOrder)) {
119+
if ((row < 0) || (row > yOrder)) {
125120
std::string errstr = "Poly2d::getCoeff - Trying to get coefficient for row " +
126-
std::to_string(row+1) + " out of " +
127-
std::to_string(azimuthOrder+1);
121+
std::to_string(row+1) + " out of " +
122+
std::to_string(yOrder+1);
128123
throw std::out_of_range(errstr);
129124
}
130-
if ((col < 0) || (col > rangeOrder)) {
131-
std::string errstr = "Poly2d::getCoeff - Trying to get coefficient for col " +
132-
std::to_string(col+1) + " out of " + std::to_string(rangeOrder+1);
125+
if ((col < 0) || (col > xOrder)) {
126+
std::string errstr = "Poly2d::getCoeff - Trying to get coefficient for col " +
127+
std::to_string(col+1) + " out of " + std::to_string(xOrder+1);
133128
throw std::out_of_range(errstr);
134129
}
135-
return coeffs[IDX1D(row,col,rangeOrder+1)];
130+
return coeffs[IDX1D(row,col,xOrder+1)];
136131
}

cxx/isce3/core/Serialization.h

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -253,12 +253,12 @@ inline void loadFromH5(isce3::io::IGroup & group, Poly2d & poly, std::string nam
253253
isce3::io::loadFromH5(group, name, poly.coeffs);
254254

255255
// Set other polynomial properties
256-
poly.rangeOrder = poly.coeffs.size() - 1;
257-
poly.azimuthOrder = 0;
258-
poly.rangeMean = 0.0;
259-
poly.azimuthMean = 0.0;
260-
poly.rangeNorm = 1.0;
261-
poly.azimuthNorm = 1.0;
256+
poly.xOrder = poly.coeffs.size() - 1;
257+
poly.yOrder = 0;
258+
poly.xMean = 0.0;
259+
poly.yMean = 0.0;
260+
poly.xNorm = 1.0;
261+
poly.yNorm = 1.0;
262262
}
263263

264264
// ------------------------------------------------------------------------

cxx/isce3/cuda/core/gpuPoly2d.cu

Lines changed: 22 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,3 @@
1-
//
2-
// Author: Liang Yu
3-
// Copyright 2018
4-
//
5-
// NOTE: gpuOrbit used as template
6-
71
#include <cuda_runtime.h>
82
#include <vector>
93
#include <isce3/core/Constants.h>
@@ -16,72 +10,71 @@ using isce3::core::Poly2d;
1610
using std::vector;
1711

1812

19-
// Advanced "copy" constructor to handle deep-copying of Poly2d data (only callable by host). Owner
20-
// member variable indicates that only the host-side copy of the gpuPoly2d can handle freeing the
13+
// Advanced "copy" constructor to handle deep-copying of Poly2d data (only callable by host). Owner
14+
// member variable indicates that only the host-side copy of the gpuPoly2d can handle freeing the
2115
// memory (device-side copy constructor for gpuPoly2d sets owner to false)
2216
__host__ gpuPoly2d::gpuPoly2d(const Poly2d &poly) :
23-
rangeOrder(poly.rangeOrder),
24-
azimuthOrder(poly.azimuthOrder),
25-
rangeMean(poly.rangeMean),
26-
azimuthMean(poly.azimuthMean),
27-
rangeNorm(poly.rangeNorm),
28-
azimuthNorm(poly.azimuthNorm),
17+
xOrder(poly.xOrder),
18+
yOrder(poly.yOrder),
19+
xMean(poly.xMean),
20+
yMean(poly.yMean),
21+
xNorm(poly.xNorm),
22+
yNorm(poly.yNorm),
2923
owner(true)
3024
{
31-
25+
3226
const int n_coeffs = poly.coeffs.size();
3327

3428
// Malloc device-side memory (this API is host-side only)
3529
cudaMalloc(&coeffs, n_coeffs*sizeof(double));
3630

37-
// Copy OrPoly2d data to device-side memory and keep device pointer in gpuOrPoly2d object. Device-side
31+
// Copy OrPoly2d data to device-side memory and keep device pointer in gpuOrPoly2d object. Device-side
3832
// copy constructor simply shallow-copies the device pointers when called
3933
cudaMemcpy(coeffs, poly.coeffs.data(), n_coeffs*sizeof(double), cudaMemcpyHostToDevice);
4034
}
4135

4236

43-
// Both the host-side and device-side copies of the gpuPoly2d will call the destructor, so we have to
44-
// implement a way of having an arbitrary copy on host OR device determine when to free the memory
37+
// Both the host-side and device-side copies of the gpuPoly2d will call the destructor, so we have to
38+
// implement a way of having an arbitrary copy on host OR device determine when to free the memory
4539
// (since only the original host-side copy should free)
4640
gpuPoly2d::~gpuPoly2d() {
4741
if (owner) {
4842
cudaFree(coeffs);
4943
}
5044
}
5145

52-
__device__ double gpuPoly2d::eval(double azi, double rng) const {
46+
__device__ double gpuPoly2d::eval(double y, double x) const {
5347

54-
double xval = (rng - rangeMean) / rangeNorm;
55-
double yval = (azi - azimuthMean) / azimuthNorm;
48+
double xval = (x - xMean) / xNorm;
49+
double yval = (y - yMean) / yNorm;
5650

5751
double scalex;
5852
double scaley = 1.;
5953
double val = 0.;
60-
for (int i=0; i<=azimuthOrder; i++,scaley*=yval) {
54+
for (int i=0; i<=yOrder; i++,scaley*=yval) {
6155
scalex = 1.;
62-
for (int j=0; j<=rangeOrder; j++,scalex*=xval) {
63-
val += scalex * scaley * coeffs[IDX1D(i,j,rangeOrder+1)];
56+
for (int j=0; j<=xOrder; j++,scalex*=xval) {
57+
val += scalex * scaley * coeffs[IDX1D(i,j,xOrder+1)];
6458
}
6559
}
6660

6761
return val;
6862
}
6963

70-
__global__ void eval_d(gpuPoly2d p, double azi, double rng, double *val)
64+
__global__ void eval_d(gpuPoly2d p, double y, double x, double *val)
7165
{
72-
*val = p.eval(azi, rng);
66+
*val = p.eval(y, x);
7367
}
7468

75-
__host__ double gpuPoly2d::eval_h(double azi, double rng)
69+
__host__ double gpuPoly2d::eval_h(double y, double x)
7670
{
7771
double *val_d;
7872
double val_h;
7973
// use unified memory?
8074
cudaMalloc((double**)&val_d, sizeof(double));
8175
dim3 grid(1), block(1);
82-
eval_d<<<grid,block>>>(*this, azi, rng, val_d);
76+
eval_d<<<grid,block>>>(*this, y, x, val_d);
8377
cudaMemcpy(&val_h, val_d, sizeof(double), cudaMemcpyDeviceToHost);
8478
cudaFree(val_d);
8579
return val_h;
8680
}
87-

0 commit comments

Comments
 (0)