In [ ]:

    

import numpy as np
import pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt

from crossing_tree import fbm, crossings, crossing_tree
from sklearn.preprocessing import StandardScaler


    

In [ ]:

    

random_state = np.random.RandomState(0x0BADC0D)


    

In [ ]:

    

gen_ = fbm(2**15+1, H=0.60, time=True)
# gen_ = fbm(2**22+1, H=0.95, time=True)
# gen_ = fbm(2**22+1, H=0.90, time=True)
gen_.initialize(random_state)


    

In [ ]:

    

T, X = gen_()


    

In [ ]:

    

print(np.std(np.diff(X)))
print(np.mean(np.abs(np.diff(X))))
print(np.median(np.abs(np.diff(X))))
print(np.percentile(np.abs(np.diff(X)), 95))
print(1.0 / np.sqrt(len(T) - 1))


    

In [ ]:

    

# scale_ = np.diff(X).std()
# scale_ = np.median(np.abs(np.diff(X)))
scale_ = np.mean(np.abs(np.diff(X))) * 1024 # * 16384


    

In [ ]:

    

# %%timeit -n 20
xi, ti = crossings(X, T, scale_, 0)
xi2, ti2 = crossings(X, T, scale_, 0.5)

print(xi.shape)
print(ti.shape)

fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)
ax.plot(T[:(1<<20)], X[:(1<<20)], "y", alpha=0.5)
ax.plot(ti[:11], xi[:11], "-k")
ax.plot(ti2[:11], xi2[:11], "-r")
ax.set_title("Scale %g"%(scale_,))
# scale_ *= 2
plt.show()


    

In [ ]:

    

from scipy.stats import ks_2samp


    

In [ ]:

    

# scale_ = np.median(np.abs(np.diff(X)))
scale_ = np.diff(X).std()
scale_ = 1.96 / np.sqrt( len(T)-1)
# scale_ = np.median(np.abs(np.diff(X)))
# scale_ = np.mean(np.abs(np.diff(X)))
for k in range(10):
    xi, ti = crossings(X, T, scale_, 0)
    xi2, ti2 = crossings(X, T, scale_, 0.25)
    scale_ *= 2

    dti_s = np.diff(ti)
    dti_s /= np.sqrt(np.mean(dti_s**2))
    
    dti2_s = np.diff(ti2)
    dti2_s /= np.sqrt(np.mean(dti2_s**2))

    fig = plt.figure(figsize=(12, 5))
    ax = fig.add_subplot(121)
    ax.hist(dti_s, bins=100, color="k", lw=0);
    ax = fig.add_subplot(122)
    ax.hist(dti2_s, bins=100, color="r", lw=0);
    plt.show()
    print(ks_2samp(dti_s, dti2_s))


    

In [ ]:

    

# scale_ = np.median(np.abs(np.diff(X)))
scale = np.diff(X).std()/2
# scale = 1.96 / np.sqrt( len(T)-1)
# scale = np.median(np.abs(np.diff(X)))
# scale = np.mean(np.abs(np.diff(X)))
origin = 0.

xi, ti, offspring, excursions, subcrossings, durations = crossing_tree(X, T, scale, origin)
print(len(offspring))


    

In [ ]:

    

def structure(offspring):
    iter_ = iter(offspring)
    try:
        value_ = next(iter_)
    except StopIteration:
        raise TypeError('reduce() of empty sequence')
    yield value_
    for index in iter_:
        value_ = value_[index]
        yield value_

list(structure(offspring))


    

In [ ]:

    

%load_ext cython


    

In [ ]:

    

%%cython
#-a
import numpy as _np
cimport numpy as np
cimport cython

from libc.math cimport isnan, fabs, floor, ceil, NAN

np.import_array()

ctypedef fused real:
    cython.floating

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
@cython.initializedcheck(False)
@cython.nonecheck(False)
def integer_xing(real[:] x, real[:] t, real scale, real origin):
    cdef np.intp_t n_samples = x.shape[0]
    cdef np.intp_t[::1] size = _np.empty(n_samples - 1, dtype=_np.int)
    cdef real[::1] first = _np.empty(n_samples - 1, dtype=_np.float)
    
    cdef np.intp_t i
    cdef real first_, last_, direction, prev_last = NAN
    cdef np.intp_t total = 0, size_
    with nogil:
        # Detect integer-level crossings, ignoring re-crossings of the same level
        for i in range(n_samples - 1):
            direction, size_ = 0.0, 0
            if x[i] < x[i+1]:
                first_, last_ = ceil((x[i] - origin) / scale), floor((x[i+1] - origin) / scale)
                direction = +1.0
            elif x[i] > x[i+1]:
                first_, last_ = floor((x[i] - origin) / scale), ceil((x[i+1] - origin) / scale)
                direction = -1.0
            if direction != 0.0:
                size_ = <int>fabs(last_ + direction - first_)
                if size_ > 0 and prev_last == first_:
                    first_ += direction
                    size_ -= 1
                if size_ > 0:
                    prev_last = last_
            first[i], size[i] = first_, size_
            total += size_

    cdef real[::1] xi = _np.empty(total, dtype=_np.float)
    cdef real[::1] ti = _np.empty(total, dtype=_np.float)

    cdef np.int_t k, j = 0
    cdef long double x_slope_, t_slope_, first_xi_, first_ti_
    with nogil:
        # Interpolate the crossing times and crossing levels
        for i in range(n_samples-1):
            size_ = size[i]
            if size_ > 0:
                x_slope_ = +scale if x[i+1] > x[i] else -scale
                t_slope_ = (t[i+1] - t[i]) / (x[i+1] - x[i])
                first_ = first[i] * scale + origin
                for k in range(size_):
                    xi[j] = first_ + x_slope_ * k
                    ti[j] = t[i] + t_slope_ * (xi[j] - x[i])
                    j += 1
## Marginally slower
#             size_ = size[i]
#             if size_ > 0:
#                 t_slope_ = (t[i+1] - t[i]) / (x[i+1] - x[i])
#                 xi[j] = first[i] * scale
#                 ti[j] = t[i] + t_slope_ * (xi[j] - x[i])
#                 j += 1
#             if size_ > 1:
#                 x_slope_ = +scale if x[i+1] > x[i] else -scale
#                 for k in range(size_ - 1):
#                     xi[j] = xi[j-1] + x_slope_
#                     ti[j] = ti[j-1] + t_slope_ * x_slope_
#                     j += 1
    return xi, ti