Skip to content

Fix barycentric interpolation on Windows #18

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 6 commits into
base: poly_refactor
Choose a base branch
from
+52 −50
Open

Fix barycentric interpolation on Windows #18

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3304a55
Use fuse types to address mix of incoming int types across platforms
ConnectedSystems Jun 30, 2021
9b297bc
specify np.array internally
ConnectedSystems Jun 30, 2021
cfe89af
specify numpy array dtype
ConnectedSystems Jun 30, 2021
6866674
cache index lookup
ConnectedSystems Jun 30, 2021
aba461a
Disable broken test (at least temporarily)
ConnectedSystems Jun 30, 2021
249a168
Merge branch 'poly_refactor' into fix-windows-barycentric
ConnectedSystems Jul 1, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions pyapprox/barycentric_interpolation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,14 +183,14 @@ def multivariate_hierarchical_barycentric_lagrange_interpolation(
try:
from pyapprox.cython.barycentric_interpolation import \
multivariate_hierarchical_barycentric_lagrange_interpolation_pyx

result = \
multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
x, fn_vals, active_dims,
active_abscissa_indices_1d.astype(np.int_),
num_abscissa_1d.astype(np.int_),
num_active_abscissa_1d.astype(np.int_),
shifts.astype(np.int_), abscissa_and_weights)
x, fn_vals, active_dims.astype(np.int64),
active_abscissa_indices_1d.astype(np.int64),
num_abscissa_1d.astype(np.int64),
num_active_abscissa_1d.astype(np.int64),
shifts.astype(np.int64), abscissa_and_weights)
if np.any(np.isnan(result)):
raise ValueError('Error values not finite')

Expand Down
79 changes: 39 additions & 40 deletions pyapprox/cython/barycentric_interpolation.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,12 +1,18 @@
# cython: infer_types=True
cimport cython

cimport numpy as np
import numpy as np


ctypedef np.double_t double_t
ctypedef np.int_t int_t
ctypedef np.int64_t int64_t
ctypedef fused int_t:
cython.integral
int
long
long long
np.int64_t
np.int32_t
np.int_t


@cython.cdivision(True) # Deactivate division by zero checking
Expand All @@ -22,7 +28,7 @@ cpdef compute_barycentric_weights_1d_pyx(np.ndarray[double_t] samples, double C_

double[:,:] weights_view = weights
double[:] samples_view = samples

weights_view[0,0] = 1.
for jj in range(1, num_samples):
for kk in range(jj):
Expand All @@ -43,45 +49,43 @@ cpdef compute_barycentric_weights_1d_pyx(np.ndarray[double_t] samples, double C_
@cython.boundscheck(False) # Deactivate bounds checking
@cython.wraparound(False) # Deactivate negative indexing.
cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
double[:,:] x, double[:,:] fn_vals, np.ndarray[int64_t] active_dims,
double[:,:] x, double[:,:] fn_vals, int_t[:] active_dims,
int_t[:,:] active_abscissa_indices_1d, int_t[:] num_abscissa_1d,
int_t[:] num_active_abscissa_1d, int_t[:] shifts,
double[:,:] abscissa_and_weights):

# active_dims is type unstable! Switches from int64 to int32

cdef:
Py_ssize_t num_act_dims_pt,ii,jj,kk,mm,cnt,dim,num_abscissa,act_dim_idx,fn_val_index,dd
bint has_inactive_abscissa, is_active_dim, done
double x_dim_k, denom, denom_d, basis

double mach_eps = np.finfo(np.float64).eps

Py_ssize_t num_pts = x.shape[1]
Py_ssize_t num_act_dims = active_dims.shape[0]

Py_ssize_t max_num_abscissa_1d = abscissa_and_weights.shape[0]//2
int_t[:] multi_index = np.empty((num_act_dims), dtype=np.int_)
np.ndarray[np.int64_t] multi_index = np.empty((num_act_dims), dtype=np.int64)

Py_ssize_t num_qoi = fn_vals.shape[1]

np.ndarray[double_t, ndim=2] result = np.empty((num_pts,num_qoi), dtype=np.float64)
double[:,:] result_view = result

# Allocate persistent memory. Each point will fill in a varying amount
# of entries. We use a view of this memory to stop reallocation for each
# of entries. We use a view of this memory to stop reallocation for each
# data point
cdef int_t[:] act_dims_pt_persistent = np.empty((num_act_dims),dtype=np.int_)
cdef int_t[:] act_dim_indices_pt_persistent = np.empty(
(num_act_dims),dtype=np.int_)
cdef np.ndarray[np.int64_t] act_dims_pt_persistent = np.empty((num_act_dims), dtype=np.int64)
cdef np.ndarray[np.int64_t] act_dim_indices_pt_persistent = np.empty((num_act_dims), dtype=np.int64)

cdef:
double[:,:] c_persistent=np.empty((num_qoi,num_act_dims),dtype=np.float64)
double[:,:] bases = np.empty(
(max_num_abscissa_1d, num_act_dims),dtype=np.float64)

for kk in range(num_pts):
# compute the active dimension of the kth point in x and the
# compute the active dimension of the kth point in x and the
# set multi_index accordingly
multi_index[:] = 0
num_act_dims_pt = 0
Expand All @@ -93,12 +97,12 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
num_abscissa = num_abscissa_1d[act_dim_idx]
x_dim_k = x[dim,kk]
for ii in range(num_abscissa):
if ( ( cnt < num_active_abscissa_1d[act_dim_idx] ) and
if ( ( cnt < num_active_abscissa_1d[act_dim_idx] ) and
( ii == active_abscissa_indices_1d[act_dim_idx][cnt] ) ):
cnt+=1
if ( abs( x_dim_k - abscissa_and_weights[2*ii,act_dim_idx] ) < mach_eps ):
is_active_dim = False
if ( ( cnt > 0 ) and
if ( ( cnt > 0 ) and
(active_abscissa_indices_1d[act_dim_idx][cnt-1]==ii)):
multi_index[act_dim_idx] = cnt-1
else:
Expand All @@ -110,7 +114,7 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
act_dim_indices_pt_persistent[num_act_dims_pt] = act_dim_idx
num_act_dims_pt+=1
# end for act_dim_idx in range(num_act_dims):

if ( has_inactive_abscissa ):
result_view[kk,:] = 0.
else:
Expand All @@ -122,22 +126,22 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
num_abscissa = num_abscissa_1d[act_dim_idx]
x_dim_k = x[dim,kk]
bases[0,dd] = abscissa_and_weights[1,act_dim_idx] /\
(x_dim_k - abscissa_and_weights[0,act_dim_idx])
(x_dim_k - abscissa_and_weights[0,act_dim_idx])
denom_d = bases[0,dd]
for ii in range(1,num_abscissa):
basis = abscissa_and_weights[2*ii+1,act_dim_idx] /\
(x_dim_k - abscissa_and_weights[2*ii,act_dim_idx])
bases[ii,dd] = basis
bases[ii,dd] = basis
denom_d += basis
# this is sometimes not a problem, just check all values are
# finite when this function is called from python

# this is sometimes not a problem, just check all values are
# finite when this function is called from python
#if ( abs(denom_d) < mach_eps ):
# raise Exception, "interpolation absacissa are not unique"
# raise Exception, "interpolation absacissa are not unique"
denom *= denom_d

# end for dd in range(num_act_dims_pt):

if ( num_act_dims_pt == 0 ):
# if point is an abscissa return the fn value at that point
# fn_val_index = np.sum(multi_index*shifts)
Expand All @@ -153,44 +157,40 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
fn_val_index = 0
for dd in range(num_act_dims):
fn_val_index += multi_index[dd]*shifts[dd]
#fn_val_index = np.sum(multi_index*shifts)

while (True):
act_dim_idx = act_dim_indices_pt_persistent[0]
for ii in range(num_active_abscissa_1d[act_dim_idx]):
fn_val_index+=shifts[act_dim_idx]*(ii-multi_index[act_dim_idx])
multi_index[act_dim_idx] = ii
basis=bases[active_abscissa_indices_1d[act_dim_idx][ii],0]
#c_persistent[:,0]+= basis * fn_vals[fn_val_index,:]
basis=bases[active_abscissa_indices_1d[act_dim_idx][ii], 0]
for mm in range(num_qoi):
c_persistent[mm,0] += basis * fn_vals[fn_val_index,mm]


for dd in range(1,num_act_dims_pt):
act_dim_idx = act_dim_indices_pt_persistent[dd]
basis = bases[active_abscissa_indices_1d[act_dim_idx][multi_index[act_dim_idx]],dd]
for mm in range(num_qoi):
c_persistent[mm,dd] += basis * c_persistent[mm,dd-1]
c_persistent[mm,dd-1] = 0.
#c_persistent[:,dd] += basis * c_persistent[:,dd-1]
#c_persistent[:,dd-1] = 0.

if (multi_index[act_dim_idx]<num_active_abscissa_1d[act_dim_idx]-1):
fn_val_index += shifts[act_dim_idx]
multi_index[act_dim_idx] += 1
multi_index[act_dim_idx] = multi_index[act_dim_idx] + 1
break
elif ( dd < num_act_dims_pt - 1 ):
fn_val_index-=shifts[act_dim_idx]*multi_index[act_dim_idx]
multi_index[act_dim_idx] = 0
else:
done = True

if ( done ):
break

for mm in range(num_qoi):
result_view[kk,mm]=c_persistent[mm,num_act_dims_pt-1]/denom
#result[kk,:] = c_persistent[:,num_act_dims_pt-1] / denom
#if np.any(np.isnan(result[kk,:])):
# #print (c_persistent [:,num_act_dims_pt-1])
# #print (denom)
# raise Exception, 'Error values not finite'

return result


Expand All @@ -199,15 +199,15 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx(
@cython.wraparound(False) # Deactivate negative indexing.
cpdef tensor_product_lagrange_interpolation_pyx(
double[:, :] x, double[:, :] fn_vals, double[:, :, :] basis_vals_1d,
int_t[:, :] active_indices, int_t[:] active_vars):
int[:, :] active_indices, int[:] active_vars):

cdef Py_ssize_t ii, jj, dd, kk
cdef Py_ssize_t nindices = active_indices.shape[1]
cdef Py_ssize_t nsamples = x.shape[1]
cdef Py_ssize_t nqoi = fn_vals.shape[1]
cdef Py_ssize_t nactive_vars = active_vars.shape[0]
cdef double basis_val = 1

values = np.zeros((nsamples, nqoi), dtype=float)
cdef double [:, :] values_view = values
for jj in range(nindices):
Expand All @@ -218,4 +218,3 @@ cpdef tensor_product_lagrange_interpolation_pyx(
for kk in range(nqoi):
values_view[ii, kk] += basis_val*fn_vals[jj, kk]
return values

11 changes: 7 additions & 4 deletions pyapprox/sparse_grid.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -534,16 +534,19 @@ def evaluate_sparse_grid_subspace(samples, subspace_index, subspace_values,
barycentric_weights_1d = []
for dd in range(num_active_sample_vars):
active_idx = active_sample_vars[dd]
abscissa_1d.append(samples_1d[active_idx][subspace_index[active_idx]])
abscissa_1d.append(samples_1d[active_idx][subspace_index[active_idx]])
abscissa_dd = abscissa_1d[dd]
interval_length = 2
if abscissa_1d[dd].shape[0] > 1:
interval_length = abscissa_1d[dd].max()-abscissa_1d[dd].min()
if abscissa_dd.shape[0] > 1:
interval_length = abscissa_dd.max() - abscissa_dd.min()

barycentric_weights_1d.append(
compute_barycentric_weights_1d(
abscissa_1d[dd], interval_length=interval_length))
abscissa_dd, interval_length=interval_length))

if num_active_sample_vars == 0:
return np.tile(subspace_values, (samples.shape[1], 1))

poly_vals = multivariate_barycentric_lagrange_interpolation(
samples, abscissa_1d, barycentric_weights_1d, subspace_values,
active_sample_vars)
Expand Down