Skip to content
Snippets Groups Projects
Commit 00b666c5 authored by FROGE Ewen's avatar FROGE Ewen
Browse files

Upload New File

parent 178efe68
Branches
No related tags found
No related merge requests found
#%%
# Import lib
import numpy as np
from scipy.ndimage import gaussian_filter1d
from scipy.signal import convolve
from MFsynthesis import synthmrw,synthMultiFractalField
import entropy.entropy as entropy # Assuming your TE library is imported
import entropy.tools as tools
import matplotlib.pyplot as plt
import os
import pickle
#%%
#Saving routine
save_dir = '/users2/local/e22froge/codes/TE_Filter/LinearShiftIncrement/Results'
os.makedirs(save_dir, exist_ok=True)
# Saving routine
def save_results(key, te_values, te_std):
""" Saves results to a file with the key components, including the rounded shift. """
shift_rounded = round(key[3], 2) # Ensure the shift is rounded to 2 decimal places
filename = f'{save_dir}/result_{key[0]}_{key[1]}_{key[2]}_shift_{shift_rounded}.pkl'
with open(filename, 'wb') as f:
pickle.dump((te_values, te_std), f)
# Loading routine
def load_results(key):
""" Loads results from a file with the key components, including the rounded shift. """
shift_rounded = round(key[3], 2) # Ensure the shift is rounded to 2 decimal places
filename = f'{save_dir}/result_{key[0]}_{key[1]}_{key[2]}_shift_{shift_rounded}.pkl'
if os.path.exists(filename):
with open(filename, 'rb') as f:
return pickle.load(f)
return None
#%%
# Define parameters, kernels and keys
entropy.set_sampling(Theiler=3,N_eff=10000,N_real=10)
# Configuration
N = 2**20
path = '/users2/local/e22froge/codes/Innovation/DataSignals/uspa_mob15_1.dat' # path to data
fid = open(path, 'r')
Modane = np.fromfile(fid, dtype='>f')[:N]
# Generate synthetic FBM signal (remove FGN if unnecessary)
fbm_signal,_ = synthmrw(dim=1, N=N, H=1/3, lambda2=0) # Placeholder for FBM
fbm_signal=fbm_signal[:,0]
fgn_signal = np.diff(fbm_signal, axis=0)
FBM_regu=synthMultiFractalField(dim=1, N=N, H=1/3, lambda2=0, regparam=[1, 2980])
FBM_regu=FBM_regu[:,0]
scales = [20, 30] # Define your scales
positions = ['causal', 'symmetric', 'anticausal']
signals=['fgn_signal', 'fbm_signal', 'Modane']
signals_regu=['fbm_signal']
max_lag_factor = 2 # We will vary lag up to 10 times the largest scale
# fbm_signal_flip=np.flip(fbm_signal)
# fgn_signal_flip=np.flip(fgn_signal)
# Modane_flip=np.flip(Modane)
# flipped_signals=['fgn_signal_flip','fbm_signal_flip','Modane_flip']
# Prepare results dictionary
results = {}
#%%
# Computation
def increment_kernel_linear_shift(scale, shift):
""" The shift is scaled between 0 and 1, where:
- shift = 0 means fully causal
-shift = 0.5 symetric
- shift = 1 means fully anticausal
"""
kernel = np.zeros(2 * scale + 1)
# Convert the shift (0 to 1) to an integer index shift based on the scale
shift_idx = int(shift * scale)
# Place +1 and -1 based on the scaled shift
kernel[scale - shift_idx] = 1 # +1 moves from kernel[scale] to kernel[0] as shift goes from 0 to 1
kernel[2*scale+1-1-shift_idx] = -1 # -1 moves from kernel[-1] to kernel[scale] as shift goes from 0 to 1
return kernel
# Now define the range of shifts you want to explore (e.g., from 0 to 1 in steps of 0.1)
shifts = np.linspace(0, 1, 11) # 11 steps from 0 to 1
# Example: Let's say we explore for different scales and compute TE for each shift
for sgn in signals_regu:
signal = eval(sgn)
for scale1 in scales:
for scale2 in scales:
if scale1 == scale2:
continue # Skip identical scales
for shift in shifts:
key = (sgn, scale1, scale2, round(shift, 2)) # Modify key to include shift rounded to 2 decimals
# Check if the result already exists
existing_result = load_results(key)
if existing_result:
results[key] = existing_result
continue # Skip computation if already saved
# Call the kernel function with the current shift for both scales
kernel1 = increment_kernel_linear_shift(scale1, shift)
kernel2 = increment_kernel_linear_shift(scale2, shift)
# Apply convolution with the shifted kernels
signal_scale1 = convolve(signal, kernel1, mode='same')
signal_scale2 = convolve(signal, kernel2, mode='same')
# Compute Transfer Entropy (TE)
max_lag = max(scale1, scale2) * max_lag_factor
te_values = []
te_std = []
for lag in range(1, max_lag + 1):
te_value = entropy.compute_TE(tools.reorder(signal_scale1), tools.reorder(signal_scale2), lag=lag, stride=lag)[0]
te_values.append(te_value)
te_std.append(entropy.get_last_info()[0])
# Save results immediately after computation
results[key] = te_values, te_std
save_results(key, te_values, te_std)
print(f'Saved result: {key}')
# %%
fig, axes = plt.subplots(6, 2, figsize=(16, 30))
fig.suptitle(f'Transfer Entropy for Different Shifts (fbm_signal)', fontsize=30)
sgn = 'fbm_signal'
scales = [20, 30]
for idx, shift in enumerate(shifts):
row, col = divmod(idx, 2) # Get the row and column for the subplot
ax = axes[row, col]
ax.set_title(f'Shift: {round(shift, 2)} - TE for Scales {scales[0]} and {scales[1]}', fontsize=18)
ax.set_xlabel('Lag', fontsize=16)
ax.set_ylabel('Transfer Entropy', fontsize=16)
key_1_to_2 = (sgn, scales[0], scales[1], round(shift, 2))
if key_1_to_2 in results:
te_values_1_to_2, te_std_1_to_2 = results[key_1_to_2]
lags_1_to_2 = range(1, len(te_values_1_to_2) + 1)
ax.errorbar(lags_1_to_2, te_values_1_to_2, yerr=te_std_1_to_2, label=f'TE: Scale {scales[0]} -> {scales[1]}', capsize=5)
key_2_to_1 = (sgn, scales[1], scales[0], round(shift, 2))
if key_2_to_1 in results:
te_values_2_to_1, te_std_2_to_1 = results[key_2_to_1]
lags_2_to_1 = range(1, len(te_values_2_to_1) + 1)
ax.errorbar(lags_2_to_1, te_values_2_to_1, yerr=te_std_2_to_1, label=f'TE: Scale {scales[1]} -> {scales[0]}', capsize=5)
ax.axvline(x=scales[0], color='k', linestyle=':')
ax.axvline(x=scales[1], color='k', linestyle=':')
scale_diff = abs(scales[0] - scales[1])
ax.axvline(x=scale_diff, color='k', linestyle=':')
ax.legend(fontsize=12)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.savefig('/users2/local/e22froge/codes/TE_Filter/LinearShiftIncrement/TE_fbm_signal_shifts.pdf', bbox_inches='tight')
plt.show()
# %%
# %%
kernel_type = 'increment'
impulse = np.zeros(100)
impulse[len(impulse) // 2] = 1
fig, axes = plt.subplots(6, 2, figsize=(16, 30))
fig.suptitle(f'Impulse Response of Increment Kernel for Different Shifts', fontsize=30)
for j, shift in enumerate(shifts):
row, col = divmod(j, 2) # Get the row and column for the subplot
ax = axes[row, col]
ax.set_title(f'Shift: {round(shift, 2)}', fontsize=18)
ax.set_xlabel('Time', fontsize=16)
ax.set_ylabel('Amplitude', fontsize=16)
for scale in scales:
kernel = increment_kernel_linear_shift(scale, shift) # Use shift instead of position
impulse_response = convolve(impulse, kernel, mode='same')
ax.plot(impulse_response, label=f'Scale={scale}')
ax.axhline(0, color='black', linestyle='--')
ax.legend(fontsize=12)
plt.tight_layout(rect=[0, 0, 1, 0.95])
plt.savefig('/users2/local/e22froge/codes/TE_Filter/LinearShiftIncrement/Impulse_responses_Increment_Shifts.pdf', bbox_inches='tight')
plt.show()
# %%
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment