OSCR

The dominance of large-scale phase dynamics in human cortex, from delta to gamma.

Code ↔ Paper

7 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 7 matches
  1. [1] § Methods › SVD for empirical fourier decomposition ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 264–281 · score 0.61 · inferior superior, anterior posterior, sensor, MEG, SF, phase
  2. [2] § Methods › Multi-scale differencing of phase ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 689–790 · score 0.60 · smallest triangle, geodesic distances, edge, bin, vectors, weight
  3. [3] § Results › Overview ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 594–624 · score 0.60 · inferior superior, anterior posterior, Wave maps, angle, vector, phase
  4. [4] § Methods › Multi-scale differencing of phase ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 308–340 · score 0.59 · approximately equilateral, sEEG, collated, hemisphere, triplet, angle
  5. [5] § Methods › Numerical methods › Phase estimation ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 159–234 · score 0.58 · logarithmically spaced, Morlet wavelets, oversampling, cycle, phase
  6. [6] § Methods › Surrogate testing ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 689–790 · score 0.54 · singular vector, sEEG, wavelength, empirical, weightings, cortical
  7. [7] § Methods › Multi-scale differencing of phase ↔ estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb, lines 13–20 · score 0.50 · inter contact distances, sEEG, irregularity, MEG, phase, SF

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 793 lines · 46 KB · GPL-3.0 · 7 matches

  1. # %%
  2. import numpy as np
  3. #import cupy as cp #uncomment this if cuda gpu available and cupy installed
  4. from matplotlib.pyplot import *
  5. import matplotlib.ticker as ticker
  6. from scipy.optimize import curve_fit
  7. from pickle import dump,load,HIGHEST_PROTOCOL
  8. from os.path import exists
  9. import warnings
  10. from glob import glob
  11. from IPython.display import clear_output
  13. # %%
  14. ###select one###
  16. ###sEEG inter-contact distances have been computed using geodesic distances on surface models of cortex###
  17. measurement = 'sEEG' #example subjects
  19. ###MEG inter-contact distances have been computed from a spheroid approximation of the sensory array###
  20. #measurement = 'MEG' #example subjects
  23. # %%
  24. #the data-files used in these examples can be downloaded from https://osf.io/48ztm/files/osfstorage#
  26. if measurement=='sEEG':
  27. file_path = 'E:\\sensor_coordinates_distance_data_files\\' #path to where you put the input files listed below
  28. output_path = 'E:\\output\\' #path where you want the output files to go (keep short, output files have long names!)
  29. subjects = ['subject15','subject22'] #example subjects
  30. geodesic_distance_files = ['subject15_depth_geodesic_distances','subject22_depth_geodesic_distances'] #file-name without '.pkl' affix
  31. triangle_files = ['subject15_depth_equilateral_triplets','subject22_depth_equilateral_triplets'] #file-name without '.pkl' affix
  32. coordinate_files = ['subject15_depth_coordinatesAP-IS-LR','subject22_depth_coordinatesAP-IS-LR'] #file-name without '.csv' affix
  33. TS_files = ['subject15_depth_timeseries','subject22_depth_timeseries'] #file-name without '.npy' affix
  34. tdel_values = [1.0,1.0] # one over the sampling frequency, in ms.
  35. frequencies = [1.0, 1.148698354997035, 1.3195079107728942, 1.5157165665103982, 1.7411011265922482,
  36. 2.0, 2.29739670999407, 2.639015821545789, 3.0314331330207964, 3.4822022531844965,
  37. 4.0, 4.59479341998814, 5.278031643091579, 6.062866266041593, 6.964404506368995,
  38. 8.0, 9.18958683997628, 10.556063286183157, 12.125732532083186, 13.92880901273799,
  39. 16.0, 18.37917367995256, 21.112126572366314, 24.25146506416638, 27.857618025475986,
  40. 32.0, 36.75834735990512, 42.22425314473263, 48.50293012833276, 55.71523605095197,
  41. 64.0, 73.51669471981025, 84.44850628946526, 97.00586025666551] #temporal frequencies to analyze
  42. min_of_range = 1.0/50.0 #m, for plotting
  43. elif measurement=='MEG':
  44. file_path = 'E:\\sensor_coordinates_distance_data_files\\' #path to where you put the input files listed below
  45. output_path = 'E:\\output\\' #path where you want the output files to go (keep short, output files have long names!)
  46. subjects = ['subjectG0','subjectP9'] #example subjects
  47. geodesic_distance_files = ['subjectG0_MEG_spheroid_distances','subjectP9_MEG_spheroid_distances'] #file-name without '.pkl' affix
  48. triangle_files = ['subjectG0_MEG_equilateral_triplets','subjectP9_MEG_equilateral_triplets'] #file-name without '.pkl' affix
  49. coordinate_files = ['subjectG0_MEG_coordinatesAP-IS-LR','subjectP9_MEG_coordinatesAP-IS-LR'] #file-name without '.csv' affix
  50. TS_files = ['subjectG0_MEG_timeseries','subjectP9_MEG_timeseries'] #file-name without '.npy' affix
  51. tdel_values = [3.2,3.2] # one over the sampling frequency, in ms.
  52. frequencies = [1.0, 1.148698354997035, 1.3195079107728942, 1.5157165665103982, 1.7411011265922482,
  53. 2.0, 2.29739670999407, 2.639015821545789, 3.0314331330207964, 3.4822022531844965,
  54. 4.0, 4.59479341998814, 5.278031643091579, 6.062866266041593, 6.964404506368995,
  55. 8.0, 9.18958683997628, 10.556063286183157, 12.125732532083186, 13.92880901273799,
  56. 16.0, 18.37917367995256, 21.112126572366314, 24.25146506416638, 27.857618025475986,
  57. 32.0] #temporal frequencies to analyze
  58. min_of_range = 1.0/20.0 #m, for plotting
  59. print(subjects)
  62. # %%
  63. #more parameters
  64. aggregate_subjects = {} #for storing the spectra and other used parametres for later analysis
  65. verbose = True #used as a global variable, set to False once everything is working ok
  67. distance_units = 1000.0 #conversion of given sensor units of coordinates to metres; here mm->m
  68. minf = min(frequencies) #minimum frequency, used to window the Morlet wavelet inputs
  69. nFrequencies = len(frequencies) #number of frequencies
  71. include_normed_power = False #if False will compute the SF spectrum of the phase only
  72. forward_waves = True #separate standing wave components into pure traveling wave components
  73. nBases = 14 #number of bases to used in SVD
  74. N_cycles = 2 #number of cycles is the Morlet wavelet to estimate time-series of phase
  76. nonskinny_threshold = np.pi/4.0 #45 degrees; used for making triangular regions between contacts; should be as close as 60 degrees as possible without reducing the number of triangles too severely
  77. smallest_triangle_size = 0.01 #units of linear m; used for binning the triangles
  78. largest_triangle_size = 0.32 #units of linear m; used for binning the triangles
  80. N_histogram_bins = 32 #for plotting
  81. if verbose==False: warnings.filterwarnings("ignore")
  84. # %%
  85. #main functions
  86. def complex_spatial_SVD(phi_cts,nBases=3,bases_only=False,use_bases=None,unit_phase_out=True): #singular value decomposition of the complex-valued phase data
  87. print('Doing spatial SVD')
  88. phi_Cs = phi_cts.reshape(-1,phi_cts.shape[-1])
  89. if use_bases is None: #generate new bases from phi
  90. u_sB,s_B,vh_BC= np.linalg.svd(phi_Cs.T,full_matrices=False) #b<B
  91. s_b = s_B[:nBases] #diagonals of sigma
  92. print('Percentage root variance',100.0*s_b/s_B.sum())
  93. bases_sb = u_sB[:,:nBases] #left singular vectors
  94. print('Percentage variance',100.0*(s_b)**2/((s_B)**2).sum())
  95. bases_sb = u_sB[:,:nBases] #left singular vectors
  96. else: bases_sb = use_bases #otherwise, use_bases is a previously generated basis set
  97. if bases_only==False: #calculate fits, right singular vectors, and reduced rank model
  98. if use_bases is None: betas_bC = (np.diag(s_b)@vh_BC[:nBases,:]) #weighted right singular vectors
  99. else: betas_bC = (phi_Cs@bases_sb.conj()).T
  100. betas_ctb = betas_bC.T.reshape(phi_cts.shape[0],-1,bases_sb.shape[1])
  101. if unit_phase_out:
  102. model_Cs = np.exp(1j*np.angle(bases_sb@betas_bC).T) #reduced rank model with unit length complex phases
  103. model_cts = model_Cs.reshape(phi_cts.shape)
  104. fit_C = (phi_Cs/model_Cs).mean(-1).real #for good fits, phi/model will point along the real axis
  105. fit_ct = fit_C.reshape(phi_cts.shape[0],phi_cts.shape[1])
  106. if use_bases is None: return bases_sb,s_b,fit_ct,betas_ctb,model_cts #number of outputs change depending on input flags
  107. else: return bases_sb,fit_ct,betas_ctb,model_cts #number of outputs change depending on input flags
  108. else:
  109. model_Cs = (bases_sb@betas_bC).T #reduced rank model
  110. model_cts = model_Cs.reshape(phi_cts.shape)
  111. if use_bases is None: return bases_sb,s_b,betas_ctb,model_cts #number of outputs change depending on input flags
  112. else: return bases_sb,betas_ctb,model_cts #number of outputs change depending on input flags
  113. else: return bases_sb #number of outputs change depending on input flags
  115. def make_SF_power(dmetres_SS,binned_triplets_av,normalization_a,sigmas_b,Bases_sb):
  116. #estimate wavelength and power for empirical data
  117. bases_avb = Bases_sb[binned_triplets_av,:]
  118. if verbose: print('bases_avb shape',bases_avb.shape)
  119. r_avb = np.absolute(bases_avb) #magnitude of contribution by each vertex (over triangles, bases)
  120. if verbose: print('r_avb min',r_avb.min(),'r_avb max',r_avb.max())
  121. estpower_ba = r_avb.mean(1).T * normalization_a[np.newaxis,:] * sigmas_b[:,np.newaxis] #weight variance explained by average vertex contribution for each triangle and normalize for number of triangles at this size
  122. if verbose: print('estpower_ba shape',estpower_ba.shape)
  123. dphidx_ba,dphi_ba = dphi_dx_from_triangles(bases_avb.transpose(2,0,1),dmetres_SS,binned_triplets_av)
  124. if verbose: print('dphidx_ba shape',dphidx_ba.shape)
  125. distance_per_radian_ba = 1.0/dphidx_ba
  126. metres_per_cycle_ba = 2.0*np.pi*distance_per_radian_ba #dmetres_SS already in metres
  127. if verbose: print('metres_per_cycle_ba shape',metres_per_cycle_ba.shape)
  128. return estpower_ba,metres_per_cycle_ba,dphidx_ba,dphi_ba
  130. def dphi_dx_from_triangles(phi_cav,dmetres_SS,triplets_av): #compute the rate of change of phase across triangles
  131. a = dmetres_SS[triplets_av[:,0],triplets_av[:,1]]
  132. b = dmetres_SS[triplets_av[:,1],triplets_av[:,2]]
  133. c = dmetres_SS[triplets_av[:,2],triplets_av[:,0]]
  134. ABC_v_xyz = coordinates_of_triangle_given_lengths(a,b,c) #nvd
  135. ABC_vd = ABC_v_xyz[:,:,:2] #don't need zero Z coordinate
  136. _ABC_vd = ABC_vd / np.linalg.norm(ABC_vd,axis=-2)[:,np.newaxis,:] #normalized area to one
  137. relative_phase = np.angle(phi_cav[:,:,0:1]/phi_cav[:,:,0:3]) #relative to first index; units of radians
  138. dphi_dx = []
  139. dphi = []
  140. for c in range(phi_cav.shape[0]):
  141. dphi_dx_c = []
  142. dphi_c = []
  143. for n in range(phi_cav.shape[1]):
  144. y = relative_phase[c,n,:,np.newaxis] #vo
  145. X = ABC_vd[n,:,:] #vd
  146. X1 = np.concatenate((X,np.ones((X.shape[0],1),float)),axis=1)
  147. b = np.linalg.inv(X1.T@X1)@(X1.T@y) #dv,vd->dd; dv,vo->d0: dd,d0->d0
  148. dphi_dx_c += [np.linalg.norm(b[:,0],axis=0)]
  149. _X = _ABC_vd[n,:,:] #vd
  150. _X1 = np.concatenate((_X,np.ones((_X.shape[0],1),float)),axis=1)
  151. _b = np.linalg.inv(_X1.T@_X1)@(_X1.T@y) #dv,vd->dd; dv,vo->d0: dd,d0->d0
  152. dphi_c += [np.linalg.norm(_b[:,0],axis=0)]
  153. dphi_dx += [np.array(dphi_dx_c)]
  154. dphi += [np.array(dphi_c)]
  155. dphi_dx = np.array(dphi_dx)
  156. dphi = np.array(dphi)
  157. return dphi_dx,dphi
  159. def wavelet(data_cts,N_cycles,frequency,tdel,minf,use_cuda=False,PHname=None): #compute the phase and power using Morlet wavelets
  160. #read in the files if they arlready exist
  161. if PHname is not None:
  162. LPname = PHname.replace('PH_','LP_')
  163. print('Looking for',PHname)
  164. print('Looking for',LPname)
  165. if PHname is not None and exists(PHname+'.npy') and exists(LPname+'.npy'):
  166. print('Reading in phase file')
  167. phase_CTS = np.load(PHname+'.npy')
  168. if verbose: print('phase_CTS min max',phase_CTS.min(),phase_CTS.max(),'phase_CTS shape',phase_CTS.shape)
  169. print('Reading in power file')
  170. logpower_CTS = np.load(LPname+'.npy')
  171. if verbose: print('logpower_CTS min max',logpower_CTS.min(),logpower_CTS.max(),'logpower_CTS shape',logpower_CTS.shape)
  172. return np.exp(1j*phase_CTS),np.power(10,logpower_CTS) #convert the phase to complex-valued unit vectors (phi), convert the logpower to power
  173. else:
  174. #data_cts is the time series data in shape cases (trials, epochs), times (samples), sensors (contacts, electrodes)
  175. #N_cycles is the number of cycles in the Morlet wavelet; typically 2 in Alexander et al.
  176. #frequency is the centre frequency, typically chosen from oversampled, logarithmically spaced frequncies
  177. #tdel is the time interval between samples e.g. 1ms
  178. #minf is the minimum frequency used in the analysis, need to ensure the final phase/power estimates all have the same number of samples
  179. print('Estimating short time-series Fourier components for temporal frequency %.2fHz'%frequency)
  180. nSamples = data_cts.shape[1]
  181. nSensors = data_cts.shape[2]
  182. N_cycles_of_samples = int(N_cycles*1000.0/(frequency*tdel))
  183. #the extra samples of data required at both start and beginning of estimated phase/power values
  184. PH_unpad_at_minf = int(N_cycles*500.0/(minf*tdel))
  185. #the number of samples of phase/power that will be generated
  186. PHSamples = nSamples - 2*PH_unpad_at_minf
  187. #half the wavelet cycles at lowest frequency minus half the wavelet cycles at this frequency
  188. PH_pad_at_f = PH_unpad_at_minf - int(N_cycles*500.0/(frequency*tdel))
  189. W = gausian_window(N_cycles_of_samples) #size of Morlet window
  190. f_indexes = []
  191. #how many sets of indexes to make? one for each sample in phase
  192. for t in range(PHSamples):
  193. #make a list of all data samples that are used in calculating each phase estimate
  194. padded_range = range(t+PH_pad_at_f,t+PH_pad_at_f+N_cycles_of_samples)
  195. f_indexes += padded_range
  196. if use_cuda: #works with cupy if available https://cupy.dev
  197. f_indexes = cp.asarray(f_indexes)
  198. power_cTs = cp.zeros((data_cts.shape[0],PHSamples,nSensors),float)
  199. phi_cTs = cp.zeros((data_cts.shape[0],PHSamples,nSensors),complex)
  200. W = cp.asarray(W)
  201. for c,data_ts in enumerate(data_cts):
  202. print('#',end='')
  203. data_ts = cp.asarray(data_ts)
  204. for s in range(nSensors):
  205. data_tw = data_ts[f_indexes,s].reshape(PHSamples,N_cycles_of_samples) #tw
  206. data_tw = data_tw - data_tw.mean(1)[:,cp.newaxis] #mean of time-series to zero
  207. t_f = 2.0*cp.pi*N_cycles*cp.arange(N_cycles_of_samples)/N_cycles_of_samples
  208. t_phi = cp.exp(1j*(t_f)) * W
  209. fourier_components = (data_tw * t_phi[cp.newaxis,:]).sum(1)
  210. phi_cTs[c,:,s] = cp.exp(1j*cp.angle(fourier_components)) #over w
  211. power_cTs[c,:,s] = cp.absolute(fourier_components)
  212. print()
  213. return cp.asnumpy(phi_cTs),cp.asnumpy(power_cTs)
  214. else:
  215. f_indexes = np.asarray(f_indexes)
  216. power_cTs = np.zeros((data_cts.shape[0],PHSamples,nSensors),float)
  217. phi_cTs = np.zeros((data_cts.shape[0],PHSamples,nSensors),complex)
  218. for c,data_ts in enumerate(data_cts):
  219. print('#',end='')
  220. for s in range(nSensors):
  221. data_tw = data_ts[f_indexes,s].reshape(PHSamples,N_cycles_of_samples) #tw
  222. data_tw = data_tw - data_tw.mean(1)[:,np.newaxis] #mean of time-series to zero
  223. t_f = 2.0*np.pi*N_cycles*np.arange(N_cycles_of_samples)/N_cycles_of_samples
  224. t_phi = np.exp(1j*(t_f)) * W
  225. fourier_components = (data_tw * t_phi[np.newaxis,:]).sum(1)
  226. phi_cTs[c,:,s] = np.exp(1j*np.angle(fourier_components)) #over w
  227. power_cTs[c,:,s] = np.absolute(fourier_components)
  228. print()
  229. if PHname is not None: #i.e. looked for save files but didn't find them
  230. print('Writing to',PHname+'.npy',end = ' ')
  231. np.save(PHname+'.npy',np.angle(phi_cTs)) #saved as the real angle
  232. print('Writing to',LPname+'.npy',end = ' ')
  233. np.save(LPname+'.npy',np.log10(power_cTs)) #save the log power
  234. return phi_cTs,power_cTs
  236. def gaus(x,a,x0,sigma): #define a Gaussian
  237. return a*np.exp(-(x-x0)**2/(2*sigma**2))
  239. def gausian_window(x): #make a Gaussian window
  240. if type(x) is int:
  241. n = x
  242. x = np.arange(n)
  243. else:
  244. n = len(x)
  245. cosine_window = 0.5*(1.0 - np.cos(2*np.pi*x/(n-1))) #approximate the window as a cosine first
  246. mean = sum(x*cosine_window)/n
  247. sigma = np.sqrt(sum(cosine_window*(x-mean)**2)/n)
  248. popt,pcov = curve_fit(gaus,x,cosine_window,p0=[1,mean,sigma]) #turn it into a Gaussian window
  249. G = gaus(x,*popt)
  250. return G
  252. # %%
  253. #helper functions
  254. def save_pkl(output_path,obj,name): #the format used for saving dictionaries, lists
  255. print('Saving',output_path+ name + '.pkl')
  256. with open(output_path+ name + '.pkl', 'wb') as f:
  257. dump(obj, f, HIGHEST_PROTOCOL)
  259. def load_pkl(output_path,name): #open previously saved dictionaries, lists
  260. print('Loading',output_path + name + '.pkl')
  261. with open(output_path + name + '.pkl', 'rb') as f:
  262. return load(f)
  264. def get_coordinates_and_distances(file_path,coordinate_file,distance_file,distance_units): # convenience function for grabbing sensor coordinates and cross-sensor distances
  265. print('Opening coordinate file',coordinate_file)
  266. coordinates_sd = np.loadtxt(file_path+coordinate_file+'.csv',skiprows=1,delimiter=',').T #native coordinate system units; A-P,I-S,L-R: anterior-posterior, inferior-superior, left-right
  267. coordinates_metres = coordinates_sd/distance_units
  268. dmetres_ss = load_pkl(file_path,distance_file)/distance_units #mm->m in this data
  269. print('dmetres_ss min',dmetres_ss[dmetres_ss>0.0].min(),'dmetres_ss max',dmetres_ss[dmetres_ss<0.5].max())
  270. if 'depth' in coordinate_file:
  271. #initialize distances between contacts
  272. dmetres_SS = np.linalg.norm(coordinates_metres[:,np.newaxis,:]-coordinates_metres[np.newaxis,:,:],axis=2)
  273. if verbose: print('dmetres_SS min',dmetres_SS.min(),'dmetres_SS max',dmetres_SS.max())
  274. print('dmetres_SS[dmetres_SS>0] min',dmetres_SS[dmetres_SS>0].min(),'dmetres_SS[dmetres_SS<0.5] max',dmetres_SS[dmetres_SS<0.5].max())
  275. #grab geodesic distances
  276. for s in range(dmetres_SS.shape[0]):
  277. for ss in range(dmetres_SS.shape[1]):
  278. dmetres_SS[s,ss] = max(dmetres_SS[s,ss],dmetres_ss[s,ss]) #take the largest of the two because geodesic can give zero for neighbouring contacts
  279. elif 'MEG' in coordinate_file:
  280. dmetres_SS = dmetres_ss
  281. return coordinates_metres,dmetres_SS
  283. def angles_of_triangle_from_lengths(dmetres_SS,triplets_av): #used to construct approximate equilateral triangles from inter-contact distances
  284. """a, b and c are lengths of the sides of a triangle"""
  285. a = dmetres_SS[triplets_av[:,0],triplets_av[:,1]]
  286. b = dmetres_SS[triplets_av[:,1],triplets_av[:,2]]
  287. c = dmetres_SS[triplets_av[:,2],triplets_av[:,0]]
  288. # Square of lengths a2, b2, c2
  289. a2 = a**2
  290. b2 = b**2
  291. c2 = c**2
  292. # From Cosine law
  293. alpha = np.arccos((b2 + c2 - a2) / (2 * b * c))
  294. beta = np.arccos((a2 + c2 - b2) / (2 * a * c))
  295. gamma = np.arccos((a2 + b2 - c2) / (2 * a * b))
  296. return alpha,beta,gamma
  299. def triangle_normal(triangles):
  300. # The cross product of two sides is a normal vector
  301. return np.cross(triangles[:,1] - triangles[:,0],
  302. triangles[:,2] - triangles[:,0], axis=1)
  304. def triangle_area(triangles):
  305. # The norm of the cross product of two sides is twice the area
  306. return np.linalg.norm(triangle_normal(triangles), axis=1) / 2
  308. def collate_nonskinny_triangles(dmetres_SS,nonskinny_threshold=np.pi/4.0): #make a list of triplets, with each triplet specifying an approximately equilateral triangle, where nonskinny_threshold is the smallest allowed angle
  309. nSensors = dmetres_SS.shape[0]
  310. if verbose: print('nSensors',nSensors)
  311. if verbose: print('dmetres_SS shape',dmetres_SS.shape)
  312. if verbose: print('dmetres_SS')
  313. if verbose: print(dmetres_SS)
  314. print('Making nonskinny triangles of contacts')
  315. unique_triplets = []
  316. for s in range(nSensors):
  317. print('.',end='')
  318. for ss in range(nSensors):
  319. for sss in range(nSensors):
  320. if s!=ss and ss!=sss and sss!=s and \
  321. dmetres_SS[s,ss]>0.0 and dmetres_SS[ss,sss]>0.0 and dmetres_SS[sss,s]>0.0 and \
  322. dmetres_SS[s,ss]<0.5 and dmetres_SS[ss,sss]<0.5 and dmetres_SS[sss,s]<0.5: #0.5m is the flag for cross-hemisphere distances in sEEG distance calculation:
  323. u = set(np.array([s,ss,sss]))
  324. if u not in unique_triplets: unique_triplets += [u]
  325. print()
  326. for u,un in enumerate(unique_triplets): unique_triplets[u] = np.sort(np.array(list(un))) #set becomes array
  327. unique_triplets = np.array(unique_triplets)
  328. print(unique_triplets.shape)
  329. alpha_a,beta_a,gamma_a = angles_of_triangle_from_lengths(dmetres_SS,unique_triplets) #ss,av
  330. angles_va = np.array([alpha_a,beta_a,gamma_a])
  331. angles_min_a = angles_va.min(0)
  332. nonskinnyity_met = angles_min_a>nonskinny_threshold
  333. nonskinny_triplets = np.array(unique_triplets)[nonskinnyity_met]
  334. if verbose:
  335. print('unique_triplets shape',unique_triplets.shape)
  336. print('nonskinnyity_met shape',nonskinnyity_met.shape)
  337. print('nonskinny_triplets shape',nonskinny_triplets.shape)
  338. print('nonskinny_triplets')
  339. print(nonskinny_triplets)
  340. return nonskinny_triplets #the indexes of contacts in each triangle
  342. def bin_triplets_linearSF(dmetres_SS,nonskinny_triplets,return_binwise=False,smallest_bin_size=0.01,largest_bin_size=0.32): #can return triangles in bins, or as a flat list ordered by triangle size
  343. print('Binning triangles')
  344. try:
  345. #calculate areas of each triangle from geodesic distances
  346. nominal_min_SF = 500.0 # cycles per metre
  347. n_bins = int(largest_bin_size/smallest_bin_size) #this is silly way to get n=32
  348. #get the triangle edges from the geodesic distances and triplets
  349. ta = nonskinny_triplets
  350. dm = dmetres_SS
  351. dm_ta01 = dm[ta[:,0],ta[:,1]]
  352. dm_ta12 = dm[ta[:,1],ta[:,2]]
  353. dm_ta20 = dm[ta[:,2],ta[:,0]]
  354. nominal_triangle_coordinates_te = coordinates_of_triangle_given_lengths(dm_ta01,dm_ta12,dm_ta20) #te: triangle x edge
  355. triplet_areas = triangle_area(nominal_triangle_coordinates_te)
  356. root_areas = np.sqrt(triplet_areas) #linear size of each triangle
  357. order_root_areas = root_areas.argsort(0)
  358. ordered_root_areas = root_areas[order_root_areas]
  359. ordered_nonskinny_triplets = nonskinny_triplets[order_root_areas]
  360. #fill the bins
  361. nonskinny_triplets_linearSF_bins = [] #triangles assigned to each bin
  362. nonskinny_triplets_a = [] #flat list of triangles assigned, ordered by triangle linear area
  363. root_areas_linearSF_bins = [] #areas of triangles assigned to each bin
  364. root_areas_a = [] #flat list of areas of triangles assigned to each bin
  365. n_per_linearSF_bins = [] #number of triangles in each bin
  366. if ordered_root_areas.shape[0]==0: return None
  367. linearSF_bin_edges = [ordered_root_areas[0]] #first two edges, because n_bins+1 edges
  368. linearSF_ranges_bins = [] #min and max of each bin
  369. old_size_cutoff = 1.0/nominal_min_SF
  370. for SF_cutoff in np.append(np.linspace(1.0/largest_bin_size,1.0/smallest_bin_size,n_bins),nominal_min_SF)[:-1][::-1]: #last value is taken first, but is already used in initialization
  371. size_cutoff = 1.0/SF_cutoff
  372. size_cutoff_boolean = ((ordered_root_areas<size_cutoff) * (ordered_root_areas>=old_size_cutoff)).astype(bool)
  373. nThese = size_cutoff_boolean.sum(0)
  374. these_nonskinny_triplets = ordered_nonskinny_triplets[size_cutoff_boolean,:]
  375. these_root_areas = ordered_root_areas[size_cutoff_boolean]
  376. if verbose:
  377. print(old_size_cutoff,size_cutoff,nThese)
  378. if nThese>0: print(these_root_areas.min(0),these_root_areas.max(0))
  379. else: print('-','-')
  380. if nThese>0:
  381. nonskinny_triplets_linearSF_bins += [list(these_nonskinny_triplets)]
  382. nonskinny_triplets_a += list(these_nonskinny_triplets)
  383. root_areas_linearSF_bins += [list(these_root_areas)]
  384. root_areas_a += list(these_root_areas)
  385. linearSF_bin_edges += [size_cutoff]
  386. linearSF_ranges_bins += [np.array([these_root_areas.min(0),these_root_areas.max(0)])] #zeroth is zero, 1st is smallest_bin_size, nth is largest_bin_size
  387. n_per_linearSF_bins += [these_nonskinny_triplets.shape[0]]
  388. else:
  389. nonskinny_triplets_linearSF_bins += [[np.array([-1,-1,-1])]]
  390. root_areas_linearSF_bins += [[-1]]
  391. linearSF_bin_edges += [-1]
  392. linearSF_ranges_bins += [[-1,-1]]
  393. n_per_linearSF_bins += [0]
  394. old_size_cutoff = 1*size_cutoff
  395. if len(root_areas_a)==0: return None
  396. root_areas_a = np.array(root_areas_a)
  397. print('root_areas shape',root_areas.shape,'ordered_root_areas shape',ordered_root_areas.shape,'root_areas_a shape',root_areas_a.shape)
  398. n_per_linearSF_bins = np.array(n_per_linearSF_bins)
  399. #some book-keeping
  400. nonzero_root_areas_linearSF_bins = []
  401. nonzero_n_per_linearSF_bins = []
  402. nonzero_nonskinny_triplets_linearSF_bins = []
  403. nonzero_linearSF_bin_edges = []
  404. nonzero_linearSF_ranges_bins = [linearSF_ranges_bins[0]]
  405. for B,n_bin in enumerate(n_per_linearSF_bins):
  406. if n_bin>0:
  407. nonzero_root_areas_linearSF_bins += [root_areas_linearSF_bins[B]]
  408. nonzero_n_per_linearSF_bins += [n_per_linearSF_bins[B]]
  409. nonzero_nonskinny_triplets_linearSF_bins += [nonskinny_triplets_linearSF_bins[B]]
  410. if len(nonzero_linearSF_bin_edges)==0: nonzero_linearSF_bin_edges += [linearSF_bin_edges[B]]
  411. nonzero_linearSF_bin_edges += [linearSF_bin_edges[B+1]] #n_bins + 1 values
  412. nonzero_linearSF_ranges_bins += [linearSF_ranges_bins[B]]
  413. nonzero_linearSF_bin_edges[-1] = ordered_root_areas[-1] #correct the last edges
  414. nonzero_linearSF_ranges_bins[-1][1] = ordered_root_areas[-1] #correct the second term of the last nonzero bin range
  415. nonzero_linearSF_bin_edges = np.array(nonzero_linearSF_bin_edges)
  416. nonzero_linearSF_ranges_bins = np.array(nonzero_linearSF_ranges_bins)
  417. nonzero_n_per_linearSF_bins = np.array(nonzero_n_per_linearSF_bins)
  418. range_a = [ordered_root_areas[0],ordered_root_areas[-1]]
  419. print('nonzero_linearSF_ranges_bins')
  420. print(nonzero_linearSF_ranges_bins)
  421. max_of_linearSF_bins = max(nonzero_n_per_linearSF_bins)
  422. normalization_linearSF_bins = -np.ones(n_per_linearSF_bins.shape,float)
  423. normalization_linearSF_a = [] #'linearSF' because its based on assumption of linearly increasing SF
  424. print('n_per_linearSF_bins sum',n_per_linearSF_bins.sum(),'nonzero_n_per_linearSF_bins sum',nonzero_n_per_linearSF_bins.sum())
  425. for B,bin_n in enumerate(n_per_linearSF_bins):
  426. if bin_n>0:
  427. normalization_linearSF_bins[B] = max_of_linearSF_bins/bin_n #will eventually multiply by normalization weight
  428. normalization_linearSF_a += bin_n*[normalization_linearSF_bins[B]] #grab copies equal to the number of triangles in bin
  429. normalization_linearSF_a = np.array(normalization_linearSF_a)
  430. print('normalization_linearSF_a shape',normalization_linearSF_a.shape)
  431. if verbose:
  432. print('linearSF_bin_edges')
  433. print(linearSF_bin_edges)
  434. print('linearSF_ranges_bins')
  435. print(linearSF_ranges_bins)
  436. print('n_per_linearSF_bins')
  437. print(n_per_linearSF_bins)
  438. print('normalization_linearSF_a')
  439. print(normalization_linearSF_a)
  440. print('nonzero_linearSF_bin_edges')
  441. print(nonzero_linearSF_bin_edges)
  442. print('nonzero_n_per_linearSF_bins')
  443. print(nonzero_n_per_linearSF_bins)
  444. plot_bin_counts_and_normalization(nonzero_linearSF_bin_edges,nonzero_n_per_linearSF_bins,root_areas_a,normalization_linearSF_a,SFscale=True,cut_final_bin=True)
  445. if return_binwise: return n_per_linearSF_bins,nonskinny_triplets_linearSF_bins,normalization_linearSF_bins,\
  446. linearSF_ranges_bins,root_areas_linearSF_bins,linearSF_bin_edges #first array is used to avoid empty bins
  447. else: return n_per_linearSF_bins,np.array(nonskinny_triplets_a),normalization_linearSF_a,range_a,root_areas_a,linearSF_bin_edges
  448. except:
  449. print('No binned triplets returned')
  450. return None
  452. def coordinates_of_triangle_given_lengths(a,b,c): #abstract triangle coordinates from geodesic distances
  453. """a, b and c are lengths of the sides of a triangle"""
  454. A = np.array([[0]*len(a),[0]*len(a),[0]*len(a)]) # coordinates of vertex A
  455. B = np.array([c,[0]*len(a),[0]*len(a)]) # coordinates of vertex B
  456. #if verbose: print(A.shape,B.shape)
  457. C_x = b * (b**2 + c**2 - a**2) / (2 * b * c)
  458. C_y = np.sqrt(b**2 - C_x**2) # square root
  459. C = np.array([C_x,C_y,[0]*len(a)]) # coordinates of vertex C
  460. return np.concatenate([[A],[B],[C]],axis=0).transpose(2,0,1) #tvd: triangle, vertex, cartesian dimension
  462. def decompose_rphi_into_forward_waves(rphi_cts,tdel,frequency): #removes the negative going component of the wave and estimates normalized velocity
  463. #tdel is 1/sampling_frequency in ms
  464. #frequency in Hertz
  465. #rphi_cts is complex valued phase, possibly with magnitude weighting
  466. print('Decompose spatial vectors of phase into pure traveling wave components')
  467. one_cycle_of_samples = int(1000.0/(tdel*frequency)) #only requires four samples over the cycle, so this int could be replace with 4
  468. rphi_Cs = rphi_cts.reshape(-1,rphi_cts.shape[2])
  469. vel_C = [] #normalized velocity, 0 is pure SW, 1 is pure TW
  470. TWf_C = [] #forward spatial component
  471. for C in range(rphi_Cs.shape[0]):
  472. rphi_Ts = np.exp(1j*np.linspace(-np.pi,np.pi,one_cycle_of_samples,endpoint=False))[:,np.newaxis] * rphi_Cs[C,:][np.newaxis,:]
  473. u_r,s_r,vt_r = np.linalg.svd(rphi_Ts.real)
  474. TWf_s = vt_r[0,:] + 1j*vt_r[1,:] #throw away the time dimension we added, only need the spatial dimension
  475. TWf_C += [TWf_s[np.newaxis]]
  476. vel_C += [np.array([s_r[1]/s_r[0]])[np.newaxis]]
  477. vel_ct = np.concatenate(vel_C,axis=0).reshape(rphi_cts.shape[0],rphi_cts.shape[1])
  478. TWf_cts = np.concatenate(TWf_C,axis=0).reshape(rphi_cts.shape[0],rphi_cts.shape[1],rphi_cts.shape[2])
  479. return vel_ct,TWf_cts
  481. # %%
  482. #plotting functions
  483. def sf_histogram(imaging_path,plot_name,output_file,list_of_estpower,list_of_wavelength,list_of_max_distances, #plot the histogram of SF
  484. max_of_range,min_of_range,N_histogram_bins,
  485. stacked=False,colourmap='RdYlBu_r',linewidth=1,ticklabelsize=10.0,axislabelsize=12.0,listwise_weights=None):
  486. print('Plotting',plot_name)
  487. if verbose: print('max_of_range',max_of_range,'min_of_range',min_of_range)
  488. if verbose: print('1.0/max_of_range',1.0/max_of_range,'1.0/min_of_range',1.0/min_of_range)
  489. cmap = cm.get_cmap(colourmap)
  490. colours = []
  491. for l in range(len(list_of_estpower)):
  492. colours += [cmap(l/len(list_of_estpower))] #avoid division by zero
  493. #colours = ['red','green','blue','magenta','cyan']
  494. fig1 = figure(figsize=(5,5),dpi=90)
  495. axs1 = []
  496. axs1 += [fig1.add_subplot(1,1,1)]
  497. axs1[-1].set_title(plot_name.replace('_',' '))
  498. axs1[-1].tick_params(axis='x',which = 'both',labelsize=ticklabelsize)
  499. axs1[-1].tick_params(axis='y',labelsize=ticklabelsize)
  500. axs1[-1].set_xscale('log')
  501. if 'sEEG' in plot_name:
  502. locs = np.append(np.arange(1.0,10.0,1.0),np.arange(10.0,100.0,10.0))
  503. locs = np.delete(locs,np.where(locs==9.0))
  504. locs = np.delete(locs,np.where(locs==40.0))
  505. locs = np.delete(locs,np.where(locs==60.0))
  506. locs = np.delete(locs,np.where(locs==70.0))
  507. locs = np.delete(locs,np.where(locs==80.0))
  508. locs = np.delete(locs,np.where(locs==90.0))
  509. elif 'MEG' in plot_name:
  510. locs = np.append(np.arange(1.0,10.0,1.0),np.arange(10.0,50.0,10.0))
  511. locs = np.delete(locs,np.where(locs==9.0))
  512. locs = np.delete(locs,np.where(locs==40.0))
  513. axs1[-1].xaxis.set_minor_locator(ticker.FixedLocator(locs))
  514. axs1[-1].xaxis.set_major_locator(ticker.NullLocator())
  515. axs1[-1].xaxis.set_minor_formatter(ticker.ScalarFormatter())
  516. axs1[-1].set_xlabel('Spatial frequency (cycles/metre)',fontsize=axislabelsize)
  517. axs1[-1].set_ylabel('Estimated power (A.U.)',fontsize=axislabelsize)
  518. if listwise_weights is None: listwise_weights = np.ones((len(list_of_estpower)),float)
  519. logscale_bins = np.power(10,np.linspace(np.log10(1.0/max_of_range),np.log10(1.0/min_of_range),N_histogram_bins))
  520. list_of_hist_power = []
  521. list_of_hist_sf = []
  522. if stacked:
  523. cumulative_hist = np.zeros((N_histogram_bins-1),float)
  524. for l,estpower,wavelength,case_weight in zip(range(len(list_of_estpower)),list_of_estpower,list_of_wavelength,listwise_weights):
  525. wavelength = np.array(wavelength)
  526. estpower = np.array(estpower)
  527. sf = 1.0/wavelength
  528. power_hist = np.histogram(sf,bins=logscale_bins,weights=estpower)
  529. w_power_hist = case_weight * power_hist[0] #case_weight is a post-hoc way to add weighting to each curve e.g. to mimic sigma
  530. cumulative_hist += w_power_hist
  531. axs1[-1].plot((power_hist[1][:-1]+power_hist[1][1:])/2.0,cumulative_hist,c=colours[l],linewidth=linewidth)
  532. for md in list_of_max_distances:
  533. axs1[-1].plot(1.0/md,0.0,'.',color='k')
  534. list_of_hist_power += [w_power_hist]
  535. list_of_hist_sf += [(power_hist[1][:-1]+power_hist[1][1:])/2.0]
  536. else:
  537. for l,estpower,wavelength,case_weight in zip(range(len(list_of_estpower)),list_of_estpower,list_of_wavelength,listwise_weights):
  538. wavelength = np.array(wavelength)
  539. estpower = np.array(estpower)
  540. sf = 1.0/wavelength
  541. power_hist = np.histogram(sf,bins=logscale_bins,weights=estpower)
  542. w_power_hist = case_weight * power_hist[0]
  543. axs1[-1].plot((power_hist[1][:-1]+power_hist[1][1:])/2.0,w_power_hist,c=colours[l],linewidth=linewidth)
  544. for md in list_of_max_distances:
  545. axs1[-1].plot(1.0/md,0.0,'.',color='k')
  546. list_of_hist_power += [w_power_hist]
  547. list_of_hist_sf += [(power_hist[1][:-1]+power_hist[1][1:])/2.0]
  548. display(fig1)
  549. fig1.savefig((imaging_path+plot_name+'_'+output_file).replace(' ','_').replace('c/m','cperm')+'.jpg',facecolor='white',transparent=False,dpi=600)
  550. fig1.savefig((imaging_path+plot_name+'_'+output_file).replace(' ','_').replace('c/m','cperm')+'.svg',facecolor='white',transparent=False,dpi=600)
  551. fig1.clf()
  552. close(fig1)
  553. del fig1
  554. return list_of_hist_power,list_of_hist_sf
  556. def plot_bin_counts_and_normalization(bin_edges,bin_counts,used_root_areas,normalization_a,SFscale=True,cut_final_bin=False): #for checking the triangle construction, binning, normalization
  557. print('Plotting bin counts')
  558. fig = figure(figsize=(5,2.5))
  559. ax = fig.add_subplot(1,1,1)
  560. if SFscale:
  561. ax.plot(1.0/bin_edges[1:],bin_counts,'4',color='g')
  562. ax.plot(1.0/bin_edges[:-1],bin_counts,'3',color='r')
  563. ax.set_xlabel('Bin edges (1/m)')
  564. ax.set_ylabel('Triangle count')
  565. else:
  566. ax.plot(bin_edges[1:],bin_counts,'3',color='g')
  567. ax.plot(bin_edges[:-1],bin_counts,'4',color='r')
  568. ax.set_xlabel('Bin edges (m)')
  569. ax.set_ylabel('Triangle count')
  570. xlims = ax.get_xlim()
  571. ax.set_xlim(xlims[0],1.2/bin_edges[1])
  572. display(fig)
  573. fig.clf()
  574. close(fig)
  575. del fig
  576. print('Plotting binwise normalization')
  577. fig = figure(figsize=(5,2.5))
  578. ax = fig.add_subplot(1,1,1)
  579. if SFscale:
  580. ax.plot(1.0/used_root_areas,1.0/normalization_a,'.',color='k',markersize=0.2) #normalization_a is divisor in SF power calculation
  581. ax.set_xlabel('Triangle SF')
  582. ax.set_ylabel('1.0 / normalization factor')
  583. else:
  584. ax.plot(used_root_areas,1.0/normalization_a,'.',color='k',markersize=0.2) #normalization_a is divisor in SF power calculation
  585. ax.set_xlabel('Triangle size')
  586. ax.set_ylabel('1.0 / normalization factor')
  587. xlims = ax.get_xlim()
  588. ax.set_xlim(xlims[0],1.2/bin_edges[1])
  589. display(fig)
  590. fig.clf()
  591. close(fig)
  592. del fig
  594. def phase_plot(bases_sb,frequency,contact_xyz,vector_name='wave map',mag=False,figout=None): #assumes complex valued bases
  595. print('Plotting phases')
  596. nBases = bases_sb.shape[1]
  597. cols = 3
  598. rows = (nBases//3)+1
  599. fig = figure(figsize=(cols*8,rows*8),dpi=200)
  600. axs = []
  601. if mag: alpha = np.absolute(bases_sb)/np.absolute(bases_sb).max()
  602. else: alpha = np.ones(bases_sb.shape,float)
  603. for b in range(nBases):
  604. axs += [fig.add_subplot(rows,cols,b+1,projection='3d')]
  605. axs[b].view_init(elev=90,azim=-90)
  606. axs[b].set_title('Frequency %.1fHz %s %d'%(frequency,vector_name,b+1))
  607. for s,xyz in enumerate(contact_xyz):
  608. color = cm.hsv(np.angle(bases_sb[s,b]) / (2.0*np.pi) + 0.5) #put angles on range 0 to 1
  609. axs[b].plot(xyz[0],xyz[1],xyz[2],'.',markersize=20.0,markerfacecolor=color,markeredgecolor=None,alpha=alpha[s,b])
  610. axs[b].plot(xyz[0],xyz[1],xyz[2],'.',markersize=20.0,markeredgecolor=color,fillstyle='none',alpha=1.0)
  611. axs[b].set_aspect('equal')
  612. axs[b].set_xlabel('Anterior-Posterior')
  613. axs[b].set_ylabel('Inferior-Superior')
  614. fig.subplots_adjust(wspace=0.0, hspace=0.0)
  615. if figout is not None:
  616. fig.savefig(figout+'.png',facecolor='white',transparent=False)
  617. fig.savefig(figout+'.svg',facecolor='white',transparent=False)
  618. else:
  619. rect = fig.patch
  620. rect.set_facecolor('white')
  621. display(fig)
  622. fig.clf()
  623. close(fig)
  624. del fig
  626. def draw_triangles(contact_xyz,binned_triplets,bin_counts,file_path,title,root_areas): #for visualization of the equilateral triangles, draws in simple cartesian coordinates
  627. print('Plotting triangles over size bins')
  628. nBinPlots = bin_counts.shape[0]
  629. cols = 3
  630. rows = (nBinPlots//3)+1
  631. fig = figure(figsize=(cols*10,rows*10),dpi=150)
  632. axs = []
  633. cumulative_bin_count = 0
  634. for br in range(nBinPlots):
  635. axs += [fig.add_subplot(rows,cols,br+1,projection='3d')]
  636. axs[br].view_init(elev=90,azim=-90)
  637. axs[br].plot(contact_xyz[:,0],contact_xyz[:,1],contact_xyz[:,2],'.',color='k')
  638. if bin_counts[br]>0:
  639. triplets = binned_triplets[cumulative_bin_count:cumulative_bin_count+bin_counts[br]]
  640. mean_root_area = root_areas[cumulative_bin_count:cumulative_bin_count+bin_counts[br]].mean(0)
  641. axs[br].set_title('Size bin %d, mean size %.3f'%(br+1,mean_root_area))
  642. for t,trip in enumerate(triplets):
  643. if triplets.shape[0]>1: colour = cm.cool(t/(triplets.shape[0])) #0 to 1
  644. else: colour = cm.cool(0.0)
  645. T = np.concatenate((trip,trip[0:1]),axis=0) #close the triangle
  646. axs[br].plot(contact_xyz[T,0],contact_xyz[T,1],contact_xyz[T,2],color=colour,alpha=1.0/np.sqrt(bin_counts[br])) #density of colour on plot increases with line length
  647. if t==0: axs[br].plot(contact_xyz[T,0],contact_xyz[T,1],contact_xyz[T,2],color=cm.cool(0.0),alpha=1.0)
  648. elif t==triplets.shape[0]-1: axs[br].plot(contact_xyz[T,0],contact_xyz[T,1],contact_xyz[T,2],color=cm.cool(1.0),alpha=1.0)
  649. axs[br].set_aspect('equal')
  650. xticks = axs[br].get_xticks()
  651. axs[br].set_xticks(xticks[::2])
  652. zticks = axs[br].get_zticks()
  653. axs[br].set_zticks(zticks[::4])
  654. cumulative_bin_count += bin_counts[br]
  655. set_axes_equal(axs[br]) # aspect ratio is 1:1:1 in data space
  656. tight_layout()
  657. fig.savefig((file_path+title).replace('c/m','cperm')+'.png')
  658. fig.savefig((file_path+title).replace('c/m','cperm')+'.svg')
  659. display(fig)
  660. fig.clf()
  661. close(fig)
  662. del fig
  664. def set_axes_equal(ax): #to draw veridical relative distances in plots
  665. """
  666. Make axes of 3D plot have equal scale so that spheres appear as spheres,
  667. cubes as cubes, etc.
  669. Input
  670. ax: a matplotlib axis, e.g., as output from plt.gca().
  671. """
  672. x_limits = ax.get_xlim3d()
  673. y_limits = ax.get_ylim3d()
  674. z_limits = ax.get_zlim3d()
  675. x_range = abs(x_limits[1] - x_limits[0])
  676. x_middle = np.mean(x_limits)
  677. y_range = abs(y_limits[1] - y_limits[0])
  678. y_middle = np.mean(y_limits)
  679. z_range = abs(z_limits[1] - z_limits[0])
  680. z_middle = np.mean(z_limits)
  681. # The plot bounding box is a sphere in the sense of the infinity
  682. # norm, hence I call half the max range the plot radius.
  683. plot_radius = 0.5*max([x_range, y_range, z_range])
  684. ax.set_xlim3d([x_middle - plot_radius, x_middle + plot_radius])
  685. ax.set_ylim3d([y_middle - plot_radius, y_middle + plot_radius])
  686. ax.set_zlim3d([z_middle - plot_radius, z_middle + plot_radius])
  689. # %%
  690. # main script
  691. binning_success = 0 #keep track of whether triangles were successfuly constructed and binned
  692. for j,subject in enumerate(subjects): #loop through the participant's files
  693. tdel = tdel_values[j]
  694. print('Doing subject',subject)
  695. #MEG, sEEG files; made using other resources not included in this code base
  696. coordinate_file = coordinate_files[j]
  697. distance_file = geodesic_distance_files[j]
  698. timeseries_file = TS_files[j]
  699. triangle_file = triangle_files[j]
  700. print(triangle_file)
  701. metres_sd,dmetres_SS = get_coordinates_and_distances(file_path,coordinate_file,distance_file,distance_units)
  702. #save relevant paramenters with subject for later analysis
  703. aggregate_subjects[subject] = {}
  704. aggregate_subjects[subject]['tdel'] = tdel
  705. aggregate_subjects[subject]['metres_sd'] = metres_sd
  706. aggregate_subjects[subject]['dmetres_SS'] = dmetres_SS
  707. nSensors = metres_sd.shape[0]
  708. #look for triplet files, make them if they don't exist (slow, taking many hours for large number of contacts, so we save them after making)
  709. if exists(file_path+triangle_file+'_nonskinny_threshold_%.4f.pkl'%nonskinny_threshold): nonskinny_triplets = load_pkl(file_path,triangle_file+'_nonskinny_threshold_%.4f'%nonskinny_threshold)
  710. else:
  711. nonskinny_triplets = collate_nonskinny_triangles(dmetres_SS,nonskinny_threshold=nonskinny_threshold)
  712. save_pkl(file_path,nonskinny_triplets,triangle_file+'_nonskinny_threshold_%.4f'%nonskinny_threshold)
  713. if verbose: print('nonskinny_triplets shape',nonskinny_triplets.shape)
  714. #put the triangles into bins of the approximately the same linear sizes
  715. binning_out = bin_triplets_linearSF(dmetres_SS,nonskinny_triplets,return_binwise=False,smallest_bin_size=smallest_triangle_size,largest_bin_size=largest_triangle_size)
  716. if binning_out is not None: #if None, failed to bin the triplets with current values for nonskinny_threshold, smallest_bin_proportion, or because small number of distances
  717. #housekeeping from the binning
  718. binning_success += 1
  719. analysis_info_string = 'nB_%d_nC_%d_nonskinny_th_%.4f_t_range_%.3f_to_%.3f'%(nBases,N_cycles,nonskinny_threshold,smallest_triangle_size,largest_triangle_size)
  720. print('analysis_info_string',analysis_info_string)
  721. n_per_linearSF_bins,binned_triplets_av,normalization_a,range_a,root_areas_a,linearSF_bin_edges = binning_out
  722. aggregate_subjects[subject]['binned_triplets_av'] = binned_triplets_av
  723. aggregate_subjects[subject]['normalization_linearSF_a'] = normalization_a #for naming consistency
  724. aggregate_subjects[subject]['range_a'] = range_a
  725. aggregate_subjects[subject]['root_areas_a'] = root_areas_a
  726. aggregate_subjects[subject]['linearSF_bin_edges'] = linearSF_bin_edges
  727. bin_counts = n_per_linearSF_bins[n_per_linearSF_bins>0]
  728. do_bins_string = ''
  729. min_triangle = range_a[0]
  730. max_triangle = range_a[-1]
  731. print('n_per_linearSF_bins',n_per_linearSF_bins)
  732. print(j,binning_success,subject)
  733. if verbose:
  734. plot_title = triangle_file+'_nonskinny_threshold_%.4f_smallest_triangle_size_%.3f_largest_triangle_size_%.3f%s'%(nonskinny_threshold,smallest_triangle_size,largest_triangle_size,do_bins_string)
  735. draw_triangles(metres_sd,binned_triplets_av,bin_counts,output_path,plot_title,root_areas_a)
  736. else: continue #try the next subject
  737. max_distance = dmetres_SS[dmetres_SS<0.5].max() #0.5 (metres) is used as the mask term since it is larger than the largest human head
  738. max_of_range = 2*max_triangle #used for plotting
  739. aggregate_subjects[subject]['max_distance'] = max_distance
  740. aggregate_subjects[subject]['max_of_range'] = max_of_range
  741. aggregate_subjects[subject]['min_of_range'] = min_of_range
  742. timeseries_cts = np.load(file_path+timeseries_file+'.npy')
  743. if verbose: print('timeseries_cts shape',timeseries_cts.shape)
  744. for frequency in frequencies:
  745. aggregate_subjects[subject][frequency] = {}
  746. list_of_estpower = []
  747. list_of_wavelength = []
  748. #estimate phase
  749. data_label = 'empirical'
  750. print('Empirical bases from SVD')
  751. PHname = "%sPH_%.4fHz_%s_%s_%.4fminf" %(file_path,frequency,subject,measurement,minf)
  752. phi_cts,pow_cts = wavelet(timeseries_cts,N_cycles,frequency,tdel,minf,PHname=PHname)
  753. if include_normed_power:
  754. data_label = 'npower'
  755. pow___s = pow_cts.reshape(-1,pow_cts.shape[2]).mean(0) #normalize power per each grey-matter contact
  756. pow_cts = pow_cts / pow___s[np.newaxis,np.newaxis,:]
  757. phi_cts = pow_cts * phi_cts #now weighted by normalized power
  758. else: data_label = 'phase'
  759. del pow_cts
  760. # and make left singular vectors
  761. bases_sb,sigmas_b,betas_ctb,reduced_cts = complex_spatial_SVD(phi_cts,nBases=nBases,bases_only=False,use_bases=None,unit_phase_out=False)
  762. if verbose:
  763. print('bases_sb shape',bases_sb.shape)
  764. print('absolute(bases_sb) min',np.absolute(bases_sb).min(),'absolute(bases_sb) max',np.absolute(bases_sb).max())
  765. phase_plot(bases_sb,frequency,metres_sd,mag=True)
  766. # remove the negative wave components
  767. if forward_waves:
  768. vel_ct,TWf_cts = decompose_rphi_into_forward_waves(bases_sb.T[:,np.newaxis,:],tdel,frequency)
  769. if verbose: print('TWf_cts shape',TWf_cts.shape)
  770. print('Normalized velocity of spatial vector of phase (SW=0,TW=1)')
  771. print(vel_ct.T)
  772. Bases_sb = TWf_cts[:,0,:].T
  773. if verbose: print('absolute(Bases_sb) min',np.absolute(Bases_sb).min(),'absolute(Bases_sb) max',np.absolute(Bases_sb).max())
  774. print('Empirical traveling wave components')
  775. if verbose: phase_plot(Bases_sb,frequency,metres_sd,mag=True)
  776. else: Bases_sb = bases_sb
  777. aggregate_subjects[subject][frequency]['Bases_sb'] = Bases_sb
  778. #estimate the SF
  779. estpower_ba,metres_per_cycle_ba,dphidx_ba,dphi_ba = make_SF_power(dmetres_SS,binned_triplets_av,normalization_a,sigmas_b,Bases_sb)
  780. list_of_estpower += [estpower_ba.reshape(-1)]
  781. list_of_wavelength += [metres_per_cycle_ba.reshape(-1)]
  782. plot_name = '%s %s SF %s TF %.2fHz'%(subject,measurement,data_label,frequency)
  783. aggregate_subjects[subject][frequency]['list_of_estpower'] = list_of_estpower
  784. aggregate_subjects[subject][frequency]['list_of_wavelength'] = list_of_wavelength
  785. save_pkl(output_path,aggregate_subjects[subject],(plot_name+'_'+analysis_info_string).replace(' ','_'))
  786. sf_histogram(output_path,plot_name,analysis_info_string,list_of_estpower,list_of_wavelength,\
  787. [max_triangle],max_of_range,min_of_range,N_histogram_bins,colourmap='cool')
  788. del aggregate_subjects[subject][frequency]
  789. clear_output(wait=True) #.pynb file can get very large, especially with verbose==True
  790. print('binning success',binning_success)
  792. # %%

estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb at commit 42c7fe1, under GPL-3.0 · at the source

Overview

Authors: David M Alexander1,2, Laura Dugué1,3
ORCID iDs: David M Alexander, Laura Dugué
  1. Université Paris Cité, CNRS, Integrative Neuroscience and Cognition Center Paris France
  2. GHU-Paris Psychiatrie et Neurosciences, Hôpital Sainte-Anne Paris France
  3. Institut Universitaire de France (IUF) Paris France
Institutions: Centre Neurosciences intégratives et Cognition (France); GHU Paris psychiatrie & neurosciences; Institut Universitaire de France (France)
Journal: eLife, volume 13, article RP100674
Dates: published online 9 April 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.7554/elife.100674 · PMID 41954599 · PMCID PMC13065328 · OpenAlex W4403758897
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), intracranial EEG (iEEG / ECoG / SEEG) (modality), human (organism)
Methods: Preprocessing, Connectivity, Statistics, Spectral & time-frequency, Evoked potentials, Machine learning, fMRI & imaging
Keywords: cortex, traveling waves, sEEG, spatial frequency, temporal frequency, large-scale, Human
MeSH: Cerebral Cortex*, Delta Rhythm*, Gamma Rhythm*, Electroencephalography, Gray Matter, Humans (* major topic)
Journal subjects: Neuroscience
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: European Research Council (European Union's Horizon 2020 research and innovation programme (grant No 852139 - Laura Dugué))
Citations: cited by 2 papers (Europe PMC); 103 references in the paper

Abstract

The organization of the phase of electrical activity in the cortex is critical to inter-site communication, but the balance of this communication across large-scale (>8 cm), macroscopic (>1 cm), and mesoscopic (1 cm to 1 mm) ranges is an open question. The spatial frequencies (i.e. the spatial scales) of cortical waves have been characterized in the gray matter for micro- and mesoscopic scales of cortex and show decreasing spatial power with increasing spatial frequency. This research, however, has been limited by the size of the measurement array, thus excluding large-scale traveling waves. Obversely, poor spatial resolution of extracranial measurements prevents incontrovertible large-scale estimates of spatial power. We estimate the spatial frequency spectrum of phase dynamics in order to quantify the uncertain large-scale range, utilizing stereotactic electroencephalogram to measure local-field potentials within the gray matter. We take advantage of the large extent of spatial coverage of the cortical sheet, and irregular sampling is offset by use of linear algebra techniques. We find the spatial power of the phase is highest at the lowest spatial frequencies (longest wavelengths), consistent with the power spectra ranges for micro- and meso-scale dynamics, but here shown up to the size of the measurement array (up to 8–16 cm). This result arises across a wide range of temporal frequencies, from the delta band (1–3 Hz) through to the high gamma range (60–100 Hz).

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.

duguelab/traveling-wave-analysis

License: GPL-3.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: 42c7fe112c99a06e05aaf31939a0eb3e218e0fcc, 6 January 2026
Languages: Python (5), Jupyter (5)
Size: 15 files, 10 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, 5 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (10 files), Matplotlib (5 files), SciPy (2 files)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
12 files
At the source: github.com/DugueLab/Traveling-wave-analysis/blob/main/estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 10 scripts, each with its path and the digest of its content;
  • 7 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The key Python routines used in this research can be found here: https://github.com/DugueLab/Traveling-wave-analysis/blob/main/estimation_of_SF_of_phase_of_cortical_activity_on_irregular_measurement_arrays.ipynb (copy archived at Alexander, 2026).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 2 authors, 7 keywords, 6 MeSH terms, 1 funder, 95 references.

Cite

This paper

Alexander, D. M., & Dugué, L. (2026). The dominance of large-scale phase dynamics in human cortex, from delta to gamma. eLife, 13, RP100674. https://doi.org/10.7554/elife.100674

BibTeX

@article{alexander2026dominance,
author = {Alexander, David M and Dugué, Laura},
title = {{The dominance of large-scale phase dynamics in human cortex, from delta to gamma}},
journal = {eLife},
year = {2026},
month = apr,
volume = {13},
pages = {RP100674},
publisher = {eLife Sciences Publications, Ltd},
issn = {2050-084X},
doi = {10.7554/elife.100674},
url = {https://doi.org/10.7554/elife.100674},
pmid = {41954599},
pmcid = {PMC13065328}
}

RIS

TY - JOUR
AU - Alexander, David M
AU - Dugué, Laura
TI - The dominance of large-scale phase dynamics in human cortex, from delta to gamma
T2 - eLife
J2 - eLife
PY - 2026
DA - 2026/04/09
VL - 13
SP - RP100674
SN - 2050-084X
PB - eLife Sciences Publications, Ltd
DO - 10.7554/elife.100674
UR - https://doi.org/10.7554/elife.100674
LA - en
ER -

CSL-JSON

{
"id": "10.7554/elife.100674",
"type": "article-journal",
"title": "The dominance of large-scale phase dynamics in human cortex, from delta to gamma",
"container-title": "eLife",
"author": [
{
"family": "Alexander",
"given": "David M"
},
{
"family": "Dugué",
"given": "Laura"
}
],
"container-title-short": "eLife",
"volume": "13",
"page": "RP100674",
"DOI": "10.7554/elife.100674",
"PMID": "41954599",
"PMCID": "PMC13065328",
"ISSN": "2050-084X",
"publisher": "eLife Sciences Publications, Ltd",
"URL": "https://doi.org/10.7554/elife.100674",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
9
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1371/journal.pcbi.1013007 [code]
Traveling waves in the human visual cortex: An MEG-EEG model-based approach
Journal: n/a
In common: SciPy, Matplotlib, NumPy, EEG, 15 references, author Laura Dugué
[2] doi:10.7554/elife.106753 [code]
Traveling waves across scales: Different mechanisms but same canonical computation?
Journal: n/a
In common: 14 references, author Laura Dugué
[3] doi:10.1038/s41467-026-71386-z [code]
Planar, spiral, and concentric traveling waves distinguish behavioral states in human memory.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, 10 references
[4] doi:10.1073/pnas.2527296123
Traveling-wave transcranial alternating current stimulation (twtACS) causally links neural timing to cognitive function.
Journal: Proceedings of the National Academy of Sciences of the United States of America
In common: intracranial EEG (iEEG / ECoG / SEEG), 8 references
[5] doi:10.3389/fncom.2026.1844662
Quantifying cortex-wide traveling brain waves of complex patterns with a graph-based algorithm.
Journal: Frontiers in computational neuroscience
In common: intracranial EEG (iEEG / ECoG / SEEG), 8 references
[6] doi:10.7554/elife.108208 [code]
Realistic coupling enables flexible macroscopic traveling waves in the mouse cortex.
Journal: eLife
In common: NumPy, 7 references
[7] doi:10.1016/j.isci.2026.116728 [code]
Awake cortex stabilizes traveling waves for global and reliable information routing.
Journal: iScience
In common: SciPy, NumPy, 4 references
[8] doi:10.1038/s41467-026-72931-6 [code]
Three parsimonious spatiotemporal patterns in cerebellum reveal individual traits in function and behavior.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, 3 references
[9] doi:10.1038/s41467-026-76011-7 [code]
Human cortex organizes dynamic co-fluctuations along the sensorimotor-association axis.
Journal: Nature communications
In common: SciPy, Matplotlib, NumPy, 3 references
[10] doi:10.1038/s42003-025-09444-3 [code]
Decoupling of neurophysiological activity from structure mirrors global microarchitectural and neuromodulatory trends.
Journal: Communications biology
In common: SciPy, Matplotlib, NumPy, 3 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.