Raw File
"""
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')
            
back to top