diff --git a/pyapprox/barycentric_interpolation.py b/pyapprox/barycentric_interpolation.py index 9ff67508..75984ab5 100644 --- a/pyapprox/barycentric_interpolation.py +++ b/pyapprox/barycentric_interpolation.py @@ -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') diff --git a/pyapprox/cython/barycentric_interpolation.pyx b/pyapprox/cython/barycentric_interpolation.pyx index f9ed581f..7b18f158 100644 --- a/pyapprox/cython/barycentric_interpolation.pyx +++ b/pyapprox/cython/barycentric_interpolation.pyx @@ -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 @@ -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): @@ -43,37 +49,35 @@ 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) @@ -81,7 +85,7 @@ cpdef multivariate_hierarchical_barycentric_lagrange_interpolation_pyx( (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 @@ -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: @@ -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: @@ -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) @@ -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] 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)