In [5]:
from __future__ import division, absolute_import
%matplotlib inline
import numpy as np
from numpy import log
import matplotlib.pyplot as plt
from IPython.html.widgets import *
# from wtmm import wtmm, perform_cwt, skeletor
from wtmm import skeletor
import scipy.ndimage
from collections import OrderedDict
from scipy import signal
from numpy import log
In [6]:
def forgery(Iterations=0, Multifractal=1):
if Multifractal:
turns=((0.25,0.5),(0.75, 0.25))
else:
turns=((0.4,0.6),(0.6, 0.4))
first_turn, second_turn = turns
ys = [0,1]
ts = [0,1]
for i in range(0, Iterations + 1):
j=0
while ts[j] < 1:
dt = ts[j+1] - ts[j]
dy = ys[j+1] - ys[j]
ts.insert(j+1, ts[j] + first_turn[0]*dt)
ts.insert(j+2, ts[j] + second_turn[0]*dt)
ys.insert(j+1, ys[j] + first_turn[1]*dy)
ys.insert(j+2, ys[j] + second_turn[1]*dy)
j += 3
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={'axisbg':'#EEEEEE', 'axisbelow':True})
ax.grid(color='w', linewidth=2, linestyle='solid')
ax.plot(ts, ys, color='b', alpha=0.4)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
return fig
def forgery_prices(Iterations=10, Multifractal=1):
if Multifractal:
turns=((0.25,0.5),(0.75, 0.25))
else:
turns=((0.4,0.6),(0.6, 0.4))
first_turn, second_turn = turns
ys = [0,1]
ts = [0,1]
for i in range(0, Iterations + 1):
j=0
while ts[j] < 1:
dt = ts[j+1] - ts[j]
dy = ys[j+1] - ys[j]
ts.insert(j+1, ts[j] + first_turn[0]*dt)
ts.insert(j+2, ts[j] + second_turn[0]*dt)
ys.insert(j+1, ys[j] + first_turn[1]*dy)
ys.insert(j+2, ys[j] + second_turn[1]*dy)
j += 3
return np.array(ts), np.array(ys)
cartoon_t, cartoon_s = forgery_prices(5, True)
print("cartoon's signal is {} in length".format(len(cartoon_s)))
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={'axisbg':'#EEEEEE', 'axisbelow':True})
ax.grid(color='w', linewidth=2, linestyle='solid')
ax.plot(cartoon_t, cartoon_s, color='b', alpha=0.4)
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
plt.show()
In [37]:
def perform_cwt(sig, width_step=0.5, max_scale=None, wavelet=signal.ricker, epsilon=0.1, order=1, mask_result=True, plot=True):
"""
Perform the continuous wavelet transform against the incoming signal. This function will normalize the signal
(to 0-1 in the y axis) for you, as well as taking the -1 * abs( log( ) ) of the matrix that is found. Literature
suggests that len/4 is a good balance for finding the bifurcations vs execution time
This will automatically create the maxima only mask of the wavelet coef matrix for you. To see the original, use
plot=True
:param sig: 1 dimensional array -- the signal to be hit with the wavelet
:param width_step: what width step to use between the min and the max
:param max_scale: the maximum scale to use. Default = len(sig)/4
:param wavelet: what wavelet to use as the mother
:param epsilon: how to score the maxima's intensity (e.g. intensity / epsilon )
:param order: how many neighbors to look at when finding the local maxima
:param plot: whether to plot the original CWT coefficient matrix as a heatmap
:return: the mask, see above
"""
if not max_scale:
max_scale = len(sig) / 4
widths = np.arange(1, max_scale, width_step)
# normalize the signal to fit in the wavelet
sig_max = sig.max()
sig_min = sig.min()
sig = (sig - (sig_min - 0.01)) / (sig_max - sig_min + 0.02)
# Run the transform
w_coefs = abs(-1 * log(abs(signal.cwt(sig, wavelet, widths))))
# Create the mask, keeping only the maxima
if mask_result:
mask = _create_w_coef_mask(w_coefs, epsilon=epsilon, order=order)
else:
mask = w_coefs
if plot:
plt.figure(figsize=(14, 10))
plt.pcolormesh(w_coefs)
plt.colorbar()
ax = plt.gca()
ax.set_ylim(ax.get_ylim()[::-1])
ax.xaxis.tick_top()
plt.show()
return mask
def _create_w_coef_mask(w_coefs, epsilon=0.1, order=1):
"""
Create a new matrix, the same shape as the wavelet coefficient one, but with zeros everywhere except for local
maxima's. Epsilon here is used for ranking the strength of the local maxima.
Assumes that the coefficient matrix coming in is already in absolute terms
:param w_coefs: wavelet coefficient matrix
:param epsilon: divided against the maxima, used for transparent ranking
:param order: how many neighboors on a given row to look at to determine maxima
:return: same shape array, see above
"""
mask = np.zeros_like(w_coefs, dtype=int)
for n, row in enumerate(w_coefs):
maxs = signal.argrelmax(row, order=order)[0]
mask[n, maxs] = row[maxs] / epsilon
return mask
In [69]:
%%time
sig = cartoon_s
sig_t = cartoon_t
w_coefs = perform_cwt(sig, width_step=0.25, max_scale=100)
In [64]:
plt.figure(figsize=(14,10))
plt.pcolormesh(w_coefs, cmap='Greys')
plt.colorbar()
ax = plt.gca()
ax.set_ylim(ax.get_ylim()[::-1])
ax.xaxis.tick_top()
plt.show()
In [75]:
data = np.copy(w_coefs)
In [74]:
bifucations = skeletor(data, smallest_scale=1, proximity=9)
plt.figure(figsize=(14,10))
for n, (k, v) in enumerate(bifucations.iteritems()):
rows, cols = zip(*v)
plt.plot(cols, rows)
ax = plt.gca()
ax.set_ylim(ax.get_ylim()[::-1])
ax.xaxis.tick_top()
plt.show()
In [ ]: