fixed the code again and added all the processing code needed for each plot into the plot function itself. ignore plot_behavior for now. thats a task for later.

Author: j2figuer

src/plotting/plot_photometry.py

@@ -3,10 +3,12 @@
 from matplotlib.widgets import Slider
 from scipy.stats import linregress
 
-def plot_raw_ratiometric_and_565_photometry_signals(xvals, pulse_times_in_mins, raw_green, raw_red, raw_405, relative_raw_signal, sampling_rate): 
+def plot_raw_photometry_signals(visits, raw_green, raw_red, raw_405, relative_raw_signal, sampling_rate): 
     """
     Plots the raw 470, 405, 565 and ratiometric 470/405 fluorescence signals.
     """
+    xvals = np.arange(0,len(raw_green))/sampling_rate/60 
+    pulse_times_in_mins = [time / 60000 for time in visits]
 
     raw = plt.figure(figsize=(16, 10))
     plt.suptitle('Raw & ratiometric 470/405 fluorescence signals', fontsize=16)
@@ -42,6 +44,7 @@ def plot_405_470_correlation(raw_405, raw_green, slope_405x470, intercep_405x470
     """
     Plots the correlation between the 405 and 470 signals.
     """
+    slope_405x470, intercep_405x470, r_value_405x470 = linregress(x=raw_405, y=raw_green)
     plt.figure()
     plt.scatter(raw_405[::5], raw_green[::5],alpha=0.1, marker='.')
     x = np.array(plt.xlim())
@@ -62,6 +65,7 @@ def plot_405_565_correlation(raw_405, raw_red, slope_405x565, intercep_405x565,
     """
     Plots the correlation between the 405 and 565 signals.
     """
+    slope_405x565, intercep_405x565, r_value_405x565 = linregress(x=raw_405, y=raw_red)
     plt.figure()
     plt.scatter(raw_405[::5], raw_red[::5],alpha=0.1, marker='.')
     x = np.array(plt.xlim())
@@ -82,7 +86,7 @@ def plot_470_565_correlation(raw_green, raw_red, slope_470x565, intercep_470x565
     """
     Plots the correlation between the 470 and 565 signals.
     """
-
+    slope_470x565, intercep_470x565, r_value_470x565 = linregress(x=raw_green, y=raw_red)
     plt.figure()
     plt.scatter(raw_green[::5], raw_red[::5],alpha=0.1, marker='.')
     x = np.array(plt.xlim())
@@ -103,7 +107,7 @@ def plot_ratio_565_correlation(ratio_highpass, red_highpass, slope_filtered, int
     """
     Plots the correlation between the 470/405 ratio and the 565 signal.
     """
-    
+    slope_filtered, intercept_filtered, r_value_filtered = linregress(x=red_highpass, y=ratio_highpass)
     plt.figure(figsize=(13, 10))
     plt.scatter(red_highpass[::5], ratio_highpass[::5],alpha=0.1, marker='.')
     x = np.array(plt.xlim())
@@ -175,10 +179,62 @@ def update(val):
     plt.show()
     return None
 
-def plot_signals_aligned_port_entry(port_entry_indices, aligned_time, mean_ratio, sem_ratio, mean_green, sem_green, mean_405, sem_405, mean_red, sem_red, total_rewarded_trials, total_omitted_trials):
+def plot_signals_aligned_port_entry(sampling_rate, time_seconds, visits, port_entry_indices, aligned_time,ratio_highpass, highpass_405, green_highpass, red_highpass, mean_ratio, sem_ratio, mean_green, sem_green, mean_405, sem_405, mean_red, sem_red, total_rewarded_trials, total_omitted_trials):
     """
     """
-
+    time_window_in_sec = 10
+    
+    samples_window = time_window_in_sec * sampling_rate
+    total_rewarded_trials = (rwd == 1).sum()
+    print("Total rewarded trials: "+str(total_rewarded_trials))
+
+    total_omitted_trials = (rwd == 0).sum()
+    print("Total omitted trials: "+str(total_omitted_trials))
+
+    # Initialize lists for aligned traces
+    aligned_ratio = []
+    aligned_green = []
+    aligned_405 = []
+    aligned_red = []
+
+#################################################################################
+    # QUESTION FOR STEPH:
+    # I need this code to run some of the plotting from this point forward.
+    # Is there a way to get this information from convert_behvaior?
+    # For now, I'll keep everything the way it is using "visinds" and "port_entry_indeces"
+    visinds = sampledata.loc[sampledata.port.notnull()].index.values
+    port_entry_indices = visinds 
+#################################################################################
+
+    # Align windows
+    for idx in visits:
+        if idx - samples_window >= 0 and idx + samples_window < len(time_seconds):  # Ensure indices are in bounds
+            aligned_ratio.append(ratio_highpass[idx - samples_window:idx + samples_window])
+            aligned_green.append(green_highpass[idx - samples_window:idx + samples_window])
+            aligned_405.append(highpass_405[idx - samples_window:idx + samples_window])
+            aligned_red.append(red_highpass[idx - samples_window:idx + samples_window])
+
+    # Convert to numpy arrays for easier averaging
+    aligned_ratio = np.array(aligned_ratio)
+    aligned_green = np.array(aligned_green)
+    aligned_405 = np.array(aligned_405)
+    aligned_red = np.array(aligned_red)
+
+    # Calculate means and SEMs
+    mean_ratio = np.mean(aligned_ratio, axis=0)
+    sem_ratio = np.std(aligned_ratio, axis=0) / np.sqrt(len(aligned_ratio))
+
+    mean_green = np.mean(aligned_green, axis=0)
+    sem_green = np.std(aligned_green, axis=0) / np.sqrt(len(aligned_green))
+
+    mean_405 = np.mean(aligned_405, axis=0)
+    sem_405 = np.std(aligned_405, axis=0) / np.sqrt(len(aligned_405))
+
+    mean_red = np.mean(aligned_red, axis=0)
+    sem_red = np.std(aligned_red, axis=0) / np.sqrt(len(aligned_red))
+
+    # Generate time axis for aligned data
+    aligned_time = np.linspace(-time_window_in_sec, time_window_in_sec, 2 * samples_window)
     # Plot averaged traces
     plt.figure(figsize=(16, 10))
     plt.title("Average Traces Aligned to Port Entry ("+str(len(port_entry_indices))+" total trials)")
@@ -208,10 +264,10 @@ def plot_signals_aligned_port_entry(port_entry_indices, aligned_time, mean_ratio
 
     return None
 
-def plot_normalized_signals(pulse_times_in_mins, green_zscored, zscored_405, red_zscored, ratio_zscored, xvals):
+def plot_normalized_signals(sampling_rate, pulse_times_in_mins, green_zscored, zscored_405, red_zscored, ratio_zscored):
     """
     """
-
+    xvals = np.arange(0,len(green_zscored))/sampling_rate/60 
     zscrd = plt.figure(figsize=(16, 10))
     plt.suptitle('Z-scored signals calculated after preprocessing raw singals \n by applying a high-pass filter at 0.001 Hz and a low-pass filter at 10 Hz')
 
@@ -245,7 +301,42 @@ def plot_normalized_signals(pulse_times_in_mins, green_zscored, zscored_405, red
 def plot_ratio_and_565_signals_aligned_port_entry(rat, date, time_window, time_window_xvals, ratio_traces_mean, ratio_traces_sem, red_traces_mean, red_traces_sem, visinds): 
     """
     """
-
+    ratio_traces = []
+    red_traces = [] 
+    time_window_in_sec = 10
+
+    time_window_xvals = np.arange(-time_window_in_sec*86,time_window_in_sec*86+1)/86
+
+    # Loop through each index in visinds
+    for i in visinds:
+        # Extract a time window of data around the current index
+        start_idx = i - 86 * time_window_in_sec
+        end_idx = i + 86 * time_window_in_sec
+        ratio_trace = sampledata.loc[start_idx:end_idx, "ratio_z_scored"].values
+        red_trace = sampledata.loc[start_idx:end_idx, "red_z_scored"].values
+
+        # Only use traces of the correct length (matching xvals)
+        if len(ratio_trace) == len(time_window_xvals):
+            ratio_traces.append(ratio_trace)
+
+        if len(red_trace) == len(time_window_xvals):
+            red_traces.append(red_trace)
+
+    # Calculate mean and SEM 
+    if ratio_traces: 
+        ratio_traces_mean = np.mean(ratio_traces, axis=0)
+        ratio_traces_sem = np.std(ratio_traces, axis=0) / np.sqrt(len(ratio_traces))
+    else:  # If no traces, create empty arrays
+        ratio_traces_rwd_mean = np.zeros_like(time_window_xvals)
+        ratio_traces_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM 
+    if red_traces: 
+        red_traces_mean = np.mean(red_traces, axis=0)
+        red_traces_sem = np.std(red_traces, axis=0) / np.sqrt(len(red_traces))
+    else:  # If no traces, create empty arrays
+        red_traces_mean = np.zeros_like(time_window_xvals)
+        red_traces_sem = np.zeros_like(time_window_xvals)
     fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, sharex=True, sharey=True, figsize=(19, 10))
 
     plt.suptitle(f'Z-scored ratiometric 470/405 & 565 fluorescence signals {rat} {date} {time_window}sec', fontsize=16)
@@ -269,9 +360,73 @@ def plot_ratio_and_565_signals_aligned_port_entry(rat, date, time_window, time_w
 
     return None
 
-def plot_ratio_and_565_signals_aligned_port_entry_separated_by_rewarded_or_omitted(rat, date, time_window_in_sec, xvals, ratio_rwd_mean, ratio_rwd_sem, ratio_om_mean, ratio_om_sem, red_rwd_mean, red_rwd_sem, red_om_mean, red_om_sem, total_rewarded_trials, total_omitted_trials):
+def plot_ratio_and_565_signals_aligned_port_entry_separated_by_rewarded_or_omitted(rat, date, visinds, time_window_in_sec, xvals, ratio_rwd_mean, ratio_rwd_sem, ratio_om_mean, ratio_om_sem, red_rwd_mean, red_rwd_sem, red_om_mean, red_om_sem, total_rewarded_trials, total_omitted_trials):
     """
     """
+    ratio_rwd_traces = []
+    ratio_om_traces = []
+    red_rwd_traces = []
+    red_om_traces = []
+
+    time_window_in_sec = 10
+
+    time_window_xvals = np.arange(-time_window_in_sec*86,time_window_in_sec*86+1)/86
+
+    # Loop through each index in visinds
+    for i in visinds:
+        # Extract a time window of data around the current index
+        start_idx = i - 86 * time_window_in_sec
+        end_idx = i + 86 * time_window_in_sec
+        ratio_trace = sampledata.loc[start_idx:end_idx, "ratio_z_scored"].values
+        red_trace = sampledata.loc[start_idx:end_idx, "red_z_scored"].values
+
+        # Only use traces of the correct length (matching xvals)
+        if len(ratio_trace) == len(time_window_xvals):
+            # Check if the trial is rewarded or omitted
+            if sampledata.loc[i, 'rwd'] == 1:
+                ratio_rwd_traces.append(ratio_trace)  # Add to rewarded traces
+            else:
+                ratio_om_traces.append(ratio_trace)  # Add to omitted traces
+
+            # Only use traces of the correct length (matching xvals)
+        if len(red_trace) == len(time_window_xvals):
+            # Check if the trial is rewarded or omitted
+            if sampledata.loc[i, 'rwd'] == 1:
+                red_rwd_traces.append(red_trace)  # Add to rewarded traces
+            else:
+                red_om_traces.append(red_trace)  # Add to omitted traces
+
+    # Calculate mean and SEM for rewarded traces
+    if ratio_rwd_traces:  # If there are any rewarded traces
+        ratio_rwd_mean = np.mean(ratio_rwd_traces, axis=0)
+        ratio_rwd_sem = np.std(ratio_rwd_traces, axis=0) / np.sqrt(len(ratio_rwd_traces))
+    else:  # If no rewarded traces, create empty arrays
+        ratio_rwd_mean = np.zeros_like(time_window_xvals)
+        ratio_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for rewarded traces
+    if red_rwd_traces:  # If there are any rewarded traces
+        red_rwd_mean = np.mean(red_rwd_traces, axis=0)
+        red_rwd_sem = np.std(red_rwd_traces, axis=0) / np.sqrt(len(red_rwd_traces))
+    else:  # If no rewarded traces, create empty arrays
+        red_rwd_mean = np.zeros_like(time_window_xvals)
+        red_rwd_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for omitted traces
+    if ratio_om_traces:  # If there are any omitted traces
+        ratio_om_mean = np.mean(ratio_om_traces, axis=0)
+        ratio_om_sem = np.std(ratio_om_traces, axis=0) / np.sqrt(len(ratio_om_traces))
+    else:  # If no omitted traces, create empty arrays
+        ratio_om_mean = np.zeros_like(time_window_xvals)
+        ratio_om_sem = np.zeros_like(time_window_xvals)
+
+    # Calculate mean and SEM for omitted traces
+    if red_om_traces:  # If there are any omitted traces
+        red_om_mean = np.mean(red_om_traces, axis=0)
+        red_om_sem = np.std(red_om_traces, axis=0) / np.sqrt(len(red_om_traces))
+    else:  # If no omitted traces, create empty arrays
+        red_om_mean = np.zeros_like(time_window_xvals)
+        red_om_sem = np.zeros_like(time_window_xvals)
 
     fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, sharex=True, sharey=True, figsize=(19, 10))