"""

d, n = T.shape

l = n - m + 1

D = np.empty((d, l), dtype=np.float64)

for i in range(d):

if np.any(~Q_subseq_isfinite[i]):

D[i, :] = np.inf

else:

D[i, :] = core.mass_absolute(Q[i], T[i], p=p)

D[i][~(T_subseq_isfinite[i])] = np.inf

T,

m,

subseq_idx,

nn_idx,

k,

include=None,

discords=False,

discretize_func=None,

n_bit=8,

p=2.0,

):

"""

T = core._preprocess(T)

core.check_window_size(m, max_size=T.shape[1], n=T.shape[1])

subseqs, _ = core.preprocess_non_normalized(T[:, subseq_idx : subseq_idx + m], m)

neighbors, _ = core.preprocess_non_normalized(T[:, nn_idx : nn_idx + m], m)

if discretize_func is None:

T_isfinite = np.isfinite(T)

T_min = T[T_isfinite].min()

T_max = T[T_isfinite].max()

discretize_func = partial(

_maamp_discretize, a_min=T_min, a_max=T_max, n_bit=n_bit

)

disc_subseqs = discretize_func(subseqs)

disc_neighbors = discretize_func(neighbors)

D = np.linalg.norm(disc_subseqs - disc_neighbors, axis=1, ord=p)

S = core._subspace(D, k, include=include, discords=discords)

"""

return (

np.round(((a - a_min) / (a_max - a_min)) * ((2**n_bit) - 1.0)).astype(np.int64)

+ 1

T,

m,

subseq_idx,

nn_idx,

include=None,

discords=False,

discretize_func=None,

n_bit=8,

p=2.0,

):

"""

T = core._preprocess(T)

core.check_window_size(m, max_size=T.shape[1], n=T.shape[1])

if discretize_func is None:

T_isfinite = np.isfinite(T)

T_min = T[T_isfinite].min()

T_max = T[T_isfinite].max()

discretize_func = partial(

_maamp_discretize, a_min=T_min, a_max=T_max, n_bit=n_bit

)

bit_sizes = np.empty(T.shape[0])

S = [None] * T.shape[0]

for k in range(T.shape[0]):

subseqs, _ = core.preprocess_non_normalized(

T[:, subseq_idx[k] : subseq_idx[k] + m], m

)

neighbors, _ = core.preprocess_non_normalized(

T[:, nn_idx[k] : nn_idx[k] + m], m

)

disc_subseqs = discretize_func(subseqs)

disc_neighbors = discretize_func(neighbors)

D = np.linalg.norm(disc_subseqs - disc_neighbors, axis=1, ord=p)

S[k] = core._subspace(D, k, include=include, discords=discords)

bit_sizes[k] = core._mdl(disc_subseqs, disc_neighbors, S[k], n_bit=n_bit)

query_idx,

T_A,

T_B,

m,

T_B_subseq_isfinite,

p=2.0,

include=None,

discords=False,

excl_zone=None,

):

"""

d, n = T_A.shape

l = n - m + 1

start_row_idx = 0

D = _multi_mass_absolute(

T_B[:, query_idx : query_idx + m],

T_A,

m,

T_B_subseq_isfinite[:, query_idx],

T_B_subseq_isfinite,

p,

)

if include is not None:

core._apply_include(D, include)

start_row_idx = include.shape[0]

if discords:

D[start_row_idx:][::-1].sort(axis=0, kind="mergesort")

else:

D[start_row_idx:].sort(axis=0, kind="mergesort")

D_prime = np.zeros(l, dtype=np.float64)

for i in range(d):

D_prime[:] = D_prime + D[i]

D[i, :] = D_prime / (i + 1)

if excl_zone is not None:

core.apply_exclusion_zone(D, query_idx, excl_zone, np.inf)

"""

T, T_subseq_isfinite = core.preprocess_non_normalized(T, m)

if T.ndim <= 1: # pragma: no cover

err = f"T is {T.ndim}-dimensional and must be at least 2-dimensional"

raise ValueError(f"{err}")

core.check_window_size(m, max_size=T.shape[1], n=T.shape[1])

if include is not None: # pragma: no cover

include = core._preprocess_include(include)

excl_zone = int(

np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM)

) # See Definition 3 and Figure 3

D = _maamp_multi_distance_profile(

query_idx, T, T, m, T_subseq_isfinite, p, include, discords, excl_zone

)

start,

T_A,

T_B,

m,

excl_zone,

T_B_subseq_isfinite,

p=2.0,

include=None,

discords=False,

):

"""

D = _maamp_multi_distance_profile(

start, T_A, T_B, m, T_B_subseq_isfinite, p, include, discords, excl_zone

)

d = T_A.shape[0]

P = np.full(d, np.inf, dtype=np.float64)

I = np.full(d, -1, dtype=np.int64)

for i in range(d):

min_index = np.argmin(D[i])

I[i] = min_index

P[i] = D[i, min_index]

if np.isinf(P[i]): # pragma nocover

I[i] = -1

"""

d = T.shape[0]

l = T.shape[1] - m + 1

p_norm = np.empty((d, l), dtype=np.float64)

p_norm_first = np.empty((d, l), dtype=np.float64)

for i in range(d):

p_norm[i] = np.power(core.mass_absolute(T[i, start : start + m], T[i], p=p), p)

p_norm_first[i] = np.power(core.mass_absolute(T[i, :m], T[i], p=p), p)

d,

k,

idx,

p_norm,

T,

m,

excl_zone,

T_A_subseq_isfinite,

T_B_subseq_isfinite,

p,

p_norm_even,

p_norm_odd,

p_norm_first,

):

"""

for i in range(d):

# Numba's prange requires incrementing a range by 1 so replace

# `for j in range(k-1,0,-1)` with its incrementing compliment

for rev_j in prange(1, k):

j = k - rev_j

# GPU Stomp Parallel Implementation with Numba

# DOI: 10.1109/ICDM.2016.0085

# See Figure 5

if idx % 2 == 0:

# Even

p_norm_even[i, j] = (

p_norm_odd[i, j - 1]

- abs(T[i, idx - 1] - T[i, j - 1]) ** p

+ abs(T[i, idx + m - 1] - T[i, j + m - 1]) ** p

)

else:

# Odd

p_norm_odd[i, j] = (

p_norm_even[i, j - 1]

- abs(T[i, idx - 1] - T[i, j - 1]) ** p

+ abs(T[i, idx + m - 1] - T[i, j + m - 1]) ** p

)

if idx % 2 == 0:

p_norm_even[i, 0] = p_norm_first[i, idx]

if not T_B_subseq_isfinite[i, idx]:

p_norm[i] = np.inf

else:

p_norm[i] = p_norm_even[i]

else:

p_norm_odd[i, 0] = p_norm_first[i, idx]

if not T_B_subseq_isfinite[i, idx]:

p_norm[i] = np.inf

else:

p_norm[i] = p_norm_odd[i]

p_norm[i][~(T_A_subseq_isfinite[i])] = np.inf

p_norm[i][p_norm[i] < config.STUMPY_P_NORM_THRESHOLD] = 0

T,

m,

range_stop,

excl_zone,

T_A_subseq_isfinite,

T_B_subseq_isfinite,

p,

p_norm,

p_norm_first,

k,

range_start=1,

include=None,

discords=False,

):

"""

p_norm_odd = p_norm.copy()

p_norm_even = p_norm.copy()

d = T.shape[0]

P = np.empty((d, range_stop - range_start), dtype=np.float64)

I = np.empty((d, range_stop - range_start), dtype=np.int64)

p_norm = np.empty((d, k), dtype=np.float64)

p_norm_prime = np.empty(k, dtype=np.float64)

start_row_idx = 0

if include is not None:

tmp_swap = np.empty((include.shape[0], k), dtype=np.float64)

restricted_indices = include[include < include.shape[0]]

unrestricted_indices = include[include >= include.shape[0]]

mask = np.ones(include.shape[0], dtype=bool)

mask[restricted_indices] = False

for idx in range(range_start, range_stop):

_compute_multi_p_norm(

d,

k,

idx,

p_norm,

T,

m,

excl_zone,

T_A_subseq_isfinite,

T_B_subseq_isfinite,

p,

p_norm_even,

p_norm_odd,

p_norm_first,

)

# `include` processing must occur here since we are dealing with distances

if include is not None:

core._apply_include(

p_norm,

include,

restricted_indices,

unrestricted_indices,

mask,

tmp_swap,

)

start_row_idx = include.shape[0]

if discords:

p_norm[start_row_idx:][::-1].sort(axis=0)

else:

p_norm[start_row_idx:].sort(axis=0)

core._compute_multi_PI(d, idx, p_norm, p_norm_prime, range_start, P, I, p)

"""

T_A = T

T_B = T_A

T_A, T_A_subseq_isfinite = core.preprocess_non_normalized(T_A, m)

T_B, T_B_subseq_isfinite = core.preprocess_non_normalized(T_B, m)

if T_A.ndim <= 1: # pragma: no cover

err = f"T is {T_A.ndim}-dimensional and must be at least 1-dimensional"

raise ValueError(f"{err}")

core.check_window_size(m, max_size=min(T_A.shape[1], T_B.shape[1]), n=T_A.shape[1])

if include is not None:

include = core._preprocess_include(include)

d, n = T_B.shape

l = n - m + 1

excl_zone = int(

np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM)

) # See Definition 3 and Figure 3

P = np.empty((d, l), dtype=np.float64)

I = np.empty((d, l), dtype=np.int64)

start = 0

stop = l

P[:, start], I[:, start] = _get_first_maamp_profile(

start,

T_A,

T_B,

m,

excl_zone,

T_B_subseq_isfinite,

p,

include,

discords,

)

p_norm, p_norm_first = _get_multi_p_norm(start, T_A, m, p=p)

P[:, start + 1 : stop], I[:, start + 1 : stop] = _maamp(

T_A,

m,

stop,

excl_zone,

T_A_subseq_isfinite,

T_B_subseq_isfinite,

p,

p_norm,

p_norm_first,

l,

start + 1,

include,

discords,

)