<jupyter_start><jupyter_text>EEG Analysis Generating EEG From Simulation Dipolemoments<jupyter_code># Import packages

import numpy as np               # array manipulation
np.seterr(divide='ignore',       # ignore warnings
          invalid='ignore')      

from scipy import stats as st    # stats
from scipy import signal as ss
import matplotlib                # plotting
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid.inset_locator import (inset_axes, InsetPosition,mark_inset) # insets for plots

# EEG Model Fitting
from fooof import FOOOF          # single PSD model
from fooof import FOOOFGroup     # group PSD model

# Import utilities to manage frequency band definitions
from fooof.bands import Bands   
from fooof.analysis import get_band_peak_fg

# Import plotting function for model parameters and components
from fooof.plts.periodic import plot_peak_fits, plot_peak_params
from fooof.plts.aperiodic import plot_aperiodic_params, plot_aperiodic_fits
from fooof.plts.spectra import plot_spectrum
from fooof.sim.gen import gen_aperiodic
from fooof.plts.spectra import plot_spectrum
from fooof.plts.annotate import plot_annotated_model
from fooof.objs.utils import average_fg<jupyter_output><empty_output><jupyter_text>Calculating EEG From 4-sphere Head Model Example (Cannot run)<jupyter_code># Load in Data

DIPOLEMOMENTS_H = [np.load(dp) for dp in <path-to-Healthy-DIPOLEMOMENTS>]
DIPOLEMOMENTS_D = [np.load(dp) for dp in <path-to-Depression-DIPOLEMOMENTS>]

DIPOLEMOMENTS = [DIPOLEMOMENTS_H,DIPOLEMOMENTS_D]

#create empty dicts to hold data
EEG_h,EEG_d = {},{}

#params for forward solution
radii = [79000., 80000., 85000., 90000.]           # pia, CSF, skull, scalp (um)
sigmas = [0.47, 1.71, 0.02, 0.41]                  # ^ conductances (S/m)
rz = np.array([0., 0., 78250.])                    # dipole location, center of L23 (um)
r2 = np.array([[0., 0., 90000]])                   # EEG sensor location (um)

EEG_args = LFPy.FourSphereVolumeConductor(radii, sigmas, r2) #set 4sphere model

for i, DIPOLEMOMENT in enumerate(DIPOLEMOMENTS):   # for each condition
    for j, run in enumerate(DIPOLEMOMENT):         # for each run in condition
        dp = np.add(np.add(run['HL23PYR'],run['HL23SST']),np.add(run['HL23PV'],run['HL23VIP'])) # add dipole of each pop
        EEG = EEG_args.calc_potential(dp, rz)      # calc EEG timeseries from summed dipolemoments

        key = str(1000+j) if i==0 else str(2000+j) # create key=SEED
        d = {str(key):{'ts_raw':EEG.flatten()}}
        EEG_h.update(d) if i==0 else EEG_d.update(d)

sf = 40000                                         # sampling frequency
nperseg_welch = sf*3                               # 3s window
noverlap_welch = perseg_welch*.5                   # 50% overlap
end_freq = 120                                     # frequency to cut power spectrum array, to save space

for e,signal in enumerate([EEG_h,EEG_d]):          # for each condition
    for name, value in signal.items():             # for each run
        for timeseries,ts in value.items():        # for each timeseries
            w = {}
            ts = ts.flatten()[80000:]              # flatten and remove transient 2s * 40000 data points / second

            # Welch
            freq_wel,ps_wel = ss.welch(ts, fs=sf, nperseg=nperseg_welch, noverlap=noverlap_welch)
            
            end_wel = int(end_freq/freq_wel[1])    # gives index of end freq
            w = {'ps_wel':ps_wel[:end_wel]}

        value.update(w)                            # update dictionary
        
        
        
# Output:

# EEG_h:{'1000':'ts_raw':[],'ps_wel':[],
#        '1001':'ts_raw':[],'ps_wel':[],
#        '1002':'ts_raw':[],'ps_wel':[],
#        '1003':'ts_raw':[],'ps_wel':[],
#        ...
#        '1059':'ts_raw':[],'ps_wel':[]}
    
# EEG_d:{'2000':'ts_raw':[],'ps_wel':[],
#        '2001':'ts_raw':[],'ps_wel':[],
#        '2002':'ts_raw':[],'ps_wel':[],
#        '2003':'ts_raw':[],'ps_wel':[],
#        ...
#        '2059':'ts_raw':[],'ps_wel':[]}<jupyter_output><empty_output><jupyter_text>Dipolemoment -> EEG<jupyter_code># Load in data

# DIPOLEMOMENTS
PN = np.load('Data/Calculating_EEG/PN_DIPOLEMOMENT.npy')
SST = np.load('Data/Calculating_EEG/SST_DIPOLEMOMENT.npy')
PV = np.load('Data/Calculating_EEG/PV_DIPOLEMOMENT.npy')
VIP = np.load('Data/Calculating_EEG/VIP_DIPOLEMOMENT.npy')


# Corresponding EEG
EEG_H = np.load('Data/Healthy/EEG_H.npy')
# Plot Dipoles and EEG
fig,axes = plt.subplots(figsize=(15,7),ncols=1,nrows=5)

x_dip = np.arange(2000,4000,.025)
x_eeg = np.arange(2000,30000+.025,.025)
names = ['Pyr','SST','PV','VIP','EEG']

# Dipoles
axes[0].plot(x_dip,PN[:,2],c='k')#,DIPOLEMOMENT['HL23PN1'][:,2][80000:],c='k')      # Pyr z dipole
axes[1].plot(x_dip,SST[:,2],c='red')    # SST z dipole
axes[2].plot(x_dip,PV[:,2],c='green')  # PV z dipole
axes[3].plot(x_dip,VIP[:,2],c='orange') # VIP z dipole

# EEG
axes[4].plot(x_eeg,EEG_H[1][80000:])
axes[4].set_ylabel('Potential (mV)',fontsize=14)

# Formatting
fs,titlefs = 14,18
[ax.set_title(name,loc='left',fontsize=fs) for ax,name in zip(axes,names)]
[ax.set_xlim(2000,4000) for ax in axes]
[ax.set_xticklabels([]) for ax in axes[:-1]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in axes]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in axes]
axes[-1].ticklabel_format(style='plain')
fig.text(x=.01,y=.4,s='Current Dipole Moment (nA/um)',rotation=90,fontsize=fs)

plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Example: Healthy Condition 1) PSD using Welch's Method for One Simulation<jupyter_code># Load in Data
#EEG_H = np.load('Data/Healthy/EEG_H.npy')

# Set welch method parameters
ts = EEG_H[0][80000:]                         # EEG timeseries after transient
sampling_frequency = 40000                    # .025 ms = .000025 = 1/40000
nperseg = sampling_frequency * 3              # 1 s window
noverlap = nperseg*.5                         # overlap of windows

freq,Pxx = ss.welch(ts,fs=sampling_frequency,nperseg=nperseg,noverlap=noverlap)

print('Frequency resolution is',freq[1],'Hz')
# Set frequency bounds
endfreq = 100                                 # end freq
endfreq_ind = np.where(freq==endfreq)[0][0]   # index of end freq

freq_trim = freq[:endfreq_ind]
Pxx_trim = Pxx[:endfreq_ind]
print('Frequency resolution is',round(freq[1],2),'Hz, index of',endfreq,'Hz in frequency vector is position',endfreq_ind)
# Plot Single run PSD
fig = plt.figure(figsize=(8,6))        
axes = plt.subplot(111)

# plot PSD in main plot
axes.plot(freq_trim,Pxx_trim,c='k')

# plot PSD in inset with log-log axes
inset = axes.inset_axes([.5,.5,.5,.5])
inset.plot(freq_trim,Pxx_trim,c='k')

# Formatting
axes.set_title('Healthy Power Spectral Density',fontsize=titlefs,y=1.05)
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in [axes,inset]] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in [axes,inset]]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in [axes,inset]]
axes.set_xlabel('Frequency (Hz)',fontsize=fs)
axes.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
inset.set_xscale('log')
inset.set_yscale('log')
inset.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
inset.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
axes.set_xticks([5,10,15,20,25,30])
axes.set_xlim(4,30)
inset.set_xlim(4,100)
inset.set_yticklabels([])
labels = [r'$\theta$',r'$\alpha$',r'$\beta$']
[axes.annotate(xy=(x,0),text=t,fontsize=fs) for x,t in zip([5.5,9.5,22],labels)]
[inset.annotate(xy=(x,.5e-16),text=t,fontsize=fs) for x,t in zip([5,8.5,22],labels)]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>2) Describing Single PSD with FOOOF<jupyter_code># Define model parameters
H1 = FOOOF(peak_width_limits = (1.5,12.0),  # Limits on possible peak width, in Hz, as (lower_bound, upper_bound)
           max_n_peaks = 3,                 # Maximum number of peaks to fit.
           min_peak_height = 0,             # Absolute threshold for detecting peaks, in units of the input data.
           peak_threshold = 2,              # Relative threshold for detecting peaks, in units of standard deviation of input data
           aperiodic_mode = 'fixed')        # Which approach to take for fitting the aperiodic component.
# Set model bounds
startfreqind = np.where(freq==1)[0][0]
endfreqind = np.where(freq==30)[0][0]

# Fit model
H1.fit(freqs=freq[startfreqind:endfreqind],          # frequencies to fit model
      power_spectrum = Pxx[startfreqind:endfreqind]) # power spectra to fit
# Plot original data with fit
H1.plot(plt_log=True) # change plt_log=False
# Plot annotated model
plot_annotated_model(H1, annotate_aperiodic=False,plt_log=True)
# Print out model fit results
print('\naperiodic params:\n')
print('    Offset       Exponent') # Offset = y-int, Exponent=slope of aperiodic
print(H1.aperiodic_params_)

print('\npeak params:\n')
print(' Center Freq     Power      Bandwidth') # Center of gaussian, peak of gaussian above aperiodic, width of gaussian
print(H1.peak_params_)

print('\nr-squared:', H1.r_squared_) # coefficient of determination, how much data variation is determined by model

print('\nfit error:', H1.error_) # error of model
# Generate Aperiodic and Periodic components of power spectra
H_AP_fit = gen_aperiodic(H1.freqs,H1.aperiodic_params_)        # generate aperiodic fit to absolute

H_P_fit = H1.power_spectrum-H_AP_fit                           # generate periodic PS by subtracting aperiodic from absolute
# Plot components and full PSD together
fig = plt.figure(figsize=(15,6))

absolute = plt.subplot(131)
aperiodic = plt.subplot(132)
periodic = plt.subplot(133)

absolute.plot(freq[:endfreq_ind],Pxx[:endfreq_ind],color='k') # plot absolute
aperiodic.plot(H1.freqs,H_AP_fit,color='k')                   # plot aperiodic
periodic.plot(H1.freqs,H_P_fit,color='k')                     # plot periodic

axes = [absolute,aperiodic,periodic]

# Formatting
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in axes] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in axes]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in axes]
[ax.set_xticks([5,10,15,20,25,30]) for ax in axes]
_=[ax.set_yticklabels([]) for ax in axes[1:]]
absolute.set_title('Absolute PSD',fontsize=titlefs)
aperiodic.set_title('Aperiodic PSD',fontsize=titlefs)
periodic.set_title('Periodic PSD',fontsize=titlefs)
[ax.set_xlabel('Frequency (Hz)',fontsize=fs) for ax in axes]
absolute.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
aperiodic.set_ylabel('Power (a.u)',fontsize=fs)
[ax.set_xlim(4,30) for ax in axes]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>3) Mean PSD for Healthy Condition<jupyter_code># Load in data
EEG_H = np.load('Data/Healthy/EEG_H.npy')

print(len(EEG_H),'timeseries loaded.')
# Generate PSDs
ts = EEG_H[0][80000:]                         # EEG timeseries after transient
sampling_frequency = 40000                    # .025 ms = .000025 = 1/40000
nperseg = sampling_frequency * 3              # 1 s window
noverlap = nperseg*.5                         # overlap of windows

H_PS = [0,0,0,0,0]                            # hold power spectra

for r,ts in enumerate(EEG_H):
    print('Calculating PSD for simulation',r)
    freq,Pxx = ss.welch(ts[80000:],fs=sampling_frequency,nperseg=nperseg,noverlap=noverlap)
    
    endfreqind = np.where(freq==100)[0][0]    # end index
    
    H_PS[r] = Pxx[:endfreqind]

freq_trim = freq[:endfreqind]
print('\nFrequency resolution is',round(freq_trim[1],2),'Hz')

H_PS = np.array(H_PS)
# Check if correct shape
H_PS.shape # 5 runs, 300 datapoints
# Calculate mean and standard deviation
H_PS_mean = np.mean(H_PS,axis=0) # mean for freq power
H_PS_sd = np.std(H_PS,axis=0)    # sd for each freq power
# Plot group-averaged PSD with mean and shaded sd

fig = plt.figure(figsize=(8,6))
ax1 = plt.subplot(111)

# Set values up for plotting
x = freq[:endfreqind]
y = H_PS_mean
lower_sd = H_PS_mean - H_PS_sd   # fill in later
upper_sd = H_PS_mean + H_PS_sd   # fill in later

# plot PSD in main plot
ax1.plot(x,y,color='k')                                    # plot mean
ax1.fill_between(x,lower_sd,upper_sd,color='k',alpha=.15)  # shade sd

# plot PSD in inset with log-log axes
inset = ax1.inset_axes([.5,.5,.5,.5])                      # inset for plot
inset.plot(x,y,color='k')                                  # same mean as above
inset.fill_between(x,lower_sd,upper_sd,color='k',alpha=.15)# same sd as above

# Formatting
ax1.set_title('Healthy Power Spectral Density',fontsize=titlefs,y=1.05)
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in [ax1,inset]] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in [ax1,inset]]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in [ax1,inset]]
ax1.set_xlabel('Frequency (Hz)',fontsize=fs)
ax1.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
inset.set_xscale('log')
inset.set_yscale('log')
inset.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
inset.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax1.set_xticks([5,10,15,20,25,30])
ax1.set_xlim(4,30)
inset.set_xlim(4,100)
inset.set_yticklabels([])
labels = [r'$\theta$',r'$\alpha$',r'$\beta$']
[ax1.annotate(xy=(x,0),text=t,fontsize=fs) for x,t in zip([5.5,9.5,22],labels)]
[inset.annotate(xy=(x,.5e-16),text=t,fontsize=fs) for x,t in zip([5,8.5,22],labels)]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>4) Describing Multiple PSDs with FOOOF (FOOOF Group Analysis)<jupyter_code># Load in Data
# EEG_H = np.load('Data/Healthy/EEG_H.npy')

# Calculate PSD for each condition
sampling_frequency = 40000                    # .025 ms = .000025 = 1/40000
nperseg = sampling_frequency * 3              # 1 s window
noverlap = nperseg*.5                         # overlap of windows

H_PS = [0,0,0,0,0]                            # hold power spectra
endfreqind = np.where(freq==100)[0][0]        # end index

for r,ts in enumerate(EEG_H):
    freq,Pxx = ss.welch(ts[80000:],fs=sampling_frequency,nperseg=nperseg,noverlap=noverlap)
    H_PS[r] = Pxx[:endfreqind]

H_PS = np.array(H_PS)

print('Frequency resolution is',round(freq[1],2))
# Fits all power spectra in a group of power spectra

H = FOOOFGroup(peak_width_limits=[2, 12.0],  
            max_n_peaks= 3,                 # aprox 1 peak / freq band
            min_peak_height=0,
            peak_threshold=2.0,
            aperiodic_mode='fixed',
            verbose = True)

startfreqind = np.where(freq==3)[0][0]
endfreqind = np.where(freq==30)[0][0]

H.fit(freqs=freq[startfreqind:endfreqind],              # frequencies to fit model
      power_spectra = H_PS[:,startfreqind:endfreqind])  # power spectra to fit     
# Print out model fit results
print('\naperiodic params:\n')
print('    Offset       Exponent') # Offset = y-int, Exponent=slope of aperiodic
print(H.get_params('aperiodic'))

print('\npeak params:\n')
print(' Center Freq     Power     Bandwidth    Sim #') # Center of gaussian, peak of gaussian above aperiodic, width of gaussian
print(H.get_params('peak_params'))
# Plot fits and group average

fig,axes = plt.subplots(figsize=(16,6),nrows=1,ncols=4,sharey=True) #axes = [subplot,subplot,subplot,subplot]

[H.get_fooof(ind=i).plot(ax=axes[i],add_legend=False) for i in range(3)] # list of first 3 fits in group

# Define frequency bands of interest
bands = Bands({'theta' : [4, 8],
               'alpha' : [8, 12],
               'beta' : [12, 30]})

h_av = average_fg(H, bands, avg_method='mean')             # average all model fits using predefined bands
h_av.plot(plot_peaks='shade',add_legend=False,ax=axes[3])  # plot group average 


# Formatting
[ax.set_yticklabels([]) for ax in axes[1:]]
[ax.set_ylabel('') for ax in axes[1:]]
_=[ax.set_title(t,fontsize=titlefs) for ax,t in zip(axes,['Sim 1','Sim 2','Sim 3','Group Average'])]
# Generate Aperiodic and Periodic fits for each simulation

# Aperiodic
H_APs = np.zeros_like(H_PS[:,startfreqind:endfreqind]) # hold aperiodic PS

for i in range(5):
    H_APs[i] = gen_aperiodic(H.freqs,H.get_params('aperiodic')[i])

# Periodic
H_Ps = np.zeros_like(H_PS[:,startfreqind:endfreqind]) # hold Periodic PS

for i in range(5):
    PS_abs = H.get_fooof(ind=i).power_spectrum
    H_Ps[i] = PS_abs - H_APs[i]

print('Absolute power spectra:',H_PS[:,startfreqind:endfreqind].shape)
print('Aperiodic power spectra:',H_APs.shape)
print('Periodic power spectra:',H_APs.shape)
# Average over simulations
H_PS_mean = np.mean(H_PS[:,startfreqind:endfreqind],axis=0)
H_PS_sd = np.std(H_PS[:,startfreqind:endfreqind],axis=0)

H_APs_mean = np.mean(H_APs,axis=0)
H_APs_sd = np.std(H_APs,axis=0)

H_Ps_mean = np.mean(H_Ps,axis=0)
H_Ps_sd = np.std(H_Ps,axis=0)

print('Mean absolute power spectra:',H_PS_mean.shape)
print('Mean Aperiodic power spectra:',H_APs_mean.shape)
print('Mean Periodic power spectra:',H_Ps_mean.shape)
# Setup values for plotting

# shared for all
x = H.freqs

# abs
H_PS_sd_l = H_PS_mean - H_PS_sd
H_PS_sd_u = H_PS_mean + H_PS_sd

# aperiodic
H_APs_sd_l = H_APs_mean - H_APs_sd
H_APs_sd_u = H_APs_mean + H_APs_sd

# periodic
H_Ps_sd_l = H_Ps_mean - H_Ps_sd
H_Ps_sd_u = H_Ps_mean + H_Ps_sd
# Plot Group Averaged Power Spectra
fig = plt.figure(figsize=(15,6))

absolute = plt.subplot(131)
aperiodic = plt.subplot(132)
periodic = plt.subplot(133)

absolute.plot(x,H_PS_mean,color='k')
absolute.fill_between(x,H_PS_sd_l,H_PS_sd_u,color='k',alpha=.15) 

aperiodic.plot(x,H_APs_mean,color='k')           
aperiodic.fill_between(x,H_APs_sd_l,H_APs_sd_u,color='k',alpha=.15) 

periodic.plot(x,H_Ps_mean,color='k')   
periodic.fill_between(x,H_Ps_sd_l,H_Ps_sd_u,color='k',alpha=.15) 



# Formatting
axes = [absolute,aperiodic,periodic]
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in axes] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in axes]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in axes]
[ax.set_xticks([5,10,15,20,25,30]) for ax in axes]
_=[ax.set_yticklabels([]) for ax in axes[1:]]
absolute.set_title('Absolute PSD',fontsize=titlefs)
aperiodic.set_title('Aperiodic PSD',fontsize=titlefs)
periodic.set_title('Periodic PSD',fontsize=titlefs)
[ax.set_xlabel('Frequency (Hz)',fontsize=fs) for ax in axes]
absolute.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
aperiodic.set_ylabel('Power (a.u)',fontsize=fs)
[ax.set_xlim(4,30) for ax in axes]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Compare Healthy and Depression Microcircuits  Task 1: Compare Healthy and Depression PSDs<jupyter_code># Load in Data

# Calculate PSD for each condition 
sampling_frequency = 40000                    # .025 ms = .000025 = 1/40000
nperseg = sampling_frequency * 3              # 1 s window
noverlap = nperseg*.5                         # overlap of windows

# Healthy

# Depression

# Check if correct shape

# Calculate mean and standard deviation for each condition

# Plot group-averaged PSD with mean and shaded sd
fig = plt.figure(figsize=(8,6))
ax1 = plt.subplot(111)                 # main plot
inset = ax1.inset_axes([.5,.5,.5,.5])  # inset with log-log axes

# Healthy
h_x = freq[:endfreqind]
h_y = H_PS_mean
h_lower_sd = H_PS_mean - H_PS_sd   # fill in later
h_upper_sd = H_PS_mean + H_PS_sd   # fill in later

ax1.plot(h_x,h_y,color='k')                  
ax1.fill_between(h_x,h_lower_sd,h_upper_sd,color='k',alpha=.15)
inset.plot(h_x,h_y,color='k')
inset.fill_between(h_x,h_lower_sd,h_upper_sd,color='k',alpha=.15)

# Depression
# plot mean depression PSD +/- sd in the same subplot as healhty
d_x = 
d_y = 
d_lower_sd = 
d_upper_sd = 

dclr = (0.8, 0.0, 0.45)

ax1.plot(d_x,d_y,color=dclr)
ax1.fill_between(d_x,d_lower_sd,d_upper_sd,color=dclr,alpha=.15)
inset.plot(d_x,d_y,color=dclr)
inset.fill_between(d_x,d_lower_sd,d_upper_sd,color=dclr,alpha=.15)

# Formatting
ax1.set_title('Healthy Power Spectral Density',fontsize=titlefs,y=1.05)
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in [ax1,inset]] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in [ax1,inset]]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in [ax1,inset]]
ax1.set_xlabel('Frequency (Hz)',fontsize=fs)
ax1.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
inset.set_xscale('log')
inset.set_yscale('log')
inset.get_xaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
inset.get_yaxis().set_major_formatter(matplotlib.ticker.ScalarFormatter())
ax1.set_xticks([5,10,15,20,25,30])
ax1.set_xlim(4,30)
inset.set_xlim(4,100)
inset.set_yticklabels([])
labels = [r'$\theta$',r'$\alpha$',r'$\beta$']
[ax1.annotate(xy=(x,0),text=t,fontsize=fs) for x,t in zip([5.5,9.5,22],labels)]
[inset.annotate(xy=(x,.5e-16),text=t,fontsize=fs) for x,t in zip([5,8.5,22],labels)]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Task 2: Compare Healthy and Depression FOOOF models<jupyter_code># Load in Data

# Calculate PSD for each condition

# Fits all power spectra in a group of power spectra

# Print out model fit results

# Plot fits and group average
fig,axes = plt.subplots(figsize=(16,6),nrows=2,ncols=4,sharey=True)

# Define frequency bands of interest
bands = Bands({'theta' : [4, 8],
               'alpha' : [8, 12],
               'beta' : [12, 30]})

# Healthy
[H.get_fooof(ind=i).plot(ax=axes[0][i],add_legend=False) for i in range(3)] 
h_av = average_fg(H, bands, avg_method='mean')                # average all model fits using predefined bands
h_av.plot(plot_peaks='shade',add_legend=False,ax=axes[0][3])  # plot group average 

# Depression
# axes for individual simulation plots are axes[1][i]
# axes for group average plot is axes[1][3]


[ax.set_yticklabels([]) for ax in axes[0][1:]]
[ax.set_ylabel('') for ax in axes[0][1:]]
[ax.set_ylabel('') for ax in axes[1][1:]]
[ax.set_xlabel('') for ax in axes[0]]
_=[ax.set_title(t,fontsize=titlefs) for ax,t in zip(axes[0],['Sim 1','Sim 2','Sim 3','Group Average'])]
# Generate Aperiodic and Periodic fits for each simulation

# Aperiodic


# Periodic


# Average over simulations

#--------------Indexing healthy, do not remove---------------
H_PS_mean = np.mean(H_PS[:,startfreqind:endfreqind],axis=0)
H_PS_sd = np.std(H_PS[:,startfreqind:endfreqind],axis=0)
#------------------------------------------------------------


# Setup values for plotting

# abs

# aperiodic

# periodic

# Plot Group Averaged Power Spectra
fig = plt.figure(figsize=(15,6))

absolute = plt.subplot(131)    # absolute subplot for healthy and depression
aperiodic = plt.subplot(132)   # aperiodic suboplot for healthy and depression
periodic = plt.subplot(133)    # periodic suboplot for healthy and depression

# Healthy
# for all - plot mean, then fill between mean+sd, mean-sd

absolute.plot(x,H_PS_mean,color='k')
absolute.fill_between(x,H_PS_sd_l,H_PS_sd_u,color='k',alpha=.15) 

aperiodic.plot(x,H_APs_mean,color='k')           
aperiodic.fill_between(x,H_APs_sd_l,H_APs_sd_u,color='k',alpha=.15) 

periodic.plot(x,H_Ps_mean,color='k')   
periodic.fill_between(x,H_Ps_sd_l,H_Ps_sd_u,color='k',alpha=.15) 


# Depression
# color = dclr


# Formatting
axes = [absolute,aperiodic,periodic]
[[ax.axvline(i,0,1, linestyle='--',color='k',alpha=.4) for ax in axes] for i in [4,8,12]]
[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in axes]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in axes]
[ax.set_xticks([5,10,15,20,25,30]) for ax in axes]
_=[ax.set_yticklabels([]) for ax in axes[1:]]
absolute.set_title('Absolute PSD',fontsize=titlefs)
aperiodic.set_title('Aperiodic PSD',fontsize=titlefs)
periodic.set_title('Periodic PSD',fontsize=titlefs)
[ax.set_xlabel('Frequency (Hz)',fontsize=fs) for ax in axes]
absolute.set_ylabel(r'Power $\left( \frac{ V^{2} }{Hz} \right)$',fontsize=fs)
aperiodic.set_ylabel('Power (a.u)',fontsize=fs)
[ax.set_xlim(4,30) for ax in axes]
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Bonus: Frequency band power change<jupyter_code># labels
colors = ['k',dclr]
labels = ['Healthy','Depression']

# Define the frequency range for our simulations
freq_range = [freq[startfreqind],freq[endfreqind]]

# Define frequency bands of interest
bands = Bands({'theta' : [4, 8],
               'alpha' : [8, 12],
               'beta' : [12, 30]})

fig = plt.figure(figsize=(18,6))
for b,band in enumerate(bands.bands.items()):
    
    h_b = get_band_peak_fg(H, band[1]) # get peaks from Healthy, defined by band limits
    d_b = 
    
    ax1 = plt.subplot(1,3,b+1) # plot peaks
    
    bs = [h_b,d_b]
    freq_range = [band[1][0]-4,band[1][1]+4] # xrange
    plot_peak_params(bs, freq_range=freq_range,labels=labels, colors=colors,ax=ax1) # plot band power
    
    ax1.set_title(band[0],fontsize=titlefs)
    ax1.set_ylim(0,1)<jupyter_output><empty_output><jupyter_text>Calculate significance of band power difference<jupyter_code># View params
H.get_params('peak_params')
# center frequency, power, bandwidth, sim #
# ttest between values

# Healthy condition
healthy = []
bandlims = [4,8]

# collect data
for peak in H.get_params('peak_params'):
    if bandlims[0]<peak[0]<bandlims[1]:
        healthy.append(peak[1])
        
# depression
depression = []


# compute ttest
t,p = st.ttest_ind(healthy,depression)

print('Ttest for peak power',bandlims[0],'-',bandlims[1],'Hz\n')
if p<0.05:
    print('p =',p,'SIG')
else:
    print('p =',p,'nonsig')<jupyter_output><empty_output><jupyter_text>Aperiodic parameter difference<jupyter_code># View params
H.get_params('aperiodic_params')
# center frequency, power, bandwidth, sim #
# Plot Params
offset_H = H.get_params('aperiodic_params')[:,0]
exponent_H = H.get_params('aperiodic_params')[:,1]

fig,axes = plt.subplots(figsize=(8,5),ncols=2,nrows=1)
axes[0].scatter(np.zeros_like(offset_H),offset_H,c='k',s=100,alpha=.5)
axes[0].scatter(np.ones_like(offset_D),offset_D,c=dclr,s=100,alpha=.5)
axes[1].scatter(np.zeros_like(exponent_H),exponent_H,c='k',s=100,alpha=.5)
axes[1].scatter(np.ones_like(exponent_D),exponent_D,c=dclr,s=100,alpha=.5)

[ax.spines[side].set_visible(False) for side in ['right','top'] for ax in axes]
[ax.tick_params(axis='both', which='major', labelsize=fs) for ax in axes]
[ax.set_xlim(-.5,1.5) for ax in axes]
[ax.set_xticks([0,1]) for ax in axes]
[ax.set_xticklabels(['Healthy','Depression']) for ax in axes]

axes[0].set_ylabel('Offset',fontsize=titlefs)
axes[1].set_ylabel('Exponent',fontsize=titlefs)
plt.suptitle('Aperiodic Parameter Comparison',fontsize=titlefs)

plt.tight_layout()
# ttest between values

# Healthy condition
offset_H = H.get_params('aperiodic_params')[:,0]
exponent_H = H.get_params('aperiodic_params')[:,1]


# Depression condition

# compute ttest
t,p = st.ttest_ind(exponent_H,exponent_D)

if p<0.05:
    print('p =',p,'SIG')
else:
    print('p =',p,'non')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15429435409521702, 'fraction_numerical': 0.032415090164841245, 'mean_word_length': 4.376288659793815, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '197987', 'license': 'MIT', 'n_tokens_mistral': 11066, 'n_tokens_neox': 10228, 'n_words': 2086, 'path': 'day4/Microcircuit_Modelling/Session_2_EEG.ipynb', 'repo_name': 'edickie/kcni-school-lessons', 'size': '27117'}