Signal basis of of ALSSM and transformed (orthonormal) space [ex116.0]

Description

Out:

[[21.]]
[[770.]]
[[50666.]]
[[3956810.]]
[[0.99255632]]
[[1.01642477]]
[[1.01222752]]
[[1.01153216]]
/Users/wbr8/homedata/projects/lmlib-branches/reto-doc/examples-teaching/mse/fig-examp/fig-examp-feature-space-bases.py:99: UserWarning: This figure was using constrained_layout==True, but that is incompatible with subplots_adjust and or tight_layout: setting constrained_layout==False.
  plt.subplots_adjust(left=0.08, right=0.96, top=0.96, bottom=0.06, hspace=0.5)

from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import matrix_power
from scipy.signal import find_peaks

import lmlib as lm
from lmlib.utils.generator import gen_wgn, gen_sine, gen_rand_walk, gen_exponential

np.random.seed(0)

# Feature space parameters
use_offset = True
draw_trajs = True
draw_clusters = True
show_dimensions_3d = [1,2,3]
use_pca = True

# Signal
K = 2000
ks = [40, 120, 310, 480, 610, 850, 920, 1100, 1250, 1320, 1550, 1700, 1910]
spikes_len = [60, 45, 50, 60, 45, 80, 60, 50, 81, 45, 60, 50, 45]
k0_sin = [0, -2, 5, 0, -2, 1, 0, 4, 1, -2, 0, 6, -2]
matches = [0, 1, 2, 0, 1, 3, 0, 2, 3, 1, 0, 2, 1]
if (draw_clusters):
    colors = ['b', 'g', 'r', 'orange']
else:
    colors = ['b', 'b', 'b', 'b']



NN = [1,0.8,.5,.2]

# Model Parameters
a = -10 #
b = 10 #
g_spike = 50
order_spike = 4

y_baseline = 0.2*gen_sine(K, K // 2) * gen_sine(K, K // 3) + 0.01 * gen_rand_walk(K, seed=0)
y_signal = np.zeros(K)

s_len = (b-a)
for ik, (k, k0, match) in enumerate(zip(ks, k0_sin, matches)):
    k0 = 0
    nn = NN[match]
    y_signal[k+a:k+b] = gen_sine(s_len, s_len / nn , 1, [k0]) * gen_exponential(s_len, ((s_len - 1) / s_len) ** 5)
y = y_signal + 0.0 * y_baseline + gen_wgn(K, sigma=0.01, seed=0)

# Defining ALSSM models
alssm_spike = lm.AlssmPoly(poly_degree=order_spike - 1, label="alssm-pulse")

# Defining segments with a left- resp. right-sided decaying window and a center segment with nearly rectangular window
segmentC = lm.Segment(a=a, b=b, direction=lm.FORWARD, g=g_spike, label='center-segment')

# Defining the final cost function (a so called composite cost = CCost)
# mapping matrix between models and segments (rows = models, columns = segments)
F = [[1]]
ccost = lm.CompositeCost((alssm_spike, ), ( segmentC, ), F)



# Filter
separam = lm.SEParam(ccost)
xs_hat = separam.filter_minimize_x(y)

V = separam.get_V()

xs_pulse = xs_hat[:, 0:order_spike]
zs = np.einsum('mn, kn->km', V, xs_pulse)

peaks, _ = find_peaks(np.convolve(y**2, 1/20*np.ones(20), 'same'), height=0.05, distance=20)

# Trajectories
trajs_pulse = lm.map_trajectories(ccost.trajectories(xs_hat[ks],  F=[[1,]], thds=0.001), ks, K, merge_ks=False, merge_seg=True)


bases = ccost.trajectories(np.eye(order_spike),  F=[[1,]], thds=0.001)
bases_V = ccost.trajectories((np.linalg.inv(np.transpose(V)))@np.eye(order_spike),  F=[[1,]], thds=0.001)


for i in np.arange(order_spike):
    print(np.transpose(bases[i][0][1])@bases[i][0][1])

for i in np.arange(order_spike):
    print(np.transpose(bases_V[i][0][1])@bases_V[i][0][1])



# Plot
fig = plt.figure(figsize=(7, 4), constrained_layout=True)
plt.subplots_adjust(left=0.08, right=0.96, top=0.96, bottom=0.06, hspace=0.5)

spec = fig.add_gridspec(3, 2)
ax0 = fig.add_subplot(spec[0, :])
ax1 = fig.add_subplot(spec[1:3, 0])
ax2 = fig.add_subplot(spec[1:3, 1])


# plot signal
ax0.plot(range(K), y, lw=0.5, c='k', label='observation')
if (draw_trajs):
    for traj, m in zip(trajs_pulse, matches):
        ax0.plot(range(K), traj, lw=1.0, c=colors[m], label='Pulse')
ax0.scatter(ks, y[ks], marker=7, c='k')
ax0.set(xlabel='k', ylabel=r'$y$')
# ax0.legend(loc=1)

# plot bases
NOF_BV = len(bases) # number of base vectors
#  ax1.plot(range(K), y, lw=0.5, c='k', label='observation')
for ib in range(NOF_BV):
   ax1.plot(bases[ib][0][0], bases[ib][0][1], '-o', linewidth = 0.5, markersize=4,  label='x: v'+str(ib))
#  ax1.scatter(bases[ib][0][0], bases[ib][0][1], lw=1.0,  label='')

ax1.legend(loc=4)

# ax1.set_ylim(ax1.get_ylim()[0]/4, ax1.get_ylim()[1]/4)
ax1.set(xlabel='i', ylabel=r'$x$')


# plot bases
NOF_BV = len(bases_V) # number of base vectors
#  ax1.plot(range(K), y, lw=0.5, c='k', label='observation')
for ib in range(NOF_BV):
   ax2.plot(bases_V[ib][0][0], bases_V[ib][0][1], '-o', linewidth = 0.5, markersize=4,  label='z: v'+str(ib))
#   ax2.scatter(bases_V[ib][0][0], bases_V[ib][0][1], lw=1.0,  label='')

ax2.legend(loc=4)


# ax2.set_ylim(ax1.get_ylim()[0]/4, ax1.get_ylim()[1]/4)
ax2.set(xlabel='i', ylabel=r'$z$')




plt.show()

Detection of QRS wave onset in ECG signal [ex902.0] Feature Space Transformation [ex115.0]