26 KiB
26 KiB
In [ ]:
import os
import numpy as np
import plotly.graph_objects as go
import plotly.express as px
from scipy.io import wavfile
from data_processing import load_all_wavs
import random
# Constants
SAMPLE_RATE = 19531
LEAD_LIMIT = 256
TIME_LIMIT_SEC = 0.1
DATA_DIR = "data"
# Utility functions for visualizations
def load_and_trim_data(data_dir, sample_rate, lead_limit, time_limit_sec):
"""Load and trim the data to specified lead and time limits."""
all_data = load_all_wavs(data_dir)
trimmed_data = []
max_samples = int(sample_rate * time_limit_sec)
for lead in all_data[:lead_limit]:
trimmed_data.append(lead[:max_samples])
return np.array(trimmed_data)
def plot_individual_leads(data, sample_rate):
"""Plot individual leads using Plotly for interactive zoom."""
fig = go.Figure()
for i, lead in enumerate(data):
fig.add_trace(go.Scatter(
x=np.linspace(0, len(lead) / sample_rate, num=len(lead)),
y=lead,
mode='lines',
name=f'Lead {i+1}'
))
fig.update_layout(
title='Individual Leads',
xaxis_title='Time [s]',
yaxis_title='Amplitude'
)
return fig
def plot_lead_correlations(data):
"""Plot correlation between leads using Plotly."""
correlations = np.corrcoef(data)
fig = px.imshow(correlations,
labels=dict(color="Correlation"),
x=[f'Lead {i+1}' for i in range(len(data))],
y=[f'Lead {i+1}' for i in range(len(data))],
color_continuous_scale='RdBu_r',
zmin=-1, zmax=1)
fig.update_layout(
title='Correlation Between Leads'
)
return fig
def plot_highly_correlated_leads(data, threshold=0.75):
"""Plot leads that have high correlation with each other using Plotly."""
correlations = np.corrcoef(data)
high_corr_pairs = np.argwhere(correlations > threshold)
grouped_pairs = {}
for (i, j) in high_corr_pairs:
if i >= j:
continue
if i not in grouped_pairs:
grouped_pairs[i] = []
grouped_pairs[i].append(j)
figs = []
for i, group in grouped_pairs.items():
fig = go.Figure()
fig.add_trace(go.Scatter(
x=np.arange(len(data[i])),
y=data[i],
mode='lines',
name=f'Lead {i+1}'
))
for j in group:
fig.add_trace(go.Scatter(
x=np.arange(len(data[j])),
y=data[j],
mode='lines',
name=f'Lead {j+1}',
line=dict(dash='dash')
))
fig.update_layout(
title=f'Highly Correlated Leads Group {i+1}',
xaxis_title='Sample',
yaxis_title='Amplitude'
)
figs.append(fig)
return figs
def plot_top_correlated_pairs(data, top_n=3):
"""Plot the top N most highly correlated pairs of leads and their correlations over time."""
correlations = np.corrcoef(data)
np.fill_diagonal(correlations, 0) # Ignore self-correlations
top_pairs = np.unravel_index(np.argsort(correlations.ravel())[-top_n:], correlations.shape)
top_pairs = list(zip(top_pairs[0], top_pairs[1]))
figs = []
for i, (lead1, lead2) in enumerate(top_pairs):
fig = go.Figure()
fig.add_trace(go.Scatter(
x=np.arange(len(data[lead1])),
y=data[lead1],
mode='lines',
name=f'Lead {lead1+1}'
))
fig.add_trace(go.Scatter(
x=np.arange(len(data[lead2])),
y=data[lead2],
mode='lines',
name=f'Lead {lead2+1}',
line=dict(dash='dash')
))
fig.update_layout(
title=f'Top Correlated Pair {i+1}: Lead {lead1+1} and Lead {lead2+1}',
xaxis_title='Sample',
yaxis_title='Amplitude'
)
figs.append(fig)
# Create correlation matrix
lead1_data = data[lead1]
lead2_data = data[lead2]
arr = np.stack((lead1_data, lead2_data), axis=1)
correlation_matrix = np.cov(arr) / np.max(lead1_data) / np.max(lead2_data)
# Replace NaN values with 0
#correlation_matrix = np.nan_to_num(correlation_matrix)
fig = px.imshow(correlation_matrix,
color_continuous_scale='RdBu_r')
fig.update_layout(
title=f'Correlation Over Time for Lead {lead1+1} and Lead {lead2+1}',
)
figs.append(fig)
return figs
# Load and trim the data
trimmed_data = load_and_trim_data(DATA_DIR, SAMPLE_RATE, LEAD_LIMIT, TIME_LIMIT_SEC)
# Plot individual leads
fig_individual_leads = plot_individual_leads(trimmed_data, SAMPLE_RATE)
fig_individual_leads.show()
# Plot lead correlations
fig_lead_correlations = plot_lead_correlations(trimmed_data)
fig_lead_correlations.show()
# Plot highly correlated leads
figs_highly_correlated_leads = plot_highly_correlated_leads(trimmed_data)
for fig in figs_highly_correlated_leads:
fig.show()
# Plot top correlated pairs and their correlations over time
#figs_top_correlated_pairs = plot_top_correlated_pairs(trimmed_data)
#for fig in figs_top_correlated_pairs:
# fig.show()
In [ ]:
import numpy as np
import plotly.graph_objects as go
from scipy.io import wavfile
from data_processing import load_all_wavs
import random
# Constants
SAMPLE_RATE = 19531
LEAD_LIMIT = 128
TIME_LIMIT_SEC = 0.1
DATA_DIR = "data"
RECONSTRUCTION_START_FRACTION = 0.5 # Fraction of time at which to start reconstruction
# Utility functions for visualizations
def load_and_trim_data(data_dir, sample_rate, lead_limit, time_limit_sec):
"""Load and trim the data to specified lead and time limits."""
all_data = load_all_wavs(data_dir)
trimmed_data = []
max_samples = int(sample_rate * time_limit_sec)
for lead in all_data[:lead_limit]:
trimmed_data.append(lead[:max_samples])
return np.array(trimmed_data)
def plot_random_lead_with_reconstruction(data, reconstruction_start_fraction, top_n=8):
"""Plot a random lead with reconstruction based on correlations with other leads."""
num_leads = data.shape[0]
random_lead_index = random.randint(0, num_leads - 1)
# Get the random lead data
random_lead_data = data[random_lead_index]
# Calculate correlations
correlations = np.corrcoef(data)
np.fill_diagonal(correlations, 0) # Ignore self-correlations
# Find the top N most correlated leads
most_correlated_indices = np.argsort(correlations[random_lead_index])[-top_n:]
# Plot the original random lead
fig = go.Figure()
time = np.linspace(0, len(random_lead_data) / SAMPLE_RATE, num=len(random_lead_data))
fig.add_trace(go.Scatter(
x=time,
y=random_lead_data,
mode='lines',
name=f'Original Lead {random_lead_index+1}'
))
# Perform reconstruction
reconstruction_start_index = int(len(random_lead_data) * reconstruction_start_fraction)
reconstructed_data = np.zeros_like(random_lead_data)
reconstructed_data[reconstruction_start_index] = random_lead_data[reconstruction_start_index-2]*0.25 + random_lead_data[reconstruction_start_index-1]*0.25 + random_lead_data[reconstruction_start_index]*0.5
for i in range(reconstruction_start_index + 1, len(random_lead_data)):
for correlated_index in most_correlated_indices:
reconstructed_data[i] += correlations[random_lead_index, correlated_index] * \
(data[correlated_index, i] - data[correlated_index, i-1])
reconstructed_data[i] += reconstructed_data[i-1]
# Plot the reconstructed lead
fig.add_trace(go.Scatter(
x=time[reconstruction_start_index:],
y=reconstructed_data[reconstruction_start_index:],
mode='lines',
name=f'Reconstructed Lead {random_lead_index+1}',
line=dict(dash='dash')
))
fig.update_layout(
title=f'Lead {random_lead_index+1} with Reconstructed Data',
xaxis_title='Time [s]',
yaxis_title='Amplitude',
legend_title_text='Data'
)
fig.show()
# Load and trim the data
trimmed_data = load_and_trim_data(DATA_DIR, SAMPLE_RATE, LEAD_LIMIT, TIME_LIMIT_SEC)
# Plot random lead with reconstruction
for i in range(8):
plot_random_lead_with_reconstruction(trimmed_data, RECONSTRUCTION_START_FRACTION)
In [ ]:
import numpy as np
import plotly.graph_objects as go
from scipy.io import wavfile
from data_processing import load_all_wavs
import random
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, RationalQuadratic
from sklearn.preprocessing import StandardScaler
# Constants
SAMPLE_RATE = 19531
LEAD_LIMIT = 128
TIME_LIMIT_SEC = 0.1
DATA_DIR = "data"
RECONSTRUCTION_START_FRACTION = 0.5 # Fraction of time at which to start reconstruction
KERNEL_TYPE = 'RBF' # Configurable kernel: 'RBF', 'Matern', 'RationalQuadratic'
ALPHA = 1e-6 # Small noise term to ensure positive definite kernel
NUM_GP_SAMPLES = 1 # Number of GP samples to draw
TOP_N = 3 # Number of most correlated leads to consider
GP_TRAINING_POINTS = 64 # Number of recent points to use for GP fitting
# Utility functions for visualizations
def load_and_trim_data(data_dir, sample_rate, lead_limit, time_limit_sec):
"""Load and trim the data to specified lead and time limits."""
all_data = load_all_wavs(data_dir)
trimmed_data = []
max_samples = int(sample_rate * time_limit_sec)
for lead in all_data[:lead_limit]:
trimmed_data.append(lead[:max_samples])
return np.array(trimmed_data)
def plot_combined_reconstruction(data, reconstruction_start_fraction, kernel_type, num_gp_samples, top_n, gp_training_points):
"""Plot a random lead with combined GP extrapolation and correlated lead contributions."""
num_leads = data.shape[0]
random_lead_index = random.randint(0, num_leads - 1)
# Get the random lead data
random_lead_data = data[random_lead_index]
# Normalize the data
scaler = StandardScaler()
random_lead_data_normalized = scaler.fit_transform(random_lead_data.reshape(-1, 1)).flatten()
# Calculate correlations
correlations = np.corrcoef(data)
np.fill_diagonal(correlations, 0) # Ignore self-correlations
# Find the top N most correlated leads
most_correlated_indices = np.argsort(correlations[random_lead_index])[-top_n:]
# Plot the original random lead
fig = go.Figure()
time = np.linspace(0, len(random_lead_data) / SAMPLE_RATE, num=len(random_lead_data))
fig.add_trace(go.Scatter(
x=time,
y=random_lead_data,
mode='lines',
name=f'Original Lead {random_lead_index+1}'
))
# Perform GP extrapolation
reconstruction_start_index = int(len(random_lead_data) * reconstruction_start_fraction)
# Fit Gaussian Process
kernel = None
if kernel_type == 'RBF':
kernel = RBF(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
elif kernel_type == 'Matern':
kernel = Matern(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
elif kernel_type == 'RationalQuadratic':
kernel = RationalQuadratic(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
gp = GaussianProcessRegressor(kernel=kernel, alpha=ALPHA)
train_time = time[reconstruction_start_index - gp_training_points:reconstruction_start_index].reshape(-1, 1)
train_data = random_lead_data_normalized[reconstruction_start_index - gp_training_points:reconstruction_start_index]
gp.fit(train_time, train_data)
extrapolated_time = time[reconstruction_start_index:].reshape(-1, 1)
extrapolated_mean_normalized, extrapolated_std_normalized = gp.predict(extrapolated_time, return_std=True)
# Inverse transform the normalized predictions
extrapolated_mean = scaler.inverse_transform(extrapolated_mean_normalized.reshape(-1, 1)).flatten()
extrapolated_std = extrapolated_std_normalized * scaler.scale_[0]
# Initialize the combined reconstruction with GP mean
combined_reconstruction = np.copy(extrapolated_mean)*0 + random_lead_data[reconstruction_start_index-1] *0.33 + random_lead_data[reconstruction_start_index]*0.66
# Calculate variances (inverse of weights)
gp_variance = extrapolated_std ** 2
# Perform reconstruction
reconstruction_start_index = int(len(random_lead_data) * reconstruction_start_fraction)
reconstructed_data = np.zeros_like(random_lead_data)
reconstructed_data[reconstruction_start_index] = random_lead_data[reconstruction_start_index-2]*0.25 + random_lead_data[reconstruction_start_index-1]*0.25 + random_lead_data[reconstruction_start_index]*0.5
for i in range(reconstruction_start_index + 1, len(random_lead_data)):
for correlated_index in most_correlated_indices:
reconstructed_data[i] += correlations[random_lead_index, correlated_index] * \
(data[correlated_index, i] - data[correlated_index, i-1])
reconstructed_data[i] += reconstructed_data[i-1]
#combined_reconstruction = reconstructed_data
combined_reconstruction = (extrapolated_mean / gp_variance + reconstructed_data[reconstruction_start_index:]) / (1 / gp_variance + 1)
# Plot the GP extrapolated mean
#fig.add_trace(go.Scatter(
# x=time[reconstruction_start_index:],
# y=extrapolated_mean,
# mode='lines',
# name=f'GP Extrapolated Lead {random_lead_index+1}',
# line=dict(dash='dot')
#))
# Plot the combined reconstruction
fig.add_trace(go.Scatter(
x=time[reconstruction_start_index:],
y=combined_reconstruction,
mode='lines',
name=f'Combined Reconstruction Lead {random_lead_index+1}',
line=dict(dash='dash')
))
fig.add_trace(go.Scatter(
x=time[reconstruction_start_index:],
y=reconstructed_data[reconstruction_start_index:],
mode='lines',
name=f'Peer Reconstruction Lead {random_lead_index+1}',
line=dict(dash='dot')
))
fig.update_layout(
title=f'Lead {random_lead_index+1} with Combined GP and Correlated Lead Contributions',
xaxis_title='Time [s]',
yaxis_title='Amplitude',
legend_title_text='Data'
)
fig.show()
# Load and trim the data
trimmed_data = load_and_trim_data(DATA_DIR, SAMPLE_RATE, LEAD_LIMIT, TIME_LIMIT_SEC)
# Plot combined reconstruction
plot_combined_reconstruction(trimmed_data, RECONSTRUCTION_START_FRACTION, KERNEL_TYPE, NUM_GP_SAMPLES, TOP_N, GP_TRAINING_POINTS)
In [ ]:
import numpy as np
import plotly.graph_objects as go
from scipy.io import wavfile
from data_processing import load_all_wavs
import random
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, Matern, RationalQuadratic
from sklearn.preprocessing import StandardScaler
# Constants
SAMPLE_RATE = 19531
LEAD_LIMIT = 128
TIME_LIMIT_SEC = 0.1
DATA_DIR = "data"
RECONSTRUCTION_START_FRACTION = 0.5 # Fraction of time at which to start reconstruction
KERNEL_TYPE = 'RBF' # Configurable kernel: 'RBF', 'Matern', 'RationalQuadratic'
ALPHA = 1e-6 # Small noise term to ensure positive definite kernel
NUM_GP_SAMPLES = 1 # Number of GP samples to draw
TOP_N = 3 # Number of most correlated leads to consider
GP_TRAINING_POINTS = 64 # Number of recent points to use for GP fitting
WINDOW_SIZE = 10 # Window size for autoregressive GP
def load_and_trim_data(data_dir, sample_rate, lead_limit, time_limit_sec):
"""Load and trim the data to specified lead and time limits."""
all_data = load_all_wavs(data_dir)
trimmed_data = []
max_samples = int(sample_rate * time_limit_sec)
for lead in all_data[:lead_limit]:
trimmed_data.append(lead[:max_samples])
return np.array(trimmed_data)
def plot_combined_reconstruction(data, reconstruction_start_fraction, kernel_type, num_gp_samples, top_n, gp_training_points, window_size):
"""Plot a random lead with combined GP extrapolation and correlated lead contributions."""
num_leads = data.shape[0]
random_lead_index = random.randint(0, num_leads - 1)
# Get the random lead data
random_lead_data = data[random_lead_index]
# Normalize the data
scaler = StandardScaler()
random_lead_data_normalized = scaler.fit_transform(random_lead_data.reshape(-1, 1)).flatten()
# Calculate correlations
correlations = np.corrcoef(data)
np.fill_diagonal(correlations, 0) # Ignore self-correlations
# Find the top N most correlated leads
most_correlated_indices = np.argsort(correlations[random_lead_index])[-top_n:]
# Plot the original random lead
fig = go.Figure()
time = np.linspace(0, len(random_lead_data) / SAMPLE_RATE, num=len(random_lead_data))
fig.add_trace(go.Scatter(
x=time,
y=random_lead_data,
mode='lines',
name=f'Original Lead {random_lead_index+1}'
))
# Determine reconstruction start index
reconstruction_start_index = int(len(random_lead_data) * reconstruction_start_fraction)
# Fit Gaussian Process on the initial training data
kernel = None
if kernel_type == 'RBF':
kernel = RBF(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
elif kernel_type == 'Matern':
kernel = Matern(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
elif kernel_type == 'RationalQuadratic':
kernel = RationalQuadratic(length_scale=1.0, length_scale_bounds=(1e-5, 1e5))
gp = GaussianProcessRegressor(kernel=kernel, alpha=ALPHA)
extrapolated_mean = np.zeros(len(random_lead_data) - reconstruction_start_index)
extrapolated_std = np.zeros(len(random_lead_data) - reconstruction_start_index)
combined_reconstruction = np.zeros(len(random_lead_data)) + random_lead_data[reconstruction_start_index]
reconstructed_data = np.zeros(len(random_lead_data))
for i in range(reconstruction_start_index, len(random_lead_data)):
# Use past samples to ensure the window is full
start_index = max(0, i - window_size)
# Prepare training data for GP
train_time = time[start_index:i].reshape(-1, 1)
train_data = combined_reconstruction[start_index:i]
train_data_normalized = scaler.transform(train_data.reshape(-1, 1)).flatten()
# Fit the GP on the current window
gp.fit(train_time, train_data_normalized)
# Predict the next value
next_time = np.array([[time[i]]])
mean, std = gp.predict(next_time, return_std=True)
# Inverse transform the normalized predictions
mean = scaler.inverse_transform(mean.reshape(-1, 1)).flatten()[0]
std = std[0] * scaler.scale_[0]
# Calculate the variance
gp_variance = std ** 2
# Calculate contributions from the most correlated leads
correlated_contribution = combined_reconstruction[i-1]
for correlated_index in most_correlated_indices:
correlated_contribution += correlations[random_lead_index, correlated_index] * (data[correlated_index, i] - data[correlated_index, i-1])
# Combine GP prediction with peer-based reconstruction
combined_value = (mean / gp_variance + correlated_contribution) / (1 / gp_variance + 1)
combined_reconstruction[i] = combined_value
# Update the reconstructed data
reconstructed_data[i] = correlated_contribution
# Store GP predictions
if i - reconstruction_start_index < len(extrapolated_mean):
extrapolated_mean[i - reconstruction_start_index] = mean
extrapolated_std[i - reconstruction_start_index] = std
# Plot the combined reconstruction
fig.add_trace(go.Scatter(
x=time[reconstruction_start_index:],
y=combined_reconstruction[reconstruction_start_index:],
mode='lines',
name=f'Combined Reconstruction Lead {random_lead_index+1}',
line=dict(dash='dash')
))
fig.add_trace(go.Scatter(
x=time[reconstruction_start_index:],
y=reconstructed_data[reconstruction_start_index:],
mode='lines',
name=f'Peer Reconstruction Lead {random_lead_index+1}',
line=dict(dash='dot')
))
fig.update_layout(
title=f'Lead {random_lead_index+1} with Combined GP and Correlated Lead Contributions',
xaxis_title='Time [s]',
yaxis_title='Amplitude',
legend_title_text='Data'
)
fig.show()
# Load and trim the data
trimmed_data = load_and_trim_data(DATA_DIR, SAMPLE_RATE, LEAD_LIMIT, TIME_LIMIT_SEC)
# Plot combined reconstruction
for _ in range(8):
plot_combined_reconstruction(trimmed_data, RECONSTRUCTION_START_FRACTION, KERNEL_TYPE, NUM_GP_SAMPLES, TOP_N, GP_TRAINING_POINTS, WINDOW_SIZE)
In [ ]: