Skip to content

Commit

Permalink
MAINT: Fix Cython/C type problems
Browse files Browse the repository at this point in the history
- Single-source Cython/C type definitions: `core._ext.types`
- Introduce semantic type aliases: ADJ, NODE, DEGREE, WEIGHT, FIELD, etc.
- Remove deprecated type aliases: `np.{int|float}`
- Replace C type aliases by Numpy type aliases
- Add missing Numpy type casts before Python/Cython boundary
- Disambiguate `timeseries._ext.numerics._twin_surrogates_{s|r}`
  (blame: 098489b & 2f6a5b9)
- Fix some broken output pointers
- Reduce redundant copying
- Vectorize some simple Python loops
- Eliminate some dead code
- Potentially related issues: #128, #141, #142, #145
  • Loading branch information
ntfrgl committed Apr 8, 2022
1 parent a6c4c83 commit 3dab5bf
Show file tree
Hide file tree
Showing 24 changed files with 753 additions and 992 deletions.
54 changes: 17 additions & 37 deletions pyunicorn/climate/_ext/numerics.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -15,57 +15,37 @@
# for complex systems science: The pyunicorn package"

cimport cython
from cpython cimport bool


import numpy as np
cimport numpy as np


BOOLTYPE = np.uint8
INTTYPE = np.int
INT8TYPE = np.int8
INT16TYPE = np.int16
INT32TYPE = np.int32
INT64TYPE = np.int64
FLOATTYPE = np.float
FLOAT32TYPE = np.float32
FLOAT64TYPE = np.float64
ctypedef np.uint8_t BOOLTYPE_t
ctypedef np.int_t INTTYPE_t
ctypedef np.int8_t INT8TYPE_t
ctypedef np.int16_t INT16TYPE_t
ctypedef np.int32_t INT32TYPE_t
ctypedef np.int64_t INT64TYPE_t
ctypedef np.float_t FLOATTYPE_t
ctypedef np.float32_t FLOAT32TYPE_t
ctypedef np.float64_t FLOAT64TYPE_t

from ...core._ext.types import FIELD, DFIELD
from ...core._ext.types cimport MASK_t, FIELD_t, DFIELD_t

cdef extern from "src_numerics.c":
void _cython_calculate_mutual_information(
float *anomaly, int n_samples, int N, int n_bins, double scaling,
double range_min, long *symbolic, long *hist, long *hist2d,
float *mi)
void _calculate_corr_fast(int m, int tmax, int *final_mask,
void _calculate_corr_fast(int m, int tmax, bint *final_mask,
float *time_series_ranked, float *spearman_rho)


# mutual_info =================================================================

def _calculate_mutual_information_cython(
np.ndarray[float, ndim=2, mode='c'] anomaly not None,
np.ndarray[FIELD_t, ndim=2, mode='c'] anomaly not None,
int n_samples, int N, int n_bins, double scaling, double range_min):

cdef:
np.ndarray[long, ndim=2, mode='c'] symbolic = \
np.zeros((N, n_samples), dtype='int64')
np.ndarray[long, ndim=2, mode='c'] hist = \
np.zeros((N, n_bins), dtype='int64')
np.ndarray[long, ndim=2, mode='c'] hist2d = \
np.zeros((n_bins, n_bins), dtype='int64')
np.ndarray[float, ndim=2, mode='c'] mi = \
np.zeros((N, N), dtype='float32')
np.ndarray[DFIELD_t, ndim=2, mode='c'] symbolic = np.zeros(
(N, n_samples), dtype=DFIELD)
np.ndarray[DFIELD_t, ndim=2, mode='c'] hist = np.zeros(
(N, n_bins), dtype=DFIELD)
np.ndarray[DFIELD_t, ndim=2, mode='c'] hist2d = np.zeros(
(n_bins, n_bins), dtype=DFIELD)
np.ndarray[FIELD_t, ndim=2, mode='c'] mi = np.zeros(
(N, N), dtype=FIELD)

_cython_calculate_mutual_information(
<float*> np.PyArray_DATA(anomaly), n_samples, N, n_bins, scaling,
Expand All @@ -79,14 +59,14 @@ def _calculate_mutual_information_cython(
# rainfall ====================================================================

def _calculate_corr(int m, int tmax,
np.ndarray[int, ndim=2, mode='c'] final_mask not None,
np.ndarray[float, ndim=2, mode='c'] time_series_ranked not None):
np.ndarray[MASK_t, ndim=2, mode='c'] final_mask not None,
np.ndarray[FIELD_t, ndim=2, mode='c'] time_series_ranked not None):

cdef np.ndarray[float, ndim=2, mode='c'] spearman_rho = \
np.zeros((m, m), dtype='float32')
cdef np.ndarray[FIELD_t, ndim=2, mode='c'] spearman_rho = np.zeros(
(m, m), dtype=FIELD)

_calculate_corr_fast(m, tmax,
<int*> np.PyArray_DATA(final_mask),
<bint*> np.PyArray_DATA(final_mask),
<float*> np.PyArray_DATA(time_series_ranked),
<float*> np.PyArray_DATA(spearman_rho))

Expand Down
4 changes: 3 additions & 1 deletion pyunicorn/climate/_ext/src_numerics.c
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@
* License: BSD (3-clause)
*/

#include <math.h>

// mutual_info ================================================================

void _cython_calculate_mutual_information(float *anomaly, int n_samples,
Expand All @@ -18,7 +20,7 @@ void _cython_calculate_mutual_information(float *anomaly, int n_samples,
in_nodes;
double norm, rescaled, hpl, hpm, plm;

double *p_anomaly;
float *p_anomaly;
float *p_mi, *p_mi2;
long *p_symbolic, *p_symbolic1, *p_symbolic2, *p_hist, *p_hist1,
*p_hist2, *p_hist2d;
Expand Down
7 changes: 3 additions & 4 deletions pyunicorn/climate/mutual_info.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
# array object and fast numerics
import numpy as np

from ..core._ext.types import to_cy, FIELD
from ._ext.numerics import _calculate_mutual_information_cython

# Import progress bar for easy progress bar handling
Expand Down Expand Up @@ -155,10 +156,8 @@ def _cython_calculate_mutual_information(self, anomaly, n_bins=32):
# using the maximum range of the whole dataset.
scaling = 1./(range_max - range_min)

anomaly = anomaly.astype(np.float32).copy(order='c')
mi = _calculate_mutual_information_cython(anomaly, n_samples, N,
n_bins, scaling,
range_min)
mi = _calculate_mutual_information_cython(
to_cy(anomaly, FIELD), n_samples, N, n_bins, scaling, range_min)

if self.silence_level <= 1:
print("Done!")
Expand Down
12 changes: 6 additions & 6 deletions pyunicorn/climate/rainfall.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,9 @@
# Import essential packages
#

import numpy as np

from ..core._ext.types import to_cy, MASK, FIELD
from ._ext.numerics import _calculate_corr

# Import cnTsonisClimateNetwork for TsonisClimateNetwork class
Expand Down Expand Up @@ -276,10 +279,7 @@ def calculate_corr(self, final_mask, anomaly):
:return: the Spearman correlation matrix.
"""
# Get rank time series
time_series_ranked = self.rank_time_series(anomaly).astype('float32')

time_series_ranked = self.rank_time_series(anomaly)
m, tmax = anomaly.shape

return _calculate_corr(m, tmax,
final_mask.astype('int32').copy(order='c'),
time_series_ranked.copy(order='c'))
return _calculate_corr(
m, tmax, to_cy(final_mask, MASK), to_cy(time_series_ranked, FIELD))
Loading

0 comments on commit 3dab5bf

Please sign in to comment.