# This file contains all the parameters. # Additional comments are provided in the header of main.py print '********' print '* Init *' print '********' import os isunix = lambda: os.name == 'posix' import matplotlib if isunix(): matplotlib.use('Agg') # no graphic link with Unix cluster for me import brian_no_units from brian import * from scipy.io import * from scipy import weave from time import time from time import localtime from customrefractoriness import * import pickle import glob set_global_preferences(useweave=True) # to compile C code globalStartTime=time()*second #************************* # COMPUTATION PARAMETERS * #************************* neuronTimeCompression = 2.0**0 # to compress all the neuronal time constants pbTimeCompression = 2.0**0 # to compress the problem and oscillation time constants get_default_clock().set_dt(.1*ms/neuronTimeCompression) # set time step # random state randState = 28 # use this to specify it # # use this for batches: #rl = glob.glob('../data/rand*.mat') #if size(rl)==0: # randState = 0 #else: # randState = int(rl[-1][12:15])+1 #savemat('../data/rand'+'%03d' % (randState)+'.mat',{'tmp':[]}) seed(randState) imposedEnd = 3*second/pbTimeCompression # imposed end time N = 2000 # number of presynaptic neurons nG = 1 # number of gmax values nR = 1 # number of ratio LTD/LTP M = nG*nR # number of postsynaptic neurons, numbered like that [ (r_0,g_0)...(r_0,g_nG),(r_1,g_0)...(r_1,g_nG),...,(r_nR,g_0)...(r_nR,g_nG)] recomputeSpikeList = True # recompute input spike trains as opposed to load appropriate mat files dumpSpikeList = False # save input spike list in mat files computeOutput = True # compute output (STDP) layer useSavedWeight = False # load previously dumped weights timeOffset = 0*second # simulation starts at t=timeOffset useReset = False # impose resets on input layer at dates specified in reset.###.mat file graph = True # graph output monitorInput = True # monitor input spikes monitorOutput = True # monitor output spikes monitorPot = True # monitor potential in output layer monitorCurrent = False # monitor potential in output layer monitorInputPot = False # monitor potential in input layer monitorRate = False # monitor rates in output layer isMonitoring = False # flag saying if currently monitoring monitorTime = (imposedEnd-6/pbTimeCompression)*second # start monitoring only at that time (to save memory) analyzePeriod = Inf # periodically launches analyze.py if not recomputeSpikeList and dumpSpikeList: print 'Warning: dumping a spike list which is not re-computed makes no sense. Setting dumpSpikeList to False' dumpSpikeList = False if not computeOutput and ( monitorOutput or monitorPot or monitorCurrent or monitorRate): print 'Warning: can not monitor output, which is not computed. Setting monitorOutput, monitorPot and monitorCurrent to False' monitorOutput = False monitorPot = False monitorCurrent = False monitorRate = False # load pattern values (the values are only useful for plotting), but the length of the vector is used to scale gmax if os.path.exists(os.path.join('..','data','realValuedPattern.'+'%03d' % (randState)+'.mat')): realValuedPattern=loadmat(os.path.join('..','data','realValuedPattern.'+'%03d' % (randState)+'.mat')) realValuedPattern=realValuedPattern['realValuedPattern'] else: realValuedPattern=zeros(round(.5*N)) #************************** # NEURON MODEL PARAMETERS * #************************** # Types of input and output neurons poissonInput = False # poisson input neurons (as opposed to LIF) conductanceOutput = False # conductance-base output neurons (as opposed to LIF). Tim: 9/2008: has never been critical so far poissonOutput = False # stochastic (Poisson) output neurons (as opposed to deterministic). Note that their differential equations are the same as the LIF, but firing is stochastic. #neurons (Dayan&Abbott 2001) refractoryPeriod = 1*ms/neuronTimeCompression R=10*Mohm if poissonOutput: vt=400 # not used (just for graph scaling) vr=0 # not used El=-450 # resting # El=-280 # resting (use that value to have more false alarms) else: vt=-54*mV # threshold vr=-60*mV # reset El=-70*mV # resting Ee=0*mV # for excitatory conductance taum = 20*ms/neuronTimeCompression # membrane time constant taue=taum/4 # synapse time constant sigma = 0.015*(vt-vr) # white Gaussian noise. Applies to input and output neurons # array of max conductance values unitaryEffect = taue/(taum-taue)*((taum/taue)**(-taue/(taum-taue))-(taum/taue)**(-taum/(taum-taue))) # this corresponds to the maximum of the kernel (exp(-t/taum) - exp(-t/taue)) taue / (taum-taue) if poissonOutput: # gmax=7*100.0/(1.0*N)*1.0/unitaryEffect*exp(log(2)*(array(nR*range(nG))-nG/2)) # appropriate for (1000/2000 afferents, LIF+reset) gmax=1.05**-10*2.5*1500/(1.0*size(realValuedPattern))*1.0/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # appropriate for (1000/2000 afferents, PLoS patten) # gmax=1.05**-16*2.5*1500/(1.0*size(realValuedPattern))*1.0/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # appropriate for (1000/2000 afferents, N1N2) else: # LIF+reset every 250ms, all spikes # gmax = 1.05**12/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=2.5% # gmax = 1.05**18/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=5% # gmax = 1.05**24/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=10% # gmax = 1.05**30/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=20% # gmax = 1.05**36/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=40% # LIF+reset every 125ms, all spikes (Tim 03/09) # gmax = 1.05**2/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=2.5% # gmax = 1.05**8/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=5% # gmax = 1.05**14/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=10% # gmax = 1.05**20/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=20% # gmax = 1.05**26/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=40% # 8Hz sin oscillations with 1-3 spikes/cycle, nearest spike # gmax = 1.05**4/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x2^-2 # gmax = 1.05**4/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x2^-1.5 # gmax = 1.05**2/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x2^-1 # gmax = 1.05**-2/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x1.0 # gmax = 1.05**4/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x2^1.25 # gmax = 1.05**4/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x2^1.5 # gmax = 1.05**-2/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x2^1.75 # gmax = 1.05**6/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x2^2.0 # gmax = 1.05**6/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x2^2.5 # 8Hz sin oscillations with 1-3 spikes/cycle, all spikes # gmax = 1.05**-16/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=2.5% # gmax = 1.05**-8/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=5% gmax = 1.05**-2/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=10% # gmax = 1.05**8/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=20% # gmax = 1.05**12/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # x=40% # gmax = 1.05**4/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x 2^-2 # gmax = 1.05**2/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x 2^-1 # gmax = 1.05**-2/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x 1.0 # gmax = 1.05**-6/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) # time compression x 2 # gmax = 1.05**-10/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x 2^2 # gmax = 1.05**-14/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x 2^2.5 # gmax = 1.05**-22/(size(realValuedPattern))*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## time compression x 2^3 # academic patterns (Matthieu Gilson) # gmax = 1.05**30/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=25% LIF # gmax = 1.05**38/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=10% LIF # gmax = 1.05**34/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=6.25% LIF # gmax = 1.05**50/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=100% of 200 afferents, LIF # gmax = 1.05**38/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(2*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=100% of 20 afferents, LIF # gmax = 1.05**30/size(realValuedPattern)*(vt-El)/unitaryEffect*exp(3*log(1.05)*(array(nR*range(nG))-nG/2)) ## x=50% LIF if conductanceOutput: gmax /= (Ee-vt) print 'gmax=' + str(gmax) #******************** # WEIGHT PARAMETERS * #******************** # initial synaptic weight are randomly picked (uniformly) between those two bounds if poissonOutput: initialWeight_min = 0 initialWeight_max = 10 else: initialWeight_min = 0*volt initialWeight_max = 2*8.6e-5*volt burstingCriterion = .5 # unplug neurons whose mean normalized synaptic weight is above this value. This allows to save memory by not computing neurons whose normalized synaptic weights all go to 1. #*************************** # INPUT CURRENT PARAMETERS * #*************************** # used only if recomputeSpikeList = True # input current (normalization with respect to raw values in input file) if poissonInput: # in this case I values are in fact firing rates Imin = 0 Imax = 100 else: # Imax = (1.05)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (LIF+reset) # Imin = (1.0)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (LIF+reset) # Imax = (1.140)*(vt-El)/R # Tim 12/08: appropriate for (8Hz sawtooth oscillations). up to 3 spikes/cycle # Imin = (1.025)*(vt-El)/R # Tim 12/08: appropriate for (8Hz sawtooth oscillations). up to 3 spikes/cycle Imax = (1.07)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (8Hz sin oscillations). 1-3 spikes/cycle Imin = (0.95)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (8Hz sin oscillations). 1-3 spikes/cycle # Imax = (1.012)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (8Hz sin oscillations). 1-3 spikes/cycle. taum=10ms # Imin = (0.936)*(vt-El)/R*ones(N) # Tim 12/08: appropriate for (8Hz sin oscillations). 1-3 spikes/cycle. taum=10ms # Imax = (1.2)*(vt-El)/R*ones(N) # # Imin = (0.95)*(vt-El)/R*ones(N) # # Imax = (vt-El)/R * (1.07)*ones(N) # partial drive # Imin = (vt-El)/R * (0.95)*ones(N) # partial drive # Imax = (0.992)*(vt-El)/R # Tim 12/08: appropriate for (8Hz sin oscillations). 1 spike/cycle # Imin = (0.95)*(vt-El)/R # Tim 12/08: appropriate for (8Hz sin oscillations). 1 spike/cycle # oscillation amplitude (use 0 not to have any). An array is used so that different amplitudes may be used for the input neurons #a = .075*(vt-El)/R * array(12*range(9))/8.0 #for i in range(9): # a[i*12:(i+1)*12] = .075*(vt-El)/R * (i/8.0) a = .075*(vt-El)/R*ones(N) # oscillation fequence oscilFreq = 8*Hz*pbTimeCompression #****************** # STDP PARAMETERS * #****************** nearestSpike = False # Tim 08/2008: not critical with low input spike rates (10Hz) tau_post=33.7*ms #(source: Bi & Poo 2001) tau_pre=16.8*ms #(source: Bi & Poo 2001) #tau_post=20*ms #(source: Song, Miller & Abbott 2000) #tau_pre=20*ms #(source: Song, Miller & Abbott 2000) a_pre= .005 # Tim 12/08: no more that .005 w_out = - 0*.005/3 # see Matt Gilson w_in=zeros(M) # array of w_in for i in range(nR): # here nR corresponds to number of w_in/w_out ratios (and not number of LTD/LTP ratios) # see Matt Gilson w_in[i*nG:(i+1)*nG] = -w_out*1.05**0*exp((i-nR/2)*.5*log(2)) a_post=zeros(M) # array of LTD/LTP ratios for i in range(nR): # LIF+reset every 250ms, all spikes # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-5*exp((i-nR/2)*log(1.05)) # x=2.5%, 5% 10% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-2*exp((i-nR/2)*log(1.05)) # x=20% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**0*exp((i-nR/2)*log(1.05)) # x=40% # LIF+reset every 125ms, all spikes # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-1*exp((i-nR/2)*1*log(1.05)) # # 8Hz sin oscillations with 1-3 spikes/cycle, nearest spike # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**3*exp((i-nR/2)*1*log(1.05)) # time compression x2^-2 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**5*exp((i-nR/2)*1*log(1.05)) # time compression x2^-1.5 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**8*exp((i-nR/2)*1*log(1.05)) # time compression x2^-1 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**9*exp((i-nR/2)*1*log(1.05)) # time compression x1.0 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**7*exp((i-nR/2)*1*log(1.05)) # time compression x2^1 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**7*exp((i-nR/2)*1*log(1.05)) ## time compression x2^1.25 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**5*exp((i-nR/2)*1*log(1.05)) # time compression x2^1.5 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**5*exp((i-nR/2)*1*log(1.05)) ## time compression x2^1.75 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**5*exp((i-nR/2)*1*log(1.05)) ## time compression x2^2 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**3*exp((i-nR/2)*1*log(1.05)) # time compression x2^2.5 # 8Hz sin oscillations with 1-3 spikes/cycle, all spikes # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**9*exp((i-nR/2)*1*log(1.05)) # x=2.5% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**9*exp((i-nR/2)*1*log(1.05)) # x=5% a_post[i*nG:(i+1)*nG] = -a_pre*1.05**8*exp((i-nR/2)*1*log(1.05)) # x=10% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**8*exp((i-nR/2)*1*log(1.05)) # x=20% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**8*exp((i-nR/2)*1*log(1.05)) # x=40% # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**3*exp((i-nR/2)*.5*log(1.05)) ## time compression x2^-2 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**7*exp((i-nR/2)*.5*log(1.05)) # time compression x2^-1 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**8*exp((i-nR/2)*.5*log(1.05)) ## time compression x1.0 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**2*exp((i-nR/2)*.5*log(1.05)) # time compression x2 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-5*exp((i-nR/2)*.5*log(1.05)) ## time compression x2^2 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-8*exp((i-nR/2)*.5*log(1.05)) ## time compression x2^2.5 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-10*exp((i-nR/2)*.5*log(1.05)) ## time compression x2^3 # academic patterns (Matthieu Gilson) # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-11*exp((i-nR/2)*1*log(1.05)) # Tim 05/09: x=25% LIF # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-7*exp((i-nR/2)*1*log(1.05)) # Tim 05/09: x=10% LIF # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-10*exp((i-nR/2)*1*log(1.05)) # Tim 06/09: x=100% of 200 or 20 afferents, LIF # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-12*exp((i-nR/2)*1*log(1.05)) # Tim 06/09: x=50% Poisson, PLoS pattern # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-9*exp((i-nR/2)*2*log(1.05)) # Tim 06/09: x=50% Poisson, N1N2 pattern # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-9*exp((i-nR/2)*1*log(1.05)) # Tim 06/09: x=50% Poisson, PLoS pattern, mu = 0.1 # a_post[i*nG:(i+1)*nG] = -a_pre*1.05**-11*exp((i-nR/2)*1*log(1.05)) # Tim 07/09: x=50% LIF mu = 0.0 # see Gutig et al. 2003. print 'normalized a_post/a_pre ratios' + str(a_post/a_pre) #*********** # FUNCIONS * #*********** def printtime(mess): t = localtime() print '%02d' % t[3] + ':' + '%02d' % t[4] + ' ' + mess def listMatFile(directory,randState): "get list of spike lists" fileList = os.listdir(directory) fileList.sort() fileList = [f for f in fileList # e.g. spikeList.200.010.mat # 01234567890123456789 if f[0:9] in ['spikeList'] and f[10:13] in ['%03d' % (randState)] and f[-4:] in ['.mat'] ] return fileList # Called whenever output neurons fire. Resets the potential and trigger STDP updates. # Includes C code, will be compiled the first time it is called # Param: # P: NeuronGroup (output neurons) # spikes: list of the indexes of the neuron that fire. def neurons_reset(P,spikes): if size(spikes): if not poissonOutput: P.v_[spikes]=vr # reset pot if nearestSpike: nspikes = size(spikes) A_pre = mirror.A_pre_ code = ''' for(int si=0;si_gmax(i)) wnew = _gmax(i); if(wnew<0) wnew = 0.0; } else { /* soft bound */ wnew = _synW(j,i)+_gmax(i)*(w_out+A_pre(j)*exp(mu*log(1-_synW(j,i)/_gmax(i)))); if(wnew>_gmax(i)) wnew = _gmax(i); if(wnew<0) wnew = 0.0; } _synW(j,i) = wnew; _alreadyPotentiated(j,i) = true; } } } ''' weave.inline(code, ['spikes', 'nspikes', 'N', '_alreadyPotentiated', '_synW', '_gmax', 'A_pre', 'mu','w_out'], compiler='gcc', type_converters=weave.converters.blitz, extra_compile_args=['-O3']) _alreadyDepressed[:,spikes] = False P.A_post_[spikes]=a_post[spikes] # reset A_post (~start a timer for future LTD) else: # all spikes nspikes = size(spikes) A_pre = mirror.A_pre_ code = ''' for(int si=0;si_gmax(i)) wnew = _gmax(i); if(wnew<0) wnew = 0.0; } else { /* soft bound */ wnew = _synW(j,i)+_gmax(i)*(w_out+A_pre(j)*exp(mu*log(1-_synW(j,i)/_gmax(i)))); if(wnew>_gmax(i)) wnew = _gmax(i); if(wnew<0) wnew = 0.0; } _synW(j,i) = wnew; } } ''' weave.inline(code, ['spikes', 'nspikes', 'N', '_synW', '_gmax', 'A_pre', 'mu','w_out'], compiler='gcc', type_converters=weave.converters.blitz, extra_compile_args=['-O3']) P.A_post_[spikes]+=a_post[spikes] # reset A_post (~start a timer for future LTD) # Called whenever input neurons fire. Resets the potentials and trigger STDP updates. # Note that mirror is a fake group that mirrors input neuron, only used for implementation issues. # Includes C code, will be compiled the first time it is called # Param: # P: NeuronGroup (input neurons) # spikes: list of the indexes of the neuron that fire. def mirror_reset(P,spikes): if size(spikes): P.v_[spikes] = 0 if nearestSpike: nspikes = size(spikes) A_post = neurons.A_post_ code = ''' for(int si=0;si_gmax(j)) wnew = _gmax(j); if(wnew<0.0) wnew = 0.0; } else { /* soft bound */ wnew = _synW(i,j)+_gmax(j)*(w_in(j)+A_post(j)*exp(mu*log(_synW(i,j)/_gmax(j)))); if(wnew>_gmax(j)) wnew = _gmax(j); if(wnew<0.0) wnew = 0.0; } _synW(i,j) = wnew; _alreadyDepressed(i,j) = true; } } } ''' weave.inline(code, ['spikes', 'nspikes', 'M', '_alreadyDepressed', '_synW', '_gmax', 'A_post', 'mu','w_in'], compiler='gcc', type_converters=weave.converters.blitz, extra_compile_args=['-O3']) _alreadyPotentiated[spikes,:]=False P.A_pre_[spikes]=a_pre # reset A_pre (~start a timer for future LTP) else: # all spikes nspikes = size(spikes) A_post = neurons.A_post_ code = ''' for(int si=0;si_gmax(j)) wnew = _gmax(j); if(wnew<0.0) wnew = 0.0; } else { /* soft bound */ wnew = _synW(i,j)+_gmax(j)*(w_in(j)+A_post(j)*exp(mu*log(_synW(i,j)/_gmax(j)))); if(wnew>_gmax(j)) wnew = _gmax(j); if(wnew<0.0) wnew = 0.0; } _synW(i,j) = wnew; } } ''' weave.inline(code, ['spikes', 'nspikes', 'M', '_synW', '_gmax', 'A_post', 'mu','w_in'], compiler='gcc', type_converters=weave.converters.blitz, extra_compile_args=['-O3']) P.A_pre_[spikes]+=a_pre # reset A_pre (~start a timer for future LTP)