eventseries.event_series¶

Provides class for event series analysis, namely event synchronization (ES) and event coincidence analysis (ECA). In addition, a method for the generation of binary event series from continuous time series data is included. When instantiating a class, data must either be passed as an event matrix (for details see below) or as a continuous time series. Using the class, an ES or ECA matrix can be calculated to generate a climate network using the EventSeriesClimateNetwork class. Both ES and ECA may be called without instantiating an object of the class. Significance levels are provided using analytic calculations using Poisson point processes as a null model (for ECA only) or a Monte Carlo approach.

Modifications vs. original pyunicorn EventSeries¶

Modification 1 (constructor):: Removed the axis-swap heuristic that transposed data when the number of variables exceeded the number of timesteps. The heuristic assumed that datasets with more variables than timesteps were transposed by mistake, which is not always true.
Modification 2 (make_event_matrix):: Replaced np.quantile with np.nanquantile so that time series containing NaN values are handled gracefully instead of raising an error.
Modification 3 (matrix loops):: Added tqdm progress bars to _ndim_event_synchronization and _ndim_event_coincidence_analysis for visual feedback during long runs.
Modification 4 (getters):: Added get_T(), get_N(), get_timestamps(), get_taumax(), get_lag() accessor methods to expose key instance attributes without name-mangling gymnastics.
Optimization 1 — Sparse binary-search ES with lag support (finite taumax):: For finite taumax, precomputes sorted inner-event arrays once per node and uses np.searchsorted to find only the event pairs within the lag-shifted taumax window. Complexity O(N²·L·taumax/T) vs O(N²·L²). Results are exact (atol < 1e-12 vs pairwise reference).
Optimization 2 — Blocked dense ES (infinite taumax, lag=0):: Falls back to a 4-D blocked-dense tensor kernel that vectorises the double-count correction loop — the main bottleneck in the original code. Block size is tunable; joblib parallelization is supported.
Optimization 3 — Vectorized ECA via cumsum + BLAS (integer timestamps):: For uniform integer timestamps and integer-valued lag, builds lag=0 smeared indicator matrices with cumsum rolling windows, shifts by lag steps, then uses a single BLAS matrix multiply. Complexity O(T·N + N²). Precision is identical to pairwise within float32 rounding (atol < 1e-6).
Optimization 4 — Optional joblib parallelization:: All kernel dispatchers accept an n_jobs keyword argument. If joblib is installed, computation is distributed across n_jobs worker processes. Falls back to a single-process loop when joblib is absent.

Contributors¶

Guruprem Bishnoi — Modifications 1–4 and Optimizations 1–4 (2026)

class pyunicorn.eventseries.event_series.EventSeries(data, *, timestamps=None, taumax=inf, lag=0.0, threshold_method=None, threshold_values=None, threshold_types=None)[source]¶

Bases: Cached

__cache_state__() → Tuple[Hashable, ...][source]¶

Hashable tuple of mutable object attributes, which will determine the instance identity for ALL cached method lookups in this class, in addition to the built-in object id(). Returning an empty tuple amounts to declaring the object immutable in general. Mutable dependencies that are specific to a method should instead be declared via @Cached.method(attrs=(…)).

NOTE:

A subclass is responsible for the consistency and cost of this state descriptor. For example, hashing a large array attribute may be circumvented by declaring it as a property, with a custom setter method that increments a dedicated mutation counter.

__init__(data, *, timestamps=None, taumax=inf, lag=0.0, threshold_method=None, threshold_values=None, threshold_types=None)[source]¶

Initialize an instance of EventSeries. Input data must be a 2D numpy array with time as the first axis and variables as the second axis. Event data is stored as an eventmatrix.

Format of eventmatrix: An eventmatrix is a 2D numpy array with the first dimension covering the timesteps and the second dimensions covering the variables. Each variable at a specific timestep is either ‘1’ if an event occured or ‘0’ if it did not, e.g. for 3 variables with 10 timesteps the eventmatrix could look like

array([[0, 1, 0],
[1, 0, 1], [0, 0, 0], [1, 0, 1], [0, 1, 0], [0, 0, 0], [1, 0, 0], [0, 0, 1], [0, 1, 0], [0, 0, 0]])

If input data is not provided as an eventmatrix, the constructor tries to generate one using the make_event_matrix method. Default keyword arguments are used in this case.

Parameters:

data (2D Numpy array [time, variables]) – Event series array or array of non-binary variable values
timestamps (1D Numpy array) – Time points of events of data. If not provided, integer values are used
taumax (float) – maximum time difference between two events to be considered synchronous. Caution: For ES, the default is np.inf because of the intrinsic dynamic coincidence interval in ES. For ECA, taumax is a parameter of the method that needs to be defined!
lag (float) – extra time lag between the event series
threshold_method (str 'quantile' or 'value' or 1D numpy array or str 'quantile' or 'value') – specifies the method for generating a binary event matrix from an array of continuous time series. Default: None
threshold_values (1D Numpy array or float) – quantile or real number determining threshold for each variable. Default: None.
threshold_types (str 'above' or 'below' or 1D list of strings 'above' or 'below') – Determines for each variable if event is below or above threshold

__str__()[source]¶: Return a string representation of the EventSeries object.

_compute_eca_pairwise(window_type)[source]¶: Pairwise ECA loop — identical to original pyunicorn. Modification 3: tqdm progress bar added.

_compute_eca_vectorized_lag(window_type)[source]¶

Optimization 3: Vectorized ECA via cumsum rolling windows + BLAS matrix multiply, extended to support any integer-valued lag.

Core idea: build the lag=0 smeared indicator matrix with a cumsum rolling window, then shift it by lag steps along the time axis. Boundary corrections mirror the pairwise formula exactly.

For lag=0 the output is bit-for-bit identical to the pairwise path within float32 rounding (atol < 1e-6 vs pairwise reference). Complexity: O(T·N + N²) vs O(N²·L²) for the pairwise approach.

_compute_es_matrix_blocked(node_data, block_size=32, n_jobs=1)[source]¶

Optimization 2: Assemble the full N×N directed ES matrix from blocked 4-D computations (infinite taumax, lag=0).

Optimization 4: n_jobs > 1 uses joblib if available.

_compute_es_matrix_sparse_lag(node_data, n_jobs=1)[source]¶

Optimization 1: Assemble the full N×N directed ES matrix using the lag-aware sparse binary-search kernel. Works for any scalar lag value.

Row-based chunking minimises joblib serialisation overhead. Optimization 4: n_jobs > 1 uses joblib if available.

_eca_coincidence_rate(eventseriesx, eventseriesy, *, window_type='symmetric', ts1=None, ts2=None)[source]¶

Event coincidence analysis: Returns the coincidence rates of two event series for both directions

Parameters:

eventseriesx (1D Numpy array) – Event series containing ‘0’s and ‘1’s
eventseriesy (1D Numpy array) – Event series containing ‘0’s and ‘1’s
ts1 (1D Numpy array) – Event time array containing time points when events of event series 1 occur, not obligatory
ts2 (1D Numpy array) – Event time array containing time points when events of event series 2 occur, not obligatory
window_type (str {'retarded', 'advanced', 'symmetric'}) – Only for ECA. Determines if precursor coincidence rate (‘advanced’), trigger coincidence rate (‘retarded’) or a general coincidence rate with the symmetric interval [-taumax, taumax] are computed (‘symmetric’). Default: ‘symmetric’

Return type:

list

Returns:

Precursor coincidence rates [XY, YX]

_empirical_percentiles(method=None, n_surr=1000, symmetrization='directed', window_type='symmetric')[source]¶

Compute p-values of event synchronisation (ES) and event coincidence analysis (ECA) using a Monte-Carlo approach. Surrogates are obtained by shuffling the event series. ES/ECA scores of the surrogate event series are computed and p-values are the empirical percentiles of the original event series compared to the ES/ECA scores of the surrogates.

Parameters:

method (str 'ES' or 'ECA') – determines if ES or ECA should be used
n_surr (int) – number of surrogates for Monte-Carlo method
symmetrization (str {'directed', 'symmetric', 'antisym', 'mean', 'max', 'min'} for ES, str {'directed', 'mean', 'max', 'min'} for ECA) – determines if and which symmetrisation method should be used for the ES/ECA score matrix
window_type (str {'retarded', 'advanced', 'symmetric'}) – Only for ECA. Determines if precursor coincidence rate (‘advanced’), trigger coincidence rate (‘retarded’) or a general coincidence rate with the symmetric interval [-taumax, taumax] are computed (‘symmetric’). Default: ‘symmetric’

Return type:

2D numpy array

Returns:

p-values of the ES/ECA scores for all

_ndim_event_coincidence_analysis(window_type='symmetric')[source]¶

Computes NxN event coincidence matrix of event coincidence rate.

Routing (Optimization 3):: uniform integer timestamps + integer lag → vectorized cumsum + BLAS otherwise → pairwise loop (original)

Parameters:: window_type (str {'retarded', 'advanced', 'symmetric'}) – Only for ECA. Determines if precursor coincidence rate (‘advanced’), trigger coincidence rate (‘retarded’) or a general coincidence rate with the symmetric interval [-taumax, taumax] are computed (‘symmetric’). Default: ‘symmetric’
Return type:: NxN numpy array where N is the number of variables of the eventmatrix
Returns:: event coincidence matrix

_ndim_event_synchronization()[source]¶

Compute NxN event synchronization matrix [i,j] with event synchronization from j to i without symmetrization.

Routing (Optimizations 1 & 2):: finite taumax, any lag → sparse binary-search (fast, exact) infinite taumax, lag=0 → blocked-dense 4-D (vectorized) infinite taumax, lag≠0 → pairwise loop (fallback)

Return type:: NxN numpy array where N is the number of variables of the eventmatrix
Returns:: event synchronization matrix

_ndim_event_synchronization_pairwise()[source]¶

Pairwise ES loop — identical to original pyunicorn. Used as fallback when taumax=inf and lag≠0.

Modification 3: tqdm progress bar added.

_precompute_es_node_data()[source]¶: Optimization 1 / 2: Precompute inner-event times and per-event dynamic tau2 for every node. Called once before the matrix loop.

static _symmetrization_antisym(matrix)[source]¶