import sys sys.path.extend([".."]) from data_utils import * def check_substr(): """check find_substr_endchars""" a = '101010100' print a print [ch for ch in find_substr_endchars(a,'1010')] ## test entropy rate def IID_entropyrate(): """equal prob of 1s and 0s IID. Entropy rate should be 1 bit per use.""" spiketrains = [[int(uniform(0,1)+0.5) for i in range(1000)]] print "Equal prob of 1s and 0s IID. Hrate =",calc_entropyrate(spiketrains,1) def markov_order1_entropyrate(): """markov process of order 1 p(0)=0.6, p(1)=0.4""" spiketrains = [] for j in range(10): mc = [1] for i in range(10000): if mc[-1]: p = uniform(0,1) if p<0.25: mc.append(1) else: mc.append(0) else: p = uniform(0,1) if p<0.5: mc.append(1) else: mc.append(0) spiketrains.append(mc) entropyrates = [] ## careful: don't go to higher orders.. machine stalls badly.. for i in range(10): entropyrates.append(calc_entropyrate(spiketrains,i)) figure(facecolor='w') plot(entropyrates,'r-,') def plot_table(rasters,rowlabels,collabels,data,cellcolours,titlestr,figfilename): ## 'plot' a table fig = figure(figsize=(8, 6), dpi=100) ax = fig.add_axes([0.14, 0.75, 0.85, 0.1]) axes_off(ax) ## loop over rasters in reverse order, as they are plotted from bottom upwards for rasteri,raster in enumerate(rasters[::-1]): raster = array(raster) ## find out indices of 1-s and plot them: rasterindices = where(raster==1)[0] ax.plot(rasterindices,zeros(len(rasterindices))+rasteri,'|k',\ markersize=20, markeredgewidth='2') # | is the marker ax.set_ylim(-0.5,rasteri+0.5) dirtItable = ax.table(cellText=data, cellColours=cellcolours, rowLoc='right',\ rowLabels=rowlabels, colLabels=collabels, colLoc='center', loc='bottom') table_props = dirtItable.properties() table_cells = table_props['child_artists'] for cell in table_cells: cell.set_height(1.35) cell.set_fontsize(18) ax.set_title(titlestr,fontsize=14) ## tight_layout() doesn't seem to work with table #fig.tight_layout() #fig.savefig(figfilename,dpi=300) def plot_dirtIrates(spiketrains1,spiketrains2,delay1,delay2,\ dirtIcutoff,filename='none.svg'): collabels = ['Order 1','2','3'] rowlabels = ['Delay 0','1','2','3'] dirtIs = [] cellcolours = [] for measure_delay in range(4): print "delay =",measure_delay dirtIorders = [] cellcoloursorders = [] for i in range(1,4): print "markovorder =",i causal_dirtI = calc_dirtinforate(spiketrains1,spiketrains2,\ i,i,measure_delay,measure_delay) oppcausal_dirtI = calc_dirtinforate(spiketrains2,spiketrains1,\ i,i,measure_delay,measure_delay) dirtIstr = ' causal = {:1.3f}\nopp causal = {:1.3f}'\ .format(causal_dirtI,oppcausal_dirtI) print dirtIstr dirtIorders.append(dirtIstr) if causal_dirtI>dirtIcutoff: cellcoloursorders.append('r') else: cellcoloursorders.append('w') dirtIs.append(dirtIorders) cellcolours.append(cellcoloursorders) titlestr = "Copy IID spike train 1 to 2 with delays = "\ +str(delay1)+" & "+str(delay2)+\ "\n Measure with markov order & common delay as below"+\ "\n causal is 1->2, opp causal is 2->1" plot_table([spiketrains1[0][0:100],spiketrains2[0][0:100]],\ rowlabels,collabels,dirtIs,cellcolours,titlestr,filename) def copycat_dirtIrate(delay=0): """Test directed information: trains1 has full causal exc effect on trains2. Copies spiketrain1 to 2 with causal delay i.e. delay+1. spiketrains and spiketrains1 is IID, p(0)=0.5.""" spiketrains = [] spiketrains1 = [] spiketrains2 = [] for j in range(10): ## need delay+1 num of copies of 1, to access mc1[-delay-1] below ## maintain same length of spike trains, hence common start length mc = [1]*(delay+1) mc1 = [1]*(delay+1) mc2 = [1]*(delay+1) for i in range(10000): ## mc2 fully depends on mc1 with delay if mc1[-delay-1]: mc2.append(1) else: mc2.append(0) ## mc1 is IID p = uniform() if p<0.5: mc1.append(1) else: mc1.append(0) ## mc is IID p = uniform() if p<0.5: mc.append(1) else: mc.append(0) spiketrains.append(mc) spiketrains1.append(mc1) spiketrains2.append(mc2) collabels = ['Order 1','2','3'] rowlabels = ['Delay 0','1','2','3'] dirtIs = [] cellcolours = [] for measure_delay in range(4): print "delay =",measure_delay dirtIorders = [] cellcoloursorders = [] for i in range(1,4): print "order =",i causal_dirtI = calc_dirtinforate(spiketrains1,spiketrains2,\ i,i,measure_delay,measure_delay) oppcausal_dirtI = calc_dirtinforate(spiketrains2,spiketrains1,\ i,i,measure_delay,measure_delay) acausal_dirtI = calc_dirtinforate(spiketrains,spiketrains2,\ i,i,measure_delay,measure_delay) dirtIstr = ' causal = {:1.3f}\nopp causal = {:1.3f}\n acausal = {:1.3f}'\ .format(causal_dirtI,oppcausal_dirtI,acausal_dirtI) print dirtIstr dirtIorders.append(dirtIstr) if causal_dirtI>0.9: cellcoloursorders.append('r') else: cellcoloursorders.append('w') dirtIs.append(dirtIorders) cellcolours.append(cellcoloursorders) titlestr = "Copy IID spike train 1 to 2 with causal delay = "+str(delay)+\ "\n Measure with markov order & common delay as below"+\ "\n causal is 1->2, opp causal is 2->1, acausal is IID->2" plot_table([spiketrains1[0][0:100],spiketrains2[0][0:100]],\ rowlabels,collabels,dirtIs,cellcolours,titlestr,'copycat_mydefn.svg') def partialcopycat_dirtIrate(delay1=0,delay2=0,inh=True): """Test directed information: trains1 has partial causal exc effect on trains2. Copies partially spiketrain1 to 2, with 1 delayed by delay1 and 2 delayed by delay2, (these are causal delays, hence delay1+1 and delay2+1) and flipping if inh.""" spiketrains1 = [] spiketrains2 = [] if inh: one=0 zer=1 else: one=1 zer=0 for j in range(10): ## need delay+1 num of copies of 1, to access mc1[-delay1-1] below ## maintain same length of spike trains, hence common start length delay = max(delay1,delay2) mc1 = [1]*(delay+1) mc2 = [1]*(delay+1) for i in range(10000): ## mc2 depends on mc1 and mc2, if mc1[-delay1-1] and mc2[-delay2-1]: # both are 1 mc2.append(one) elif mc1[-delay1-1] or mc2[-delay2-1]: # elif either is 1 p = uniform(0,1) if p<0.75: mc2.append(one) else: mc2.append(zer) else: # both are zero mc2.append(zer) ## mc1 is IID p = uniform() if p<0.5: mc1.append(1) else: mc1.append(0) spiketrains1.append(mc1) spiketrains2.append(mc2) plot_dirtIrates(spiketrains1,spiketrains2,delay1,delay2,0.1,'partialcopycat_mydefn.svg') def bicausal_dirtIrate(delay1=0,delay2=0,inh=True): """Test directed information: bidirectional exc/inh effect.""" spiketrains1 = [] spiketrains2 = [] if inh: one=0 zer=1 else: one=1 zer=0 for j in range(10): ## need delay+1 num of copies of 1, to access mc1[-delay1-1] below ## maintain same length of spike trains, hence common start length delay = max(delay1,delay2) mc1 = [1]*(delay+1) mc2 = [1]*(delay+1) for i in range(10000): ## mc1 and mc2 depend on mc1 and mc2, but asymmetrically ## delay is on both here! if mc1[-delay1-1] and mc2[-delay2-1]: # both are 1 p = uniform(0,1) if p<0.85: mc2.append(one) else: mc2.append(zer) p = uniform(0,1) if p<0.65: mc1.append(one) else: mc1.append(zer) elif mc1[-delay1-1] or mc2[-delay2-1]: # elif either is 1 p = uniform(0,1) if p<0.6: mc2.append(one) else: mc2.append(zer) p = uniform(0,1) if p<0.5: mc1.append(one) else: mc1.append(zer) else: # both are zero p = uniform(0,1) if p<0.8: mc2.append(zer) else: mc2.append(one) p = uniform(0,1) if p<0.6: mc1.append(zer) else: mc1.append(one) spiketrains1.append(mc1) spiketrains2.append(mc2) plot_dirtIrates(spiketrains1,spiketrains2,delay1,delay2,0.01,'bicausal_mydefn.svg') if __name__ == "__main__": #copycat_dirtIrate() partialcopycat_dirtIrate(delay1=1,delay2=2,inh=False) partialcopycat_dirtIrate(delay1=1,delay2=2,inh=True) #bicausal_dirtIrate(delay1=1,delay2=2,inh=True) show()