Skip to content

Commit

Permalink
cythonized twin_surrogates
Browse files Browse the repository at this point in the history
  • Loading branch information
wbarfuss committed Sep 27, 2015
1 parent c71720a commit 098489b
Show file tree
Hide file tree
Showing 2 changed files with 89 additions and 83 deletions.
108 changes: 86 additions & 22 deletions pyunicorn/timeseries/numerics.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,13 @@ ctypedef np.float32_t FLOAT32TYPE_t
cdef extern from "stdlib.h":
double drand48()

cdef extern from "stdlib.h":
double srand48()

cdef extern from "time.h":
double time()


# recurrence plot==============================================================

def _embed_time_series(
Expand Down Expand Up @@ -376,33 +383,90 @@ def _twins(
int min_dist, int N, np.ndarray[INT8TYPE_t, ndim=2] R,
np.ndarray[INTTYPE_t, ndim=1] nR, twins):

cdef int j, k, l
cdef int j, k, l

twins.append([])

for j in xrange(N):
twins.append([])
twins_j = twins[j]

# Respect a minimal temporal spacing between twins to avoid false
# twins du to th higher sample density in phase space along the
# trajectory
for k in xrange(j - min_dist):
# Continue only if both samples have the same number of
# neighbors and more than just one neighbor (themselves)
if nR[j] == nR[k] and nR[j] != 1:
l = 0

while R[j, l] == R[k, l]:
l = l + 1

# If l is equal to the length of the time series at
# this point, j and k are twins
if l == N:
# And the twins to the twin list
twins_k = twins[k]

twins_j.append(k)
twins_k.append(j)

break


for j in xrange(N):
twins.append([])
twins_j = twins[j]
def _twin_surrogates(
int n_surrogates, int N, int dim, twins,
np.ndarray[FLOAT32TYPE_t, ndim=2] embedding,
np.ndarray[FLOATTYPE_t, ndim=3] surrogates):

# Respect a minimal temporal spacing between twins to avoid false
# twins du to th higher sample density in phase space along the
# trajectory
for k in xrange(j - min_dist):
# Continue only if both samples have the same number of
# neighbors and more than just one neighbor (themselves)
if nR[j] == nR[k] and nR[j] != 1:
l = 0
cdef int i, j, k, l, new_k, n_twins, rand

while R[j, l] == R[k, l]:
l = l + 1
# Initialize random number generator
# srand48(time(0)) -> does not work in cython somehow ?!?!?

for i in xrange(n_surrogates):
# Randomly choose a starting point in the original trajectory
k = int(floor(drand48() * N))

j = 0

while j < N:
# Assign state vector of surrogate trajectory
for l in xrange(dim):
surrogates[i, j, l] = embedding[k, l]

# Get the list of twins of state vector k in the original time
# series
twins_k = twins[k]

# Get the number of twins of k
n_twins = len(twins_k)

# If k has no twins, go to the next sample k+1, If k has twins at
# m, choose among m+1 and k+1 with equal probability
if n_twins == 0:
k += 1
else:
# Generate a random integer between 0 and n_twins
rand = int(floor(drand48() * (n_twins + 1)))

# If rand = n_twings go to smple k+1, otherwise jump to the
# future of one of the twins
if rand == n_twins:
k += 1
else:
k = twins_k[rand]
k += 1

# If l is equal to the length of the time series at
# this point, j and k are twins
if l == N:
# And the twins to the twin list
twins_k = twins[k]
# If the new k >= n_time, choose a new random starting point in the
# original time series
if k >= N:
while True:
new_k = int(floor(drand48() * N))
if new_k != k:
break

twins_j.append(k)
twins_k.append(j)
k = new_k

break
j += 1
64 changes: 3 additions & 61 deletions pyunicorn/timeseries/recurrence_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@
_bootstrap_distance_matrix_euclidean, _bootstrap_distance_matrix_supremum,\
_diagline_dist_norqa_missingvalues, _diagline_dist_norqa, \
_diagline_dist_rqa_missingvalues, _diagline_dist_rqa, _rejection_sampling,\
_white_vertline_dist, _twins
_white_vertline_dist, _twins, _twin_surrogates

#
# Class definitions
Expand Down Expand Up @@ -1562,64 +1562,6 @@ def twin_surrogates(self, n_surrogates=1, min_dist=7):
# Initialize
surrogates = np.empty((n_surrogates, N, dim))

code = r"""
int i, j, k, l, new_k, n_twins, rand;
// Initialize random number generator
srand48(time(0));
for (i = 0; i < n_surrogates; i++) {
// Randomly choose a starting point in the original trajectory
k = floor(drand48() * N);
j = 0;
while (j < N) {
// Assign state vector of surrogate trajectory
for (l = 0; l < dim; l++) {
surrogates(i,j,l) = embedding(k,l);
}
// Get the list of twins of state vector k in the original
// time series
py::list twins_k = PyList_GetItem(twins,k);
// Get the number of twins of k
n_twins = PyList_Size(twins_k);
// If k has no twins, go to the next sample k+1. If k has
// twins at m, choose among m+1 and k+1 with equal probability
if (n_twins == 0)
k++;
else {
// Generate a random integer between 0 and n_twins
rand = floor(drand48() * (n_twins + 1));
// If rand = n_twins go to sample k+1, otherwise jump to
// the future of one of the twins
if (rand == n_twins)
k++;
else {
k = twins_k[rand];
k++;
}
}
// If the new k >= n_time, choose a new random starting point
// in the original time series
if (k >= N) {
do {
new_k = floor(drand48() * N);
}
while (k == new_k);
k = new_k;
}
j++;
}
}
"""
weave_inline(locals(), code,
['n_surrogates', 'N', 'dim', 'twins', 'embedding',
'surrogates'])
_twin_surrogates(n_surrogates, N, dim, twins, embedding,
surrogates)
return surrogates

0 comments on commit 098489b

Please sign in to comment.