#!/usr/bin/python -i
import matplotlib
matplotlib.rcParams['pdf.fonttype'] = 42
mfs=22
matplotlib.rc('xtick', labelsize=mfs)
matplotlib.rc('ytick', labelsize=mfs)
matplotlib.rc('axes', labelsize=mfs)
from pylab import *
import numpy as np
import re
import os
np.set_printoptions( threshold=999999999999999)
if (len(sys.argv) < 2):
datadir = os.popen("ls -t ./data | head -1").read().strip();
else:
datadir = sys.argv[1]
datadir = datadir.replace('./data/', '');
datadir = datadir.replace('data/', '');
IMODE=1
if (len(sys.argv) > 2):
IMODE =0
ACTIVE_CUTOFF = 10
#N400B.20.I10.i6.P1.p1.T30.S1980.s
p = re.match(r'N(\d+).B(\d+).I(\d+).i(\d+).P(\d+).p(\d+).T(\d+).S(\d+).(\w+)', datadir)
print datadir
#os.system("cp constructs.cpp ./data/%s/"%( datadir))
STIMDURATION = 1800 + 100 + 100 # Taken from constructs.cpp
NTOTAL = int(p.group(1)) #inh + pyr neurons
NBRANCHES = int(p.group(2))
NINPUTS = int(p.group(3))
NPERINPUT = int(p.group(4))
NPATTERNS = int(p.group(5))
INTERSTIM = int(p.group(7))
SUFFIX = p.group(9)
NPYR = int(0.8*NTOTAL)
PYR_IDS = range(0 , NPYR)
IN_IDS = range(NPYR, NTOTAL)
def getoverlaps(patarr, threshold):
ta = patarr > threshold
cors = np.zeros((ta.shape[0], ta.shape[0]))
s = ta[i,:].sum() + ta[j,:].sum();
for i in range(ta.shape[0]):
for j in range(ta.shape[0]):
cors[i][j] = 2.0*float((ta[i,:] & ta[j,:]).sum()) / float(ta[i,:].sum() + ta[j,:].sum())
return cors
spcoeff = None
def getrates(period, patternid=0, rates=1):
if (period == 'pre'):
tstart = 0
tend = tstart + NPATTERNS*STIMDURATION
elif (period == 'during'):
tstart = NPATTERNS*STIMDURATION*0
tend = tstart + NPATTERNS*STIMDURATION
elif (period == 'post'):
tstart = NPATTERNS*STIMDURATION*1
tend = tstart + NPATTERNS*STIMDURATION
elif (period=='single'):
tstart = NPATTERNS*STIMDURATION*3
tend = tstart + NINPUTS*STIMDURATION
if (patternid>=0 and period != 'single'):
tstart = tstart + STIMDURATION*patternid
tend = tstart + STIMDURATION
if (rates==1): return averaged[:, tstart:tend]
else: return raster[:, tstart:tend]
def loadtimings(filename, tduration):
ff = open(filename, 'r')
fdata = ff.readlines()
sx = len(fdata)
sy = tduration;
raster = np.zeros( (sx, sy) );
nid=0
for l in fdata:
ar = np.fromstring(l, sep=' ' , dtype=int)
raster[nid, ar] = 1
raster[nid,0] =0 # XXX bug
nid += 1
return raster
def plotdistancesspikes():
global spcoeff
global spcounts
ry = STIMDURATION*NINPUTS + (NINPUTS*(NINPUTS-1))*STIMDURATION
NMEMS = NINPUTS
raster = loadtimings("./data/%s/spikes.dat"%(datadir), ry)
windowsize=5.
WND = np.repeat(1.0, windowsize) / windowsize
#averaged = np.zeros(raster.shape)
counts = np.zeros((NMEMS,NMEMS))
tot =0
for i in range(NMEMS):
for j in range(NMEMS):
if (i < j):
tstart = NPATTERNS*STIMDURATION + STIMDURATION*(tot)
tend = tstart +STIMDURATION
r = raster[:, tstart:tend]
sums = np.sum(r)
tot += 1
counts[i][j] = sums
figure()
imshow(raster , interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
#savefig("./data/%s/raster.png"%(datadir))
figure()
imshow(counts , interpolation='nearest', aspect='auto',cmap='jet')
#colorbar()
#tight_layout()
savefig("./data/%s/distances.png"%(datadir))
return
def plotspikes():
global spcoeff
global spcounts
ry = STIMDURATION*NPATTERNS*2 + 0*NINPUTS*STIMDURATION
raster = loadtimings("./data/%s/spikes.dat"%(datadir), ry)
windowsize=50.
WND = np.repeat(1.0, windowsize) / windowsize
averaged = np.zeros(raster.shape)
for nid in range(raster.shape[0]):
av = np.convolve(raster[nid, :], WND)
averaged[nid, : ] = av[windowsize/2:-windowsize/2+1]
figure()
imshow(averaged *(500.), interpolation='nearest', aspect='auto', cmap='jet')
colorbar()
#tight_layout()
#savefig("./data/%s/raster.png"%(datadir));
spcounts = np.zeros( (len(PYR_IDS), NPATTERNS) )
spcoeff = np.zeros( (len(PYR_IDS), NPATTERNS) )
cors = np.zeros( (NPATTERNS,len(PYR_IDS)) )
patresp = np.zeros((len(PYR_IDS), NPATTERNS))
for n in range(NPATTERNS):
tstart = NPATTERNS*STIMDURATION*1 + STIMDURATION*n
tend = tstart +STIMDURATION
r = raster[:, tstart:tend]
av = averaged[:, tstart:tend]
sums = np.sum(r,1)
corrs = np.nan_to_num(np.corrcoef(av)) # discard NaN values
for nid in PYR_IDS:
#s =0;
#s += cc[nid, NTOTAL + n*NPERINPUT+inpid];
#cors[n, nid] = s / len(PYR_IDS)
stt = sums[nid]
spcounts[nid, n] = stt
#for inpnid in range(NPERINPUT):
#print "%d %d %d\n"%( NTOTAL, NPERINPUT, NTOTAL + n*NPERINPUT +inpnid)
#spcoeff[nid, n] += corrs[nid, NTOTAL + n*NPERINPUT +inpnid]/NPERINPUT
spcounts /= 2.0
fircorr = np.zeros((len(PYR_IDS), NPATTERNS))
if (0):
for nid in [1, 2, 3]:
nn = np.zeros((NPATTERNS, STIMDURATION))
for sp in range(NPATTERNS):
tstart = NPATTERNS*STIMDURATION*1 + STIMDURATION*sp
tend = tstart +STIMDURATION
nn[sp, :] = raster[nid, tstart:tend]
gg = np.nan_to_num(np.corrcoef(nn))
figure()
imshow(gg , interpolation='nearest', aspect='auto', cmap='jet')
colorbar()
return;
spcountspre = np.zeros( (len(PYR_IDS), NPATTERNS) )
for n in range(NPATTERNS):
tstart = NPATTERNS*STIMDURATION*0 + STIMDURATION*n
tend = tstart +STIMDURATION
r = raster[:, tstart:tend]
#av = averaged[:, tstart:tend]
sums = np.sum(r,1)
for nid in PYR_IDS:
#s =0;
#s += cc[nid, NTOTAL + n*NPERINPUT+inpid];
#cors[n, nid] = s / len(PYR_IDS)
spcountspre[nid, n] = sums[nid]
spcoeff = spcoeff.transpose()
"""
figure()
title("Average firing rate per pattern recall")
ylabel("number of spikes PRE")
xlabel("#pattern")
imshow(spcountspre, interpolation='nearest', aspect='auto',cmap='jet')
np.save("./data/%s/spikespre"%(datadir), spcountspre);
colorbar()
"""
figure()
title("Average firing rate (Hz) during pattern recall")
xlabel("Memory No.")
ylabel("Pyramidal Neuron No.")
imshow(spcounts , interpolation='nearest', aspect='auto',cmap='jet')
np.save("./data/%s/spikespost"%(datadir), spcounts);
savefig("./data/%s/spikespost.pdf"%(datadir))
colorbar()
figure()
if (NPATTERNS>1):
pp = np.corrcoef((spcounts).transpose())
xlabel("Memory #")
ylabel("Memory #")
title("Firing rates Correlation")
imshow(pp.transpose() , interpolation='nearest', aspect='auto',cmap='jet', vmin=-0.2, vmax=1.)
colorbar()
np.save("./data/%s/spcountscorr"%(datadir), pp.transpose());
xticks(arange(NPATTERNS))
figure()
activeonly = np.corrcoef((spcounts>ACTIVE_CUTOFF).transpose())
xlabel("Memory #")
ylabel("Memory #")
title("Population overlap")
imshow(activeonly , interpolation='nearest', aspect='auto',cmap='jet', vmin=-0.2, vmax=1.)
colorbar()
#np.save("./data/%s/spcountscorr"%(datadir), pp.transpose());
xticks(arange(NPATTERNS))
figure()
ylabel("% neurons")
xlabel("# mem")
title("Percentage of neurons > %d spikes"%(ACTIVE_CUTOFF))
spc = np.zeros(NPATTERNS)
for n in range(NPATTERNS):
spc[n] = (spcounts[:,n] > ACTIVE_CUTOFF).sum()
vals = 100*spc/NPYR
print "VV"
print vals
bar( arange(NPATTERNS), vals )
#ylim((0,60))
np.save("./data/%s/above_%d_dt"%(datadir, ACTIVE_CUTOFF), vals)
savefig("./data/o_%s_%d.pdf"%( SUFFIX, INTERSTIM));
xticks(arange(NPATTERNS))
figure()
title("avg firing rate per pattern (interval=%d min)"%(INTERSTIM))
ylabel("avg firing rate (Hz)")
xlabel("# mem")
nbins = 20
pp = zeros( (NPATTERNS, nbins));
perr = np.std(spcounts.transpose(), 1)
vv = np.sum(spcounts.transpose(), 1)/(NPYR*1.8) ;
bar(arange(NPATTERNS), vv)
ylim((0,50))
np.save("./data/%s/avgrate"%(datadir), vv)
xticks(arange(NPATTERNS))
# hist = np.histogram(st[p, :], bins=nbins)
# pp[p, :] = hist[0]
#imshow(pp , interpolation='nearest', aspect='auto',cmap='jet')
#colorbar()
"""
figure()
title("spikes overlap per mem")
ylabel("#mem")
xlabel("#mem")
if 0:
pp = getoverlaps(spcounts.transpose(), 20)
bar([0], [pp[0,1]])
ylim((0.,1.))
else:
if (NPATTERNS > 1):
#pp = np.corrcoef((spcounts).transpose())
pp = getoverlaps(spcounts.transpose(), ACTIVE_CUTOFF)
pp = pp.transpose();
imshow(pp , interpolation='nearest', aspect='auto',cmap='jet', vmin=0, vmax=1.)
np.save("./data/%s/spoverlap"%(datadir), pp)
colorbar()
xticks(arange(NPATTERNS))
yticks(arange(NPATTERNS))
#tight_layout()
savefig("./data/%s/firing.png"%(datadir));
"""
#figure()
#title("Dist. of frequences during recall")
#hist(spcounts.flatten()/1.800, bins=100)
"""
figure()
title("#patterns / neuron > 30 spikes")
xlabel("#patterns")
ylabel("#neurons")
n, b, l = hist( (spcounts>ACTIVE_CUTOFF).sum(axis=1), bins=12, range =(0,12))
print n
print b
np.save("./data/%s/pat_neuron"%(datadir), n)
np.save("./data/%s/pat_neuron_dt"%(datadir), (spcounts>ACTIVE_CUTOFF).sum(axis=1))
figure(figsize=(15,15))
subplot(221)
title("#firing correlations coefficients")
pp = np.zeros( (len(PYR_IDS), NPATTERNS) )
ylabel("#pattern");
xlabel("#neuron");
imshow(spcoeff.transpose(), interpolation='nearest', aspect='auto',cmap='jet', vmin=-.2, vmax=1)
colorbar()
subplot(222)
title("#corr. coeff.")
ylabel("#pattern");
xlabel("#pattern");
if (NPATTERNS>1):
imshow(np.corrcoef(spcoeff), interpolation='nearest', aspect='auto',cmap='jet', vmin=-.2, vmax=1)
colorbar()
subplot(223)
title("recruited neurons")
ylabel("% recruited pyr. neurons");
xlabel("#pattern");
n1 = []
n2 = []
n3 = []
for i in range(NPATTERNS):
n1.append((spcoeff[i,:]>0.2).sum()/float(NPYR))
n2.append((spcoeff[i,:]>0.3).sum()/float(NPYR))
n3.append((spcoeff[i,:]>0.5).sum()/float(NPYR))
#plot(n1, label='0.2', marker='o')
plot(n2, label='0.3', marker='o')
#plot(n3, label='0.6', marker='o')
legend()
subplot(224)
hst = []
thresh = 0.3
for i in range(NPYR):
hst.append((spcoeff[:,i]>thresh).sum())
title('# of patterns per neuron (thresh=%f)'%(thresh))
hist(hst, bins=20)
xlabel('# of patterns')
ylabel('# of neurons')
#tight_layout()
savefig("./data/%s/firingcorr.png"%(datadir), dpi=200);
"""
def plotsyns():
global brhits
ss = np.loadtxt("./data/%s/synstate.dat"%(datadir))
brweights= np.zeros((NINPUTS, NPYR*NBRANCHES))
nrnweights= np.zeros((NINPUTS, NPYR))
brhits = np.zeros((NPYR*NBRANCHES))
brpatterns = np.zeros((NBRANCHES*NPYR, NINPUTS))
learnedhits= np.zeros(NPYR*NBRANCHES)
unlearnedhits= np.zeros(NPYR*NBRANCHES)
mem_branches = np.zeros((NPATTERNS, NPYR*NBRANCHES));
for l in ss:
bid = l[1]
nrnid = l[2]
inputid = l[4]
srcnrnid=l[3]
w = l[6]
if (nrnid < NPYR and inputid >= 0):
if (w > 0.999):
brhits[bid] += 1
brpatterns[bid, inputid] = 1
if (spcounts[nrnid, inputid] > ACTIVE_CUTOFF):
learnedhits[bid] += 1
else:
unlearnedhits[bid] += 1
if (w > 1.0):
mem_branches[inputid, bid] += 1
brweights[inputid, bid] += w
nrnweights[inputid, nrnid] += w
wmax = nrnweights.max()
brmax = brweights.max()
nrnweights /= wmax
brweights /= brmax
#figure()
#title("Number of potentiated synapses per branch");
#hist(mem_branches[0, :], range=(1,19), bins=19)
#m = np.average(mem_branches, axis=1)
#st = np.std(mem_branches, axis=1)
#figure()
#bar(range(0,len(m)), m, yerr=st)
#figure()
#title("Potentiated synapses per branch")
#imshow(mem_branches.transpose(), interpolation='nearest', aspect='auto',cmap='jet')
#colorbar()
if (NPATTERNS>1):
figure()
title("Correlation of Potentiated synapses per neuron")
pp = np.corrcoef(nrnweights)
imshow(pp , interpolation='nearest', aspect='auto',cmap='jet');
colorbar()
figure()
title("Correlation of Potentiated synapses per branch")
pp = np.corrcoef(mem_branches)
imshow(pp , interpolation='nearest', aspect='auto',cmap='jet', vmin=-0.2, vmax=1.)
colorbar()
def plotdspikes():
ry = STIMDURATION*NPATTERNS*2
filename = ("./data/%s/branchspikes.dat"%(datadir))
ff = open(filename, 'r')
fdata = ff.readlines()
sx = len(fdata)
branch_dspikes = np.zeros((NPATTERNS, NBRANCHES*NPYR))
bid =0
for l in fdata:
ar = np.fromstring(l, sep=' ' , dtype=int)
for k in range(NPATTERNS):
print bid
print (len(ar))
branch_dspikes[k, bid] += ar[NPATTERNS + k]
bid += 1
if (bid >= NPYR*NBRANCHES): break
branch_dspikes = branch_dspikes > 1
figure()
xlabel("#pattern")
ylabel("#branch")
title("dend. spikes per memory")
imshow(branch_dspikes.transpose() , interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
figure()
title("corr. coeff ")
ylabel("#pattern")
xlabel("#pattern")
pp = np.corrcoef(branch_dspikes)
imshow(pp , interpolation='nearest', aspect='auto',cmap='jet', vmin=-0.2, vmax=1.)
colorbar()
#savefig("./data/%s/dspikes.png"%(datadir));
"""
figure(figsize=(15,6))
subplot(121)
ylabel("#pattern")
xlabel("#pyr neyron")
title("dend. spike rate per neuron")
imshow(nrnspcounts , interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
subplot(122)
title("corr. coeff ")
ylabel("#pattern")
xlabel("#pattern")
pp = np.corrcoef(nrnspcounts)
imshow(pp , interpolation='nearest', aspect='auto',cmap='jet', vmin=-.2, vmax=1.)
colorbar()
#tight_layout()
savefig("./data/%s/nrndspikes.png"%(datadir));
"""
def plotproteins():
""" ss = np.loadtxt("./data/%s/branchproteins.dat"%(datadir))
figure(figsize=(15,6))
ylabel("#branch")
xlabel("protein")
title("protein levels per branch ")
imshow( ss, interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
savefig("./data/%s/proteins.png"%(datadir));
"""
ss = np.loadtxt("./data/%s/synstate.dat"%(datadir))
brweights= np.zeros((NPATTERNS, NPYR*NBRANCHES))
nrnweights= np.zeros((NPATTERNS, NPYR))
synToPattern = []
for l in ss:
bid = l[1]
nrnid = l[2]
patternid = l[4]
srcnrnid=l[3]
w = l[6]
synToPattern.append(patternid)
synsort = []
for n in range(NPATTERNS):
k =0
for l in synToPattern:
if (n == synToPattern[int(l)]):
synsort.append( k)
k+= 1
ss = np.loadtxt("./data/%s/tags.dat"%(datadir))
print len(synsort)
figure(figsize=(15,6))
ylabel("#syn")
xlabel("tag")
title("tags per synapse ")
ab = ss[ss[:,0].argsort()]
imshow(ab[:,1:], interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
savefig("./data/%s/tags.png"%(datadir));
ss = np.loadtxt("./data/%s/nrnprotein.dat"%(datadir))
figure(figsize=(15,6))
ylabel("#nrn")
xlabel("protein")
title("proteins per neuron")
imshow( ss, interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
savefig("./data/%s/nrnprotein.png"%(datadir));
ss = np.loadtxt("./data/%s/weighthistory.dat"%(datadir))
figure(figsize=(15,6))
ylabel("#synapse")
xlabel("weight")
title("weight per synapse")
ab = ss[ss[:,0].argsort()]
imshow( ab[:,1:], interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
savefig("./data/%s/nrnweights.png"%(datadir));
def plotbstrengths():
ss = np.loadtxt("./data/%s/branchstrengths.dat"%(datadir))
figure(figsize=(15,6))
#subplot(121)
ylabel("#branch")
xlabel("branch strength")
title("branch strengths ")
imshow( ss, interpolation='nearest', aspect='auto',cmap='jet')
colorbar()
#subplot(122)
#title("corr. coeff ")
#ylabel("#branch")
#xlabel("#branch")
#pp = np.corrcoef(ss)
#imshow( pp , interpolation='nearest', aspect='auto',cmap='jet')
#colorbar()
tight_layout()
savefig("./data/%s/bstrengths.png"%(datadir));
def ansims():
for sim in range(1,6):
for i in range(0,NRUNS):
ddir = "N%d.B%d.I%d.i%d.P%d.p%d.T%d.S%d.sn_sim%d"%(NTOTAL, NBRANCHES, NINPUTS, NPERINPUT, NPATTERNS, NPERPATTERN, INTERSTIM, RSEED+i, sim)
dd = np.load("./data/%s/spcountscorr.npy"%(ddir))
p += dd
if (0):
NRUNS = 10
p = np.zeros((NPATTERNS, NPATTERNS))
for i in range(0,NRUNS):
ddir = "N%d.B%d.I%d.i%d.P%d.p%d.T%d.S%d.%s"%(NTOTAL, NBRANCHES, NINPUTS, NPERINPUT, NPATTERNS, NPERPATTERN, INTERSTIM, RSEED+i, suff)
dd = np.load("./data/%s/spcountscorr.npy"%(ddir))
p += dd
p /= NRUNS
figure()
imshow(p , interpolation='nearest', aspect='auto',cmap='jet')
else:
if (IMODE): ion()
plotspikes()
plotsyns()
#plotdspikes()
#plotproteins()
#plotbstrengths()
if (IMODE): show()
