2
2
3
3
import cython
4
4
from cython import Py_ssize_t
5
+ from cython.parallel import prange
5
6
6
7
from libc.stdlib cimport malloc, free
7
8
from libc.string cimport memmove
8
9
from libc.math cimport fabs, sqrt
10
+ from cpython cimport bool
9
11
10
12
import numpy as np
11
13
cimport numpy as cnp
@@ -230,14 +232,15 @@ def kth_smallest(numeric[:] a, Py_ssize_t k) -> numeric:
230
232
231
233
@cython.boundscheck(False )
232
234
@cython.wraparound(False )
233
- def nancorr(ndarray[ float64_t , ndim = 2 ] mat, bint cov = 0 , minp = None ):
235
+ def nancorr(float64_t[:, : ] mat , bint cov = 0 , minp = None , bool parallel = False ):
234
236
cdef:
235
237
Py_ssize_t i, j, xi, yi, N, K
236
238
bint minpv
237
- ndarray[ float64_t, ndim = 2 ] result
238
- ndarray[ uint8_t, ndim = 2 ] mask
239
+ float64_t[:, : ] result
240
+ uint8_t[:, : ] mask
239
241
int64_t nobs = 0
240
242
float64_t vx, vy, sumx, sumy, sumxx, sumyy, meanx, meany, divisor
243
+ int64_t blah = 0
241
244
242
245
N, K = (< object > mat).shape
243
246
@@ -249,44 +252,82 @@ def nancorr(ndarray[float64_t, ndim=2] mat, bint cov=0, minp=None):
249
252
result = np.empty((K, K), dtype = np.float64)
250
253
mask = np.isfinite(mat).view(np.uint8)
251
254
252
- with nogil:
253
- for xi in range (K):
254
- for yi in range (xi + 1 ):
255
- nobs = sumxx = sumyy = sumx = sumy = 0
256
- for i in range (N):
257
- if mask[i, xi] and mask[i, yi]:
258
- vx = mat[i, xi]
259
- vy = mat[i, yi]
260
- nobs += 1
261
- sumx += vx
262
- sumy += vy
255
+ if parallel:
256
+ with nogil:
257
+ for xi in prange(K, schedule = ' dynamic' ):
258
+ nancorr_single_row(mat, N, K, result, xi, mask, minpv, cov)
259
+ else :
260
+ with nogil:
261
+ for xi in range (K):
262
+ nancorr_single_row(mat, N, K, result, xi, mask, minpv, cov)
263
263
264
- if nobs < minpv:
265
- result[xi, yi] = result[yi, xi] = NaN
266
- else :
267
- meanx = sumx / nobs
268
- meany = sumy / nobs
264
+ return np.asarray(result)
269
265
270
- # now the cov numerator
271
- sumx = 0
272
266
273
- for i in range (N):
274
- if mask[i, xi] and mask[i, yi]:
275
- vx = mat[i, xi] - meanx
276
- vy = mat[i, yi] - meany
267
+ @ cython.boundscheck (False )
268
+ @ cython.wraparound (False )
269
+ cdef void nancorr_single_row(float64_t[:, :] mat,
270
+ Py_ssize_t N,
271
+ Py_ssize_t K,
272
+ float64_t[:, :] result,
273
+ Py_ssize_t xi,
274
+ uint8_t[:, :] mask,
275
+ bint minpv,
276
+ bint cov = 0 ) nogil:
277
+ for yi in range (xi + 1 ):
278
+ nancorr_single(mat, N, K, result, xi, yi, mask, minpv, cov)
277
279
278
- sumx += vx * vy
279
- sumxx += vx * vx
280
- sumyy += vy * vy
281
280
282
- divisor = (nobs - 1.0 ) if cov else sqrt(sumxx * sumyy)
281
+ @ cython.boundscheck (False )
282
+ @ cython.wraparound (False )
283
+ cdef void nancorr_single(float64_t[:, :] mat,
284
+ Py_ssize_t N,
285
+ Py_ssize_t K,
286
+ float64_t[:, :] result,
287
+ Py_ssize_t xi,
288
+ Py_ssize_t yi,
289
+ uint8_t[:, :] mask,
290
+ bint minpv,
291
+ bint cov = 0 ) nogil:
292
+ cdef:
293
+ Py_ssize_t i, j
294
+ int64_t nobs = 0
295
+ float64_t vx, vy, sumx, sumy, sumxx, sumyy, meanx, meany, divisor
283
296
284
- if divisor != 0 :
285
- result[xi, yi] = result[yi, xi] = sumx / divisor
286
- else :
287
- result[xi, yi] = result[yi, xi] = NaN
297
+ nobs = sumxx = sumyy = sumx = sumy = 0
298
+ for i in range (N):
299
+ if mask[i, xi] and mask[i, yi]:
300
+ vx = mat[i, xi]
301
+ vy = mat[i, yi]
302
+ nobs += 1
303
+ sumx += vx
304
+ sumy += vy
305
+
306
+ if nobs < minpv:
307
+ result[xi, yi] = result[yi, xi] = NaN
308
+ else :
309
+ meanx = sumx / nobs
310
+ meany = sumy / nobs
311
+
312
+ # now the cov numerator
313
+ sumx = 0
314
+
315
+ for i in range (N):
316
+ if mask[i, xi] and mask[i, yi]:
317
+ vx = mat[i, xi] - meanx
318
+ vy = mat[i, yi] - meany
319
+
320
+ sumx += vx * vy
321
+ sumxx += vx * vx
322
+ sumyy += vy * vy
323
+
324
+ divisor = (nobs - 1.0 ) if cov else sqrt(sumxx * sumyy)
325
+
326
+ if divisor != 0 :
327
+ result[xi, yi] = result[yi, xi] = sumx / divisor
328
+ else :
329
+ result[xi, yi] = result[yi, xi] = NaN
288
330
289
- return result
290
331
291
332
# ----------------------------------------------------------------------
292
333
# Pairwise Spearman correlation
0 commit comments