calculate synthetic ramps given a series of time points and coefficients of unit steps at each time. \n convolution ramps

synthetic_ramp_method.ipynb

synthetic_ramp_method.py

# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: -all
#     formats: ipynb,py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.16.7
#   kernelspec:
#     display_name: main
#     language: python
#     name: python3
# ---

# %% [markdown]
# ## synthetic ramp method
#
# calculate synthetic ramps using convolution for a given series of time points:
# * define matrix $A$ with the coefficients $a_{i,j}$ of the unit steps in $n$ different series $i$ for each time $t_j$.
# * define matrix $\sigma(t-t_j)*\sigma(t)=(t-t_j)\cdot \sigma(t-t_j)$ which is lower triangular with permutations of $(t-t_j)\cdot \sigma(t-t_j)$
#     * example with $t_0<t_1$: $(t_0-t_1)\cdot \sigma(t_0-t_1)=0$, and $(t_1-t_0)\cdot \sigma(t_1-t_0)=(t_1-t_0)$ due to step definition
# * the matrix $\sigma(t-t_j)\cdot A$ has then the slopes $\frac{dy_i}{dt}$ at times $t_j$ and $(t-t_j)\cdot \sigma(t-t_j)\cdot A$ has the synthetic ramps by convolution.
# $$
# \underbrace{\left[\begin{array}{cccc}
# 1 & 0 & 0 & ...\\
# 1 & 1 & 0 & ...\\
# 1 & 1 & 1 & ...\\
# ... & ... & ... & ...
# \end{array}\right]}_{\sigma(t-t_j)}\cdot 
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,1} & a_{2,1} & ...\\
# a_{1,2} & a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{A} = 
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,0}+a_{1,1} & a_{2,0}+a_{2,1} & ...\\
# a_{1,0}+a_{1,1}+a_{1,2} & a_{2,0}+a_{2,1}a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{\mathrm{slopes}}
# $$
#
# $$
# \underbrace{\left[\begin{array}{cccc}
# t_0-t_0 & 0 & 0 & ...\\
# t_1-t_0 & t_1-t_1 & 0 & ...\\
# t_2-t_0 & t_2-t_1 & t_2-t_2 & ...\\
# ... & ... & ... & ...
# \end{array}\right]}_{(t-t_j)\cdot \sigma(t-t_j)}\cdot 
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,1} & a_{2,1} & ...\\
# a_{1,2} & a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{A} = 
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0}\cdot(t_0-t_0) & a_{2,0}\cdot(t_0-t_0) & ...\\
# a_{1,0}\cdot(t_1-t_0)+a_{1,1}\cdot(t_1-t_1) & a_{2,0}\cdot(t_1-t_0)+a_{2,1}\cdot(t_1-t_1) & ...\\
# a_{1,0}\cdot(t_2-t_0)+a_{1,1}\cdot(t_2-t_1)+a_{1,2}\cdot(t_2-t_2) & a_{2,0}\cdot(t_2-t_0)+a_{2,1}\cdot(t_2-t_1)+a_{2,2}\cdot(t_2-t_2) & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{\mathrm{synth\ curves}}
# $$
#
# example below with three series $i=\{1,2,3\}$

# %% [markdown]
#
# ### filtering to delay/advance signals according to one of the signals
#
# problem description: 
# * based on one of the signals as key signal: $y_a$
# * when $y_a$ increases, other signals should be delayed by a fixed time $\delta t$
# * when $y_a$ decreases, all other signals should be advanced by a fixed time $\delta t$
#
# solution: use filtering to obtain adapted $\sigma(t-t_j\pm\delta t)*\sigma(t)=\sigma(t-t_j\pm\delta t)\cdot (t-t_j\pm\delta t)$ , where 
# * $\pm\delta t$ is negative (delaying $\sigma$) whenever the slope (and delayed slope) of $y_a$ is positive
# * $\pm\delta t$ is positive (advancing $\sigma$) whenever the slope (and advanced slope) of $y_a$ is negative
#
# This ensures positive steady states are reached including delay, or negative steady states are reached including advancement.
#
# Method:
# 1. extend times vector in every time $t_j$ by the next time offsets by time-shift $\pm \delta t$, thereby triplicating its size
# >(TODO: save space by extending only based on actual derivatives of $y_a$)
# 2. calculate slopes for extended times vector as $\sigma(t-t_j)\cdot A$, delayed slopes as $\sigma(t-t_j-\delta t)\cdot A$ and advanced slopes $\sigma(t-t_j+\delta t)\cdot A$
# 3. based on slopes obtained, calculate elementwise the matrices $\sigma(t-t_j\pm\delta t)$ and $(t-t_j\pm\delta t)\cdot\sigma(t-t_j\pm\delta t)$, determining the sign of $\pm \delta t$ by each corresponding slope
# 4. filtered slopes are $\sigma(t-t_j\pm\delta t)\cdot A$ and filtered ramps are $(t-t_j\pm\delta t)\cdot\sigma(t-t_j\pm\delta t)\cdot A$ as above

# %%
from numpy import array,linspace,zeros,pi,exp
from matplotlib import pyplot as plt

steps_times=array([
    [0,0,0,0],
    [3,0,-1,0],
    [4,0,1,0],
    [6,3/5,0,0],
    [8,0,2,3/5],
    [8.5,0,-2,0],
    [11,-3/5,0,0],
    [14,-1.5/2,0,-3/5],
    [16,-1.5/2,0.5,-3/4],
    [17,1.5,-0.5,-3/4],
    [21,0,0,+3/2],
])
t_shift=1
times=steps_times[:,0]
steps=steps_times[:,1:]
initial_vals=array([2,1,2])

eye=array([[1 if times[i]-times[j]==0 else 0 for j in range(times.shape[0])] for i in range(times.shape[0])])
sigma=array([[1 if times[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times.shape[0])])
sigma_t=array([[(times[i]-times[j]) if times[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times.shape[0])])

slopes_min=sigma.dot(steps)
synth_ramps_min=initial_vals+sigma_t.dot(steps)

times_extended_min=array(list(set(list(times)+[x for x in list(times-t_shift)+list(times+t_shift)]))) # include all unique points shifted by +/- t_shift
times_extended_min.sort()
eye=array([[1 if times_extended_min[i]-times[j]==0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
eye_shifted=array([[1 if times_extended_min[i]-times[j]-t_shift==0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
sigma=array([[1 if times_extended_min[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
sigma_t=array([[(times_extended_min[i]-times[j]) if times_extended_min[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
sigma_shifted=array([[1 if times_extended_min[i]-times[j]-t_shift>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])

slopes=sigma.dot(steps)
synth_ramps=initial_vals+sigma_t.dot(steps)
shifted_steps=eye_shifted.dot(steps)
shifted_slopes=sigma_shifted.dot(steps)

eye_shifted_neg=array([[1 if times_extended_min[i]-times[j]+t_shift==0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
sigma_shifted_neg=array([[1 if times_extended_min[i]-times[j]+t_shift>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended_min.shape[0])])
shifted_steps_neg=eye_shifted_neg.dot(steps)
shifted_slopes_neg=sigma_shifted_neg.dot(steps)

sigma_shifted_func=zeros([times_extended_min.shape[0],times.shape[0]])
sigma_t_shifted_func=zeros([times_extended_min.shape[0],times.shape[0]])
eye_shifted_func=zeros([times_extended_min.shape[0],times.shape[0]])
for i in range(times_extended_min.shape[0]):
    for j in range(times.shape[0]):
        t_shift_func=0
        if (shifted_slopes[i,0]>0) | (slopes[i,0]>0):
            t_shift_func=+t_shift
        elif (slopes[i,0]<0) | (shifted_slopes_neg[i,0]<0):
            t_shift_func=-t_shift
        sigma_shifted_func[i,j]=1 if (times_extended_min[i]-times[j]-t_shift_func>=0) else 0
        eye_shifted_func[i,j]=1 if (times_extended_min[i]-times[j]-t_shift_func==0) else 0 # doesn't work to produce steps to integrate, some duplication of steps
        sigma_t_shifted_func[i,j]=(times_extended_min[i]-times[j]-t_shift_func) if (times_extended_min[i]-times[j]-t_shift_func>=0) else 0

#shifted_steps_func=eye_shifted_func.dot(steps) # doesn't work to produce steps to integrate, some duplication of steps
shifted_slopes_func=sigma_shifted_func.dot(steps)
shifted_synth_ramps=initial_vals+sigma_t_shifted_func.dot(steps)

eye_ext=array([[1 if times_extended_min[i]-times_extended_min[j]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended_min.shape[0])])
eye_ext_shift1=array([[1 if times_extended_min[i-1]-times_extended_min[j]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended_min.shape[0])])
#eye_ext_shift1=concatenate([[zeros(times_extended_min.shape[0])],eye_ext_shift1])
sigma_ext=array([[1 if times_extended_min[i]-times_extended_min[j]>=0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended_min.shape[0])])
sigma_t_ext=array([[(times_extended_min[i]-times_extended_min[j]) if times_extended_min[i]-times_extended_min[j]>=0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended_min.shape[0])])

shifted_steps_func=(eye_ext-eye_ext_shift1).dot(sigma_shifted_func.dot(steps))
reintegrated_ramps_from_steps=initial_vals+sigma_t_ext.dot((eye_ext-eye_ext_shift1).dot(sigma_shifted_func.dot(steps)))

idx_output = (shifted_steps_func!=0).any(axis=1)

times_extended=array(list(set(linspace(times.min(),times.max(),1000).tolist()+times_extended_min.tolist())))
times_extended.sort()

sigma_extended=array([[1 if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
sigma_t_extended=array([[(times_extended[i]-times[j]) if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
slopes_extended=sigma_extended.dot(steps)

sigma_extended_shifted=array([[1 if times_extended[i]-times[j]-t_shift>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
sigma_extended_shifted_neg=array([[1 if times_extended[i]-times[j]+t_shift>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])

slopes_extended_shifted=sigma_extended_shifted.dot(steps)
slopes_extended_shifted_neg=sigma_extended_shifted_neg.dot(steps)

eye=array([[1 if times_extended[i]-times[j]==0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
sigma=array([[1 if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
slopes_extended_shifted_func=sigma.dot(steps)


# increase size of obtained arrays from times_extended_min to times_extended
eye_ext=array([[1 if times_extended[i]-times_extended_min[j]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended.shape[0])])
eye_ext_shift1=array([[1 if times_extended[i]-times_extended_min[j-1]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended.shape[0])])
sigma_ext=array([[1 if times_extended[i]-times_extended[j]>=0 else 0 for j in range(times_extended.shape[0])] for i in range(times_extended.shape[0])])
sigma_t_ext=array([[(times_extended[i]-times_extended[j]) if times_extended[i]-times_extended[j]>=0 else 0 for j in range(times_extended.shape[0])] for i in range(times_extended.shape[0])])

shifted_steps_ext_func=(eye_ext_shift1-eye_ext).dot(sigma_shifted_func.dot(steps))
shifted_slopes_ext_func=sigma_ext.dot(shifted_steps_ext_func)
synth_ramps_ext_func=initial_vals+sigma_t_ext.dot(shifted_steps_ext_func)

print('lower triangular (t-tj)sigma(t-tj)) \n',sigma_t)
fig,ax_list=plt.subplots(nrows=5,ncols=1,height_ratios=[1,1/2,1/2,1/2,1/2],constrained_layout=True,sharex=True,figsize=[9,9])
for j in range(steps.shape[1]):
    l1,=ax_list[0].plot(times,synth_ramps_min[:,j],marker=['o','<','s'][j],label=f'series i={j+1}')
    ax_list[0].plot(times_extended_min[idx_output],shifted_synth_ramps[idx_output,j],'-.',marker=['o','<','s'][j],color=l1.get_color(),fillstyle='none',label=f'series i={j+1} (filtered i=1)')
    markerline, stemlines, baseline=ax_list[1].stem(times,steps[:,j],markerfmt=['o','<','s'][j],linefmt=['blue','orange','green'][j],label=f'series i={j+1}')
    stemlines.set_linestyle('-');markerline.set_markerfacecolor('none')
    markerline, stemlines, baseline=ax_list[3].stem(times_extended_min[idx_output],shifted_steps_func[idx_output,j],markerfmt=['o','<','s'][j],linefmt=['blue','orange','green'][j],label=f'series i={j+1}')
    stemlines.set_linestyle('-');markerline.set_markerfacecolor('none')
    ax_list[2].plot(times_extended,slopes_extended[:,j],linestyle=['-','--','-.'][j],label=f'series i={j+1}')

ax_list[0].plot(times_extended_min[idx_output],reintegrated_ramps_from_steps[idx_output,0],':',marker='none',color='red',fillstyle='none',label=f'series i={0+1} reconstructed')
ax_list[4].plot(times_extended,slopes_extended[:,0],['-','--','-.'][0],label=f'series i={0+1}, slope')
ax_list[4].plot(times_extended,slopes_extended_shifted[:,0],['-','--','-.'][1],label=f'series i={0+1}, delay')
ax_list[4].plot(times_extended,slopes_extended_shifted_neg[:,0],['-','--','-.'][2],label=f'series i={0+1}, adv')
ax_list[4].plot(times_extended+t_shift,shifted_slopes_ext_func[:,0],['-','--','-.',':'][3],label=f'series i={0+1}, func')
for j in range(len(times)):
    ax_list[0].axvline(times[j],linestyle='--',color=[0.85,0.85,0.85],zorder=-1)
    ax_list[1].axvline(times[j],linestyle='--',color=[0.85,0.85,0.85],zorder=-1)
    ax_list[2].axvline(times[j],linestyle='--',color=[0.85,0.85,0.85],zorder=-1)
for j in times_extended_min[idx_output]:
    ax_list[4].axvline(j,linestyle='--',color=[0.6,0.6,0.6],zorder=-1)
    ax_list[3].axvline(j,linestyle='--',color=[0.6,0.6,0.6],zorder=-1)
    ax_list[0].axvline(j,linestyle='--',color=[0.6,0.6,0.6],zorder=-1)

ax_list[1].plot()

ax_list[-1].set_xlabel('time $t_j$')
ax_list[0].set_ylabel('$y_i$')
ax_list[1].set_ylabel('$a_{i,j}$')
ax_list[2].set_ylabel(r'$\frac{dy_i}{dt}$')
ax_list[3].set_ylabel('$a_{i,j}$ filtered')
ax_list[0].legend(loc='center left',bbox_to_anchor=[1.01,0.5]);
ax_list[1].legend(loc='center left',bbox_to_anchor=[1.01,0.5])
ax_list[2].legend(loc='center left',bbox_to_anchor=[1.01,0.5])
ax_list[3].legend(loc='center left',bbox_to_anchor=[1.01,0.5])
ax_list[4].legend(loc='center left',bbox_to_anchor=[1.01,0.5]);

# %%
# %matplotlib widget
sigma=array([[1 if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
eye=array([[1 if times_extended[i]-times_extended_min[j]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended.shape[0])])
eye_shifted1=array([[1 if times_extended[i]-times_extended_min[j-1]==0 else 0 for j in range(times_extended_min.shape[0])] for i in range(times_extended.shape[0])])
#eye_shifted1=concatenate([[zeros(times_extended.shape[0])],eye_shifted1.T]).T

(eye_shifted1-eye)
plt.figure();plt.plot(times_extended,sigma.dot(steps))
#plt.plot(times_extended,eye.dot(sigma_shifted_func.dot(steps)))
#plt.plot(times_extended_min,sigma_shifted_func.dot(steps))
#plt.plot(times_extended,eye.dot(sigma_shifted_func.dot(steps)))
sigma=array([[1 if times_extended[i]-times_extended[j]>=0 else 0 for j in range(times_extended.shape[0])] for i in range(times_extended.shape[0])])
sigma_t=array([[(times_extended[i]-times_extended[j]) if times_extended[i]-times_extended[j]>=0 else 0 for j in range(times_extended.shape[0])] for i in range(times_extended.shape[0])])

shifted_steps_ext_func=(eye_shifted1-eye).dot(sigma_shifted_func.dot(steps)) # concatenate([[[0,0,0]],(eye_shifted1-eye).dot(sigma_shifted_func.dot(steps))[:-1,:]])#
shifted_slopes_ext_func=sigma.dot(shifted_steps_ext_func)
synth_ramps_ext_func=initial_vals+sigma_t.dot(shifted_steps_ext_func)
for j in [0]:
    markerline, stemlines, baseline=plt.stem(times_extended,shifted_steps_ext_func[:,j],linefmt=['blue','orange','green'][j])
    stemlines.set_linestyle('-');markerline.set_markerfacecolor('none')
    plt.plot(times_extended,shifted_slopes_ext_func[:,j])
    plt.plot(times_extended,synth_ramps_ext_func[:,j])

# %% [markdown]
# ## low-pass filter

# %%
n_max=120*2
t=linspace(0,2*pi*5*4,n_max)
T=2*pi*5*4/(n_max-1)
t_shift=T
omega0=1/(30*T)
beta=exp(-omega0*T)

steps=array([[-1 if j==30 or j==50 else 0] for j in range(n_max)])-array([[-1 if j==30+1 or j==50+1 else 0] for j in range(n_max)])
sigma=array([[1 if t[i]-t[j]>=0 else 0 for j in range(t.shape[0])] for i in range(t.shape[0])])
sigma_t=array([[(t[i]-t[j]) if t[i]-t[j]>=0 else 0 for j in range(t.shape[0])] for i in range(t.shape[0])])

synth_ramps=sigma_t.dot(steps)
slopes=sigma.dot(steps)


sigma_t_shifted=array([[(t[i]-t[j]-t_shift) if t[i]-t[j]-t_shift>=0 else 0 for j in range(t.shape[0])] for i in range(t.shape[0])])
synth_ramps_shifted=sigma_t_shifted.dot(steps)

synth_ramps_filtered=beta*synth_ramps_shifted+(1-beta)*synth_ramps
#synth_ramps_filtered=synth_ramps*(1-exp(-omega0*t.reshape([n_max,1])))


fig,ax_list=plt.subplots(nrows=3,ncols=1,constrained_layout=True,sharex=True,figsize=(7,7))
for j in range(steps.shape[1]):
    l1,=ax_list[0].plot(t,synth_ramps[:,j],marker=['o','<','s'][j],label=f'series i={j+1}')
    ax_list[0].plot(t,synth_ramps_shifted[:,j],'-.',marker=['o','<','s'][j],fillstyle='none',label=f'series i={j+1} (shifted i=1)')
    ax_list[0].plot(t,synth_ramps_filtered[:,j],'-.',marker=['o','<','s'][j],fillstyle='none',label=f'series i={j+1} (filtered i=1)')
    markerline, stemlines, baseline=ax_list[1].stem(t,steps[:,j],markerfmt=['o','<','s'][j],linefmt=['blue','orange','green'][j],label=f'series i={j+1}')
    stemlines.set_linestyle('-');markerline.set_markerfacecolor('none')
    ax_list[2].plot(t,slopes[:,j],linestyle=['-','--','-.'][j],label=f'series i={j+1}')


ax_list[-1].set_xlabel('time $t_j$')
ax_list[0].set_ylabel('$y_i$')
ax_list[1].set_ylabel('$a_{i,j}$')
ax_list[2].set_ylabel(r'$\frac{dy_i}{dt}$')
ax_list[0].legend(loc='center left',bbox_to_anchor=[1.01,0.5]);
ax_list[1].legend(loc='center left',bbox_to_anchor=[1.01,0.5])
ax_list[2].legend(loc='center left',bbox_to_anchor=[1.01,0.5])