Skip to content
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
Next Next commit
lint cython
  • Loading branch information
MuellerSeb committed Apr 10, 2024
commit 773f0bd3beec85abadf0ee50167be87ffc28cb44
79 changes: 35 additions & 44 deletions src/pykrige/lib/cok.pyx
Original file line number Diff line number Diff line change
@@ -1,9 +1,7 @@
#cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
# -*- coding: utf-8 -*-
# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
import numpy as np

cimport numpy as np
from libc.math cimport sqrt

import scipy.linalg

Expand All @@ -13,12 +11,14 @@ from scipy.linalg.cython_lapack cimport dgesv
from .variogram_models cimport get_variogram_model


cpdef _c_exec_loop(double [:, ::1] a_all,
double [:, ::1] bd_all,
char [::1] mask,
long n,
dict pars):
cdef long i, j, k
cpdef _c_exec_loop(
double [:, ::1] a_all,
double [:, ::1] bd_all,
char [::1] mask,
long n,
dict pars
):
cdef long i, k

npt = bd_all.shape[0]

Expand All @@ -34,13 +34,11 @@ cpdef _c_exec_loop(double [:, ::1] a_all,

nb = n + 1


c_variogram_function = get_variogram_model(pars['variogram_function'].__name__)

cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters'])

cdef double [::1,:] a_inv = np.asfortranarray(np.empty_like(a_all))

cdef double [::1, :] a_inv = np.asfortranarray(np.empty_like(a_all))

if pars['pseudo_inv']:
if pars['pseudo_inv_type'] == "pinv":
Expand All @@ -54,8 +52,8 @@ cpdef _c_exec_loop(double [:, ::1] a_all,
else:
a_inv = np.asfortranarray(scipy.linalg.inv(a_all))


for i in range(npt): # same thing as range(npt) if mask is not defined, otherwise take the non masked elements
# same thing as range(npt) if mask is not defined, otherwise take the non masked elements
for i in range(npt):
if mask[i]:
continue
bd = bd_all[i]
Expand All @@ -69,22 +67,20 @@ cpdef _c_exec_loop(double [:, ::1] a_all,
if pars['exact_values']:
check_b_vect(n, bd, b, eps)


# Do the BLAS matrix-vector multiplication call
dgemv(
# # Compute y := alpha*A*x + beta*y
'N', # char *trans, # {'T','C'}: o(A)=A'; {'N'}: o(A)=A
&nb, # int *m, # Rows of A (prior to transpose from *trans)
&nb, # int *n, # Columns of A / min(len(x))
&alpha, # np.float64_t *alpha, # Scalar multiple
&(a_inv[0,0]), # np.float64_t *a, # Matrix A: mxn
&nb, # int *lda, # The size of the first dimension of A (in memory)
&(b[0]), # np.float64_t *x, # Vector x, min(len(x)) = n
&inc, # int *incx, # The increment between elements of x (usually 1)
&beta, # np.float64_t *beta, # Scalar multiple
&(x[0]), # np.float64_t *y, # Vector y, min(len(y)) = m
&inc # int *incy # The increment between elements of y (usually 1)
)
dgemv( # Compute y := alpha*A*x + beta*y
'N', # char *trans, # {'T','C'}: o(A)=A'; {'N'}: o(A)=A
&nb, # int *m, # Rows of A (prior to transpose from *trans)
&nb, # int *n, # Columns of A / min(len(x))
&alpha, # np.float64_t *alpha, # Scalar multiple
&(a_inv[0, 0]), # np.float64_t *a, # Matrix A: mxn
&nb, # int *lda, # The size of the first dimension of A (in memory)
&(b[0]), # np.float64_t *x, # Vector x, min(len(x)) = n
&inc, # int *incx, # The increment between elements of x (usually 1)
&beta, # np.float64_t *beta, # Scalar multiple
&(x[0]), # np.float64_t *y, # Vector y, min(len(y)) = m
&inc # int *incy # The increment between elements of y (usually 1)
)

z_tmp = 0.0
ss_tmp = 0.0
Expand All @@ -99,12 +95,14 @@ cpdef _c_exec_loop(double [:, ::1] a_all,

return zvalues.base, sigmasq.base

cpdef _c_exec_loop_moving_window(double [:, ::1] a_all,
double [:, ::1] bd_all,
char [::1] mask,
long [:,::1] bd_idx,
long n_max,
dict pars):
cpdef _c_exec_loop_moving_window(
double [:, ::1] a_all,
double [:, ::1] bd_all,
char [::1] mask,
long [:, ::1] bd_idx,
long n_max,
dict pars
):
cdef long i, j, k, p_i, p_j, npt, n

npt = bd_all.shape[0]
Expand All @@ -121,19 +119,14 @@ cpdef _c_exec_loop_moving_window(double [:, ::1] a_all,
cdef double z_tmp, ss_tmp, eps=pars['eps']
cdef int nb, nrhs=1, info
cdef int [::1] ipiv = np.zeros(n_max+1, dtype='int32')
cdef double alpha=1.0, beta=0.0
cdef long [::1] bd_idx_sel

nb = n + 1



c_variogram_function = get_variogram_model(pars['variogram_function'].__name__)

cdef double [::1] variogram_model_parameters = np.asarray(pars['variogram_model_parameters'])



for i in range(npt):
if mask[i]:
continue
Expand All @@ -142,7 +135,6 @@ cpdef _c_exec_loop_moving_window(double [:, ::1] a_all,

bd_idx_sel = bd_idx[i, :]


for k in range(n):
p_i = bd_idx_sel[k]
for j in range(n):
Expand All @@ -167,8 +159,8 @@ cpdef _c_exec_loop_moving_window(double [:, ::1] a_all,
x[k] = b[k]

# higher level (and slower) call to do the same thing as dgesv below
#a2D = a_selection.base[:nb*nb].reshape((nb, nb))
#x = scipy.linalg.solve(a2D, b.base[:nb])
# a2D = a_selection.base[:nb*nb].reshape((nb, nb))
# x = scipy.linalg.solve(a2D, b.base[:nb])

dgesv(
&nb,
Expand All @@ -186,7 +178,6 @@ cpdef _c_exec_loop_moving_window(double [:, ::1] a_all,
elif info < 0:
raise ValueError('Wrong arguments')


z_tmp = 0.0
ss_tmp = 0.0
for k in range(n):
Expand Down
2 changes: 1 addition & 1 deletion src/pykrige/lib/variogram_models.pxd
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@

# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
ctypedef void (*variogram_model_t)(double [::1], long, double [::1], double [::1])

cdef variogram_model_t get_variogram_model(function_name)
31 changes: 18 additions & 13 deletions src/pykrige/lib/variogram_models.pyx
Original file line number Diff line number Diff line change
@@ -1,13 +1,8 @@
#cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
# -*- coding: utf-8 -*-
import numpy as np
cimport numpy as np
# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True
from libc.math cimport exp

# copied from variogram_model.py



# copied from variogram_model.py
cdef variogram_model_t get_variogram_model(name):
cdef variogram_model_t c_func

Expand All @@ -27,7 +22,9 @@ cdef variogram_model_t get_variogram_model(name):
return c_func


cdef void _c_linear_variogram_model(double [::1] params, long n, double [::1] dist, double[::1] out) nogil:
cdef void _c_linear_variogram_model(
double [::1] params, long n, double [::1] dist, double[::1] out
) noexcept nogil:
cdef long k
cdef double a, b
a = params[0]
Expand All @@ -36,17 +33,21 @@ cdef void _c_linear_variogram_model(double [::1] params, long n, double [::1] d
out[k] = a*dist[k] + b


cdef void _c_power_variogram_model(double [::1] params, long n, double [::1] dist, double[::1] out) nogil:
cdef void _c_power_variogram_model(
double [::1] params, long n, double [::1] dist, double[::1] out
) noexcept nogil:
cdef long k
cdef double a, b, c
a = params[0]
b = params[1]
c = params[2]
for k in range(n):
out[k] = a*(dist[k]**b) + c
out[k] = a*(dist[k]**b) + c


cdef void _c_gaussian_variogram_model(double [::1] params, long n, double [::1] dist, double [::1] out) nogil:
cdef void _c_gaussian_variogram_model(
double [::1] params, long n, double [::1] dist, double [::1] out
) noexcept nogil:
cdef long k
cdef double a, b, c
a = params[0]
Expand All @@ -56,7 +57,9 @@ cdef void _c_gaussian_variogram_model(double [::1] params, long n, double [::1]
out[k] = a*(1 - exp(-(dist[k]/(b*4.0/7.0))**2)) + c


cdef void _c_exponential_variogram_model(double [::1] params, long n, double[::1] dist, double[::1] out) nogil:
cdef void _c_exponential_variogram_model(
double [::1] params, long n, double[::1] dist, double[::1] out
) noexcept nogil:
cdef long k
cdef double a, b, c
a = params[0]
Expand All @@ -66,7 +69,9 @@ cdef void _c_exponential_variogram_model(double [::1] params, long n, double[::1
out[k] = a*(1 - exp(-dist[k]/(b/3.0))) + c


cdef void _c_spherical_variogram_model(double [::1] params, long n, double[::1] dist, double[::1] out) nogil:
cdef void _c_spherical_variogram_model(
double [::1] params, long n, double[::1] dist, double[::1] out
) noexcept nogil:
cdef long k
cdef double a, b, c
a = params[0]
Expand Down