STiMCON_Fig6.py
"""
Created on Fri May 1 10:36:03 2020
@author: sanoev
Shows how acoustic time and stimulus time is not the same in STiMCON (Figure 6)
"""
import math
import numpy as np
import STiMCON_core
import STiMCON_sen
import matplotlib.pyplot as plt
import ColorScheme as cs
cmaps = cs.CCcolormap()
bfi = cs.baseFigInfo()
#%% basic information
fs = 1000;
Freq = 4.0;
#%% create the language model constrainst:
# model is: I do/I eat/I go
LMnames = np.array(['a','b','c','d','e'])
feedbackmat = np.zeros([5,5])
feedbackmat[0] = [0, 0.0, 0.8, 0.8, 0.0]
feedbackmat[1] = [0.0, 0, 0, 0,0]
feedbackmat[2] = [0.8, 0, 0, 0,0]
feedbackmat[3] = [0.0, 0, 0, 0,0]
feedbackmat[4] = [0.8, 0, 0, 0,0]
Nnodes = len(feedbackmat)
#%% define the parameters of the model
parameters = {"Nnodes": Nnodes,
"OsFreq": Freq,
"OsAmp": 1,
"OsOffset": 0.25*math.pi,
"activation_threshold": 1,
"feedbackmat": feedbackmat,
"feedbackinf": 1.5,
"feedbackdecay": 0.01,
"feedbackdelay": int(0.9/Freq*fs),
"latinhibstrength": 0,
"selfexitation": 0,
"Inhib": -0.2,
"fs": fs,
'LMnames': LMnames}
#%% define the parameters for the sensory input
stimpara = {'word_duration': int(0.5/Freq*fs),
'onsetdelay': int(0.5/Freq*fs),
'Nnodes': Nnodes}
#% adjust OsOffset based on the onsetdelay and word_duration:
peak = (stimpara['word_duration']+stimpara['onsetdelay'])/fs
parameters['OsOffset'] = peak*Freq*(2*math.pi)
## set all the parameters
senObj = STiMCON_sen.modelSen(stimpara,parameters)
#%% things for fft calculations
def calcFFT(tofft, fs = 1000, pad = 0):
tofft = (tofft-np.mean(tofft))
pad = int(0*fs)
tofft = np.concatenate((np.zeros(pad),tofft,np.zeros(pad)),axis=0)
N = tofft.shape[-1]
Han = np.hanning(N)
tofft = tofft*Han
freq = np.fft.fftfreq(N, 1/fs)
FourSen = np.fft.fft(tofft)
return FourSen, freq
#%% create itterations of sensory input at different temporal offsets
import numpy.matlib
timeoff = np.linspace(-0.1,0.1,19)
thres = np.linspace(0.2,1.5,15)
FI = np.arange(0,100)
stimnames = ['low low','high high','low_high','high_low']
SEN = np.zeros([len(stimnames),len(timeoff),len(thres),len(FI)])
ACT = np.zeros([len(stimnames),len(timeoff),len(thres),len(FI)])
SUPRA = np.zeros([len(stimnames),len(timeoff),len(thres),len(FI)])
FIRSTSUPRA = np.zeros([len(stimnames),len(timeoff),len(thres),len(FI)])
SenType = ''
for stim in range(4):
if stim == 0: # high-high
v = 1
elif stim == 1: # low-low
v = 2
elif stim == 2: # high low
v= 3
elif stim == 3: # low high
v= 4
Ust = np.array([0,v])
st = numpy.matlib.repmat(Ust,1,5)[0]
for th in range(len(thres)):
for it in range(len(timeoff)):
stimtime= np.linspace(0,(len(st)-1)/Freq,len(st))*fs
j = np.arange(1,len(stimtime),2)
stimtime[j] = stimtime[j] + timeoff[it]*fs
if SenType == 'linear':
seninput = {'stim_ord': st,
'intensity': np.zeros(len(st))+thres[th],
'stim_time': stimtime,
'tot_length': (len(st)+1)/Freq*fs}
sensory_input = senObj.create_stim_vartimethres(seninput)
else:
seninput = {'stim_ord': st,
'intensity': np.zeros(len(st))+thres[th],
'stim_time': stimtime,
'tot_length': (len(st)+1)/Freq*fs}
sensory_input = senObj.create_stim_vartimethresG(seninput,GausW=3)
SenType = ''
STiMCON_var = STiMCON_core.modelPara(parameters)
out = STiMCON_var.runsingle(sensory_input)
# ignore the first and last big (only sustained bit..)
I = np.arange(500,len(sensory_input[0])-500)
F,freq = calcFFT(np.mean(sensory_input[:,I],0))
Fact,freq = calcFFT(np.mean(out['activation'][Ust,:][:,I],0)) # overall activation
FsupraSen,freq = calcFFT(np.mean(out['spiketimes'][:,I]==2,0)) # overall activation
idx= np.abs(freq-Freq).argmin()
SEN[stim,it,th,:] = abs(F[FI])
ACT[stim,it,th,:] = abs(Fact[FI])
SUPRA[stim,it,th,:]= abs(FsupraSen[FI])
curnodac = 10
tempFA = np.zeros(out['spiketimes'].shape)
for T in range(out['spiketimes'].shape[1]):
if sum(out['spiketimes'][:,T]==2) > 0:
idx2 = np.argmax(out['spiketimes'][:,T]==2)
if idx2 != curnodac:
tempFA[idx2,T] = 1
curnodac = idx2
FIRSTSUPRA[stim,it,th,:] = abs(calcFFT(np.mean(tempFA[:,I],0))[0])[FI]
#%% now plot
plot = ''
import copy
axs = list()
fig = plt.figure(constrained_layout=True, figsize = (bfi.figsize.Col2/3*2,6))
gs = fig.add_gridspec(20, 5)
ths=np.arange(0,15)
for stim in range(4):
D = SUPRA[stim,:,:,idx]
if plot == 'SEN':
D = copy.deepcopy(SEN[stim,:,:,idx]).transpose()**2
D = D/1000
else:
D = copy.deepcopy(ACT[stim,:,:,idx]).transpose()**2
D = D/100000
rowmean = D.mean(axis=0)
rowstd = D.std(axis=0)
axs.append(fig.add_subplot(gs[stim*5:stim*5+5,0:3]))
pos = axs[stim].imshow(D, aspect = 'auto', origin = 'lower',cmap = 'OrRd')
axs[stim].set_title(stimnames[stim])
x = np.linspace(timeoff[0]*1000,timeoff[-1]*1000,5).astype(int)
axs[stim].set_xticks(np.linspace(0,len(timeoff),len(x)))
y = np.around(np.linspace(thres[0],thres[-1],5),1)
axs[stim].set_yticks(np.linspace(0,len(thres),len(y)))
axs[stim].set_yticklabels(y)
cbar = fig.colorbar(pos, ax=axs[stim])
cbar.set_label('power')
if stim == 3:
axs[stim].set_xlabel('odd word offset (ms)')
axs[stim].set_xticklabels(x)
else:
axs[stim].tick_params(axis='x', which='both', bottom=True, labelbottom=False)
# power
axs[stim].set_ylabel('stimulus intensity')
if plot != 'SEN':
cbar.mappable.set_clim(2,4.8)
#% slice at 0.8 and 1.0
ths = [9, 6]
for it in range(len(ths)):
th = len(thres)-ths[it]
axs.append(fig.add_subplot(gs[it*4:it*4+4,3:6]))
for stim in range(4):
if plot == 'SEN':
D = copy.deepcopy(SEN[stim,:,th,idx])**(2)
else:
D = copy.deepcopy(ACT[stim,:,th,idx])**(2)
tit = 'power at intensity ' + '{:.1f}'.format(thres[th])
axs[-1].plot(timeoff*1000, D)
axs[-1].set_title(tit)
if plot != 'SEN':
axs[-1].set_ylabel('power ($x10^{5}$)')
axs[-1].set_xticks(np.arange(-100,150,50))
axs[-1].set_yticks([3e5,4e5])
axs[-1].set_yticklabels([3,4])
else:
axs[-1].set_ylabel('power ($x10^{3}$)')
axs[-1].set_yticks([1e3,2e3])
axs[-1].set_yticklabels([1,2])
if it == 0:
axs[-1].tick_params(axis='x', which='both', bottom=True, labelbottom=False)
else:
axs[-1].set_xlabel('odd word offset (ms)')
to = [4,9,14]
th = len(thres)-ths[1]
for it in range(len(to)):
axs.append(fig.add_subplot(gs[(it+3)*3+2:(it+3)*3+3+2,3:6]))
for stim in range(4):
if plot == 'SEN':
D = copy.deepcopy(SEN[stim,to[it],th,:])**(20)
else:
D = copy.deepcopy(ACT[stim,to[it],th,:])**(20)
D = D/1e56
axs[-1].plot(freq[FI[1:-80]], D[1:-80])
axs[-1].set_ylabel('$ampl^{10}(x10^{56})$')
axs[-1].set_title('odd word offset ' '{:.1f}'.format(timeoff[to[it]]*1000))
axs[-1].set_xticks(np.arange(2,10,2))
if it < len(to)-1:
axs[-1].tick_params(axis='x', which='both', bottom=True, labelbottom=False)
else:
axs[-1].set_xlabel('frequency (Hz)')
plt.show()
#fig.savefig(figloc + 'PredRhythm' + SenType + plot + '.pdf', format='pdf')