createVisualizations.py

import numpy as np
#import matplotlib
#matplotlib.use("Agg")
import matplotlib.pyplot as plt
from datetime import datetime

import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits import mplot3d
from sklearn.decomposition import FastICA
import pdb
from scipy.ndimage import gaussian_filter1d
import cv2
import pandas as pd
import scipy
#from pylab import *
# import tifffile as tiff
import matplotlib as mpl
from matplotlib.colors import Normalize
import matplotlib.cm as cm
from matplotlib import rcParams
import matplotlib.gridspec as gridspec
import itertools
import statistics
import matplotlib.colors as colors
from matplotlib.ticker import MultipleLocator
import matplotlib.gridspec as gridspec
import matplotlib.cm as cm
from collections import OrderedDict
#import sima
#import sima.motion
#import sima.segment
from scipy.signal import find_peaks
from scipy.stats.stats import pearsonr
from scipy.interpolate import interp1d
#from mtspec import mt_coherence
from scipy import stats
import matplotlib.ticker as ticker
from scipy.signal import butter, filtfilt, convolve
import pywt
from cycler import cycler
from sklearn.preprocessing import MinMaxScaler
from sklearn.preprocessing import scale
from sklearn.preprocessing import MaxAbsScaler
import tools.dataAnalysis as dataAnalysis
from tools.pyqtgraph.configfile import *
from scipy.ndimage import uniform_filter1d
import scipy.ndimage as ndimage
from scipy.signal import resample
from scipy import signal
import seaborn as sns
from itertools import cycle
params= OrderedDict([
    ('videoParameter',{
        'dpi': 500,
        }),
    ('projectionInTimeParameters', {
        'horizontalCuts' : [5,8,9,10,12] ,
        'verticalCut' : 50. ,
        'stimStart': 5. ,
        'stimLength' : 0.2 ,
        'fitStart' : 5. ,
        'baseLinePeriod' : 5. ,
        'threeDAspectRatio' : 3,
        'stimulationBarLocation' : -0.1,
        }),
    ('caEphysParameters', {
        'leaveOut' : 0.1,
        }),
    ])

class createVisualizations:

    ##########################################################################################
    def __init__(self, figureDir, mouse):

        self.mouse = mouse
        self.figureDirectory = figureDir
        if not os.path.isdir(self.figureDirectory):
            os.system('mkdir %s' % self.figureDirectory)

        configFile = self.figureDirectory + '%s.config' % mouse
        if os.path.isfile(configFile):
            self.config = readConfigFile(configFile)
        else:
            self.config = params
            writeConfigFile(self.config, configFile)

        # self.pawID = ['FL', 'FR', 'HL', 'HR']

    ##########################################################################################
    def determineFileName(self, reco, what=None, date=None):
        if (what is None) and (date is None):
            ff = self.figureDirectory + '%s' % (reco)
        elif date is None:
            ff = self.figureDirectory + '%s_%s' % (reco, what)
        else:
            ff = self.figureDirectory + '%s_%s_%s' % (date, reco, what)
        return ff
    ##########################################################################################
    def layoutOfPanel(self, ax,xLabel=None,yLabel=None,Leg=None,xyInvisible=[False,False]):


        ax.spines['top'].set_visible(False)
        ax.spines['right'].set_visible(False)
        ax.grid(False)
        #
        if xyInvisible[0]:
            ax.spines['bottom'].set_visible(False)
            ax.xaxis.set_visible(False)
        else:
            ax.spines['bottom'].set_position(('outward', 10))
            ax.xaxis.set_ticks_position('bottom')
        #
        if xyInvisible[1]:
            ax.spines['left'].set_visible(False)
            ax.yaxis.set_visible(False)
        else:
            ax.spines['left'].set_position(('outward', 10))
            ax.yaxis.set_ticks_position('left')


        if xLabel != None :
            ax.set_xlabel(xLabel)

        if yLabel != None :
            ax.set_ylabel(yLabel)

        if Leg != None :
            ax.legend(loc=Leg[0], frameon=False)
            if len(Leg)>1 :
                legend = ax.get_legend()  # plt.gca().get_legend()
                ltext = legend.get_texts()
                plt.setp(ltext, fontsize=Leg[1])


    def visualizeSignals(self, rawSignals, denoisedSignals):
        """
        Visualize the raw and denoised fluorescence signals for each ROI.

        Parameters:
        - rawSignals: 2D numpy array of the raw signals.
        - denoisedSignals: 2D numpy array of the denoised signals.
        """
        # Set the color cycle to tab10
        plt.rc('axes', prop_cycle=(cycler('color', plt.cm.tab10.colors)))

        numRois = rawSignals.shape[0]
        fig, axes = plt.subplots(numRois, 1, figsize=(15, numRois * 2))

        if numRois == 1:
            axes = [axes]  # Make sure axes is iterable

        for i, ax in enumerate(axes):
            # Plot raw signal
            ax.plot(rawSignals[i, :], label=f'Raw Signal ROI {i + 1}')
            # Plot denoised signal
            ax.plot(denoisedSignals[i, :], label=f'Denoised Signal ROI {i + 1}', linestyle='--')
            ax.legend()
            ax.set_title(f'ROI {i + 1}')
            ax.set_xlabel('Time')
            ax.set_ylabel('Fluorescence Intensity')

        plt.tight_layout()
        plt.show()

    def plotReconstructedComponents(self, rawSignals, reconstructions, n_components=5):
        """
        Plot the signal reconstructed from each principal component on separate lines for each ROI.

        Parameters:
        - rawSignals: 2D numpy array of the raw signals.
        - reconstructions: A list of 2D numpy arrays, each containing the signal
                           reconstructed from each principal component.
        - n_components: The number of principal components to plot.
        """
        numRois = rawSignals.shape[0]
        colors = plt.cm.viridis(np.linspace(0, 1, n_components))

        # Create subplots for each ROI
        fig, axes = plt.subplots(numRois, 1, figsize=(15, numRois * 3))

        if numRois == 1:
            axes = [axes]  # Make sure axes is iterable

        for roi_idx, ax in enumerate(axes):
            ax.plot(rawSignals[roi_idx, :], label='Raw Signal', color='black', linewidth=2)

            for comp_idx in range(n_components):
                ax.plot(reconstructions[comp_idx][roi_idx, :], label=f'PC {comp_idx + 1}', color=colors[comp_idx])

            ax.set_title(f'ROI {roi_idx + 1}')
            ax.set_xlabel('Time')
            ax.set_ylabel('Fluorescence Intensity')
            ax.legend()

        plt.tight_layout()
        plt.show()


    def plotLeastNoisyComponent(self, rawSignals, reconstructions, pca, n_components=5):
        """
        Plot the least noisy principal component for each ROI.

        Parameters:
        - rawSignals: 2D numpy array of the raw signals.
        - reconstructions: A list of 2D numpy arrays, each containing the signal
                           reconstructed from each principal component.
        - pca: The PCA object after fitting.
        - n_components: The number of principal components to plot.
        """
        numRois = rawSignals.shape[0]

        # Find the index of the least noisy component (the one with the smallest eigenvalue)
        least_noisy_index = np.argmax(pca.explained_variance_)

        # Create subplots for each ROI
        fig, axes = plt.subplots(numRois, 1, figsize=(15, numRois * 3))

        if numRois == 1:
            axes = [axes]  # Make sure axes is iterable

        for roi_idx, ax in enumerate(axes):
            ax.plot(rawSignals[roi_idx, :], label='Raw Signal', color='black', linewidth=2)
            ax.plot(reconstructions[1][roi_idx, :], label=f'Least Noisy PC', color='red')

            ax.set_title(f'ROI {roi_idx + 1}')
            ax.set_xlabel('Time')
            ax.set_ylabel('Fluorescence Intensity')
            ax.legend()

        plt.tight_layout()
        plt.show()

    def createEMGLaserStimPawMovementFigure(self, date, rec, EMG, swingStanceD, pawPos):
        # figure #################################
        fig_width = 20  # width in inches
        fig_height = 16  # height in inches
        fig_size = [fig_width, fig_height]
        params = {'axes.labelsize': 12, 'axes.titlesize': 12, 'font.size': 11, 'xtick.labelsize': 11, 'ytick.labelsize': 11, 'figure.figsize': fig_size, 'savefig.dpi': 600,
                  'axes.linewidth': 1.3, 'ytick.major.size': 4,  # major tick size in points
                  'xtick.major.size': 4  # major tick size in points
                  # 'edgecolor' : None
                  # 'xtick.major.size' : 2,
                  # 'ytick.major.size' : 2,
                  }
        rcParams.update(params)

        # set sans-serif font to Arial
        rcParams['font.sans-serif'] = 'Arial'

        # create figure instance
        fig = plt.figure()

        # define sub-panel grid and possibly width and height ratios
        gs = gridspec.GridSpec(3, 1,  # ,
                               # width_ratios=[1.2,1]
                               # height_ratios=[1, 1]
                               )

        # define vertical and horizontal spacing between panels
        gs.update(wspace=0.3, hspace=0.3)

        # possibly change outer margins of the figure
        plt.subplots_adjust(left=0.05, right=0.96, top=0.96, bottom=0.05)

        # sub-panel enumerations
        plt.figtext(0.06, 0.98, '%s   %s   %s ' % (self.mouse, date, rec), clip_on=False, color='black', size=14)

        # plt.figtext(0.06, 0.92, 'A',clip_on=False,color='black', weight='bold',size=22)
        def normalize_data(data):
            """Normalizes the input data to the range of 0 to 1."""
            data_min = np.min(data)
            data_max = np.max(data)
            data_range = data_max - data_min
            normalized_data = (data - data_min) / data_range
            return normalized_data

        chan = ['chan0', 'chan1']
        chanId = ['biceps', 'triceps']
        EMG_col = ['C6', 'C8']
        gssub1 = gridspec.GridSpecFromSubplotSpec(1, 2, subplot_spec=gs[2], hspace=0.4)
        ax4 = plt.subplot(gssub1[0])
        ax4b = plt.subplot(gssub1[1])
        for c in range(2):
            print(c)
            current_key = f'current_{chan[c]}'
            binWidth = 1.E-3  # in sec
            EMGcountwindow = 0.05

            col = ['C0', 'C1', 'C2', 'C3']
            paws = ['FL', 'FR', 'HL', 'HR']
            spacing = [100, 200, 500, 500]
            gssub = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=gs[0], hspace=0.4)
            ax0 = plt.subplot(gssub[0])
            pawMax = [0]
            pawMin = [0]
            i = 0

            # ax0 = plt.subplot(gssub[i])
            pawMin.append(np.min(pawPos[i][:, 1]))
            pawMax.append(np.max(pawPos[i][:, 1]))
            # ax0.plot(pawPos[i][:,0],pawPos[i][:,1],c=col[i],label=paws[i])
            # ax0.plot(EMGTimes,EMGcurrent*300)
            # plt.plot(EMGTimes,currentHP)
            # plt.show()
            idxSwings = swingStanceD['swingP'][i][1]
            indecisiveSteps = swingStanceD['swingP'][i][3]
            recTimes = swingStanceD['forFit'][i][2]
            highpassfreq = 300.
            dt = np.mean(EMG['time'][1:] - EMG['time'][:-1])
            rate = 1. / dt
            print(current_key)
            currentEMG = dataAnalysis.butter_highpass(EMG[current_key], rate, highpassfreq, order=4)
            # Define filter parameters
            fs = 20000  # sampling frequency
            cutoff = 300  # cutoff frequency for low-pass filter (Hz)
            order = 4  # order of the filter

            # Create filter coefficients
            # b, a = butter(order, cutoff / (fs / 2), btype='low')

            b, a = butter(order, cutoff / (fs / 2), 'highpass')
            # Apply filter to EMG signal
            emg_filtered = filtfilt(b, a, EMG[current_key])

            # Smooth the filtered EMG signal using a moving average filter
            window_size = 1  # Number of samples to average over
            EMG_smooth = np.convolve(emg_filtered, np.ones(window_size) / window_size, mode='same')

            # Apply a low-pass filter with a cutoff frequency of 500 Hz
            b, a = signal.butter(4, 500 / (fs / 2), 'lowpass')
            EMG_filtered = signal.filtfilt(b, a, EMG[current_key])

            # Rectify the filtered EMG signal
            EMG_rectified = np.abs(EMG_filtered)

            # Apply a moving average filter with a window length of 100 ms
            window_length = int(fs * 0.1)
            window = np.ones(window_length) / float(window_length)
            EMG_smooth = np.convolve(EMG_rectified, window, mode='same')

            # Scale the smoothed rectified signal by a constant factor

            EMG_norm = normalize_data(EMG_smooth)
            envelope = gaussian_filter1d(EMG_norm, sigma=500)

            t_emg = EMG['time']  # time points for EMG data
            time30s = t_emg <= 30
            # if c == 0:
            #     np.savetxt('biceps_enveloppe.txt', envelope[time30s])
            #     np.savetxt('biceps_raw.txt', EMG[current_key][time30s])
            #     np.savetxt('biceps_time.txt', t_emg[time30s])
            pos_interp = np.interp(t_emg, pawPos[i][:, 0], pawPos[i][:, 1])
            start = 30
            end = 45
            ax1 = plt.subplot(gssub[c + 1])
            N = ['reflect', 'nearest', 'mirror', 'wrap', 'constant']

            ax0.plot(t_emg, pos_interp, label='FL x-pos')
            # ax0.plot(pawPos[i][:,0], uniform_filter1d(newEMGatPawPosTime, size=m))
            # ax1.plot(t_emg, uniform_filter1d(currentEMG, mode='nearest', size=1), alpha=0.5,color=EMG_col[c])
            ax1.plot(t_emg, currentEMG, alpha=0.4, color=EMG_col[c])
            # ax1.plot(t_emg, EMG_denoised, alpha=0.2,color=EMG_col[c])
            ax1.set_xlim(start, end)
            ax0.set_xlim(start, end)
            self.layoutOfPanel(ax0, xLabel='time (s)', yLabel='FL x-pos', xyInvisible=[True, False], Leg=[1, 9])
            self.layoutOfPanel(ax1, xLabel='time (s)', yLabel=f'EMG {chanId[c]} (mV)', xyInvisible=[(True if c == 0 else False), False], Leg=[1, 9])
            gssub1a = gridspec.GridSpecFromSubplotSpec(2, 1, subplot_spec=gs[1], hspace=0.4)
            if c == 0:
                ax2 = plt.subplot(gssub1a[0])
            else:
                ax2 = plt.subplot(gssub1a[1])
            ax3 = ax2.twinx()
            ax3.plot(t_emg, pos_interp, label='FL x-pos', alpha=0.2)
            # pdb.set_trace()
            peaksEMG, _ = find_peaks(envelope)
            peaksPaw, _ = find_peaks(pos_interp)
            valleysEMG, _ = find_peaks(-envelope)
            valleysPaw, _ = find_peaks(-pos_interp)

            ax2.plot(t_emg, (envelope), alpha=0.8, color=EMG_col[c], label=f'EMG {chanId[c]}', )
            ax2.plot(t_emg, (envelope), alpha=0.8, color=EMG_col[c], label=f'EMG {chanId[c]}', )
            self.layoutOfPanel(ax2, xLabel='time (s)', yLabel=f'norm EMG envelope', xyInvisible=[False, False], Leg=[1, 9])
            ax3.spines['left'].set_visible(False)
            ax3.spines['top'].set_visible(False)
            ax3.spines['bottom'].set_visible(False)
            ax2.set_xlim(start, end)
            beforeAfter = [0.4, 0.4]  # in sec
            dt = 0.01
            tbins = np.linspace(-beforeAfter[0], beforeAfter[1], int((beforeAfter[0] + beforeAfter[1]) / dt) + 1, endpoint=True)
            emg_swing_start = []
            emg_stance_start = []

            for n in range(len(idxSwings) - 1):

                idxStart = np.argmin(np.abs(pawPos[0][:, 0] - recTimes[idxSwings[n][0]]))
                idxEnd = np.argmin(np.abs(pawPos[0][:, 0] - recTimes[idxSwings[n][1]]))
                idxStartNext = np.argmin(np.abs(pawPos[0][:, 0] - recTimes[idxSwings[n + 1][0]]))
                start_time = pawPos[0][:, 0][idxStart]
                end_time = pawPos[0][:, 0][idxEnd]
                startNext = pawPos[0][:, 0][idxStartNext]

                idx_start_time = np.argmin(np.abs(t_emg - start_time))
                idx_end_time = np.argmin(np.abs(t_emg - end_time))
                idx_startNext = np.argmin(np.abs(t_emg - startNext))

                maskSwingStart = (t_emg > (t_emg[idx_start_time] - beforeAfter[0])) & (t_emg < (t_emg[idx_start_time] + beforeAfter[1]))
                maskStanceStart = (t_emg > (t_emg[idx_end_time] - beforeAfter[0])) & (t_emg < (t_emg[idx_end_time] + beforeAfter[1]))

                emg_swing_start.append(envelope[maskSwingStart])
                emg_stance_start.append(envelope[maskStanceStart])

                ax0.plot(t_emg[idx_start_time:idx_end_time], pos_interp[idx_start_time:idx_end_time], alpha=0.2, color='red', label='swing period')
                ax1.axvspan(t_emg[idx_start_time], t_emg[idx_end_time], color='C0', alpha=0.05)
                if n == 0:
                    ax0.legend(frameon=False, loc='upper right')
                    ax1.legend(frameon=False, loc='upper right')
            # # pdb.set_trace()
            # cnt, edges = np.histogram(np.concatenate(emg_swing_start).ravel(),
            #                           bins=tbins)
            # timePoints = (edges[1:] + edges[:-1]) / 2
            # swing_centered_emg=cnt/len(emg_swing_start)
            #
            # cnt2, edges2 = np.histogram(np.concatenate(emg_stance_start).ravel(),
            #                           bins=tbins)
            # timePoints2 = (edges2[1:] + edges2[:-1]) / 2
            # stance_centered_emg=cnt2/len(emg_stance_start)
            #

            max_len = max(len(a) for a in emg_swing_start)
            emg_swing_start_padded = [np.pad(a, (0, max_len - len(a)), mode='constant') for a in emg_swing_start]
            avg_swing_centered_emg = np.mean(emg_swing_start_padded, axis=0)
            max_len_stance = max(len(a) for a in emg_stance_start)
            emg_stance_start_padded = [np.pad(a, (0, max_len_stance - len(a)), mode='constant') for a in emg_stance_start]
            avg_stance_centered_emg = np.mean(emg_stance_start_padded, axis=0)
            sem_stance = stats.sem(emg_stance_start_padded, axis=0)
            sem_swing = stats.sem(emg_swing_start_padded, axis=0)

            time_swing_avg = np.linspace(-beforeAfter[0], beforeAfter[1], len(avg_swing_centered_emg))
            time_stance_avg = np.linspace(-beforeAfter[0], beforeAfter[1], len(avg_stance_centered_emg))

            ax4.plot(time_swing_avg, avg_swing_centered_emg, color=EMG_col[c], label=(f'EMG {chanId[c]} n={len(idxSwings)} swings' if c == 0 else f'EMG {chanId[c]}'), lw=2, alpha=0.6)
            ax4b.plot(time_stance_avg, avg_stance_centered_emg, color=EMG_col[c], label=(f'EMG {chanId[c]}'), lw=2)
            ax4.fill_between(time_swing_avg, avg_swing_centered_emg - sem_swing, avg_swing_centered_emg + sem_swing, color=EMG_col[c], alpha=0.3)
            ax4b.fill_between(time_stance_avg, avg_stance_centered_emg - sem_stance, avg_stance_centered_emg + sem_stance, color=EMG_col[c], alpha=0.5)
            # for n in range(len(idxSwings) - 1):
            #     time_swing = np.linspace(-beforeAfter[0], beforeAfter[1], len(emg_swing_start[n]))
            #     time_stance = np.linspace(-beforeAfter[0], beforeAfter[1], len(emg_stance_start[n]))
            #     ax4.plot(time_swing,emg_swing_start[n],color=EMG_col[c], label=(f'EMG {chanId[c]}' if n==0 else None), alpha=0.2)
            #
            #     ax4b.plot(time_stance,emg_stance_start[n], color=EMG_col[c], label=(f'EMG {chanId[c]}' if n==0 else None), alpha=0.2)
            #     if n==0:
            #         ax4.legend(frameon=False, loc='upper right')
            #         ax4b.legend(frameon=False, loc='upper right')
            ax4.axvline(0, ls='--', color='grey', lw=1, alpha=0.1)
            ax4b.axvline(0, ls='--', color='grey', lw=1, alpha=0.1)
            ax4.set_xlim(-0.3, 0.3)
            ax4b.set_xlim(-0.3, 0.3)
            ax4.set_ylim(0.1, 0.6)
            ax4b.set_ylim(0.1, 0.6)
            # ax4.step(timePoints,swing_centered_emg,color=EMG_col[c], label=f'EMG {chanId[c]}', where='mid')
            # ax4b.step(timePoints2, stance_centered_emg, color=EMG_col[c], label=f'EMG {chanId[c]}', where='mid')
            self.layoutOfPanel(ax4, xLabel='time centered on swing (s)', yLabel=f' normalized EMG (avg.)', xyInvisible=[False, False], Leg=[1, 9])
            self.layoutOfPanel(ax4b, xLabel='time centered on stance (s)', yLabel=f' normalized EMG (avg.)', xyInvisible=[False, False], Leg=[1, 9])
        rec = rec.replace('/', '-')
        fname = self.determineFileName(rec, what='EMG', date=date)

        # plt.savefig('caTriggeredAverages_%s.pdf' % caTriggeredAverages[nDays][0])  # define vertical and horizontal spacing between panels  # plt.show()
        # plt.savefig(fname + '.png')
        plt.savefig(fname + '.pdf')

    def createFiberOptoFigure(self, trainData):
        # figure #################################
        fig_width = 20  # width in inches
        fig_height = 16  # height in inches
        fig_size = [fig_width, fig_height]
        params = {'axes.labelsize': 12, 'axes.titlesize': 12, 'font.size': 11, 'xtick.labelsize': 11, 'ytick.labelsize': 11, 'figure.figsize': fig_size, 'savefig.dpi': 600,
                  'axes.linewidth': 1.3, 'ytick.major.size': 4,  # major tick size in points
                  'xtick.major.size': 4  # major tick size in points
                  # 'edgecolor' : None
                  # 'xtick.major.size' : 2,
                  # 'ytick.major.size' : 2,
                  }
        rcParams.update(params)

        # set sans-serif font to Arial
        rcParams['font.sans-serif'] = 'Arial'

        # create figure instance
        fig = plt.figure()

        # define sub-panel grid and possibly width and height ratios
        gs = gridspec.GridSpec(1, 1,  # ,
                               # width_ratios=[1.2,1]
                               # height_ratios=[1, 1]
                               )

        # define vertical and horizontal spacing between panels
        gs.update(wspace=0.3, hspace=0.3)

        # possibly change outer margins of the figure
        plt.subplots_adjust(left=0.05, right=0.96, top=0.96, bottom=0.05)

        # sub-panel enumerations
        # plt.figtext(0.06, 0.98, '%s   %s   %s ' % (self.mouse, date, rec), clip_on=False, color='black', size=14)
        gssub = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=gs[0], hspace=0.4)
        ax0 = plt.subplot(gssub[0])
        ax1= ax0.twinx()

        for train_number, data in trainData.items():
            # pdb.set_trace()
            if train_number=='0':
                time=data['Time']
                avg_signal=data['AverageSignal']
                trigger=data['Trigger']
                sem_signal=data['semSignal']
                sns.lineplot(data=data, x='Time', y='AverageSignal', ax=ax0, color='C2', lw=2, markers=True)
                sns.lineplot(data=data, x='Time', y='Trigger', ax=ax1, color='C3', lw=0.2, markers=True)
                ax1.fill(time, trigger, alpha=0.3, color='C3')
                ax0.fill_between(time, avg_signal - sem_signal, avg_signal +sem_signal,
                                 color='C2', alpha=0.3)
                self.layoutOfPanel(ax0, xLabel='time (s)', yLabel='Delta F/F', xyInvisible=[False, False], Leg=[1, 9])
                self.layoutOfPanel(ax1, xLabel='time (s)', yLabel='Delta F/F', xyInvisible=[False, True], Leg=[1, 9])
                plt.show()

    # plt.figtext(0.06, 0.92, 'A',clip_on=False,color='black', weight='bold',size=22)
    #
    # fname = self.determineFileName(rec, what='EMG', date=date)
    #
    #
    # # plt.savefig(fname + '.png')
    # plt.savefig(fname + '.pdf')

    def createFiberOptoFigureTrials(self, combined_data):
        # figure #################################
        fig_width = 20  # width in inches
        fig_height = 16  # height in inches
        fig_size = [fig_width, fig_height]
        params = {'axes.labelsize': 12, 'axes.titlesize': 12, 'font.size': 11, 'xtick.labelsize': 11, 'ytick.labelsize': 11, 'figure.figsize': fig_size, 'savefig.dpi': 600,
                  'axes.linewidth': 1.3, 'ytick.major.size': 4,  # major tick size in points
                  'xtick.major.size': 4  # major tick size in points
                  # 'edgecolor' : None
                  # 'xtick.major.size' : 2,
                  # 'ytick.major.size' : 2,
                  }
        rcParams.update(params)

        # set sans-serif font to Arial
        rcParams['font.sans-serif'] = 'Arial'

        # create figure instance
        fig = plt.figure()

        # define sub-panel grid and possibly width and height ratios
        gs = gridspec.GridSpec(1, 1,  # ,
                               # width_ratios=[1.2,1]
                               # height_ratios=[1, 1]
                               )

        # define vertical and horizontal spacing between panels
        gs.update(wspace=0.3, hspace=0.3)

        # possibly change outer margins of the figure
        plt.subplots_adjust(left=0.05, right=0.96, top=0.96, bottom=0.05)

        # Sample structure for combined_data
        # combined_data = pd.DataFrame({
        #     'Trial': ['trial_0043', 'trial_0043', 'trial_0044', 'trial_0044', ...],
        #     'Time': [...],
        #     'Signal': [...],
        #     'Trigger': [...],
        #     ...
        # })

        # Function to plot a single trial
        # Function to plot a single trial with specified colors
        def plot_trial(ax0, ax1, trial_data, trial_name, avg_signal_color, trigger_color):
            time = trial_data['Time']
            avg_signal = trial_data['AverageSignal']
            trigger = trial_data['Trigger']
            sem_signal = trial_data['semSignal']

            sns.lineplot(data=trial_data, x='Time', y='AverageSignal', ax=ax0, color=avg_signal_color, lw=2,
                         markers=True, label=f'Trial {trial_name}')
            # sns.lineplot(data=trial_data, x='Time', y='Trigger', ax=ax1, color=trigger_color, lw=0.2, markers=True)
            ax1.fill_between(time, trigger, alpha=0.1, color=trigger_color)
            ax0.fill_between(time, avg_signal - sem_signal, avg_signal + sem_signal, color=avg_signal_color, alpha=0.2)

        # Main plotting
        gs = gridspec.GridSpec(1, 1)
        gssub = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=gs[0], hspace=0.4)
        ax0 = plt.subplot(gssub[0])
        ax1 = ax0.twinx()

        # Get color cycle from 'tab10' colormap
        color_cycle = cycle(plt.cm.tab10.colors)
        allTrialsName= combined_data['Trial'].unique()
        allTrialsNames = '_'.join(combined_data['Trial'].unique())
        # Plot each trial with its respective colors
        for trial_name in combined_data['Trial'].unique():
            trial_data = combined_data[combined_data['Trial'] == trial_name]
            avg_signal_color, trigger_color = next(color_cycle), next(color_cycle)
            plot_trial(ax0, ax1, trial_data, trial_name, avg_signal_color, trigger_color)

        # Setting layout and labels
        ax0.set_xlabel('Time (s)')
        ax0.set_ylabel('Delta F/F')
        ax1.set_ylabel('Trigger Signal')
        self.layoutOfPanel(ax0, xLabel='time (s)', yLabel='Delta F/F', xyInvisible=[False, False], Leg=[1, 9])
        self.layoutOfPanel(ax1, xLabel='time (s)', yLabel='Delta F/F', xyInvisible=[False, True], Leg=[1, 9])
        from datetime import datetime
        # Get today's date
        today = datetime.today()

        # Format the date as a string
        date_string = today.strftime("%Y-%m-%d")

        fname = self.determineFileName(allTrialsNames, what='dligth_Cb-VTA-NAc_batch1', date=date_string)
        print('saved as %s.pdf' % fname)

        # plt.savefig(fname + '.png')
        plt.savefig(fname + '.pdf')


    def plot_optogenetic_response(self,trainData, train_number, winIdxs, durationTrain):
        fig = plt.figure()

        # Define sub-panel grid
        gs = gridspec.GridSpec(1, 1)
        gssub = gridspec.GridSpecFromSubplotSpec(1, 1, subplot_spec=gs[0], hspace=0.4)

        ax0 = plt.subplot(gssub[0])
        ax1 = ax0.twinx()

        # Plotting the signal and trigger
        sns.lineplot(x='Time', y='Signal', data=trainData[train_number],
                     label=f'Trigger at {winIdxs[train_number]["start"]}s, Duration: {durationTrain}s', ax=ax0)
        sns.lineplot(x='Time', y='Trigger', data=trainData[train_number],
                     label=f'Trigger at {winIdxs[train_number]["start"]}s, Duration: {durationTrain}s', ax=ax1, color='red')

        # Setting plot title and labels
        plt.title('Response to Optogenetic Stimulation')
        plt.xlabel('Time Relative to Trigger (s)')
        plt.ylabel('Signal')

        # Mark trigger start
        plt.axvline(x=0, color='black', linestyle='--')
        plt.legend()

        plt.tight_layout()
        plt.show()


    def plot_time_signal(self,time, signal, title='Time vs. Signal', xlabel='Time', ylabel='Signal'):
        plt.figure()
        sns.lineplot(x=time, y=signal)

        plt.title(title)
        plt.xlabel(xlabel)
        plt.ylabel(ylabel)

        plt.tight_layout()
        plt.show()

    import matplotlib.pyplot as plt
    import seaborn as sns


    # Example usage
    # plot_signal_and_trigger(time_data, signal_data, trigger_data, trigger_time
    def plot_fiber_photometry(self,aniDic):
        """
        Plot fiber photometry data over time including signal, reference signal, delta F/F, and trigger events.

        Parameters:
        - aniDic: Dictionary containing 'Signal', 'RefSignal', 'Time', 'Trigger', and 'dFF' keys.
        """
        fig = plt.figure()

        # Define sub-panel grid
        gs = gridspec.GridSpec(1, 1)
        gssub = gridspec.GridSpecFromSubplotSpec(3, 1, subplot_spec=gs[0], hspace=0.4)

        ax0 = plt.subplot(gssub[0])
        ax1 = plt.subplot(gssub[1])
        ax2 = ax1.twinx()


        # Extract data from the dictionary
        time = aniDic['Time']
        signal = aniDic['Signal']
        ref_signal = aniDic['refSignal']
        dFF = aniDic['dFF']
        trigger = aniDic['Trigger']
        triggerTime = aniDic['TriggerTime']
        # Create figure and axis


        # Plot each component on a separate subplot

        ax0.plot(time, signal, label='Signal')
        ax0.plot(time, aniDic['Ypred'], label='Signal', c='red', alpha=0.1)
        ax0.set_title('Signal')
        ax0.set_ylabel('Fluorescence')
        # plt.show()
        # pdb.set_trace()
        if ref_signal is not None:
            ax0.plot(time, ref_signal, label='Ref Signal', color='orange')


        ax1.plot(time, dFF, label='ΔF/F', color='green')
        ax1.set_title('Delta F/F')
        ax1.set_ylabel('ΔF/F')

        ax2.plot(triggerTime, trigger, label='Trigger', color='red')
        ax2.set_title('Trigger')
        ax2.set_xlabel('Time (s)')
        ax2.set_ylabel('Trigger State')
        self.layoutOfPanel(ax0, xLabel='time (s)', yLabel='raw Fluorescence', xyInvisible=[False, False], Leg=[1, 9])
        self.layoutOfPanel(ax1, xLabel='time (s)', yLabel='Delta F/F', xyInvisible=[False, False], Leg=[1, 9])
        self.layoutOfPanel(ax2, xLabel='', yLabel='trigger', xyInvisible=[True, False], Leg=[1, 9])
        # plt.tight_layout()

        # Show plot
        plt.show()