From a5411f6ad22549861550c3d1f1ae4b6704d61065 Mon Sep 17 00:00:00 2001 From: Christos Kotsalos Date: Wed, 6 Jul 2022 13:54:00 +0200 Subject: [PATCH 1/3] update model for python3 and modeldb-ci support --- analysis.py | 235 +++++++------ arm.py | 145 ++++---- armGraphs.py | 613 ++++++++++++++++---------------- arminterface.py | 916 ++++++++++++++++++++++++------------------------ dummyArm.py | 138 ++++---- error.py | 92 +++-- evol_islands.py | 124 +++---- izhi.py | 8 +- main.py | 16 +- mosinit.py | 2 +- network.py | 298 ++++++++-------- nsloc.py | 2 +- server.py | 652 +++++++++++++++++----------------- shared.py | 61 ++-- stimuli.py | 14 +- 15 files changed, 1646 insertions(+), 1670 deletions(-) diff --git a/analysis.py b/analysis.py index 51ceb24..b101e5a 100644 --- a/analysis.py +++ b/analysis.py @@ -5,7 +5,7 @@ Version: 2014July9 """ -from pylab import pcolor, nonzero, mean, histogram, arange, bar, vstack,scatter, figure, hold, isscalar, gca, unique, subplot, axes, shape, imshow, colorbar, plot, xlabel, ylabel, title, xlim, ylim, clim, show, zeros, legend, savefig, cm, specgram, get_cmap, psd +from pylab import pcolor, nonzero, mean, histogram, arange, bar, vstack,scatter, figure, isscalar, gca, unique, subplot, axes, shape, imshow, colorbar, plot, xlabel, ylabel, title, xlim, ylim, clim, show, zeros, legend, savefig, cm, specgram, get_cmap, psd from scipy.io import loadmat from scipy import loadtxt, size, array, linspace, ceil from datetime import datetime @@ -20,33 +20,33 @@ ## Create colormap def bicolormap(gap=0.1,mingreen=0.2,redbluemix=0.5,epsilon=0.01): - from matplotlib.colors import LinearSegmentedColormap as makecolormap - - mng=mingreen; # Minimum amount of green to add into the colors - mix=redbluemix; # How much red to mix with the blue an vice versa - eps=epsilon; # How much of the center of the colormap to make gray - omg=1-gap # omg = one minus gap - - cdict = {'red': ((0.00000, 0.0, 0.0), - (0.5-eps, mix, omg), - (0.50000, omg, omg), - (0.5+eps, omg, 1.0), - (1.00000, 1.0, 1.0)), - - 'green': ((0.00000, mng, mng), - (0.5-eps, omg, omg), - (0.50000, omg, omg), - (0.5+eps, omg, omg), - (1.00000, mng, mng)), - - 'blue': ((0.00000, 1.0, 1.0), - (0.5-eps, 1.0, omg), - (0.50000, omg, omg), - (0.5+eps, omg, mix), - (1.00000, 0.0, 0.0))} - cmap = makecolormap('bicolormap',cdict,256) - - return cmap + from matplotlib.colors import LinearSegmentedColormap as makecolormap + + mng=mingreen; # Minimum amount of green to add into the colors + mix=redbluemix; # How much red to mix with the blue an vice versa + eps=epsilon; # How much of the center of the colormap to make gray + omg=1-gap # omg = one minus gap + + cdict = {'red': ((0.00000, 0.0, 0.0), + (0.5-eps, mix, omg), + (0.50000, omg, omg), + (0.5+eps, omg, 1.0), + (1.00000, 1.0, 1.0)), + + 'green': ((0.00000, mng, mng), + (0.5-eps, omg, omg), + (0.50000, omg, omg), + (0.5+eps, omg, omg), + (1.00000, mng, mng)), + + 'blue': ((0.00000, 1.0, 1.0), + (0.5-eps, 1.0, omg), + (0.50000, omg, omg), + (0.5+eps, omg, mix), + (1.00000, 0.0, 0.0))} + cmap = makecolormap('bicolormap',cdict,256) + + return cmap ## Raster plot def plotraster(filename=None): # allspiketimes, allspikecells, EorI, ncells, connspercell, backgroundweight, firingrate, duration): # Define a function for plotting a raster @@ -54,14 +54,14 @@ def plotraster(filename=None): # allspiketimes, allspikecells, EorI, ncells, con EorIcolors = array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise cellcolors = EorIcolors[array(s.EorI)[array(s.allspikecells,dtype=int)]] # Set each cell to be either orange or turquoise figure() # Open a new figure - scatter(s.allspiketimes,s.allspikecells,10,cellcolors,linewidths=0.5,marker='|') # Create raster + scatter(s.allspiketimes,s.allspikecells,10,cellcolors,linewidths=0.5,marker='|') # Create raster xlabel('Time (ms)') ylabel('Cell ID') title('cells=%i syns/cell=%0.1f noise=%0.1f rate=%0.1f Hz' % (s.ncells,s.connspercell,s.backgroundweight[0],s.firingrate),fontsize=12) xlim(0,s.duration) ylim(0,s.ncells) plottime = time()-plotstart # See how long it took - print(' Done; time = %0.1f s' % plottime) + print((' Done; time = %0.1f s' % plottime)) if filename: savefig(filename) #show() @@ -80,23 +80,23 @@ def plotPETH(): ylabel('Spikes/bin') ylim(0,s.scale*binsize*2) h=axes() - h.set_xticks(range(0,len(binedges),len(binedges)/10 )) + h.set_xticks(list(range(0,len(binedges),len(binedges)/10))) h.set_xticklabels(binedges[0:-1:len(binedges)/10].astype(int)) legend(s.popnames) ## Plot power spectra density def plotpsd(): - colorspsd=array([[0.42,0.67,0.84],[0.42,0.83,0.59],[0.90,0.76,0.00],[0.90,0.32,0.00],[0.34,0.67,0.67],[0.42,0.82,0.83],[0.90,0.59,0.00],[0.33,0.67,0.47],[1.00,0.85,0.00],[0.71,0.82,0.41],[0.57,0.67,0.33],[1.00,0.38,0.60],[0.5,0.2,0.0],[0.0,0.2,0.5]]) + colorspsd=array([[0.42,0.67,0.84],[0.42,0.83,0.59],[0.90,0.76,0.00],[0.90,0.32,0.00],[0.34,0.67,0.67],[0.42,0.82,0.83],[0.90,0.59,0.00],[0.33,0.67,0.47],[1.00,0.85,0.00],[0.71,0.82,0.41],[0.57,0.67,0.33],[1.00,0.38,0.60],[0.5,0.2,0.0],[0.0,0.2,0.5]]) - lfpv=[[] for c in range(len(s.lfppops))] + lfpv=[[] for c in range(len(s.lfppops))] # Get last modified .mat file if no input and plot for c in range(len(s.lfppops)): - lfpv[c] = s.lfps[:,c] + lfpv[c] = s.lfps[:,c] lfptot = sum(lfpv) - + # plot pops separately plotPops = 0 - if plotPops: + if plotPops: figure() # Open a new figure for p in range(len(s.lfppops)): psd(lfpv[p],Fs=200, linewidth= 2,color=colorspsd[p]) @@ -137,14 +137,14 @@ def plotconn(): # Plot grid lines - hold(True) + #hold(True) for pop in range(s.npops): plot(array([0,s.npops])-0.5,array([pop,pop])-0.5,'-',c=(0.7,0.7,0.7)) plot(array([pop,pop])-0.5,array([0,s.npops])-0.5,'-',c=(0.7,0.7,0.7)) # Make pretty - h.set_xticks(range(s.npops)) - h.set_yticks(range(s.npops)) + h.set_xticks(list(range(s.npops))) + h.set_yticks(list(range(s.npops))) h.set_xticklabels(s.popnames) h.set_yticklabels(s.popnames) h.xaxis.set_ticks_position('top') @@ -158,42 +158,42 @@ def plotconn(): ## Plot weight changes def plotweightchanges(filename=None): if s.usestdp: - # create plot - figh = figure(figsize=(1.2*8,1.2*6)) - figh.subplots_adjust(left=0.02) # Less space on left - figh.subplots_adjust(right=0.98) # Less space on right - figh.subplots_adjust(top=0.96) # Less space on bottom - figh.subplots_adjust(bottom=0.02) # Less space on bottom - figh.subplots_adjust(wspace=0) # More space between - figh.subplots_adjust(hspace=0) # More space between - h = axes() - - # create data matrix + # create plot + figh = figure(figsize=(1.2*8,1.2*6)) + figh.subplots_adjust(left=0.02) # Less space on left + figh.subplots_adjust(right=0.98) # Less space on right + figh.subplots_adjust(top=0.96) # Less space on bottom + figh.subplots_adjust(bottom=0.02) # Less space on bottom + figh.subplots_adjust(wspace=0) # More space between + figh.subplots_adjust(hspace=0) # More space between + h = axes() + + # create data matrix wcs = [x[-1][-1] for x in s.allweightchanges] # absolute final weight - wcs = [x[-1][-1]-x[0][-1] for x in s.allweightchanges] # absolute weight change - pre,post,recep = zip(*[(x[0],x[1],x[2]) for x in s.allstdpconndata]) - ncells = int(max(max(pre),max(post))+1) - wcmat = zeros([ncells, ncells]) + wcs = [x[-1][-1]-x[0][-1] for x in s.allweightchanges] # absolute weight change + pre,post,recep = list(zip(*[(x[0],x[1],x[2]) for x in s.allstdpconndata])) + ncells = int(max(max(pre),max(post))+1) + wcmat = zeros([ncells, ncells]) - for iwc,ipre,ipost,irecep in zip(wcs,pre,post,recep): + for iwc,ipre,ipost,irecep in zip(wcs,pre,post,recep): wcmat[int(ipre),int(ipost)] = iwc *(-1 if irecep>=2 else 1) - # plot - imshow(wcmat,interpolation='nearest',cmap=bicolormap(gap=0,mingreen=0.2,redbluemix=0.1,epsilon=0.01)) - xlabel('post-synaptic cell id') - ylabel('pre-synaptic cell id') - h.set_xticks(s.popGidStart) - h.set_yticks(s.popGidStart) - h.set_xticklabels(s.popnames) - h.set_yticklabels(s.popnames) - h.xaxis.set_ticks_position('top') - xlim(-0.5,ncells-0.5) - ylim(ncells-0.5,-0.5) - clim(-abs(wcmat).max(),abs(wcmat).max()) - colorbar() + # plot + imshow(wcmat,interpolation='nearest',cmap=bicolormap(gap=0,mingreen=0.2,redbluemix=0.1,epsilon=0.01)) + xlabel('post-synaptic cell id') + ylabel('pre-synaptic cell id') + h.set_xticks(s.popGidStart) + h.set_yticks(s.popGidStart) + h.set_xticklabels(s.popnames) + h.set_yticklabels(s.popnames) + h.xaxis.set_ticks_position('top') + xlim(-0.5,ncells-0.5) + ylim(ncells-0.5,-0.5) + clim(-abs(wcmat).max(),abs(wcmat).max()) + colorbar() if filename: savefig(filename) - #show() + #show() changeOverTime = 0 if changeOverTime: @@ -221,7 +221,7 @@ def plotweightchanges(filename=None): ylabel('Synaptic connection id') colorbar() #show() - + ## plot motor subpopulations connectivity changes @@ -232,17 +232,17 @@ def plotmotorpopchanges(): Ewpost = [] EwpreSum = [] EwpostSum = [] - if showInh: + if showInh: Iwpre = [] Iwpost = [] IwpreSum = [] - IwpostSum = [] + IwpostSum = [] for imus in range(len(s.motorCmdCellRange)): Ewpre.append([x[0][-1] for (icon,x) in enumerate(s.allweightchanges) if s.allstdpconndata[icon][1] in s.motorCmdCellRange[imus]]) Ewpost.append([x[-1][-1] for (icon,x) in enumerate(s.allweightchanges) if s.allstdpconndata[icon][1] in s.motorCmdCellRange[imus]]) EwpreSum.append(sum(Ewpre[imus])) EwpostSum.append(sum(Ewpost[imus])) - + if showInh: motorInhCellRange = s.motorCmdCellRange[imus] - s.popGidStart[s.EDSC] + s.popGidStart[s.IDSC] @@ -251,17 +251,17 @@ def plotmotorpopchanges(): IwpreSum.append(sum(Iwpre[imus])) IwpostSum.append(sum(Iwpost[imus])) - print '\ninitial E weights: ',EwpreSum - print 'final E weigths: ',EwpostSum - print 'absolute E difference: ',array(EwpostSum) - array(EwpreSum) - print 'relative E difference: ',(array(EwpostSum) - array(EwpreSum)) / array(EwpreSum) + print('\ninitial E weights: ',EwpreSum) + print('final E weigths: ',EwpostSum) + print('absolute E difference: ',array(EwpostSum) - array(EwpreSum)) + print('relative E difference: ',(array(EwpostSum) - array(EwpreSum)) / array(EwpreSum)) if showInh: - print '\ninitial I weights: ',IwpreSum - print 'final I weigths: ',IwpostSum - print 'absolute I difference: ',array(IwpostSum) - array(IwpreSum) - print 'relative I difference: ',(array(IwpostSum) - array(IwpreSum)) / array(IwpreSum) - + print('\ninitial I weights: ',IwpreSum) + print('final I weigths: ',IwpostSum) + print('absolute I difference: ',array(IwpostSum) - array(IwpreSum)) + print('relative I difference: ',(array(IwpostSum) - array(IwpreSum)) / array(IwpreSum)) + # plot figh = figure(figsize=(1.2*8,1.2*6)) @@ -300,7 +300,7 @@ def plotmotorpopchanges(): ax2.set_xticks(ind+width/2) ax2.set_xticklabels( ('shext','shflex','elext','elflex') ) ax2.grid() - + ## plot 3d architecture: def plot3darch(): @@ -318,9 +318,9 @@ def plot3darch(): ylocs=[3,2,1] zlocs=[0.1,0.5,1.2] ax.scatter(xlocs,ylocs, zlocs, s=10, c=zlocs, edgecolors='none',cmap = 'jet_r' , linewidths=0.0, alpha=1, marker='o') - azim = 40 + azim = 40 elev = 60 - ax.view_init(elev, azim) + ax.view_init(elev, azim) #xlim(min(s.xlocs),max(s.xlocs)) #ylim(min(s.ylocs),max(s.ylocs)) #ax.set_zlim(min(s.zlocs),max(s.zlocs)) @@ -337,7 +337,7 @@ def plot3darch(): def errorfill(x, y, yerr, lw=1, elinewidth=1, color=None, alpha_fill=0.2, ax=None): ax = ax if ax is not None else gca() if color is None: - color = ax._get_lines.color_cycle.next() + color = next(ax._get_lines.color_cycle) if isscalar(yerr) or len(yerr) == len(y): ymin = y - yerr ymax = y + yerr @@ -347,7 +347,7 @@ def errorfill(x, y, yerr, lw=1, elinewidth=1, color=None, alpha_fill=0.2, ax=Non ax.fill_between(x, ymax, ymin, color=color, lw= elinewidth, alpha=alpha_fill) #%% function to obtain unique list of lists -def uniqueList(seq): +def uniqueList(seq): seen = {} result = [] indices = [] @@ -358,8 +358,8 @@ def uniqueList(seq): result.append(item) indices.append(index) return result,indices - -#%% function to read data + +#%% function to read data def loadData(folder, islands, dataFrom): #%% Load data from files if islands > 1: @@ -367,40 +367,40 @@ def loadData(folder, islands, dataFrom): ind_cands_isl=[] ind_fits_isl=[] ind_cs_isl=[] - + stat_gens_isl=[] # statistics.csv for islands stat_worstfits_isl=[] stat_bestfits_isl=[] stat_avgfits_isl=[] stat_stdfits_isl=[] - + fits_sort_isl=[] #sorted data - gens_sort_isl=[] + gens_sort_isl=[] cands_sort_isl=[] params_sort_isl=[] - + for island in range(islands): ind_gens=[] # individuals data ind_cands=[] ind_fits=[] ind_cs=[] - + eval_gens=[] # error files for each evaluation eval_cands=[] eval_fits=[] eval_params=[] - - stat_gens=[] # statistics.csv + + stat_gens=[] # statistics.csv stat_worstfits=[] stat_bestfits=[] stat_avgfits=[] stat_stdfits=[] - + if islands > 0: folderFinal = folder+"_island_"+str(island) - else: + else: folderFinal = folder - + with open('../data/%s/individuals.csv'% (folderFinal)) as f: # read individuals.csv reader=csv.reader(f) for row in reader: @@ -409,7 +409,7 @@ def loadData(folder, islands, dataFrom): ind_fits.append(float(row[2])) cs = [float(row[i].replace("[","").replace("]","")) for i in range(3,len(row))] ind_cs.append(cs) - + with open('../data/%s/statistics.csv'% (folderFinal)) as f: # read statistics.csv reader=csv.reader(f) for row in reader: @@ -418,30 +418,30 @@ def loadData(folder, islands, dataFrom): stat_bestfits.append(float(row[3])) stat_avgfits.append(float(row[4])) stat_stdfits.append(float(row[6])) - + # unique generation number (sometimes repeated due to rerunning in hpc) stat_gens, stat_gens_indices = unique(stat_gens,1) # unique individuals - stat_worstfits, stat_bestfits, stat_avgfits, stat_stdfits = zip(*[[stat_worstfits[i], stat_bestfits[i], stat_avgfits[i], stat_stdfits[i]] for i in stat_gens_indices]) - - if dataFrom == 'fitness': + stat_worstfits, stat_bestfits, stat_avgfits, stat_stdfits = list(zip(*[[stat_worstfits[i], stat_bestfits[i], stat_avgfits[i], stat_stdfits[i]] for i in stat_gens_indices])) + + if dataFrom == 'fitness': for igen in range(max(ind_gens)): # read error files from evaluations for ican in range(max(ind_cands)): try: - f=open('../data/%s/gen_%d_cand_%d_error'%(folderFinal, igen,ican)); + f=open('../data/%s/gen_%d_cand_%d_error'%(folderFinal, igen,ican)); eval_fits.append(pickle.load(f)) - f=open('../data/%s/gen_%d_cand_%d_params'%(folderFinal, igen,ican)); + f=open('../data/%s/gen_%d_cand_%d_params'%(folderFinal, igen,ican)); eval_params.append(pickle.load(f)) eval_gens.append(igen) eval_cands.append(ican) except: - pass + pass #eval_fits.append(0.15) #eval_params.append([]) - - # find x corresponding to smallest error from function evaluations + + # find x corresponding to smallest error from function evaluations if dataFrom == 'fitness': #fits_sort, fits_sort_indices, fits_sort_origind = unique(eval_fits, True, True) - fits_sort_indices = sorted(range(len(eval_fits)), key=lambda k: eval_fits[k]) + fits_sort_indices = sorted(list(range(len(eval_fits))), key=lambda k: eval_fits[k]) fits_sort = [eval_fits[i] for i in fits_sort_indices] gens_sort = [eval_gens[i] for i in fits_sort_indices] cands_sort = [eval_cands[i] for i in fits_sort_indices] @@ -452,33 +452,32 @@ def loadData(folder, islands, dataFrom): fits_unique = [ind_fits[i] for i in unique_indices] gens_unique = [ind_gens[i] for i in unique_indices] cands_unique = [ind_cands[i] for i in unique_indices] - - sort_indices = sorted(range(len(fits_unique)), key=lambda k: fits_unique[k]) # sort fits + + sort_indices = sorted(list(range(len(fits_unique))), key=lambda k: fits_unique[k]) # sort fits fits_sort = [fits_unique[i] for i in sort_indices] gens_sort = [gens_unique[i] for i in sort_indices] cands_sort = [cands_unique[i] for i in sort_indices] params_sort = [params_unique[i] for i in sort_indices] - + # if multiple islands, save data for each if islands > 1: ind_gens_isl.append(ind_gens) # individuals data for islands ind_cands_isl.append(ind_cands) ind_fits_isl.append(ind_fits) ind_cs_isl.append(ind_cs) - + stat_gens_isl.append(stat_gens) # statistics.csv for islands stat_worstfits_isl.append(stat_worstfits) stat_bestfits_isl.append(stat_bestfits) stat_avgfits_isl.append(stat_avgfits) stat_stdfits_isl.append(stat_stdfits) - + fits_sort_isl.append(fits_sort) #sorted data - gens_sort_isl.append(gens_sort) + gens_sort_isl.append(gens_sort) cands_sort_isl.append(cands_sort) params_sort_isl.append(params_sort) - + if islands > 1: return ind_gens_isl, ind_cands_isl, ind_fits_isl, ind_cs_isl, stat_gens_isl, \ stat_worstfits_isl, stat_bestfits_isl, stat_avgfits_isl, stat_stdfits_isl, \ fits_sort_isl, gens_sort_isl, cands_sort_isl, params_sort_isl - diff --git a/arm.py b/arm.py index 28b95c1..2e764e1 100644 --- a/arm.py +++ b/arm.py @@ -4,7 +4,7 @@ Code to connect a virtual arm to the M1 model Multipe arm options depending on self.type: - 'randomInput': the object provides random values matching the format (used for test/debugging) -- 'dummyArm': simple kinematic arm implemented in python +- 'dummyArm': simple kinematic arm implemented in python - 'musculoskeletal': realistic musculoskeletal arm implemented in C++ and interfaced via pipes Adapted from arm.hoc in arm2dms @@ -24,7 +24,7 @@ [SH,EL] = [X,Y] = [0,1] [SH_EXT, SH_FLEX, EL_EXT, EL_FLEX] = [0,1,2,3] - + class Arm: #%% init def __init__(self, type, anim, graphs): # initialize variables @@ -43,19 +43,19 @@ def pos2angles(self, armPos, armLen): y = armPos[1] l1 = armLen[0] l2 = armLen[1] - elang = abs(2*arctan(sqrt(((l1 + l2)**2 - (x**2 + y**2))/((x**2 + y**2) - (l1 - l2)**2)))); - phi = arctan2(y,x); + elang = abs(2*arctan(sqrt(((l1 + l2)**2 - (x**2 + y**2))/((x**2 + y**2) - (l1 - l2)**2)))); + phi = arctan2(y,x); psi = arctan2(l2 * sin(elang), l1 + (l2 * cos(elang))); - shang = phi - psi; + shang = phi - psi; return [shang,elang] # convert joint angles to cartesian position - def angles2pos(self, armAng, armLen): + def angles2pos(self, armAng, armLen): elbowPosx = armLen[0] * cos(armAng[0]) # end of elbow elbowPosy = armLen[0] * sin(armAng[0]) wristPosx = elbowPosx + armLen[1] * cos(+armAng[0]+armAng[1]) # wrist=arm position wristPosy = elbowPosy + armLen[1] * sin(+armAng[0]+armAng[1]) - return [wristPosx,wristPosy] + return [wristPosx,wristPosy] #%% setTargetByID def setTargetByID(self, id, startAng, targetDist, armLen): @@ -114,12 +114,12 @@ def resetArm(self, s, t): self.error = 0 # error signal (eg. difference between ) self.critic = 0 # critic signal (1=reward; -1=punishment) self.initArmMovement = self.initArmMovement + s.testTime - + def setPMdInput(self, s): for c in range(0, s.popnumbers[s.PMd], 2): - try: - gid = s.popGidStart[s.PMd] + c # find gid of PMd + try: + gid = s.popGidStart[s.PMd] + c # find gid of PMd lid = s.gidDic[gid] # find local index corresponding to gid if (gid in s.targetPMdInputs[s.targetid]): # if gid in PMdinputs for this target #print "yes:",gid @@ -140,22 +140,22 @@ def RLcritic(self, t): diff = self.error - mean(self.errorAll[-RLsteps:-1]) # difference between error at t and at t-(RLdt/dt) eg. t-50/5 = t-10steps else: diff = 0 - if diff < -self.minRLerror: # if error negative: LTP + if diff < -self.minRLerror: # if error negative: LTP self.critic = 1 # return critic signal to model.py so can update synaptic weights elif diff > self.minRLerror: # if error positive: LTD self.critic = -1 else: # if difference not significant: no weight change self.critic = 0 - else: # if + else: # if self.critic = 0 return self.critic #%% plot joint angles def plotTraj(self, filename): - fig = figure() + fig = figure() l = 1.1*sum(self.armLen) ax = fig.add_subplot(111, autoscale_on=False, xlim=(-l/2, +l), ylim=(-l/2, +l)) # create subplot - posX, posY = zip(*[self.angles2pos([x[SH],x[EL]], self.armLen) for x in self.angAll]) + posX, posY = list(zip(*[self.angles2pos([x[SH],x[EL]], self.armLen) for x in self.angAll])) ax.plot(posX, posY, 'r') targ = Circle((self.targetPos),0.04, color='g', fill=False) # target ax.add_artist(targ) @@ -163,12 +163,12 @@ def plotTraj(self, filename): ax.set_title('X-Y Hand trajectory') xlabel('x') ylabel('y') - print 'saving '+filename + print('saving '+filename) fig.savefig(filename) #%% plot joint angles def plotAngs(self): - fig = figure() + fig = figure() ax = fig.add_subplot(111) # create subplot sh = [x[SH] for x in self.angAll] el = [x[EL] for x in self.angAll] @@ -176,8 +176,8 @@ def plotAngs(self): ax.plot(el, 'b', label='elbow') shTarg = self.pos2angles(self.targetPos, self.armLen)[0] elTarg = self.pos2angles(self.targetPos, self.armLen)[1] - ax.plot(range(0,len(sh)), [shTarg] * len(sh), 'r:', label='sh target') - ax.plot(range(0,len(el)), [elTarg] * len(el), 'b:', label='el target') + ax.plot(list(range(0,len(sh))), [shTarg] * len(sh), 'r:', label='sh target') + ax.plot(list(range(0,len(el))), [elTarg] * len(el), 'b:', label='el target') ax.set_title('Joint angles') xlabel('time') ylabel('angle') @@ -185,7 +185,7 @@ def plotAngs(self): #%% plot motor commands def plotMotorCmds(self): - fig = figure() + fig = figure() ax = fig.add_subplot(111) # create subplot shext = [x[SH_EXT] for x in self.motorCmdAll] elext = [x[EL_EXT] for x in self.motorCmdAll] @@ -202,7 +202,7 @@ def plotMotorCmds(self): #%% plot RL critic signal and error def plotRL(self): - fig = figure() + fig = figure() ax = fig.add_subplot(111) # create subplot ax.plot(self.errorAll, 'r', label='error') ax.plot((array(self.criticAll)+1.0) * max(self.errorAll) / 2.0, 'b', label='RL critic') @@ -214,12 +214,12 @@ def plotRL(self): ################################ ### SETUP ################################ - def setup(self, s):#, nduration, loopstep, RLinterval, pc, scale, popnumbers, p): + def setup(self, s):#, nduration, loopstep, RLinterval, pc, scale, popnumbers, p): self.duration = s.duration#/1000.0 # duration in msec - self.interval = s.loopstep#/1000.0 # interval between arm updates in ,sec + self.interval = s.loopstep#/1000.0 # interval between arm updates in ,sec self.RLinterval = s.RLinterval # interval between RL updates in msec self.minRLerror = s.minRLerror # minimum error change for RL (m) - self.armLen = s.armLen # elbow - shoulder from MSM;radioulnar - elbow from MSM; + self.armLen = s.armLen # elbow - shoulder from MSM;radioulnar - elbow from MSM; self.handPos = [0,0] # keeps track of hand (end-effector) x,y position self.handPosAll = [] # list with all handPos self.handVel = [0,0] # keeps track of hand (end-effector) x,y velocity @@ -267,11 +267,11 @@ def setup(self, s):#, nduration, loopstep, RLinterval, pc, scale, popnumbers, p) currentPval += angInterval - # initialize dummy or musculoskeletal arm - if s.rank == 0: - if self.type == 'dummyArm': + # initialize dummy or musculoskeletal arm + if s.rank == 0: + if self.type == 'dummyArm': self.setupDummyArm() # setup dummyArm (eg. graph animation) - elif self.type == 'musculoskeletal': + elif self.type == 'musculoskeletal': damping = 1 # damping of muscles (.osim parameter) shExtGain = 2 # gain factor to multiply force of shoulder extensor muscles (.osim parameter) shFlexGain = 1 # gain factor to multiply force of shoulder flexor muscles (.osim parameter) @@ -279,18 +279,18 @@ def setup(self, s):#, nduration, loopstep, RLinterval, pc, scale, popnumbers, p) elFlexGain = 0.8 # gain factor to multiply force of elbox flexor muscles (.osim parameter) # call function to initialize virtual arm params and run virtual arm C++ executable arminterface.setup(self.duration/1000.0, self.interval, self.startAng[SH], self.startAng[EL], self.targetPos[X], self.targetPos[Y], damping, shExtGain, shFlexGain, elExtGain, elFlexGain) - + # set PMd inputs - if s.PMdinput == 'targetSplit': + if s.PMdinput == 'targetSplit': self.setPMdInput(s) # set PMd inputs - ################################ - ### RUN ################################ - def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs=[]): + ### RUN + ################################ + def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs=[]): # Append to list the the value of relevant variables for this time step (only worker0) - if s.rank == 0: + if s.rank == 0: self.handPosAll.append(list(self.handPos)) # list with all handPos self.handVelAll.append(list(self.handVel)) # list with all handVel self.angAll.append(list(self.ang)) # list with all ang @@ -303,7 +303,7 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= ############################ # dummyArm or musculoskeletal: gather spikes for motor command ############################ - if self.type == 'dummyArm' or self.type == 'musculoskeletal': + if self.type == 'dummyArm' or self.type == 'musculoskeletal': ## Exploratory movements if s.explorMovs and t-s.timeoflastexplor >= self.randDur: # if time to update exploratory movement seed(s.id32('%d'%(int(t)+s.randseed))) # init seed @@ -313,7 +313,7 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= s.timeoflastexplor = t if s.explorMovs == 1: # add random noise to EDSC+IDSC population - IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) + IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) IDSCgidsInverse = [j for j in range(s.popGidStart[s.IDSC],s.popGidEnd[s.IDSC]+1) if j not in list(IDSCgids)] for (i,x) in enumerate(s.backgroundsources): # for all background input netstims if s.backgroundgid[i] in s.motorCmdCellRange[self.randMus] or s.backgroundgid[i] in list(IDSCgids): # if connected to chosen muscle cells @@ -322,12 +322,12 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= elif s.backgroundgid[i] in [y for z in range(s.nMuscles) if z != self.randMus for y in s.motorCmdCellRange[z]+IDSCgidsInverse]: x.interval = s.backgroundrateMin**-1*1e3 # otherwise set to minimum #if s.rank==0: print 'exploratory movement, randMus',self.randMus,' strength:',self.randMul,' duration:', self.randDur - #print 'IDSC cells:', IDSCgids + #print 'IDSC cells:', IDSCgids #print 'IDSC inverse cells:', IDSCgidsInverse - + elif s.explorMovs == 2: # add random noise to EB5 population - self.targetCells = range(s.popGidStart[s.EB5], s.popGidEnd[s.EB5]+1) # EB5 cell gids + self.targetCells = list(range(s.popGidStart[s.EB5], s.popGidEnd[s.EB5]+1)) # EB5 cell gids self.randNumCells = randint(1, int(s.explorCellsFraction*len(self.targetCells))) # num of cells to stimumales self.randCells = sample(self.targetCells, int(self.randNumCells)) # select random gids for (i,x) in enumerate(s.backgroundsources): # for all input background netstims @@ -335,34 +335,34 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= x.interval = s.backgroundrateExplor**-1*1e3 # increase input elif s.backgroundgid[i] in [y for y in self.targetCells if y not in self.randCells]: # for the rest of cells x.interval = s.backgroundrateMin**-1*1e3 # set to normal level - #if s.rank==0: - #print 'Nodes:', s.rank,' - exploratory movement, numcells:',self.randNumCells,' strength:',self.randMul,' duration:', self.randDur, 'cells:', self.randCells + #if s.rank==0: + #print 'Nodes:', s.rank,' - exploratory movement, numcells:',self.randNumCells,' strength:',self.randMul,' duration:', self.randDur, 'cells:', self.randCells ## Reset arm and set target after every trial -start from center etc - if s.trialReset and t-s.timeoflastreset > s.testTime: + if s.trialReset and t-s.timeoflastreset > s.testTime: self.resetArm(s, t) s.targetid = s.trialTargets[self.trial] # set target based on trial number - self.targetPos = self.setTargetByID(s.targetid, self.startAng, self.targetDist, self.armLen) + self.targetPos = self.setTargetByID(s.targetid, self.startAng, self.targetDist, self.armLen) if s.PMdinput == 'targetSplit': self.setPMdInput(s) ## Only move after initial period - avoids initial transitory spiking period (NSLOC sync spikes), and allows for variables with history to clear # can be justified as preparatory period (eg. watiing for go cue) if t > self.initArmMovement: - ## Gather spikes #### from all vectors to then calculate motor command + ## Gather spikes #### from all vectors to then calculate motor command for i in range(s.nMuscles): cmdVecs = [array(s.hostspikevecs[x]) for x in range(s.cellsperhost) if (s.gidVec[x] in s.motorCmdCellRange[i])] self.motorCmd[i] = sum([len(x[(x < t) * (x > t-self.cmdtimewin)]) for x in cmdVecs]) s.pc.allreduce(self.vec.from_python([self.motorCmd[i]]), 1) # sum - self.motorCmd[i] = self.vec.to_python()[0] + self.motorCmd[i] = self.vec.to_python()[0] # else: # for i in range(s.nMuscles): # stimulate all muscles equivalently so arm doesnt move # self.motorCmd[i] = 0.2 * self.cmdmaxrate - ## Calculate final motor command - if s.rank==0: - self.motorCmd = [x / self.cmdmaxrate for x in self.motorCmd] # normalize motor command + ## Calculate final motor command + if s.rank==0: + self.motorCmd = [x / self.cmdmaxrate for x in self.motorCmd] # normalize motor command if s.antagInh: # antagonist inhibition if self.motorCmd[SH_EXT] > self.motorCmd[SH_FLEX]: # sh ext > sh flex self.motorCmd[SH_FLEX] = self.motorCmd[SH_FLEX]**2 / self.motorCmd[SH_EXT] / s.antagInh @@ -372,7 +372,7 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= self.motorCmd[EL_FLEX] = self.motorCmd[EL_FLEX]**2 / self.motorCmd[EL_EXT] / s.antagInh elif self.motorCmd[EL_EXT] < self.motorCmd[EL_FLEX]: # el ext > el flex self.motorCmd[EL_EXT] = self.motorCmd[EL_EXT]**2 / self.motorCmd[EL_FLEX] / s.antagInh - + ############################ # ALL arms: Send motor command to virtual arm; receive new position; update proprioceptive population (ASC) @@ -387,29 +387,29 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= dataReceived = [self.ang[SH], self.ang[EL]] if not dataReceived or dataReceived==[-3,3]: # if error receiving packet dataReceived = [self.ang[SH], self.ang[EL]] # use previous packet - print 'Missed packet at t=%.2f',t + print('Missed packet at t=%.2f',t) elif self.type == 'dummyArm': # DUMMYARM dataReceived = self.runDummyArm(self.motorCmd) # run dummyArm elif self.type == 'randomOutput': # RANDOMOUTPUT - dataReceived = [0,0] - dataReceived[0] = uniform(self.minPval, self.maxPval) # generate 2 random values - dataReceived[1] = uniform(self.minPval, self.maxPval) - # broadcast dataReceived to other workers. + dataReceived = [0,0] + dataReceived[0] = uniform(self.minPval, self.maxPval) # generate 2 random values + dataReceived[1] = uniform(self.minPval, self.maxPval) + # broadcast dataReceived to other workers. n = s.pc.broadcast(self.vec.from_python(dataReceived), 0) # convert python list to hoc vector for broadcast data received from arm else: # other workers n = s.pc.broadcast(self.vec, 0) # receive shoulder and elbow angles from worker0 so can compare with cells in this worker - dataReceived = self.vec.to_python() + dataReceived = self.vec.to_python() if self.type == 'musculoskeletal': [self.ang[SH], self.ang[EL]] = dataReceived - self.handPos = self.angles2pos(self.ang, self.armLen) - self.angVel[SH] = self.angVel[EL] = 0 + self.handPos = self.angles2pos(self.ang, self.armLen) + self.angVel[SH] = self.angVel[EL] = 0 else: [self.ang[SH], self.ang[EL], self.angVel[SH], self.angVel[EL], self.handPos[SH], self.handPos[EL]] = dataReceived # map data received to shoulder and elbow angles - #[self.ang[SH], self.ang[EL], self.angVel[SH], self.angVel[EL]] = dataReceived # map data received to shoulder and elbow angles - + #[self.ang[SH], self.ang[EL], self.angVel[SH], self.angVel[EL]] = dataReceived # map data received to shoulder and elbow angles + #### Update proprio pop ASC for c in range(0, self.numPcells, 2): - try: + try: id = s.gidDic[self.pStart + c] # find local index corresponding to gid if (self.ang[SH] >= self.prange[c,0] and self.ang[SH] < self.prange[c,1]): # in angle in range -> high firing rate s.cells[id].interval=1000/self.maxPrate # interval in ms as a function of rate @@ -417,35 +417,35 @@ def run(self, t, s): #pc, cells, gidVec, gidDic, cellsperhost=[], hostspikevecs= s.cells[id].interval=1000/self.minPrate # interval in ms as a function of rate except: pass # local index corresponding to gid not found in this node - try: + try: id = s.gidDic[self.pStart + c + 1] # find local index corresponding to gid if (self.ang[EL] >= self.prange[c+1,0] and self.ang[EL] < self.prange[c+1,1]): # in angle in range -> high firing rate s.cells[id].interval=1000/self.maxPrate # interval in ms as a function of rate else: # if angle not in range -> low firing rate - s.cells[id].interval=1000/self.minPrate # interval in ms as a function of rate + s.cells[id].interval=1000/self.minPrate # interval in ms as a function of rate except: pass # local index corresponding to gid not found in this node - #### Calculate error between hand and target for interval between RL updates + #### Calculate error between hand and target for interval between RL updates if s.rank == 0 and self.initArmMovement: # do not update between trials #print 't=%.2f, xpos=%.2f'%(t,self.targetPos[X]) self.error = sqrt((self.handPos[X] - self.targetPos[X])**2 + (self.handPos[Y] - self.targetPos[Y])**2) - + return self.critic ################################ ### CLOSE ################################ - def close(self, s): + def close(self, s): if self.type == 'randomOutput': print('\nClosing random output virtual arm...') if self.type == 'dummyArm': if s.explorMovs == 1: # remove explor movs related noise to cells for imus in range(s.nMuscles): - IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) + IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) for (i,x) in enumerate(s.backgroundsources): if s.backgroundgid[i] in s.motorCmdCellRange[imus]+list(IDSCgids): x.interval = 0.0001**-1*1e3 @@ -460,21 +460,21 @@ def close(self, s): self.initArmMovement = int(s.initArmMovement) if s.rank == 0: - print('\nClosing dummy virtual arm ...') + print('\nClosing dummy virtual arm ...') if self.anim: ioff() # turn interactive mode off - close(self.fig) # close arm animation graph + close(self.fig) # close arm animation graph if self.graphs: # plot graphs self.plotTraj(s.filename) #self.plotAngs() #self.plotMotorCmds() #self.plotRL() - + if self.type == 'musculoskeletal': if s.explorMovs == 1: # remove explor movs related noise to cells for imus in range(s.nMuscles): - IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) + IDSCgids = array(s.motorCmdCellRange[self.randMus]) - int(s.popGidStart[s.EDSC]) + int(s.popGidStart[s.IDSC]) for (i,x) in enumerate(s.backgroundsources): if s.backgroundgid[i] in s.motorCmdCellRange[imus]+list(IDSCgids): x.interval = 0.0001**-1*1e3 @@ -489,7 +489,7 @@ def close(self, s): self.initArmMovement = int(s.initArmMovement) if s.rank == 0: - print('\nClosing dummy virtual arm ...') + print('\nClosing dummy virtual arm ...') arminterface.closeSavePlot(self.duration/1000.0, self.interval) if self.graphs: # plot graphs @@ -497,8 +497,3 @@ def close(self, s): #self.plotAngs() self.plotMotorCmds() self.plotRL() - - - - - diff --git a/armGraphs.py b/armGraphs.py index 44a899a..bb53c86 100644 --- a/armGraphs.py +++ b/armGraphs.py @@ -1,11 +1,11 @@ # armGraphs.py -- Python code for the following purposes: # (modified from testArm.py) # *Read joint angles, muscle exc, muscle activation and force from .pnt files -# *Plot for each muscle: exc, act, force, muscle lengths (3 fibres+avg) +# *Plot for each muscle: exc, act, force, muscle lengths (3 fibres+avg) # *Plot joint angles and trajectory over time # Last update: 6/17/13 (salvadord) -import sys #for exit +import sys #for exit import struct import time import matplotlib @@ -15,7 +15,7 @@ import matplotlib.animation as animation from mpl_toolkits.mplot3d import Axes3D from pylab import figure, show -from numpy import * +from numpy import * import csv #import os @@ -27,322 +27,322 @@ #secLength = 6.01 # length of simulation in seconds (first 10 ms used to set start pos) #msecInterval = 10.0 # interval at which packets are sent in miliseconds #n = int(secLength*1000/msecInterval) # calculate number of samples -#toSamples = 1/msecInterval*1000 +#toSamples = 1/msecInterval*1000 numMuscles = 4 # number of muscles = shFlex (PECT), shExt (DELT), elFlex (BIC), elExt (TRIC) #muscleNames = ["Shoulder ext (post Deltoid+Infraspinatus)", "Shoulder flex (Pectoralis+ant Deltoid)", "Elbow ext (Triceps)", "Elbow flex (Biceps+Brachialis)"] muscleNames = ["Shoulder ext", "Shoulder flex", "Elbow ext", "Elbow flex"] -numMusBranches = 18 # number of muscle branches -numJoints = 2 # number of joints (DOFs) +numMusBranches = 18 # number of muscle branches +numJoints = 2 # number of joints (DOFs) armLen = [0.4634 - 0.173, 0.7169 - 0.4634] # elbow - shoulder from MSM;radioulnar - elbow from MSM; useJointPos = 0 # use joint positions vs joint angles for 2d arm animation showBranches = 0 # include muscle branches in graphs -verbose = 0 # whether to show output on screen +verbose = 0 # whether to show output on screen ############################### # Function definition ############################### - -# function to read the result .pnt file containing muscle excitation and force + +# function to read the result .pnt file containing muscle excitation and force def readPntFiles(msmFolder, pntFile, secLength, msecInterval): - - n = int(secLength*1000/msecInterval) # calculate number of samples - - ################################### - # read joint position - - # open pnt file and read data - fileName = msmFolder+"SUNY_arm_2DOFs_horizon_static_coordinate_status.pnt" - try: - #f = open(fileName, "r" ) - jointData = [] - with open(fileName, "r") as f: - #next(f) - for line in f: - print line - jointData.append([float(i) for i in line.split()]) - jointData = array(jointData).transpose() - jointData = jointData[:, len(jointData[0])-n:] # make number of rows equal to n (packets received) - - # create arrays to store read data (joints include shoulder, elbow and wrist * 3 coords (xyz)) - jointPosSeq = zeros(((numJoints+1)*3, n)) - - # assign activation and force to output arrays - # file has format: time,ground_thorax_xyz,sternoclavicular_xyz,acromioclavicular_xyz,shoulder_xyz,elbow_xyz,radioulnar_xyz,radius_hand_xyz - - #jointPosSeq[0:6,:] = jointData[10:16,:] - #jointPosSeq[6:9,:] = jointData[19:22,:] - - jointPosSeq[0,:] = jointData[10,:] - jointPosSeq[1,:] = jointData[12,:] - jointPosSeq[2,:] = jointData[13,:] - jointPosSeq[3,:] = jointData[15,:] - jointPosSeq[4,:] = jointData[19,:] - jointPosSeq[5,:] = jointData[21,:] - except: - jointPosSeq=[] - if verbose: - print "coordinate pnt file not available" - - ######################################################## - # Read muscle activation and force - - # open pnt file and read data - #fileName = msmFolder+"SUNY_arm_2DOFs_horizon_static_muscle_status.pnt" - fileName = pntFile - #f = open(fileName, "r" ) - muscleData = [] - with open(fileName, "r") as f: - #next(f) - for line in f: - muscleData.append([float(i) for i in line.split()]) - muscleData = array(muscleData).transpose() # transpose - muscleData = muscleData[:, len(muscleData[0])-n:] # make number of rows equal to n (packets received) - - # create arrays to store read data - musExcSeq = zeros((numMusBranches, n)) - musActSeq = zeros((numMusBranches, n)) - musForcesSeq = zeros((numMusBranches, n)) - - # assign activation and force to output arrays - # file has format: time DELT1_excitation DELT1_activation DELT1_force DELT2_excitation DELT2_activation DELT2_force ... - musExcSeq = muscleData[1::3 , :] - musActSeq = muscleData[2::3 , :] - musForcesSeq = muscleData[3::3 , :] - - return jointPosSeq, musExcSeq, musActSeq, musForcesSeq - -# Plot for each muscle: exc, act, force, muscle lengths (3 fibres+avg), joint angles and trajectory over time + + n = int(secLength*1000/msecInterval) # calculate number of samples + + ################################### + # read joint position + + # open pnt file and read data + fileName = msmFolder+"SUNY_arm_2DOFs_horizon_static_coordinate_status.pnt" + try: + #f = open(fileName, "r" ) + jointData = [] + with open(fileName, "r") as f: + #next(f) + for line in f: + print(line) + jointData.append([float(i) for i in line.split()]) + jointData = array(jointData).transpose() + jointData = jointData[:, len(jointData[0])-n:] # make number of rows equal to n (packets received) + + # create arrays to store read data (joints include shoulder, elbow and wrist * 3 coords (xyz)) + jointPosSeq = zeros(((numJoints+1)*3, n)) + + # assign activation and force to output arrays + # file has format: time,ground_thorax_xyz,sternoclavicular_xyz,acromioclavicular_xyz,shoulder_xyz,elbow_xyz,radioulnar_xyz,radius_hand_xyz + + #jointPosSeq[0:6,:] = jointData[10:16,:] + #jointPosSeq[6:9,:] = jointData[19:22,:] + + jointPosSeq[0,:] = jointData[10,:] + jointPosSeq[1,:] = jointData[12,:] + jointPosSeq[2,:] = jointData[13,:] + jointPosSeq[3,:] = jointData[15,:] + jointPosSeq[4,:] = jointData[19,:] + jointPosSeq[5,:] = jointData[21,:] + except: + jointPosSeq=[] + if verbose: + print("coordinate pnt file not available") + + ######################################################## + # Read muscle activation and force + + # open pnt file and read data + #fileName = msmFolder+"SUNY_arm_2DOFs_horizon_static_muscle_status.pnt" + fileName = pntFile + #f = open(fileName, "r" ) + muscleData = [] + with open(fileName, "r") as f: + #next(f) + for line in f: + muscleData.append([float(i) for i in line.split()]) + muscleData = array(muscleData).transpose() # transpose + muscleData = muscleData[:, len(muscleData[0])-n:] # make number of rows equal to n (packets received) + + # create arrays to store read data + musExcSeq = zeros((numMusBranches, n)) + musActSeq = zeros((numMusBranches, n)) + musForcesSeq = zeros((numMusBranches, n)) + + # assign activation and force to output arrays + # file has format: time DELT1_excitation DELT1_activation DELT1_force DELT2_excitation DELT2_activation DELT2_force ... + musExcSeq = muscleData[1::3 , :] + musActSeq = muscleData[2::3 , :] + musForcesSeq = muscleData[3::3 , :] + + return jointPosSeq, musExcSeq, musActSeq, musForcesSeq + +# Plot for each muscle: exc, act, force, muscle lengths (3 fibres+avg), joint angles and trajectory over time def plotGraphs(jointPosSeq, jointAnglesSeq, musLengthsSeq, musExcSeq, musActSeq, musForcesSeq, t1, t2, msecInterval, armAnimation, saveGraphs, saveName): - # graph parameters - toDegrees = 360/(2*pi) - toSamples = 1/msecInterval*1000 - legFont = 11 - linWidth = 3 - color1 = "red" - color2 = "blue" - color3 = "green" - color4 = "purple" - line1 = "-" - line2 = ":" - gridOn = 1 - - # plot with 6 subplots for joint angles, trajectory and each of muscles - fig1 = figure() - - # time variables - T = arange(t1, t2, msecInterval/1000.0); - t1Samples = t1*toSamples - t2Samples = t2*toSamples - - if armAnimation: - speedFactor = 1 # speed of animation - - #create figure and axes - figAnim = figure() - axAnim = figAnim.add_subplot(111) - axAnim.grid(gridOn); - # set axes size to fit arm - axAnim.set_xlim(-armLen[0]-armLen[1]-0.1, armLen[0]+armLen[1]+0.1) - axAnim.set_ylim(-armLen[0]-armLen[1]-0.1, armLen[0]+armLen[1]+0.1) - - ########################### - # 2D arm movement animation - armImages=[] - for t in arange(t1Samples,t2Samples): - # Update the arm position based on jointPosSeq - if useJointPos: - # use x = z (pnt file) ; y = x (pnt file) - shoulderPosx = jointPosSeq[1, t] - shoulderPosy = jointPosSeq[0, t] - elbowPosx = jointPosSeq[3, t] - elbowPosy = jointPosSeq[2, t] - wristPosx = jointPosSeq[5, t] - wristPosy = jointPosSeq[4, t] - - # update jointAnglesSeq based on pos !!! - - # Update the arm position based on jointAnglesSeq - else: - armAng = jointAnglesSeq[:,t] - shoulderPosx = 0 - shoulderPosy = 0 - elbowPosx = armLen[0] * cos(armAng[0]) # end of elbow - elbowPosy = armLen[0] * sin(armAng[0]) - wristPosx = elbowPosx + armLen[1] * cos(+armAng[0]+armAng[1]) # wrist=arm position - wristPosy = elbowPosy + armLen[1] * sin(+armAng[0]+armAng[1]) - - # create - armLine1 = lines.Line2D([0, elbowPosx-shoulderPosx], [0, elbowPosy-shoulderPosy], color=color1, linestyle=line1, linewidth=linWidth) - armLine2 = lines.Line2D([elbowPosx-shoulderPosx, wristPosx-shoulderPosx], [elbowPosy-shoulderPosy, wristPosy-shoulderPosy], color=color2, linestyle=line1, linewidth=linWidth) - axAnim.add_line(armLine1) - axAnim.add_line(armLine2) - #label = plt.legend(armLine1, str(t/toSamples) ) - label = text.Text(x=0, y=0.5, text="time = "+str(t/toSamples), weight='bold' ) - axAnim.add_artist(label) - armImages.append([armLine1, armLine2, label]) - - # add blank frames - blankFrames = int(1*toSamples) - for t in range(blankFrames): - # Update the arm position - armLine1 = lines.Line2D([0, 0], [0, 0], color=color1, linestyle=line1, linewidth=linWidth) - armLine2 = lines.Line2D([0,0 ], [0,0], color=color2, linestyle=line1, linewidth=linWidth) - axAnim.add_line(armLine1) - axAnim.add_line(armLine2) - armImages.append([armLine1, armLine2]) - - # generate animation - armAnim = animation.ArtistAnimation(figAnim, armImages, interval=msecInterval/speedFactor, repeat_delay=500, blit=True) - - ########################### - # Plot joint angles vs time - ax = fig1.add_subplot(321) - #T = arange(0, secLength, msecInterval/1000.0); - - T=T[:len(jointAnglesSeq[0,t1Samples:t2Samples])] - ax.plot(T,jointAnglesSeq[0,t1Samples:t2Samples]*toDegrees,color=color1,linestyle=line1, linewidth=linWidth, label="shoulder") - ax.plot(T,jointAnglesSeq[1,t1Samples:t2Samples]*toDegrees,color=color2,linestyle=line1, linewidth=linWidth, label="elbow") - - ax.set_ylabel('angle (deg)', fontsize = legFont) - ax.set_xlabel('time (sec)', fontsize = legFont) - ax.set_title('Joint angles') - ax.legend(loc='upper center', bbox_to_anchor=(1.0, 1.0), borderaxespad=0., prop={'size':legFont}) - #ax.set_xlim([t1, t2]) - #ax.set_ylim(bmmYlims_sh) - ax.grid(gridOn) - - ############################ - # Plot x-y pos vs time - ax = fig1.add_subplot(322) - - # calculate x and y trajectories based on angles - if useJointPos: - xTraj = jointPosSeq[4, t1Samples:t2Samples] - yTraj = jointPosSeq[5, t1Samples:t2Samples] - else: - xTraj = armLen[0]*cos(jointAnglesSeq[0,t1Samples:t2Samples])+armLen[1]*cos(jointAnglesSeq[1,t1Samples:t2Samples]) - yTraj = armLen[0]*sin(jointAnglesSeq[0,t1Samples:t2Samples])+armLen[0]*sin(jointAnglesSeq[1,t1Samples:t2Samples]) - - #ax.plot(xTraj, yTraj,color=color2,linestyle=line1, linewidth=linWidth) - ax.plot(T, xTraj,color=color1,linestyle=line1, linewidth=linWidth, label="x") - ax.plot(T, yTraj,color=color2,linestyle=line1, linewidth=linWidth, label="y") - - ax.set_ylabel('position (m)', fontsize = legFont) - #ax.set_xlabel('x position (m)', fontsize = legFont) - ax.set_xlabel('time (sec)', fontsize = legFont) - ax.set_title('X-Y trajectory') - ax.legend(loc='upper center', bbox_to_anchor=(1.0, 1.0), borderaxespad=0., prop={'size':legFont}) - #ax.set_xlim([t1, t2]) - #ax.set_ylim(bmmYlims_sh) - ax.grid(gridOn) - - ############################ - # Plot excitation, activation and force for each muscle - - # calculate normalized force (activation and length already normalized) - #musActSeqNorm = musActSeq/musActSeq[:,t1Samples:t2Samples].max() - musForcesSeqNorm = musForcesSeq/musForcesSeq[:,t1Samples:t2Samples].max() - #musLengthsSeqNorm = musLengthsSeq/musLengthsSeq[:,t1Samples:t2Samples].max() - - # Note arrangement of muscle branches in data arrays: - #DELT1(0) DELT2(1) DELT3(2) Infraspinatus(3) Latissimus_dorsi_1(4) Latissimus_dorsi_2(5) Latissimus_dorsi_3(6) - #Teres_minor(7) PECM1(8) PECM2(9) PECM3(10) Coracobrachialis(11) TRIlong(12) TRIlat(13) TRImed(14) BIClong(15) BICshort(16) BRA(17) - # is different from muscle groups: - # Sh ext = DELT3, Infraspinatus, Latissimus_dorsi_1-3, Teres_minor - # Sh flex = PECM1, PECM2, PECM3, DELT1, Coracobrachialis - # El ext = TRIlong, TRIlat, TRImed - # El flex = BIClong, BICshort, BRA - shext=[2,3,4,5,6,7] - shflex=[0,8,9,10,11] - elext=[12,13,14] - elflex=[15,16,17] - musList=[shext, shflex,elext,elflex] - - for iMus in range(numMuscles): - ax = fig1.add_subplot(3,2,iMus+3) - # set number of muscle branches - assume each node has 3 branches (treat the Brachialis as a branch of Biceps=elbow flexor) - #iBranches = 3 - - ### Excitation and Activation #### - # equivalent for all branches of same muscle group - offset = 2 # use offset 3 because only DELT3 is used (order of muscle branches doesn't correspond muscle groups!) - ax.plot(T, musExcSeq[musList[iMus][offset],t1Samples:t2Samples],color=color1,linestyle=line1, linewidth=linWidth, label="excitation") - ax.plot(T, musActSeq[musList[iMus][offset],t1Samples:t2Samples],color=color2,linestyle=line1, linewidth=linWidth, label="activation") - - # for show branches plot individual branches and mean value for force and length - if showBranches: - ### Force and Length ### - for iBranch in range(len(musList[iMus])): - ax.plot(T, musForcesSeqNorm[musList[iMus][iBranch],t1Samples:t2Samples],color=color3,linestyle=line2, linewidth=linWidth-1) - ax.plot(T, musLengthsSeq[musList[iMus][iBranch],t1Samples:t2Samples],color=color4,linestyle=line2, linewidth=linWidth-1) - ax.plot(T, musForcesSeqNorm[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color3,linestyle=line1, linewidth=linWidth, label="force (mean)") - ax.plot(T, musLengthsSeq[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color4,linestyle=line1, linewidth=linWidth, label="length (mean)") - # for NOT show branches show mean value for force and single value for length - else: - ### Force ### - # For shoulder extensor group show only posterior Deltoid, branch 3 (DELT3) or Infraspinatus (INFSP) - # branch 2 (DELT2 = lateral deltoid) also available but currently not included in shoulder extensor group - if iMus == 0: - offset = 2 # DELT3 - #offset = 12 # INFSP - ax.plot(T, musForcesSeqNorm[musList[iMus][offset],t1Samples:t2Samples],color=color3,linestyle=line1, linewidth=linWidth, label="force") - # For rest of muscles use mean value of all branches - else: - offset=0 - ax.plot(T, musForcesSeqNorm[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color3,linestyle=line1, linewidth=linWidth, label="force") - - ### Length #### - # show length only of one muscle indicated by the index 'offset' - DELT3, PECM1, BIClong, TRIlong - maxLength = 0.20 - ax.plot(T, musLengthsSeq[musList[iMus][offset],t1Samples:t2Samples]/maxLength,color=color4,linestyle=line1, linewidth=linWidth, label="length") - - # show branche label - if (showBranches): - ax.plot(-1, -1,color=color3,linestyle=line2, linewidth=linWidth, label="force (branches)") - ax.plot(-1, -1,color=color4,linestyle=line2, linewidth=linWidth, label="length (branches)") - - # axis properties - ax.set_ylabel('normalized value', fontsize = legFont) - ax.set_ylim([0,1]) - ax.set_xlim([t1,t2]) - ax.set_xlabel('time (sec)', fontsize = legFont) - ax.set_title(muscleNames[iMus]) - if iMus==3: - ax.legend(loc='upper center', bbox_to_anchor=(-0.2, 1.8), borderaxespad=0., prop={'size':legFont}) - ax.grid(gridOn) - - # show graphs - fig1.tight_layout() - show() - - # save graphs using startPos and pattern in filename - if saveGraphs: - saveFolder = 'gif/' - fig1.savefig(saveFolder+saveName, bbox_inches=0) - - #if armAnimation: - #armAnim.save('test.mp4') - #armAnim.save(saveFolder+saveName+'.mp4',writer = writer) - - -# run single test (udp transfer, read files, plot graphs) + # graph parameters + toDegrees = 360/(2*pi) + toSamples = 1/msecInterval*1000 + legFont = 11 + linWidth = 3 + color1 = "red" + color2 = "blue" + color3 = "green" + color4 = "purple" + line1 = "-" + line2 = ":" + gridOn = 1 + + # plot with 6 subplots for joint angles, trajectory and each of muscles + fig1 = figure() + + # time variables + T = arange(t1, t2, msecInterval/1000.0); + t1Samples = t1*toSamples + t2Samples = t2*toSamples + + if armAnimation: + speedFactor = 1 # speed of animation + + #create figure and axes + figAnim = figure() + axAnim = figAnim.add_subplot(111) + axAnim.grid(gridOn); + # set axes size to fit arm + axAnim.set_xlim(-armLen[0]-armLen[1]-0.1, armLen[0]+armLen[1]+0.1) + axAnim.set_ylim(-armLen[0]-armLen[1]-0.1, armLen[0]+armLen[1]+0.1) + + ########################### + # 2D arm movement animation + armImages=[] + for t in arange(t1Samples,t2Samples): + # Update the arm position based on jointPosSeq + if useJointPos: + # use x = z (pnt file) ; y = x (pnt file) + shoulderPosx = jointPosSeq[1, t] + shoulderPosy = jointPosSeq[0, t] + elbowPosx = jointPosSeq[3, t] + elbowPosy = jointPosSeq[2, t] + wristPosx = jointPosSeq[5, t] + wristPosy = jointPosSeq[4, t] + + # update jointAnglesSeq based on pos !!! + + # Update the arm position based on jointAnglesSeq + else: + armAng = jointAnglesSeq[:,t] + shoulderPosx = 0 + shoulderPosy = 0 + elbowPosx = armLen[0] * cos(armAng[0]) # end of elbow + elbowPosy = armLen[0] * sin(armAng[0]) + wristPosx = elbowPosx + armLen[1] * cos(+armAng[0]+armAng[1]) # wrist=arm position + wristPosy = elbowPosy + armLen[1] * sin(+armAng[0]+armAng[1]) + + # create + armLine1 = lines.Line2D([0, elbowPosx-shoulderPosx], [0, elbowPosy-shoulderPosy], color=color1, linestyle=line1, linewidth=linWidth) + armLine2 = lines.Line2D([elbowPosx-shoulderPosx, wristPosx-shoulderPosx], [elbowPosy-shoulderPosy, wristPosy-shoulderPosy], color=color2, linestyle=line1, linewidth=linWidth) + axAnim.add_line(armLine1) + axAnim.add_line(armLine2) + #label = plt.legend(armLine1, str(t/toSamples) ) + label = text.Text(x=0, y=0.5, text="time = "+str(t/toSamples), weight='bold' ) + axAnim.add_artist(label) + armImages.append([armLine1, armLine2, label]) + + # add blank frames + blankFrames = int(1*toSamples) + for t in range(blankFrames): + # Update the arm position + armLine1 = lines.Line2D([0, 0], [0, 0], color=color1, linestyle=line1, linewidth=linWidth) + armLine2 = lines.Line2D([0,0 ], [0,0], color=color2, linestyle=line1, linewidth=linWidth) + axAnim.add_line(armLine1) + axAnim.add_line(armLine2) + armImages.append([armLine1, armLine2]) + + # generate animation + armAnim = animation.ArtistAnimation(figAnim, armImages, interval=msecInterval/speedFactor, repeat_delay=500, blit=True) + + ########################### + # Plot joint angles vs time + ax = fig1.add_subplot(321) + #T = arange(0, secLength, msecInterval/1000.0); + + T=T[:len(jointAnglesSeq[0,t1Samples:t2Samples])] + ax.plot(T,jointAnglesSeq[0,t1Samples:t2Samples]*toDegrees,color=color1,linestyle=line1, linewidth=linWidth, label="shoulder") + ax.plot(T,jointAnglesSeq[1,t1Samples:t2Samples]*toDegrees,color=color2,linestyle=line1, linewidth=linWidth, label="elbow") + + ax.set_ylabel('angle (deg)', fontsize = legFont) + ax.set_xlabel('time (sec)', fontsize = legFont) + ax.set_title('Joint angles') + ax.legend(loc='upper center', bbox_to_anchor=(1.0, 1.0), borderaxespad=0., prop={'size':legFont}) + #ax.set_xlim([t1, t2]) +#ax.set_ylim(bmmYlims_sh) + ax.grid(gridOn) + + ############################ + # Plot x-y pos vs time + ax = fig1.add_subplot(322) + + # calculate x and y trajectories based on angles + if useJointPos: + xTraj = jointPosSeq[4, t1Samples:t2Samples] + yTraj = jointPosSeq[5, t1Samples:t2Samples] + else: + xTraj = armLen[0]*cos(jointAnglesSeq[0,t1Samples:t2Samples])+armLen[1]*cos(jointAnglesSeq[1,t1Samples:t2Samples]) + yTraj = armLen[0]*sin(jointAnglesSeq[0,t1Samples:t2Samples])+armLen[0]*sin(jointAnglesSeq[1,t1Samples:t2Samples]) + + #ax.plot(xTraj, yTraj,color=color2,linestyle=line1, linewidth=linWidth) + ax.plot(T, xTraj,color=color1,linestyle=line1, linewidth=linWidth, label="x") + ax.plot(T, yTraj,color=color2,linestyle=line1, linewidth=linWidth, label="y") + + ax.set_ylabel('position (m)', fontsize = legFont) + #ax.set_xlabel('x position (m)', fontsize = legFont) + ax.set_xlabel('time (sec)', fontsize = legFont) + ax.set_title('X-Y trajectory') + ax.legend(loc='upper center', bbox_to_anchor=(1.0, 1.0), borderaxespad=0., prop={'size':legFont}) + #ax.set_xlim([t1, t2]) +#ax.set_ylim(bmmYlims_sh) + ax.grid(gridOn) + + ############################ + # Plot excitation, activation and force for each muscle + + # calculate normalized force (activation and length already normalized) + #musActSeqNorm = musActSeq/musActSeq[:,t1Samples:t2Samples].max() + musForcesSeqNorm = musForcesSeq/musForcesSeq[:,t1Samples:t2Samples].max() + #musLengthsSeqNorm = musLengthsSeq/musLengthsSeq[:,t1Samples:t2Samples].max() + + # Note arrangement of muscle branches in data arrays: + #DELT1(0) DELT2(1) DELT3(2) Infraspinatus(3) Latissimus_dorsi_1(4) Latissimus_dorsi_2(5) Latissimus_dorsi_3(6) + #Teres_minor(7) PECM1(8) PECM2(9) PECM3(10) Coracobrachialis(11) TRIlong(12) TRIlat(13) TRImed(14) BIClong(15) BICshort(16) BRA(17) + # is different from muscle groups: + # Sh ext = DELT3, Infraspinatus, Latissimus_dorsi_1-3, Teres_minor + # Sh flex = PECM1, PECM2, PECM3, DELT1, Coracobrachialis + # El ext = TRIlong, TRIlat, TRImed + # El flex = BIClong, BICshort, BRA + shext=[2,3,4,5,6,7] + shflex=[0,8,9,10,11] + elext=[12,13,14] + elflex=[15,16,17] + musList=[shext, shflex,elext,elflex] + + for iMus in range(numMuscles): + ax = fig1.add_subplot(3,2,iMus+3) + # set number of muscle branches - assume each node has 3 branches (treat the Brachialis as a branch of Biceps=elbow flexor) + #iBranches = 3 + + ### Excitation and Activation #### + # equivalent for all branches of same muscle group + offset = 2 # use offset 3 because only DELT3 is used (order of muscle branches doesn't correspond muscle groups!) + ax.plot(T, musExcSeq[musList[iMus][offset],t1Samples:t2Samples],color=color1,linestyle=line1, linewidth=linWidth, label="excitation") + ax.plot(T, musActSeq[musList[iMus][offset],t1Samples:t2Samples],color=color2,linestyle=line1, linewidth=linWidth, label="activation") + + # for show branches plot individual branches and mean value for force and length + if showBranches: + ### Force and Length ### + for iBranch in range(len(musList[iMus])): + ax.plot(T, musForcesSeqNorm[musList[iMus][iBranch],t1Samples:t2Samples],color=color3,linestyle=line2, linewidth=linWidth-1) + ax.plot(T, musLengthsSeq[musList[iMus][iBranch],t1Samples:t2Samples],color=color4,linestyle=line2, linewidth=linWidth-1) + ax.plot(T, musForcesSeqNorm[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color3,linestyle=line1, linewidth=linWidth, label="force (mean)") + ax.plot(T, musLengthsSeq[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color4,linestyle=line1, linewidth=linWidth, label="length (mean)") + # for NOT show branches show mean value for force and single value for length + else: + ### Force ### + # For shoulder extensor group show only posterior Deltoid, branch 3 (DELT3) or Infraspinatus (INFSP) + # branch 2 (DELT2 = lateral deltoid) also available but currently not included in shoulder extensor group + if iMus == 0: + offset = 2 # DELT3 + #offset = 12 # INFSP + ax.plot(T, musForcesSeqNorm[musList[iMus][offset],t1Samples:t2Samples],color=color3,linestyle=line1, linewidth=linWidth, label="force") + # For rest of muscles use mean value of all branches + else: + offset=0 + ax.plot(T, musForcesSeqNorm[musList[iMus],t1Samples:t2Samples].mean(axis=0),color=color3,linestyle=line1, linewidth=linWidth, label="force") + + ### Length #### + # show length only of one muscle indicated by the index 'offset' - DELT3, PECM1, BIClong, TRIlong + maxLength = 0.20 + ax.plot(T, musLengthsSeq[musList[iMus][offset],t1Samples:t2Samples]/maxLength,color=color4,linestyle=line1, linewidth=linWidth, label="length") + + # show branche label + if (showBranches): + ax.plot(-1, -1,color=color3,linestyle=line2, linewidth=linWidth, label="force (branches)") + ax.plot(-1, -1,color=color4,linestyle=line2, linewidth=linWidth, label="length (branches)") + + # axis properties + ax.set_ylabel('normalized value', fontsize = legFont) + ax.set_ylim([0,1]) + ax.set_xlim([t1,t2]) + ax.set_xlabel('time (sec)', fontsize = legFont) + ax.set_title(muscleNames[iMus]) + if iMus==3: + ax.legend(loc='upper center', bbox_to_anchor=(-0.2, 1.8), borderaxespad=0., prop={'size':legFont}) + ax.grid(gridOn) + + # show graphs + fig1.tight_layout() + show() + + # save graphs using startPos and pattern in filename + if saveGraphs: + saveFolder = 'gif/' + fig1.savefig(saveFolder+saveName, bbox_inches=0) + + #if armAnimation: + #armAnim.save('test.mp4') + #armAnim.save(saveFolder+saveName+'.mp4',writer = writer) + + +# run single test (udp transfer, read files, plot graphs) def readAndPlot(jointAnglesSeq, musLengthsSeq, msmFolder, armAnimation, saveGraphs, saveName, timeRange, msecInterval): - # Sim parameters - #armAnimation = 1 # show 2D arm animation - #saveGraphs = 1 # save graph and animation - - # define time interval to display - #timeInterval = [0.1, 30] - - # Send muscle excitations to MSM and receive joint angles and muscle lengths - #jointAnglesSeq, musLengthsSeq = sendAndReceiveMsmData(initJointAngles, musExcSeq, readSimFromFile) - - # Read data from .pnt files - jointPosSeq,musExcSeq, musActSeq, musForcesSeq = readPntFiles(msmFolder, timeRange[1], msecInterval) - - # Plot results (last 2 arguments = initial and end times in seconds) - plotGraphs(jointPosSeq, jointAnglesSeq, musLengthsSeq, musExcSeq, musActSeq, musForcesSeq, timeRange[0], timeRange[1], msecInterval, armAnimation, saveGraphs, saveName) + # Sim parameters + #armAnimation = 1 # show 2D arm animation + #saveGraphs = 1 # save graph and animation + + # define time interval to display + #timeInterval = [0.1, 30] + + # Send muscle excitations to MSM and receive joint angles and muscle lengths + #jointAnglesSeq, musLengthsSeq = sendAndReceiveMsmData(initJointAngles, musExcSeq, readSimFromFile) + + # Read data from .pnt files + jointPosSeq,musExcSeq, musActSeq, musForcesSeq = readPntFiles(msmFolder, timeRange[1], msecInterval) + + # Plot results (last 2 arguments = initial and end times in seconds) + plotGraphs(jointPosSeq, jointAnglesSeq, musLengthsSeq, musExcSeq, musActSeq, musForcesSeq, timeRange[0], timeRange[1], msecInterval, armAnimation, saveGraphs, saveName) ############################## # Main script @@ -351,11 +351,10 @@ def readAndPlot(jointAnglesSeq, musLengthsSeq, msmFolder, armAnimation, saveGrap jointAnglesSeq = zeros((numJoints, n)) musLengthsSeq = zeros((numMusBranches, n)) armAnimation = 1 # show 2D arm animation -saveGraphs = 1 # save graph and animation +saveGraphs = 1 # save graph and animation timeInterval = [0.1, 30] saveName='temp' msmFolder = "/home/salvadord/Documents/ISB/Models_linux/msarm/source/test/" readAndPlot(jointAnglesSeq, musLengthsSeq, msmFolder, armAnimation, saveGraphs, saveName, timeInterval) ''' - diff --git a/arminterface.py b/arminterface.py index 80a350e..2b2a513 100644 --- a/arminterface.py +++ b/arminterface.py @@ -1,20 +1,20 @@ -# arminterface_pipe.py -- Python code for interfacing the sim with +# arminterface_pipe.py -- Python code for interfacing the sim with # the musculoskeletal arm using pipes (code simplified from arminterface.py) -# +# # Last update: 8/29/14 (salvadordura@gmail.com) -from socket import * -import select # for socket reading -import struct # to pack messages for socket comms -from time import * +from socket import * +import select # for socket reading +import struct # to pack messages for socket comms +from time import * import numpy import pickle -import subprocess # to tun the muscskel arm code +import subprocess # to tun the muscskel arm code import os # to kill subprocess import signal # to kill subprocess import xml.etree.ElementTree as ET import atexit -from cStringIO import StringIO +from io import StringIO import random import os.path @@ -36,21 +36,21 @@ msmFeedback = 'both' # Flag to choose whether to send joint position ('pos'); joint velocity ('vel') to msm; or joint velocity calculated from joint positions ('vel2') -wamForwardType = 'pos' +wamForwardType = 'pos' # MSM returns length of all branches per muscle - decide which to use as proprioceptive info for BMM (0 = mean); one value per muscle muscleLengthBranch = [3,1,1,1] -# Flag to run MSM from Python +# Flag to run MSM from Python msmRun = 1 #msmFolder = "/home/salvadord/Documents/ISB/Models_linux/msarm/source/test/" msmFolder = "msarm/" #update to use os.getcwd() -# .osim file with arm model description +# .osim file with arm model description # will be copied with timestamp in 'osimFile' to enable multiple instances running simultaneously osimOriginal = msmFolder + "SUNY_arm_2DOFs_horizon.osim" osimFile = osimOriginal - + # Flag to show MSM animation msmAnim = 0 @@ -58,7 +58,7 @@ msmGraphs = 0 #if msmGraphs or saveDataMuscles: -import armGraphs # to plot muscskel arm graphs +import armGraphs # to plot muscskel arm graphs class PipeReader(object): @@ -77,481 +77,479 @@ def readlines(self): # Initialize and setup the sockets, data format and other variables # def setup(secLength, msecInterval, shInit, elInit, targetx, targety, damping, shExtGain, shFlexGain, elExtGain, elFlexGain): - # input variables - global verbose - global savedDataSent - global savedDataReceived - - # output variables - global msmCommand - global msmPipe - global jointAnglesSeq - global musLengthsSeq - global packetID - global jointAngleID - global muscleLengthID - global time1 - global xmlFile - global osimFile - global msmSockMaxIter - global proc_stdout - global armReady - global inputbuffer - global pntFile - - - # Packet ID numbers (0 is first packet, with sim start time) - packetID = 0 - jointAngleID = -1 - muscleLengthID = -1 - - - # if plot MSM graphs initialize required arrays (numJoints, numMusBranches and n declared in armGraphs.py) - if msmGraphs: - # Create arrays to store received data - n = int(secLength*1000/msecInterval) # calculate number of samples - #n = n - 1 - jointAnglesSeq = numpy.zeros((armGraphs.numJoints, n)) - musLengthsSeq = numpy.zeros((armGraphs.numMusBranches, n)) - - - # Run MSM simulation asynchronously (i.e. Python goes on while MSM is running) - if msmRun: - # Set paths to run MSM from Python - msmCommand = msmFolder+"Run" # runs MSM sim - if msmAnim: - xmlOriginal = msmFolder+"SUNY_arm_horizon_fwd_10ms.xml" # xml file with sim parametersd - with 3D visualization - else: - xmlOriginal = msmFolder+"SUNY_arm_horizon_fwd_no_visual_10ms.xml" # xml file with sim parameters - no visualization - - # make copy of .xml and .osim to enable running multiple isntances simultaneously - while 1: - random.seed() - rand=int(random.random()*100000000) - xmlTemp = "xml_temp_"+str(rand)+".xml" - xmlFile = msmFolder+xmlTemp - if not os.path.isfile(xmlFile): # if file doesn't exist, exit loop; otherwise try different filename - break - - osimTemp = "osim_temp_"+str(rand)+".osim" - osimFile = msmFolder+osimTemp - cp_xml_str='cp %s %s' % (xmlOriginal, xmlFile) - cp_osim_str = 'cp %s %s' % (osimOriginal, osimFile) - os.system(cp_xml_str) - os.system(cp_osim_str) - - # generate name of temporary pnt file to store muscle data - pntTemp = "muscleData_temp_"+str(rand)+".pnt" - pntFile = msmFolder+pntTemp - - # set temporary .pnt file name in temporary .xml file - msmSetPntFileName(xmlFile, pntFile) - - # set temporary .osim file name in temporary .xml file - msmSetOsimFileName(xmlFile, osimTemp) - - # set msm duration = sim duration via XML - msmSetDurationXML(secLength, xmlFile) - - # set initial joint angles via XML - msmSetJointAnglesXML(shInit,elInit) - - # set damping - msmSetDampingXML(damping) - - # set max isometric force - shMuscleList = ['DELT1', 'DELT2', 'DELT3', 'Infraspinatus', 'Latissimus_dorsi_1', 'Latissimus_dorsi_2', 'Latissimus_dorsi_3','Teres_minor', 'PECM1', 'PECM2', 'PECM3', 'Coracobrachialis'] # only shoulder muscles - shMuscleOriginalValues = [shFlexGain*1142.60, shFlexGain*1142.60, shExtGain*259.88, shExtGain*1210.84, shExtGain*389.10, shExtGain*389.10, 281.66, shExtGain*354.25, 364.41, 515.41, 390.55, 242.46] - - elMuscleList = ['TRIlong', 'TRIlat', 'TRImed', 'BIClong', 'BICshort', 'BRA'] # only elbow muscles - elMuscleOriginalValues = [elExtGain*798.52, elExtGain*624.30, elExtGain*624.30, elFlexGain*624.30, elFlexGain*435.56, elFlexGain*987.26] - - for i in range(len(shMuscleList)): - msmSetMaxIsometricForceXML(shMuscleList[i], shMuscleOriginalValues[i]*3.0) #*3 - - for i in range(len(elMuscleList)): - msmSetMaxIsometricForceXML(elMuscleList[i], elMuscleOriginalValues[i]*1.0) - - # set target location - msmSetTargetPositionXML(targetx,targety) - - # run MSM and create pipe to read output - msmPipe = subprocess.Popen([msmCommand, xmlFile], preexec_fn=os.setsid, stdin=subprocess.PIPE, stdout=subprocess.PIPE) - proc_stdout = PipeReader(msmPipe.stdout.fileno()) - - # initialize timeInterval - time1 = time() - - # Flag to ensure first pacekt is sent once the virtual arm executable is ready - armReady = 0 - - # Buffer to store data during pipe communication - inputbuffer = StringIO() - - # save data - if saveDataExchanged: - savedDataSent = [] - savedDataReceived = [] + # input variables + global verbose + global savedDataSent + global savedDataReceived + + # output variables + global msmCommand + global msmPipe + global jointAnglesSeq + global musLengthsSeq + global packetID + global jointAngleID + global muscleLengthID + global time1 + global xmlFile + global osimFile + global msmSockMaxIter + global proc_stdout + global armReady + global inputbuffer + global pntFile + + + # Packet ID numbers (0 is first packet, with sim start time) + packetID = 0 + jointAngleID = -1 + muscleLengthID = -1 + + + # if plot MSM graphs initialize required arrays (numJoints, numMusBranches and n declared in armGraphs.py) + if msmGraphs: + # Create arrays to store received data + n = int(secLength*1000/msecInterval) # calculate number of samples + #n = n - 1 + jointAnglesSeq = numpy.zeros((armGraphs.numJoints, n)) + musLengthsSeq = numpy.zeros((armGraphs.numMusBranches, n)) + + + # Run MSM simulation asynchronously (i.e. Python goes on while MSM is running) + if msmRun: + # Set paths to run MSM from Python + msmCommand = msmFolder+"Run" # runs MSM sim + if msmAnim: + xmlOriginal = msmFolder+"SUNY_arm_horizon_fwd_10ms.xml" # xml file with sim parametersd - with 3D visualization + else: + xmlOriginal = msmFolder+"SUNY_arm_horizon_fwd_no_visual_10ms.xml" # xml file with sim parameters - no visualization + + # make copy of .xml and .osim to enable running multiple isntances simultaneously + while 1: + random.seed() + rand=int(random.random()*100000000) + xmlTemp = "xml_temp_"+str(rand)+".xml" + xmlFile = msmFolder+xmlTemp + if not os.path.isfile(xmlFile): # if file doesn't exist, exit loop; otherwise try different filename + break + + osimTemp = "osim_temp_"+str(rand)+".osim" + osimFile = msmFolder+osimTemp + cp_xml_str='cp %s %s' % (xmlOriginal, xmlFile) + cp_osim_str = 'cp %s %s' % (osimOriginal, osimFile) + os.system(cp_xml_str) + os.system(cp_osim_str) + + # generate name of temporary pnt file to store muscle data + pntTemp = "muscleData_temp_"+str(rand)+".pnt" + pntFile = msmFolder+pntTemp + + # set temporary .pnt file name in temporary .xml file + msmSetPntFileName(xmlFile, pntFile) + + # set temporary .osim file name in temporary .xml file + msmSetOsimFileName(xmlFile, osimTemp) + + # set msm duration = sim duration via XML + msmSetDurationXML(secLength, xmlFile) + + # set initial joint angles via XML + msmSetJointAnglesXML(shInit,elInit) + + # set damping + msmSetDampingXML(damping) + + # set max isometric force + shMuscleList = ['DELT1', 'DELT2', 'DELT3', 'Infraspinatus', 'Latissimus_dorsi_1', 'Latissimus_dorsi_2', 'Latissimus_dorsi_3','Teres_minor', 'PECM1', 'PECM2', 'PECM3', 'Coracobrachialis'] # only shoulder muscles + shMuscleOriginalValues = [shFlexGain*1142.60, shFlexGain*1142.60, shExtGain*259.88, shExtGain*1210.84, shExtGain*389.10, shExtGain*389.10, 281.66, shExtGain*354.25, 364.41, 515.41, 390.55, 242.46] + + elMuscleList = ['TRIlong', 'TRIlat', 'TRImed', 'BIClong', 'BICshort', 'BRA'] # only elbow muscles + elMuscleOriginalValues = [elExtGain*798.52, elExtGain*624.30, elExtGain*624.30, elFlexGain*624.30, elFlexGain*435.56, elFlexGain*987.26] + + for i in range(len(shMuscleList)): + msmSetMaxIsometricForceXML(shMuscleList[i], shMuscleOriginalValues[i]*3.0) #*3 + + for i in range(len(elMuscleList)): + msmSetMaxIsometricForceXML(elMuscleList[i], elMuscleOriginalValues[i]*1.0) + + # set target location + msmSetTargetPositionXML(targetx,targety) + + # run MSM and create pipe to read output + msmPipe = subprocess.Popen([msmCommand, xmlFile], preexec_fn=os.setsid, stdin=subprocess.PIPE, stdout=subprocess.PIPE) + proc_stdout = PipeReader(msmPipe.stdout.fileno()) + + # initialize timeInterval + time1 = time() + + # Flag to ensure first pacekt is sent once the virtual arm executable is ready + armReady = 0 + + # Buffer to store data during pipe communication + inputbuffer = StringIO() + + # save data + if saveDataExchanged: + savedDataSent = [] + savedDataReceived = [] def sendAndReceiveDataPackets(simtime, msecInterval, data1, data2, data3, data4): - # input variables - global packetID - global muscleLengthID - global jointAngleID - global verbose - global anglesReceived - global savedDataSent - global savedDataReceived - global jointAnglesSeq - global musLengthsSeq - global muscleLengthBranch - global wamForwardType - global packetID - global jointAngleID - global time1 - global msmReady - global proc_stdout - global armReady - global inputbuffer - - # concatenate input arguments into a list - data = [data1, data2,data3,data4] - - # Increase packet ID. - packetID += 1.0 - - ##################### - # Send packets - ##################### - - # Send packets to MSM - - while not armReady: # Ensure virtual arm is ready to receive - ready, _, _ = select.select([proc_stdout], [], [], 0.0) - if not ready: - continue - lines = ready[0].readlines() - if verbose: - print lines - if 'READY TO RUN' in lines: - armReady=1 - - musclesExcSend = data - - try: - msmPipe.stdin.write(str(musclesExcSend)[1:-1]) - msmPipe.stdin.write("\n") - if verbose: - print "\nWriting to MSM pipe: packetID=%f"%(packetID) - print(str(musclesExcSend)) - except: - print "timeout while sending packet to msarm" - dataReceived = [-3]*numJoints # error code indicating virtual arm pipe is not available - dataReceived2 = [0]*numMuscles - - ##################### - # Receive packets - ##################### - - # Receive packets from MSM - dataReceived = [] - dataReceived2 = [] - dataReceivedPacked = [] - dataReceivedPacked2 = [] - numJoints = 2 - numMuscles = 18 - - while len(dataReceived) == 0 or len(dataReceived2) == 0: - buffline = inputbuffer.readline() # read next line in buffer - if buffline == '': # if buffer empty, read more lines from pipe - ready, _, _ = select.select([proc_stdout], [], [], 1.0) - if ready: - inputbuffer = StringIO(ready[0].readlines()) - buffline = inputbuffer.readline() - if 'Totoal' in buffline: # if end of input from virtual arm, stop - break - tmp = [float(x) for x in buffline.split()] # split line into space delimited floats - - if verbose: - print "read from msm pipe: " + str(len(buffline)) - print tmp - - if len(tmp) == numJoints: - dataReceived = tmp - elif len(tmp) == numMuscles: - dataReceived2 = tmp - - if len(dataReceived) == 0 and packetID >= int(simtime/msecInterval): - dataReceived = [-3]*numJoints # error code in case missing last packet - dataReceived2 = [0]*numMuscles - if verbose: - print "last packet: returning -3 to prevent error" - - if dataReceived2 == []: dataReceived2 = [0]*numMuscles # if last message read and dataReceived2 empty, fill with 0s - - # store data in array required to plot msm graphs (armGraphs.py) - if msmGraphs: - muscleLengthID += 1 - musLengthsSeq[:, muscleLengthID] = numpy.array(dataReceived2) - - # print received data - if verbose: - print "Received packet %f from MSM: (%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f)" % \ - (packetID, dataReceived2[0],dataReceived2[1],dataReceived2[2],dataReceived2[3],dataReceived2[4],dataReceived2[5],dataReceived2[6],dataReceived2[7],dataReceived2[8],dataReceived2[9],dataReceived2[10]) - - # store the mean value of the length of branches of each muscle in a new variable - #DELT1(0) DELT2(1) DELT3(2) Infraspinatus(3) Latissimus_dorsi_1(4) Latissimus_dorsi_2(5) Latissimus_dorsi_3(6) - #Teres_minor(7) PECM1(8) PECM2(9) PECM3(10) Coracobrachialis(11) TRIlong(12) TRIlat(13) TRImed(14) BIClong(15) BICshort(16) BRA(17) - - # Deltoid muscle length = mean(DELT1 DELT2 DELT3) - if muscleLengthBranch[0] == 0: - lengthsReceived[0] = numpy.mean(dataReceived2[0:3]) - # Deltoid muscle length = length of branch number 'muscleLengthBranch' - else: - lengthsReceived[0] = dataReceived2[muscleLengthBranch[0]-1] - - # Pectoralis muscle length = mean(PECM1 PECM2 PECM3) - if muscleLengthBranch[1] == 0: - lengthsReceived[1] = numpy.mean(dataReceived2[8:11]) - # Pectoralis muscle length = length of branch number 'muscleLengthBranch' - else: - lengthsReceived[1] = dataReceived2[muscleLengthBranch[1]+7] - - # Triceps muscle length = mean(TRIlong TRIlat TRImed) - if muscleLengthBranch[2] == 0: - lengthsReceived[2] = numpy.mean(dataReceived2[12:15]) - # Triceps muscle length = length of branch number 'muscleLengthBranch' - else: - lengthsReceived[2] = dataReceived2[muscleLengthBranch[2]+11] - - # Biceps/Brachialis muscle length = mean(BIClong BICshort BRA) - if muscleLengthBranch[3] == 0: - lengthsReceived[3] = numpy.mean(dataReceived2[15:18]) - # Biceps/Brachialis muscle length = length of branch number 'muscleLengthBranch' - else: - lengthsReceived[3] = dataReceived2[muscleLengthBranch[3]+14] - - # save received data - if saveDataExchanged: - savedDataReceived.append([simtime, getCurrTime(), lengthsReceived[0], lengthsReceived[1], lengthsReceived[2], lengthsReceived[3]]) - - ###################### - # Receive joint angles - - # store data in array required to plot msm graphs (armGraphs.py) - if msmGraphs: - jointAngleID += 1 - jointAnglesSeq[:, jointAngleID] = numpy.array(dataReceived[0:2]) - - # print received data - if verbose: - print "Received packet %f from MSM: (%.3f,%.3f)" % (packetID, dataReceived[0],dataReceived[1]) - - # if receiving joint angles, store the received shoulder and elbow angles in a new variable - anglesReceived = dataReceived - - # save received data - if saveDataExchanged: - savedDataReceived.append([simtime, getCurrTime(), anglesReceived[0], anglesReceived[1]]) - - return dataReceived + # input variables + global packetID + global muscleLengthID + global jointAngleID + global verbose + global anglesReceived + global savedDataSent + global savedDataReceived + global jointAnglesSeq + global musLengthsSeq + global muscleLengthBranch + global wamForwardType + global packetID + global jointAngleID + global time1 + global msmReady + global proc_stdout + global armReady + global inputbuffer + + # concatenate input arguments into a list + data = [data1, data2,data3,data4] + + # Increase packet ID. + packetID += 1.0 + + ##################### + # Send packets + ##################### + + # Send packets to MSM + + while not armReady: # Ensure virtual arm is ready to receive + ready, _, _ = select.select([proc_stdout], [], [], 0.0) + if not ready: + continue + lines = ready[0].readlines() + if verbose: + print(lines) + if 'READY TO RUN' in lines: + armReady=1 + + musclesExcSend = data + + try: + msmPipe.stdin.write(str(musclesExcSend)[1:-1]) + msmPipe.stdin.write("\n") + if verbose: + print("\nWriting to MSM pipe: packetID=%f"%(packetID)) + print((str(musclesExcSend))) + except: + print("timeout while sending packet to msarm") + dataReceived = [-3]*numJoints # error code indicating virtual arm pipe is not available + dataReceived2 = [0]*numMuscles + + ##################### + # Receive packets + ##################### + + # Receive packets from MSM + dataReceived = [] + dataReceived2 = [] + dataReceivedPacked = [] + dataReceivedPacked2 = [] + numJoints = 2 + numMuscles = 18 + + while len(dataReceived) == 0 or len(dataReceived2) == 0: + buffline = inputbuffer.readline() # read next line in buffer + if buffline == '': # if buffer empty, read more lines from pipe + ready, _, _ = select.select([proc_stdout], [], [], 1.0) + if ready: + inputbuffer = StringIO(ready[0].readlines()) + buffline = inputbuffer.readline() + if 'Totoal' in buffline: # if end of input from virtual arm, stop + break + tmp = [float(x) for x in buffline.split()] # split line into space delimited floats + + if verbose: + print("read from msm pipe: " + str(len(buffline))) + print(tmp) + + if len(tmp) == numJoints: + dataReceived = tmp + elif len(tmp) == numMuscles: + dataReceived2 = tmp + + if len(dataReceived) == 0 and packetID >= int(simtime/msecInterval): + dataReceived = [-3]*numJoints # error code in case missing last packet + dataReceived2 = [0]*numMuscles + if verbose: + print("last packet: returning -3 to prevent error") + + if dataReceived2 == []: dataReceived2 = [0]*numMuscles # if last message read and dataReceived2 empty, fill with 0s + + # store data in array required to plot msm graphs (armGraphs.py) + if msmGraphs: + muscleLengthID += 1 + musLengthsSeq[:, muscleLengthID] = numpy.array(dataReceived2) + + # print received data + if verbose: + print("Received packet %f from MSM: (%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f)" % \ +(packetID, dataReceived2[0],dataReceived2[1],dataReceived2[2],dataReceived2[3],dataReceived2[4],dataReceived2[5],dataReceived2[6],dataReceived2[7],dataReceived2[8],dataReceived2[9],dataReceived2[10])) + + # store the mean value of the length of branches of each muscle in a new variable + #DELT1(0) DELT2(1) DELT3(2) Infraspinatus(3) Latissimus_dorsi_1(4) Latissimus_dorsi_2(5) Latissimus_dorsi_3(6) + #Teres_minor(7) PECM1(8) PECM2(9) PECM3(10) Coracobrachialis(11) TRIlong(12) TRIlat(13) TRImed(14) BIClong(15) BICshort(16) BRA(17) + + # Deltoid muscle length = mean(DELT1 DELT2 DELT3) + if muscleLengthBranch[0] == 0: + lengthsReceived[0] = numpy.mean(dataReceived2[0:3]) + # Deltoid muscle length = length of branch number 'muscleLengthBranch' + else: + lengthsReceived[0] = dataReceived2[muscleLengthBranch[0]-1] + + # Pectoralis muscle length = mean(PECM1 PECM2 PECM3) + if muscleLengthBranch[1] == 0: + lengthsReceived[1] = numpy.mean(dataReceived2[8:11]) + # Pectoralis muscle length = length of branch number 'muscleLengthBranch' + else: + lengthsReceived[1] = dataReceived2[muscleLengthBranch[1]+7] + + # Triceps muscle length = mean(TRIlong TRIlat TRImed) + if muscleLengthBranch[2] == 0: + lengthsReceived[2] = numpy.mean(dataReceived2[12:15]) + # Triceps muscle length = length of branch number 'muscleLengthBranch' + else: + lengthsReceived[2] = dataReceived2[muscleLengthBranch[2]+11] + + # Biceps/Brachialis muscle length = mean(BIClong BICshort BRA) + if muscleLengthBranch[3] == 0: + lengthsReceived[3] = numpy.mean(dataReceived2[15:18]) + # Biceps/Brachialis muscle length = length of branch number 'muscleLengthBranch' + else: + lengthsReceived[3] = dataReceived2[muscleLengthBranch[3]+14] + + # save received data + if saveDataExchanged: + savedDataReceived.append([simtime, getCurrTime(), lengthsReceived[0], lengthsReceived[1], lengthsReceived[2], lengthsReceived[3]]) + + ###################### + # Receive joint angles + + # store data in array required to plot msm graphs (armGraphs.py) + if msmGraphs: + jointAngleID += 1 + jointAnglesSeq[:, jointAngleID] = numpy.array(dataReceived[0:2]) + + # print received data + if verbose: + print("Received packet %f from MSM: (%.3f,%.3f)" % (packetID, dataReceived[0],dataReceived[1])) + + # if receiving joint angles, store the received shoulder and elbow angles in a new variable + anglesReceived = dataReceived + + # save received data + if saveDataExchanged: + savedDataReceived.append([simtime, getCurrTime(), anglesReceived[0], anglesReceived[1]]) + + return dataReceived # Function to close sockets, save data and plot graphs def closeSavePlot(secLength, msecInterval, filestem=''): - global saveDataMuscles - global savedDataReceived - global savedDataSent - global jointAnglesSeq - global musLengthsSeq - global msmPipe - global msmRun - global xmlFile - global osimFile - global pntFile - - - # delete temporal copy of xml and osim files - rm_str = 'rm %s %s' % (xmlFile, osimFile) - os.system(rm_str) - - # close msm pipe - if msmRun: - os.killpg(msmPipe.pid, signal.SIGTERM) - - if verbose: - print "Sockets closed" - - # save data exchanged via udp to file - if saveDataExchanged: - f = open( 'data_BMM.txt', 'w' ) - for item in savedDataSent: - f.write("%s\n" % item) - for item in savedDataReceived: - f.write("%s\n" % item) - f.close() - if verbose: - print 'Data saved to file' - - # save muscle data to file - if saveDataMuscles: - msmFolder = '' # data files saved locally -# # Read data from .pnt files - jointPosSeq,musExcSeq, musActSeq, musForcesSeq = armGraphs.readPntFiles(msmFolder, pntFile, secLength, msecInterval) - with open("%s-muscles.p"%(filestem),'w') as f: - pickle.dump([musExcSeq, musActSeq, musForcesSeq], f) - - - # plot graphs - if msmGraphs: - saveName = 'msmGraphs.png' - saveGraphs = saveDataExchanged #set func argument to save graphs - armAnimation = 0# msmGraphs #set func argument to show arm animation - timeRange = [0.1, secLength] - msmFolder = '' # data files saved locally - armGraphs.readAndPlot(jointAnglesSeq, musLengthsSeq, msmFolder, armAnimation, saveGraphs, saveName, timeRange, msecInterval) - - # delete temporal pnt file - rm_str = 'rm %s' % (pntFile) - os.system(rm_str) + global saveDataMuscles + global savedDataReceived + global savedDataSent + global jointAnglesSeq + global musLengthsSeq + global msmPipe + global msmRun + global xmlFile + global osimFile + global pntFile + + + # delete temporal copy of xml and osim files + rm_str = 'rm %s %s' % (xmlFile, osimFile) + os.system(rm_str) + + # close msm pipe + if msmRun: + os.killpg(msmPipe.pid, signal.SIGTERM) + + if verbose: + print("Sockets closed") + + # save data exchanged via udp to file + if saveDataExchanged: + f = open( 'data_BMM.txt', 'w' ) + for item in savedDataSent: + f.write("%s\n" % item) + for item in savedDataReceived: + f.write("%s\n" % item) + f.close() + if verbose: + print('Data saved to file') + + # save muscle data to file + if saveDataMuscles: + msmFolder = '' # data files saved locally +# # Read data from .pnt files + jointPosSeq,musExcSeq, musActSeq, musForcesSeq = armGraphs.readPntFiles(msmFolder, pntFile, secLength, msecInterval) + with open("%s-muscles.p"%(filestem),'w') as f: + pickle.dump([musExcSeq, musActSeq, musForcesSeq], f) + + + # plot graphs + if msmGraphs: + saveName = 'msmGraphs.png' + saveGraphs = saveDataExchanged #set func argument to save graphs + armAnimation = 0# msmGraphs #set func argument to show arm animation + timeRange = [0.1, secLength] + msmFolder = '' # data files saved locally + armGraphs.readAndPlot(jointAnglesSeq, musLengthsSeq, msmFolder, armAnimation, saveGraphs, saveName, timeRange, msecInterval) + + # delete temporal pnt file + rm_str = 'rm %s' % (pntFile) + os.system(rm_str) # Time function def getCurrTime(): - return time()*1000 + return time()*1000 # # Functions to return data received via udp (called from NEURON) # def getAngleReceived(joint): - return anglesReceived[int(joint)] + return anglesReceived[int(joint)] def getLengthReceived(muscle): - return lengthsReceived[int(muscle)] + return lengthsReceived[int(muscle)] def getPacketLoss(): - if packetID==-1: - return 1 - else: - return 0 + if packetID==-1: + return 1 + else: + return 0 # Function to set joint angles by modifying the XML .osim file- also runs the MSM program in pipe def msmSetJointAnglesXML(shAng,elAng): - global simtime - global msmCommand - - # read .osim file - tree = ET.parse(osimFile) - root = tree.getroot() - - # set shoulder angle initial value - for shoulder_flex_tag in root.iter('Coordinate'): - if shoulder_flex_tag.get('name') == 'arm_flex': - break - - init_shoulder_tag = shoulder_flex_tag.find('initial_value') - init_shoulder_tag.text = str(shAng); - - # set elbow angle initial value - for elbow_flex_tag in root.iter('Coordinate'): - if elbow_flex_tag.get('name') == 'elbow_flex': - break - - init_elbow_tag = elbow_flex_tag.find('initial_value') - init_elbow_tag.text = str(elAng); - - # save xml file - tree.write(osimFile) - + global simtime + global msmCommand + + # read .osim file + tree = ET.parse(osimFile) + root = tree.getroot() + + # set shoulder angle initial value + for shoulder_flex_tag in root.iter('Coordinate'): + if shoulder_flex_tag.get('name') == 'arm_flex': + break + + init_shoulder_tag = shoulder_flex_tag.find('initial_value') + init_shoulder_tag.text = str(shAng); + + # set elbow angle initial value + for elbow_flex_tag in root.iter('Coordinate'): + if elbow_flex_tag.get('name') == 'elbow_flex': + break + + init_elbow_tag = elbow_flex_tag.find('initial_value') + init_elbow_tag.text = str(elAng); + + # save xml file + tree.write(osimFile) + def msmSetOsimFileName(xmlFile, osimTemp): - # read .osim file - - tree = ET.parse(xmlFile) - root = tree.getroot() - - # set osim file name - for tmp in root.iter('OsimFile'): - tmp.set('name', osimTemp) - - # save xml file - tree.write(xmlFile) + # read .osim file + + tree = ET.parse(xmlFile) + root = tree.getroot() + + # set osim file name + for tmp in root.iter('OsimFile'): + tmp.set('name', osimTemp) + + # save xml file + tree.write(xmlFile) def msmSetPntFileName(xmlFile, pntFile): - # read .osim file - - tree = ET.parse(xmlFile) - root = tree.getroot() - - # set osim file name - for tmp in root.iter('PntOutput'): - tmp.set('name', pntFile) - break # exit after first appearance - - # save xml file - tree.write(xmlFile) + # read .osim file + + tree = ET.parse(xmlFile) + root = tree.getroot() + + # set osim file name + for tmp in root.iter('PntOutput'): + tmp.set('name', pntFile) + break # exit after first appearance + + # save xml file + tree.write(xmlFile) # set duration of arm sim in .xml file def msmSetDurationXML(secLength, xmlFile): - # read .osim file - tree = ET.parse(xmlFile) - root = tree.getroot() - - # set end time - endTime = root.find('EndTime') - endTime.text = str(secLength) - - # save xml file - tree.write(xmlFile) - + # read .osim file + tree = ET.parse(xmlFile) + root = tree.getroot() + + # set end time + endTime = root.find('EndTime') + endTime.text = str(secLength) + + # save xml file + tree.write(xmlFile) + def msmSetDampingXML(damping): - # read .osim file - tree = ET.parse(osimFile) - root = tree.getroot() - - # set damping value - for damping_tag in root.iter('damping'): - damping_tag.text = str(damping) - - # save xml file - tree.write(osimFile) - + # read .osim file + tree = ET.parse(osimFile) + root = tree.getroot() + + # set damping value + for damping_tag in root.iter('damping'): + damping_tag.text = str(damping) + + # save xml file + tree.write(osimFile) + def msmSetTargetPositionXML(targetx, targety): - # read .osim file - tree = ET.parse(osimFile) - root = tree.getroot() - - for Body_tag in root.iter('Body'): - if Body_tag.get('name') == 'ground': - break - transform_tag = Body_tag.find('VisibleObject') - transform_tag = transform_tag.find('GeometrySet') - transform_tag = transform_tag.find('objects') - transform_tag = transform_tag.find('DisplayGeometry') - transform_tag = transform_tag.find('transform') - - transform_tag.text = str(0)+' '+str(0)+' '+str(0)+' '+str(-targetx-0.06)+' '+str(-0.05)+' '+str(-targety+0.09) # y postive = x-axis right - - # save xml file - tree.write(osimFile) - - + # read .osim file + tree = ET.parse(osimFile) + root = tree.getroot() + + for Body_tag in root.iter('Body'): + if Body_tag.get('name') == 'ground': + break + transform_tag = Body_tag.find('VisibleObject') + transform_tag = transform_tag.find('GeometrySet') + transform_tag = transform_tag.find('objects') + transform_tag = transform_tag.find('DisplayGeometry') + transform_tag = transform_tag.find('transform') + + transform_tag.text = str(0)+' '+str(0)+' '+str(0)+' '+str(-targetx-0.06)+' '+str(-0.05)+' '+str(-targety+0.09) # y postive = x-axis right + + # save xml file + tree.write(osimFile) + + def msmSetMaxIsometricForceXML(muscleName, value): - # read .osim file - tree = ET.parse(osimFile) - root = tree.getroot() - - # set shoulder angle initial value - for muscle_tag in root.iter('Schutte1993Muscle'): - if muscle_tag.get('name') == muscleName: - break - - muscleForce_tag = muscle_tag.find('max_isometric_force') - #print muscleForce_tag.text - muscleForce_tag.text = str(value); - - # save xml file - tree.write(osimFile) + # read .osim file + tree = ET.parse(osimFile) + root = tree.getroot() -def saveEMG(): - global saveDataMuscles - saveDataMuscles = 1 + # set shoulder angle initial value + for muscle_tag in root.iter('Schutte1993Muscle'): + if muscle_tag.get('name') == muscleName: + break + + muscleForce_tag = muscle_tag.find('max_isometric_force') + #print muscleForce_tag.text + muscleForce_tag.text = str(value); + # save xml file + tree.write(osimFile) +def saveEMG(): + global saveDataMuscles + saveDataMuscles = 1 diff --git a/dummyArm.py b/dummyArm.py index 01e053f..e1538db 100644 --- a/dummyArm.py +++ b/dummyArm.py @@ -1,68 +1,68 @@ # dummyArm.py -- Python code for interfacing the sim with a virtual arm -# +# # Last update: 07/21/14 (salvadordura@gmail.com) -import matplotlib +import matplotlib matplotlib.use('TkAgg') # for visualization -from socket import * -import select # for socket reading -import struct # to pack messages for socket comms +from socket import * +import select # for socket reading +import struct # to pack messages for socket comms import numpy from pylab import figure, show, ion, pause -# Initialize and setup the sockets and data format +# Initialize and setup the sockets and data format def setup(): - sendMsgFormat = "dd" # messages contain 2 doubles = joint angles of shoulder and elbow - receiveMsgFormat = "dd" # messages contain 2 doubles = velocities - receiveMsgSize = struct.calcsize(receiveMsgFormat) - localHostIP = "127.0.0.1"#"localhost"#"192.168.1.2"# # set IP for local connection - - print "Setting up connection..." # setup connection to model - sockReceive = socket(AF_INET, SOCK_DGRAM) # create sockets - hostPortReceive = 31000 # set port for receiving packets - sockReceive.bind(('', hostPortReceive)) # bind to port - sockReceive.setblocking(1) # Set blocking/non-blocking mode - print ("Created UDP socket to receive packets from NEURON model; binded to port %d"%hostPortReceive) - - sockSend = socket(AF_INET, SOCK_DGRAM) # connect to socket for sending packets - hostPortSend = 32000 # set port for sending packets - sockSend.connect((localHostIP, hostPortSend)) - sockSend.setblocking(1) # Set blocking/non-blocking mode - print ("Created UDP socket to send packets to NEURON model; socket connected to IP %s, port %d" % (localHostIP, hostPortSend)) - - return sendMsgFormat, receiveMsgFormat, sockReceive, sockSend + sendMsgFormat = "dd" # messages contain 2 doubles = joint angles of shoulder and elbow + receiveMsgFormat = "dd" # messages contain 2 doubles = velocities + receiveMsgSize = struct.calcsize(receiveMsgFormat) + localHostIP = "127.0.0.1"#"localhost"#"192.168.1.2"# # set IP for local connection + + print("Setting up connection...") # setup connection to model + sockReceive = socket(AF_INET, SOCK_DGRAM) # create sockets + hostPortReceive = 31000 # set port for receiving packets + sockReceive.bind(('', hostPortReceive)) # bind to port + sockReceive.setblocking(1) # Set blocking/non-blocking mode + print(("Created UDP socket to receive packets from NEURON model; binded to port %d"%hostPortReceive)) + + sockSend = socket(AF_INET, SOCK_DGRAM) # connect to socket for sending packets + hostPortSend = 32000 # set port for sending packets + sockSend.connect((localHostIP, hostPortSend)) + sockSend.setblocking(1) # Set blocking/non-blocking mode + print(("Created UDP socket to send packets to NEURON model; socket connected to IP %s, port %d" % (localHostIP, hostPortSend))) + + return sendMsgFormat, receiveMsgFormat, sockReceive, sockSend # Send and receive packets to/from virtual arm def sendAndReceivePackets(dataSend, sendMsgFormat, receiveMsgFormat, sockReceive, sockSend): - dataReceived = [0,0] - try: - receiveMsgSize = struct.calcsize(receiveMsgFormat) - dataReceivedPacked = sockReceive.recv(receiveMsgSize) # read packet from socket - if len(dataReceivedPacked) == receiveMsgSize: - dataReceived = struct.unpack(receiveMsgFormat, dataReceivedPacked) - print "Received packet from model: (%.2f,%.2f)" % (dataReceived[0],dataReceived[1]) - except: - print "Error receiving packet" - - inputready, outputready, e = select.select([] ,[sockSend],[], 0.0) # check if other side ready to receive - if len(outputready)>0: - try: - sent = sockSend.send(struct.pack(sendMsgFormat, dataSend[0], dataSend[1])) # send packet - print "Sent packet to virtual arm: (%.2f, %.2f)" % (dataSend[0], dataSend[1]) - except: - print "Sending socket ready but error sending packet" - else: - print "Sending socket not ready" - return dataReceived + dataReceived = [0,0] + try: + receiveMsgSize = struct.calcsize(receiveMsgFormat) + dataReceivedPacked = sockReceive.recv(receiveMsgSize) # read packet from socket + if len(dataReceivedPacked) == receiveMsgSize: + dataReceived = struct.unpack(receiveMsgFormat, dataReceivedPacked) + print("Received packet from model: (%.2f,%.2f)" % (dataReceived[0],dataReceived[1])) + except: + print("Error receiving packet") + + inputready, outputready, e = select.select([] ,[sockSend],[], 0.0) # check if other side ready to receive + if len(outputready)>0: + try: + sent = sockSend.send(struct.pack(sendMsgFormat, dataSend[0], dataSend[1])) # send packet + print("Sent packet to virtual arm: (%.2f, %.2f)" % (dataSend[0], dataSend[1])) + except: + print("Sending socket ready but error sending packet") + else: + print("Sending socket not ready") + return dataReceived # Main code for simple virtual arm duration = 4 # sec interval = 0.010 # time between packets (sec) -L1 = 1.0 # arm segment 1 length +L1 = 1.0 # arm segment 1 length L2 = 0.8 # arm segment 2 length -shang = numpy.pi/2 # shoulder angle (rad) -elang = numpy.pi/2 # elbow angle (rad) +shang = numpy.pi/2 # shoulder angle (rad) +elang = numpy.pi/2 # elbow angle (rad) shvel = 0 # shoulder velocity (rad/s) elvel = 0 # elbow velocity (rad/s) friction = 0.1 # friction coefficient @@ -75,29 +75,21 @@ def sendAndReceivePackets(dataSend, sendMsgFormat, receiveMsgFormat, sockReceiv ax.grid() line, = ax.plot([], [], 'o-', lw=2) -raw_input("Press Enter to continue...") +input("Press Enter to continue...") for i in numpy.arange(0, duration, interval): - shang = (shang + shvel * interval) % (2*numpy.pi) # update shoulder angle - elang = (elang + elvel * interval) % (2*numpy.pi)# update elbow angle - if shang<0: shang = 2*numpy.pi + shang - if elang<0: elang = 2*numpy.pi + shang - shpos = [L1*numpy.sin(shang), L1*numpy.cos(shang)] # calculate shoulder x-y pos - elpos = [L2*numpy.sin(shang+elang) + shpos[0], L2*numpy.cos(shang+elang) + shpos[1]] # calculate elbow x-y pos - - dataSend = [shang, elang] # set data to send - dataReceived = sendAndReceivePackets(dataSend, sendMsgFormat, receiveMsgFormat, sockReceive, sockSend) # send and receive data - shvel = shvel+dataReceived[0] - (friction * shvel)# update velocities based on incoming commands (accelerations) and friction - elvel = elvel+dataReceived[1] - (friction * elvel) - - line.set_data([0, shpos[0], elpos[0]], [0, shpos[1], elpos[1]]) # update line in figure - ax.set_title('Time = %.1f ms, shoulder: pos=%.2f rad, vel=%.2f, acc=%.2f ; elbow: pos = %.2f rad, vel = %.2f, acc=%.2f' % (float(i)*1000, shang, shvel, dataReceived[0] - (friction * shvel), elang, elvel, dataReceived[1] - (friction * elvel) ), fontsize=10) - show() - pause(0.0001) # pause so that the figure refreshes at every time step - - - - - - - - + shang = (shang + shvel * interval) % (2*numpy.pi) # update shoulder angle + elang = (elang + elvel * interval) % (2*numpy.pi)# update elbow angle + if shang<0: shang = 2*numpy.pi + shang + if elang<0: elang = 2*numpy.pi + shang + shpos = [L1*numpy.sin(shang), L1*numpy.cos(shang)] # calculate shoulder x-y pos + elpos = [L2*numpy.sin(shang+elang) + shpos[0], L2*numpy.cos(shang+elang) + shpos[1]] # calculate elbow x-y pos + + dataSend = [shang, elang] # set data to send + dataReceived = sendAndReceivePackets(dataSend, sendMsgFormat, receiveMsgFormat, sockReceive, sockSend) # send and receive data + shvel = shvel+dataReceived[0] - (friction * shvel)# update velocities based on incoming commands (accelerations) and friction + elvel = elvel+dataReceived[1] - (friction * elvel) + + line.set_data([0, shpos[0], elpos[0]], [0, shpos[1], elpos[1]]) # update line in figure + ax.set_title('Time = %.1f ms, shoulder: pos=%.2f rad, vel=%.2f, acc=%.2f ; elbow: pos = %.2f rad, vel = %.2f, acc=%.2f' % (float(i)*1000, shang, shvel, dataReceived[0] - (friction * shvel), elang, elvel, dataReceived[1] - (friction * elvel) ), fontsize=10) + show() + pause(0.0001) # pause so that the figure refreshes at every time step diff --git a/error.py b/error.py index 8f7e86c..0d380e3 100644 --- a/error.py +++ b/error.py @@ -4,49 +4,49 @@ def getTarget0err(f): - stdin,stdout = os.popen2("tail -n 4 "+f) - stdin.close() - lines = stdout.readlines(); stdout.close() - line = lines[0].split() - return float(line[5]) + stdin,stdout = os.popen2("tail -n 4 "+f) + stdin.close() + lines = stdout.readlines(); stdout.close() + line = lines[0].split() + return float(line[5]) def errorFromShell(filestem, readTarget0ErrFromRunFile, maxGens, maxCands, maxTargets): minVal = 1 minValtarg0 = 1 for igen in range(maxGens): - for icand in range(maxCands): - error = [] - for itarget in range(maxTargets): - try: - errfilename = '%s/gen_%s_cand_%d_target_%d_error' % (filestem,igen,icand,itarget) - runfilename = '%s/gen_%s_cand_%d.run' % (filestem,igen,icand) - with open(errfilename, 'r') as f: - error.append(pickle.load(f)) - print 'gen=%d, cand=%d, target=%d: error = %f' % (igen, icand, itarget, error[-1]) - except: - print('not found file: stem=%s, gen=%d, cand=%d, targ=%d' % (filestem,igen,icand,itarget)) - avgErr = pylab.mean(error) - if readTarget0ErrFromRunFile: - try: - target0err = getTarget0err(runfilename) - print "error target 0: %.2f \n" % (target0err) - if target0err < minValtarg0: - minValtarg0 = target0err - minCandtarg0 = icand - minGentarg0 = igen - except: - print "error reading file:", runfilename + for icand in range(maxCands): + error = [] + for itarget in range(maxTargets): + try: + errfilename = '%s/gen_%s_cand_%d_target_%d_error' % (filestem,igen,icand,itarget) + runfilename = '%s/gen_%s_cand_%d.run' % (filestem,igen,icand) + with open(errfilename, 'r') as f: + error.append(pickle.load(f)) + print('gen=%d, cand=%d, target=%d: error = %f' % (igen, icand, itarget, error[-1])) + except: + print(('not found file: stem=%s, gen=%d, cand=%d, targ=%d' % (filestem,igen,icand,itarget))) + avgErr = pylab.mean(error) + if readTarget0ErrFromRunFile: + try: + target0err = getTarget0err(runfilename) + print("error target 0: %.2f \n" % (target0err)) + if target0err < minValtarg0: + minValtarg0 = target0err + minCandtarg0 = icand + minGentarg0 = igen + except: + print("error reading file:", runfilename) - print "avg error: %.2f \n" % (avgErr) - if avgErr < minVal: - minVal = avgErr - minCand = icand - minGen = igen + print("avg error: %.2f \n" % (avgErr)) + if avgErr < minVal: + minVal = avgErr + minCand = icand + minGen = igen - print "Min error = %f ; gen = %d ; cand = %d \n" % (minVal, minGen, minCand) - if readTarget0ErrFromRunFile: - print "Min error (target 0) = %f ; gen = %d ; cand = %d \n" % (minValtarg0, minGentarg0, minCandtarg0) + print("Min error = %f ; gen = %d ; cand = %d \n" % (minVal, minGen, minCand)) + if readTarget0ErrFromRunFile: + print("Min error (target 0) = %f ; gen = %d ; cand = %d \n" % (minValtarg0, minGentarg0, minCandtarg0)) def errorFromPickle(filestem, maxGens, maxCands): from collections import OrderedDict @@ -54,27 +54,27 @@ def errorFromPickle(filestem, maxGens, maxCands): minVals = {'error0_pre': [1,0,0], 'error1_pre': [1,0,0], 'error_pre': [1,0,0], 'errord_pre': [1,0,0], \ 'error0_post': [1,0,0], 'error1_post': [1,0,0], 'error_post': [1,0,0], 'errord_post': [1,0,0], \ 'error0_lesion': [1,0,0], 'error1_lesion': [1,0,0], 'error_pre': [1,0,0], 'errord_tot': [1,0,0], 'error_fitness': [1,0,0]} - + for igen in range(maxGens): for icand in range(maxCands): error = [] try: errfilename = '%s/gen_%d_cand_%d_target_0_error' % (filestem,igen,icand) with open(errfilename, 'r') as f: - error = pickle.load(f) - #print 'gen=%d, cand=%d, errord_tot = %f, error_fitness = %f' % (igen, icand, error['errord_tot'], error['error_fitness']) - for key,val in error.iteritems(): - if val < minVals[key][0]: + error = pickle.load(f) + #print 'gen=%d, cand=%d, errord_tot = %f, error_fitness = %f' % (igen, icand, error['errord_tot'], error['error_fitness']) + for key,val in error.items(): + if val < minVals[key][0]: minVals[key][0] = val minVals[key][1] = igen minVals[key][2] = icand except: - print('not found file: stem=%s, gen=%d, cand=%d' % (filestem,igen,icand)) + print(('not found file: stem=%s, gen=%d, cand=%d' % (filestem,igen,icand))) + + for key,val in minVals.items(): + print("%s = %f ; gen = %d ; cand = %d \n" % (key, val[0], val[1], val[2])) + - for key,val in minVals.iteritems(): - print "%s = %f ; gen = %d ; cand = %d \n" % (key, val[0], val[1], val[2]) - - # set params filestem = '../data/15aug29_evolutionStrategy' @@ -88,5 +88,3 @@ def errorFromPickle(filestem, maxGens, maxCands): # call function to obtain errors from pickle file errorFromPickle(filestem, maxGens, maxCands) - - diff --git a/evol_islands.py b/evol_islands.py index 0e826de..70d8007 100644 --- a/evol_islands.py +++ b/evol_islands.py @@ -14,14 +14,14 @@ from popen2 import popen2 import pickle import multiprocessing -import Queue +import queue import subprocess ngen = -1 #global variable keeping number of generations ############################################################################### ### Simulation options -############################################################################### +############################################################################### evolAlgorithm = 'evolutionStrategyCross' #'diffEvolution' # 'evolutionStrategy' #'krichmarCustom' #'genetic'#'particleSwarm100'#'estimationDist' #'diffEvolution' # 'evolutionStrategy' # 'krichmarCustom', 'genetic' simdatadir = '../data/16jul26_'+evolAlgorithm # folder to save sim results @@ -31,7 +31,7 @@ max_migrants = 1 # migration_interval = 5 pop_size = 10 # population size per island -num_elites = pop_size/10 # num of top individuals kept each generation - maybe set to pop_size/10? +num_elites = pop_size/10 # num of top individuals kept each generation - maybe set to pop_size/10? max_generations = 2000 max_evaluations = max_generations * num_islands * pop_size targets_eval = [0] # center-out reaching target to evaluate @@ -48,15 +48,15 @@ pNames.append('backgroundrate'); pRanges.append([50,150]) pNames.append('explorMovsFactor'); pRanges.append([0.1,5]) pNames.append('cmdmaxrate'); pRanges.append([500,2000]) -pNames.append('PMdconnweight'); pRanges.append([0.5,4]) -pNames.append('PMdconnprob'); pRanges.append([1,8]) +pNames.append('PMdconnweight'); pRanges.append([0.5,4]) +pNames.append('PMdconnprob'); pRanges.append([1,8]) num_inputs = len(pNames) # evol specific params -mutation_rate = 0.4 # only for custom EC +mutation_rate = 0.4 # only for custom EC crossover_rate = 0.2 # % of children with crossover -ux_bias = round(0.1*num_inputs) +ux_bias = round(0.1*num_inputs) # Set bounds and allowed ranges for params @@ -65,11 +65,11 @@ def bound_params(candidate, args = []): for i,p in enumerate(candidate): cBound.append(max(min(p, max(pRanges[i])), min(pRanges[i]))) - # need to be integer + # need to be integer cBound[0] = round(max(min(candidate[0], max(pRanges[0])), min(pRanges[0]))/1000.0)*1000.0 # round to 1000.0 #cBound[1] = round(max(min(candidate[1], max(pRanges[1])), min(pRanges[1]))) #cBound[10] = round(max(min(candidate[10], max(pRanges[10])), min(pRanges[10]))) - + # fixed values from list #param14 = min(param14_range, key=lambda x:abs(x-c[13])) @@ -79,14 +79,14 @@ def bound_params(candidate, args = []): ############################################################################### ### Generate new set of random values for params -############################################################################### +############################################################################### def generate_rastrigin(random, args): size = args.get('num_inputs', 10) paramsRand = [] for iparam in range(len(pNames)): paramsRand.append(random.uniform(min(pRanges[iparam]),max(pRanges[iparam]))) - # need to be integer + # need to be integer paramsRand[0] = round(paramsRand[0]/1000.0)*1000.0 #paramsRand[1] = round(paramsRand[1]) #paramsRand[10] = round(paramsRand[10]) @@ -99,17 +99,17 @@ def generate_rastrigin(random, args): ############################################################################### ### Observer -############################################################################### +############################################################################### def my_observer(population, num_generations, num_evaluations, args): #ngen=num_generations best = max(population) - print('{0:6} -- {1} : {2}'.format(num_generations, - best.fitness, - str(best.candidate))) + print(('{0:6} -- {1} : {2}'.format(num_generations, + best.fitness, + str(best.candidate)))) ############################################################################### ### Custom mutator (nonuniform taking into account bounds) -############################################################################### +############################################################################### @mutator def nonuniform_bounds_mutation(random, candidate, args): """Return the mutants produced by nonuniform mutation on the candidates. @@ -117,11 +117,11 @@ def nonuniform_bounds_mutation(random, candidate, args): random -- the random number generator object candidate -- the candidate solution args -- a dictionary of keyword arguments - Required keyword arguments in args: - Optional keyword arguments in args: + Required keyword arguments in args: + Optional keyword arguments in args: - *mutation_strength* -- the strength of the mutation, where higher values correspond to greater variation (default 1) - + """ #bounder = args['_ec'].bounder #num_gens = args['_ec'].num_generations @@ -141,12 +141,12 @@ def nonuniform_bounds_mutation(random, candidate, args): ############################################################################### ### Parallel evaluation -############################################################################### +############################################################################### def parallel_evaluation_pbs(candidates, args): global ngen, targets_eval simdatadir = args.get('simdatadir') # load params ngen += 1 # increase number of generations - maxiter_wait=args.get('maxiter_wait',2000) # + maxiter_wait=args.get('maxiter_wait',2000) # default_error=args.get('default_error',0.3) #run pbs jobs @@ -155,8 +155,8 @@ def parallel_evaluation_pbs(candidates, args): for i, c in enumerate(candidates): outfilestem=simdatadir+"/gen_"+str(ngen)+"_cand_"+str(i) # set filename - for itarget in targets_eval: - with open('%s_params'% (outfilestem), 'w') as f: # save current candidate params to file + for itarget in targets_eval: + with open('%s_params'% (outfilestem), 'w') as f: # save current candidate params to file pickle.dump(c, f) command = 'mpirun -machinefile %s/nodes%d -np %d nrniv -python -mpi main.py outfilestem="%s" targetid=%d'%(simdatadir, i+1, numproc, outfilestem, itarget) # set command to run @@ -183,7 +183,7 @@ def parallel_evaluation_pbs(candidates, args): #SBATCH -e %s.err# Name of stderr output file #SBATCH --partition=compute # submit to the 'large' queue for jobs > 256 nodes #SBATCH -J %s # Job name -#SBATCH -t %s # Run time (hh:mm:ss) +#SBATCH -t %s # Run time (hh:mm:ss) #SBATCH --mail-user=%s #SBATCH --mail-type=%s #SBATCH -A %s # Allocation name to charge job against @@ -225,14 +225,14 @@ def parallel_evaluation_pbs(candidates, args): cd '/home/salvadord/m1ms/sim/' -""" % (job_name, job_name, walltime, email, mailType, project, nodes, coresPerNode, simdatadir, +""" % (job_name, job_name, walltime, email, mailType, project, nodes, coresPerNode, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir, simdatadir) - print job_string # print sbatch script + print(job_string) # print sbatch script batchfile = '%s/gen_%d.sbatch'%(simdatadir, ngen) with open(batchfile, 'w') as text_file: @@ -252,10 +252,10 @@ def parallel_evaluation_pbs(candidates, args): jobs_completed=0 while jobs_completed < total_jobs: #print outfilestem - print str(jobs_completed)+" / "+str(total_jobs)+" jobs completed" + print(str(jobs_completed)+" / "+str(total_jobs)+" jobs completed") unfinished = [[(i,j) for j,y in enumerate(x) if y is None] for i, x in enumerate(targetFitness)] unfinished = [item for sublist in unfinished for item in sublist] - print "unfinished:"+str(unfinished) + print("unfinished:"+str(unfinished)) for (icand,itarget) in unfinished: # load error from file try: @@ -264,29 +264,29 @@ def parallel_evaluation_pbs(candidates, args): errorDic=pickle.load(f) targetFitness[icand][itarget] = errorDic['errorFitness'] jobs_completed+=1 - print "icand:",icand," itarget:",itarget," error: "+str(errorDic['errorFitness']) + print("icand:",icand," itarget:",itarget," error: "+str(errorDic['errorFitness'])) except: pass #print "Waiting for job: "+str(i)+" ... iteration:"+str(num_iters[i]) num_iters+=1 - if num_iters>=maxiter_wait: #or (num_iters>maxiter_wait/2 and jobs_completed>(0.95*total_jobs)): - print "max iterations reached -- remaining jobs set to default error" + if num_iters>=maxiter_wait: #or (num_iters>maxiter_wait/2 and jobs_completed>(0.95*total_jobs)): + print("max iterations reached -- remaining jobs set to default error") for (icand,itarget) in unfinished: targetFitness[icand][itarget] = default_error jobs_completed+=1 sleep(2) # sleep 2 seconds before checking agains - print targetFitness + print(targetFitness) try: fitness = [mean(x) for x in targetFitness] - except: + except: fitness = [default_error for x in range(len(candidates))] - print 'fitness:',fitness + print('fitness:',fitness) return fitness ############################################################################### ### Multiprocessing Migration -############################################################################### +############################################################################### class MultiprocessingMigratorNoBlock(object): """Migrate among processes on the same machine. remove lock @@ -296,7 +296,7 @@ def __init__(self, max_migrants=1, migration_interval=10): self.migration_interval = migration_interval self.migrants = multiprocessing.Queue(self.max_migrants) self.__name__ = self.__class__.__name__ - + def __call__(self, random, population, args): # only migrate every migrationInterval generations if (args["_ec"].num_generations % self.migration_interval)==0: @@ -308,12 +308,12 @@ def __call__(self, random, population, args): if evaluate_migrant: fit = args["_ec"].evaluator([migrant.candidate], args) migrant.fitness = fit[0] - args["_ec"].num_evaluations += 1 - except Queue.Empty: + args["_ec"].num_evaluations += 1 + except queue.Empty: pass try: self.migrants.put(old_migrant, block=False) - except Queue.Full: + except queue.Full: pass return population @@ -345,7 +345,7 @@ def setInitial(simdatadir): # set global variable to track number of gens to initial_gen ngen = initial_gen - print initial_gen, initial_cs, initial_fit + print(initial_gen, initial_cs, initial_fit) return initial_gen, initial_cs, initial_fit @@ -353,14 +353,14 @@ def setInitial(simdatadir): ### Create islands ############################################################################### def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluations, max_generations, \ - num_inputs, mutation_rate, crossover_rate, ux_bias, pop_size, num_elites): + num_inputs, mutation_rate, crossover_rate, ux_bias, pop_size, num_elites): global num_islands - # create folder - if num_islands > 1: + # create folder + if num_islands > 1: simdatadir = simdatadir+'_island_'+str(i) mdir_str='mkdir %s' % (simdatadir) - os.system(mdir_str) + os.system(mdir_str) # if individuals.csv already exists, continue from last generation if os.path.isfile(simdatadir+'/individuals.csv'): # disabled by adding '!!' @@ -379,28 +379,28 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat seedfile.write('{0}'.format(my_seed)) seedfile.close() prng = Random() - prng.seed(my_seed) + prng.seed(my_seed) # Custom algorithm based on Krichmar's params if evolAlgorithm == 'customEvol': - # a real-valued optimization algo- rithm called Evolution Strategies (De Jong, 2002) + # a real-valued optimization algo- rithm called Evolution Strategies (De Jong, 2002) # was used with deterministic tournament selection, weak-elitism replacement, 40% Gaussian mutation and 50% - # crossover. Weak-elitism ensures the overall fitness monotonically increases each generation by replacing the - # worst fitness individual of the offspring population with the best fitness individual of the parent population. + # crossover. Weak-elitism ensures the overall fitness monotonically increases each generation by replacing the + # worst fitness individual of the offspring population with the best fitness individual of the parent population. ea = inspyred.ec.EvolutionaryComputation(prng) ea.selector = inspyred.ec.selectors.tournament_selection - ea.variator = [inspyred.ec.variators.uniform_crossover, nonuniform_bounds_mutation] + ea.variator = [inspyred.ec.variators.uniform_crossover, nonuniform_bounds_mutation] #inspyred.ec.variators.gaussian_mutation] ea.replacer = inspyred.ec.replacers.generational_replacement#inspyred.ec.replacers.plus_replacement #inspyred.ec.replacers.truncation_replacement (with num_selected=50) ea.terminator = inspyred.ec.terminators.generation_termination ea.observer = [inspyred.ec.observers.stats_observer, inspyred.ec.observers.file_observer] - final_pop = ea.evolve(generator=generate_rastrigin, + final_pop = ea.evolve(generator=generate_rastrigin, evaluator=parallel_evaluation_pbs, - pop_size=pop_size, + pop_size=pop_size, bounder=bound_params, maximize=False, max_evaluations=max_evaluations, @@ -418,7 +418,7 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat initial_gen=initial_gen, initial_cs=initial_cs, initial_fit=initial_fit) - + # Genetic elif evolAlgorithm == 'genetic': ea = inspyred.ec.GA(prng) @@ -444,7 +444,7 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat if num_islands > 1: ea.migrator = mp_migrator ea.terminator = inspyred.ec.terminators.generation_termination ea.observer = [inspyred.ec.observers.stats_observer, inspyred.ec.observers.file_observer] - final_pop = ea.evolve(generator=generate_rastrigin, + final_pop = ea.evolve(generator=generate_rastrigin, evaluator=parallel_evaluation_pbs, pop_size=pop_size, bounder=bound_params, @@ -468,7 +468,7 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat ea.variator = [inspyred.ec.variators.uniform_crossover, ea._internal_variation] ea.terminator = inspyred.ec.terminators.generation_termination ea.observer = [inspyred.ec.observers.stats_observer, inspyred.ec.observers.file_observer] - final_pop = ea.evolve(generator=generate_rastrigin, + final_pop = ea.evolve(generator=generate_rastrigin, evaluator=parallel_evaluation_pbs, pop_size=pop_size, bounder=bound_params, @@ -485,7 +485,7 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat initial_cs=initial_cs, initial_fit=initial_fit) - + # Simulated Annealing elif evolAlgorithm == 'simulatedAnnealing': @@ -590,8 +590,8 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat - best = max(final_pop) - print('Best Solution: \n{0}'.format(str(best))) + best = max(final_pop) + print(('Best Solution: \n{0}'.format(str(best)))) return ea @@ -600,10 +600,10 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat ### Main - logging, island model params, launch multiprocessing ############################################################################### if __name__ == '__main__': - # create folder + # create folder mdir_str='mkdir -p %s' % (simdatadir) - os.system(mdir_str) - + os.system(mdir_str) + # debug info logger = logging.getLogger('inspyred.ec') logger.setLevel(logging.DEBUG) @@ -611,7 +611,7 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat file_handler.setLevel(logging.DEBUG) formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') file_handler.setFormatter(formatter) - logger.addHandler(file_handler) + logger.addHandler(file_handler) # run single population or multiple islands rand_seed = int(time()) @@ -626,4 +626,4 @@ def create_island(rand_seed, island_number, mp_migrator, simdatadir, max_evaluat p.start() jobs.append(p) for j in jobs: - j.join() \ No newline at end of file + j.join() diff --git a/izhi.py b/izhi.py index 2a334ca..7d99702 100644 --- a/izhi.py +++ b/izhi.py @@ -1,7 +1,7 @@ """ IZHI -Python wrappers for the different celltypes of Izhikevich neuron. +Python wrappers for the different celltypes of Izhikevich neuron. Equations and parameter values taken from Izhikevich EM (2007). @@ -10,7 +10,7 @@ Equation for synaptic inputs taken from Izhikevich EM, Edelman GM (2008). - "Large-scale model of mammalian thalamocortical systems." + "Large-scale model of mammalian thalamocortical systems." PNAS 105(9) 3593-3598. Cell types available are based on Izhikevich, 2007 book: @@ -48,7 +48,7 @@ def createcell(section, C, k, vr, vt, vpeak, a, b, c, d, celltype, cellid): cell.d = d cell.celltype = celltype # Set cell celltype (used for setting celltype-specific dynamics) cell.cellid = cellid # Cell ID for keeping track which cell this is - cell.t0 = .0 + #cell.t0 = .0 return cell ## Cell types based on Izhikevich, 2007 book @@ -113,4 +113,4 @@ def thalamocortical(section, C=200, k=0.5, vr=-60, vt=-50, vpeak=40, a=0.1, b=15 ## Thalamic reticular nucleus cell def reticular(section, C=40, k=0.25, vr=-65, vt=-45, vpeak=0, a=0.015, b=10, c=-55, d=50, celltype=5, cellid=-1): cell = createcell(section, C, k, vr, vt, vpeak, a, b, c, d, celltype, cellid) - return cell \ No newline at end of file + return cell diff --git a/main.py b/main.py index 75a8497..c4d6d5f 100644 --- a/main.py +++ b/main.py @@ -18,18 +18,18 @@ if len(arg)==2: if hasattr(s,arg[0]) or hasattr(s,harg[1]): # Check that variable exists if arg[0] == 'outfilestem': - exec('s.'+arg[0]+'="'+arg[1]+'"') # Actually set variable - if s.rank==0: # messages only come from Master - print(' Setting %s=%s' %(arg[0],arg[1])) + exec('s.'+arg[0]+'="'+arg[1]+'"') # Actually set variable + if s.rank==0: # messages only come from Master + print((' Setting %s=%s' %(arg[0],arg[1]))) else: - exec('s.'+argv) # Actually set variable - if s.rank==0: # messages only come from Master - print(' Setting %s=%r' %(arg[0],eval(arg[1]))) - else: + exec('s.'+argv) # Actually set variable + if s.rank==0: # messages only come from Master + print((' Setting %s=%r' %(arg[0],eval(arg[1])))) + else: sys.tracebacklimit=0 raise Exception('Reading args from commandline: Variable "%s" not found' % arg[0]) elif argv=='-mpi': ismpi = True - else: pass # ignore -python + else: pass # ignore -python ############################################################################### diff --git a/mosinit.py b/mosinit.py index 79f0958..4b1f975 100644 --- a/mosinit.py +++ b/mosinit.py @@ -1 +1 @@ -execfile("main.py") +exec(compile(open("main.py", "rb").read(), "main.py", 'exec')) diff --git a/network.py b/network.py index e8a1ac6..9eaa711 100644 --- a/network.py +++ b/network.py @@ -28,7 +28,7 @@ from pylab import seed, rand, sqrt, exp, transpose, ceil, concatenate, array, zeros, ones, vstack, show, disp, mean, inf, concatenate, unique, delete from time import time, sleep from datetime import datetime -from scipy.io import savemat, loadmat +from scipy.io import savemat, loadmat import pickle from neuron import h, init, run # Import NEURON @@ -36,6 +36,7 @@ import analysis from arm import Arm # Class with arm methods and variables +import numpy as np ############################################################################### ### Sequences of commands to run full model @@ -44,16 +45,16 @@ # training and testing to 2 targets manually def runTrainTest2targets(): - # optimized values for musculoskeletal arm (here using dummy arm for demo purposes) + # optimized values for musculoskeletal arm (here using dummy arm for demo purposes) s.targetid = 0 s.trainTime = 2000 # 85000 # using 2 sec for demo purposes - s.stdpwin = 48.5 - s.eligwin = 117.8 + s.stdpwin = 48.5 + s.eligwin = 117.8 s.RLfactor = 0.01 #6 s.RLinterval = 76.8 - s.backgroundrate = s.backgroundrateTest = 134.5 + s.backgroundrate = s.backgroundrateTest = 134.5 s.backgroundrateExplor = 5 - s.cmdmaxrate = 528.8 + s.cmdmaxrate = 528.8 s.PMdconnweight = 1.0 s.PMdconnprob = 2.4 @@ -62,7 +63,7 @@ def runTrainTest2targets(): s.numTrials = ceil(s.trainTime/1000) s.trialTargets = [i%2 for i in range(int(s.numTrials+1))] # set target for each trial s.targetid=s.trialTargets[0] - + verystart=time() # store initial time s.plotraster = 1 # set plotting params @@ -75,7 +76,7 @@ def runTrainTest2targets(): s.armMinimalSave = 1 # save only arm related data # initialize network - createNetwork() + createNetwork() addStimulation() addBackground() @@ -101,7 +102,7 @@ def runTrainTest2targets(): s.explorMovs = 0 # disable exploratory movements s.duration = s.testTime # testing time s.armMinimalSave = 0 # save only arm related data - + s.targetid = 0 setupSim() runSim() @@ -111,8 +112,8 @@ def runTrainTest2targets(): if s.rank == 0: # save error to file error0 = mean(s.arm.errorAll) - print 'Target error for target ',s.targetid,' is:', error0 - s.arm.plotTraj(s.outfilestem+'_t0.png') + print('Target error for target ',s.targetid,' is:', error0) + s.arm.plotTraj(s.outfilestem+'_t0.png') # test target 1 s.targetid = 1 @@ -124,8 +125,8 @@ def runTrainTest2targets(): if s.rank == 0: # save error to file error1 = mean(s.arm.errorAll) - print 'Target error for target 0=', error0, '; target 1=', error1 - s.arm.plotTraj(s.outfilestem+'_t1.png') + print('Target error for target 0=', error0, '; target 1=', error1) + s.arm.plotTraj(s.outfilestem+'_t1.png') errorMean = (error0+error1)/2 errorFitness = errorMean + abs(error0-error1) # fitness penalizes difference between target errors @@ -134,18 +135,18 @@ def runTrainTest2targets(): errorDic['error1'] = error1 errorDic['meanError'] = errorMean errorDic['errorFitness'] = errorFitness - - print 'Mean error = %.4f ; Mean error + difference (fitness) = %.4f'%(errorMean, errorFitness) + + print('Mean error = %.4f ; Mean error + difference (fitness) = %.4f'%(errorMean, errorFitness)) s.targetid = 0 # so saves to correct file name (error of both targets saved to single file ending in target_0_error) - with open('%s_target_%d_error'% (s.outfilestem,s.targetid), 'w') as f: # save avg error over targets to outfilestem + with open('%s_target_%d_error'% (s.outfilestem,s.targetid), 'wb') as f: # save avg error over targets to outfilestem pickle.dump(errorDic, f) ## Wrapping up s.pc.runworker() # MPI: Start simulations running on each host s.pc.done() # MPI: Close MPI totaltime = time()-verystart # See how long it took in total - print('\nDone; total time = %0.1f s.' % totaltime) + print(('\nDone; total time = %0.1f s.' % totaltime)) if (s.plotraster==False and s.plotconn==False and s.plotweightchanges==False): h.quit() # Quit extra processes, or everything if plotting wasn't requested (since assume non-interactive) @@ -162,7 +163,7 @@ def runTrainTest2targetsOptim(): s.numTrials = ceil(s.trainTime/1000) s.trialTargets = [i%2 for i in range(int(s.numTrials+1))] # set target for each trial s.targetid=s.trialTargets[0] - + verystart=time() # store initial time s.plotraster = 0 # set plotting params @@ -183,7 +184,7 @@ def runTrainTest2targetsOptim(): s.timebetweensaves = s.trainTime - 1000 # initialize network - createNetwork() + createNetwork() addStimulation() addBackground() @@ -194,7 +195,7 @@ def runTrainTest2targetsOptim(): #saveData() #plotData() if s.rank == 0: # save png of traj - s.arm.plotTraj(s.outfilestem+'_train.png') # save traj fig to file + s.arm.plotTraj(s.outfilestem+'_train.png') # save traj fig to file analysis.plotweightchanges(s.outfilestem+'_train_weights.png') test = 1 @@ -209,7 +210,7 @@ def runTrainTest2targetsOptim(): s.explorMovs = 0 # disable exploratory movements s.duration = s.testTime # testing time s.armMinimalSave = 0 # save only arm related data - + s.targetid = 0 setupSim() runSim() @@ -219,8 +220,8 @@ def runTrainTest2targetsOptim(): if s.rank == 0: # save error to file error0 = mean(s.arm.errorAll) - print 'Target error for target ',s.targetid,' is:', error0 - s.arm.plotTraj(s.outfilestem+'_t0.png') + print('Target error for target ',s.targetid,' is:', error0) + s.arm.plotTraj(s.outfilestem+'_t0.png') analysis.plotraster(s.outfilestem+'_t0_raster.png') # test target 1 @@ -233,8 +234,8 @@ def runTrainTest2targetsOptim(): if s.rank == 0: # save error to file error1 = mean(s.arm.errorAll) - print 'Target error for target 0=', error0, '; target 1=', error1 - s.arm.plotTraj(s.outfilestem+'_t1.png') + print('Target error for target 0=', error0, '; target 1=', error1) + s.arm.plotTraj(s.outfilestem+'_t1.png') analysis.plotraster(s.outfilestem+'_t1_raster.png') errorMean = (error0+error1)/2 @@ -244,18 +245,18 @@ def runTrainTest2targetsOptim(): errorDic['error1'] = error1 errorDic['meanError'] = errorMean errorDic['errorFitness'] = errorFitness - - print 'Mean error = %.4f ; Mean error + difference (fitness) = %.4f'%(errorMean, errorFitness) + + print('Mean error = %.4f ; Mean error + difference (fitness) = %.4f'%(errorMean, errorFitness)) s.targetid = 0 # so saves to correct file name (error of both targets saved to single file ending in target_0_error) - with open('%s_target_%d_error'% (s.outfilestem,s.targetid), 'w') as f: # save avg error over targets to outfilestem + with open('%s_target_%d_error'% (s.outfilestem,s.targetid), 'wb') as f: # save avg error over targets to outfilestem pickle.dump(errorDic, f) ## Wrapping up s.pc.runworker() # MPI: Start simulations running on each host s.pc.done() # MPI: Close MPI totaltime = time()-verystart # See how long it took in total - print('\nDone; total time = %0.1f s.' % totaltime) + print(('\nDone; total time = %0.1f s.' % totaltime)) if (s.plotraster==False and s.plotconn==False and s.plotweightchanges==False): h.quit() # Quit extra processes, or everything if plotting wasn't requested (since assume non-interactive) @@ -265,15 +266,15 @@ def runTrainTest2targetsOptim(): def createNetwork(): ## Print diagnostic information - if s.rank==0: print("\nCreating simulation of %i cells for %0.1f s on %i hosts..." % (sum(s.popnumbers),s.duration/1000.,s.nhosts)) + if s.rank==0: print(("\nCreating simulation of %i cells for %0.1f s on %i hosts..." % (sum(s.popnumbers),s.duration/1000.,s.nhosts))) s.pc.barrier() - + ## Create empty data structures s.cells=[] # Create empty list for storing cells s.dummies=[] # Create empty list for storing fake sections s.gidVec=[] # Empty list for storing GIDs (index = local id; value = gid) s.gidDic = {} # Empyt dict for storing GIDs (key = gid; value = local id) -- ~x6 faster than gidVec.index() - + ## Set cell types celltypes=[] @@ -294,8 +295,8 @@ def createNetwork(): s.xlocs = s.modelsize*rand(s.ncells) # Create random x locations s.ylocs = s.modelsize*rand(s.ncells) # Create random y locations s.zlocs = rand(s.ncells) # Create random z locations - for c in range(s.ncells): - s.zlocs[c] = s.corticalthick * (s.zlocs[c]*(s.popyfrac[s.cellpops[c]][1]-s.popyfrac[s.cellpops[c]][0]) + s.popyfrac[s.cellpops[c]][0]) # calculate based on yfrac for population and corticalthick + for c in range(s.ncells): + s.zlocs[c] = s.corticalthick * (s.zlocs[c]*(s.popyfrac[s.cellpops[c]][1]-s.popyfrac[s.cellpops[c]][0]) + s.popyfrac[s.cellpops[c]][0]) # calculate based on yfrac for population and corticalthick ## Actually create the cells @@ -303,7 +304,7 @@ def createNetwork(): s.hostspikevecs = [] # Empty list for storing host-specific spike vectors s.cellsperhost = 0 if s.PMdinput == 'Plexon': ninnclDic = len(s.innclDic) # number of PMd created in this worker - for c in xrange(int(s.rank), s.ncells, s.nhosts): + for c in range(int(s.rank), s.ncells, s.nhosts): s.dummies.append(h.Section()) # Create fake sections gid = c if s.cellnames[gid] == 'PMd': @@ -323,16 +324,16 @@ def createNetwork(): cell = celltypes[gid](cellid = gid) # create an NSLOC cell.number = s.backgroundnumber cell.interval = s.backgroundrateMin**-1*1e3 - + elif s.cellnames[gid] == 'ASC': - cell = celltypes[gid](cellid = gid) #create an NSLOC - else: - if s.cellclasses[gid]==3: + cell = celltypes[gid](cellid = gid) #create an NSLOC + else: + if s.cellclasses[gid]==3: cell = s.fastspiking(s.dummies[s.cellsperhost], vt=-47, cellid=gid) # Don't use LTS cell, but instead a FS cell with a low threshold - else: + else: cell = celltypes[gid](s.dummies[s.cellsperhost], cellid=gid) # Create a new cell of the appropriate type (celltypes[gid]) and store it #if s.verbose>0: s.cells[-1].useverbose(s.verbose, s.filename+'los.txt') # Turn on diagnostic to file - s.cells.append(cell) + s.cells.append(cell) s.gidVec.append(gid) # index = local id; value = global id s.gidDic[gid] = s.cellsperhost # key = global id; value = local id -- used to get local id because gid.index() too slow! s.pc.set_gid2node(gid, s.rank) @@ -344,31 +345,31 @@ def createNetwork(): s.spikerecorders.append(spikerecorder) s.pc.cell(gid, s.spikerecorders[s.cellsperhost]) s.cellsperhost += 1 # contain cell numbers per host including PMd and P - print(' Number of cells on node %i: %i ' % (s.rank,len(s.cells))) + print((' Number of cells on node %i: %i ' % (s.rank,len(s.cells)))) s.pc.barrier() ## Calculate motor command cell ranges so can be used for EDSC and IDSC connectivity - nCells = s.motorCmdEndCell - s.motorCmdStartCell + nCells = s.motorCmdEndCell - s.motorCmdStartCell s.motorCmdCellRange = [] for i in range(s.nMuscles): - s.motorCmdCellRange.append(range(s.motorCmdStartCell + (nCells/s.nMuscles)*i, s.motorCmdStartCell + (nCells/s.nMuscles)*i + (nCells/s.nMuscles) )) # cells used to for shoulder motor command - + s.motorCmdCellRange.append(list(np.arange(s.motorCmdStartCell + (nCells/s.nMuscles)*i, s.motorCmdStartCell + (nCells/s.nMuscles)*i + (nCells/s.nMuscles),dtype=int))) # cells used to for shoulder motor command + ## Calculate distances and probabilities - if s.rank==0: print('Calculating connection probabilities (est. time: %i s)...' % (s.performance*s.cellsperhost**2/3e4)) + if s.rank==0: print(('Calculating connection probabilities (est. time: %i s)...' % (s.performance*s.cellsperhost**2/3e4))) conncalcstart = s.time() # See how long connecting the cells takes s.nconnpars = 5 # Connection parameters: pre- and post- cell ID, weight, distances, delays s.conndata = [[] for i in range(s.nconnpars)] # List for storing connections nPostCells = 0 EDSCpre = [] # to keep track of EB5->EDSC connection and replicate in EB5->IDSC for c in range(s.cellsperhost): # Loop over all postsynaptic cells on this host (has to be postsynaptic because of gid_connect) - gid = s.gidVec[c] # Increment global identifier + gid = s.gidVec[c] # Increment global identifier if s.cellnames[gid] == 'PMd' or s.cellnames[gid] == 'ASC': # There are no presynaptic connections for PMd or ASC. continue nPostCells += 1 - if s.toroidal: + if s.toroidal: xpath=(abs(s.xlocs-s.xlocs[gid]))**2 xpath2=(s.modelsize-abs(s.xlocs-s.xlocs[gid]))**2 xpath[xpath2ER5 conn (full conn) PMdId = (gid % s.server.numPMd) + s.ncells - s.server.numPMd #CHECK THIS! allconnprobs[PMdId] = s.connprobs[s.PMd,s.ER5] # to make this connected to ER5 allrands[PMdId] = 0 # to make this connect to ER5 - distances[PMdId] = 300 # to make delay 5 in conndata[3] + distances[PMdId] = 300 # to make delay 5 in conndata[3] makethisconnection = allconnprobs>allrands # Perform test to see whether or not this connection should be made preids = array(makethisconnection.nonzero()[0],dtype='int') # Return True elements of that array for presynaptic cell IDs - if s.PMdinput == 'targetSplit' and s.cellnames[gid] == 'ER5': # PMds 0-47 -> ER5 0-47 ; PMds 48-95 -> ER5 48-95 + if s.PMdinput == 'targetSplit' and s.cellnames[gid] == 'ER5': # PMds 0-47 -> ER5 0-47 ; PMds 48-95 -> ER5 48-95 if gid < s.popGidStart[s.ER5] + s.popnumbers[s.ER5]/2: - prePMd = [(x - s.popGidStart[s.ER5])%(s.popnumbers[s.PMd]/2) + s.popGidStart[s.PMd] for x in range(gid, gid+1)] # input from 2 PMds + prePMd = [(x - s.popGidStart[s.ER5])%(s.popnumbers[s.PMd]/2) + s.popGidStart[s.PMd] for x in range(gid, gid+1)] # input from 2 PMds else: - prePMd = [(x - s.popGidStart[s.ER5])%(s.popnumbers[s.PMd]/2) + s.popGidStart[s.PMd] + s.popnumbers[s.PMd]/2 for x in range(gid, gid+1)] # input from 2 PMds - if array(prePMd).all() < s.popGidEnd[s.PMd]: + prePMd = [(x - s.popGidStart[s.ER5])%(s.popnumbers[s.PMd]/2) + s.popGidStart[s.PMd] + s.popnumbers[s.PMd]/2 for x in range(gid, gid+1)] # input from 2 PMds + if array(prePMd).all() < s.popGidEnd[s.PMd]: #print 'prePMd=%d to ER5=%d:'%(prePMd[0],gid) preids = concatenate([preids, prePMd]) if s.cellnames[gid] == 'EDSC': # save EDSC presyn cells to replicate in IDSC, and add inputs from IDSC EDSCpre.append(array(preids)) # save EDSC presyn cells before adding IDSC input invPops = [1, 0, 3, 2] # each postsyn ESDC cell will receive input from all the antagonistic muscle IDSCs IDSCpre = [s.motorCmdCellRange[invPops[i]] - s.popGidStart[s.EDSC] + s.popGidStart[s.IDSC] for i in range(s.nMuscles) if gid in s.motorCmdCellRange[i]][0] - preids = concatenate([preids, IDSCpre]) # add IDSC presynaptic input to EDSC + preids = concatenate([preids, IDSCpre]) # add IDSC presynaptic input to EDSC elif s.cellnames[gid] == 'IDSC': # use same presyn cells as for EDSC (antagonistic inhibition) preids = array(EDSCpre.pop(0)) postids = array(gid+zeros(len(preids)),dtype='int') # Post-synaptic cell IDs @@ -423,12 +424,12 @@ def createNetwork(): for pp in range(s.nconnpars): s.conndata[pp] = array(concatenate([s.conndata[pp][c] for c in range(nPostCells)])) # Turn pre- and post- cell IDs lists into vectors s.nconnections = len(s.conndata[0]) # Find out how many connections we're going to make conncalctime = time()-conncalcstart # See how long it took - if s.rank==0: print(' Done; time = %0.1f s' % conncalctime) + if s.rank==0: print((' Done; time = %0.1f s' % conncalctime)) # set plastic connections based on plasConnsType (from evol alg) if s.plastConnsType == 0: - s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC]] # only spinal cord + s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC]] # only spinal cord elif s.plastConnsType == 1: s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC], [s.ER2,s.ER5], [s.ER5,s.EB5], [s.ER2,s.EB5], [s.ER5,s.ER2]] # + L2-L5 elif s.plastConnsType == 2: @@ -439,7 +440,7 @@ def createNetwork(): [s.ER5,s.ER6], [s.ER6,s.ER5], [s.ER6,s.EB5], \ [s.ER2,s.IL2], [s.ER2,s.IF2], [s.ER5,s.IL5], [s.ER5,s.IF5], [s.EB5,s.IL5], [s.EB5,s.IF5]] # + Inh # same with additional plasticity between PMd->L5A - elif s.plastConnsType == 4: + elif s.plastConnsType == 4: s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC], [s.PMd,s.ER5]] # only spinal cord + pmd elif s.plastConnsType == 5: s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC], [s.PMd,s.ER5], # spinal cord + pmd @@ -449,7 +450,7 @@ def createNetwork(): [s.ER2,s.ER5], [s.ER5,s.EB5], [s.ER2,s.EB5], [s.ER5,s.ER2], # + L2-L5 [s.ER5,s.ER6], [s.ER6,s.ER5], [s.ER6,s.EB5]] # + L6 elif s.plastConnsType == 7: - s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC], [s.PMd,s.ER5], # spinal cord + pmd + s.plastConns = [[s.ASC,s.ER2], [s.EB5,s.EDSC], [s.EB5,s.IDSC], [s.PMd,s.ER5], # spinal cord + pmd [s.ER2,s.ER5], [s.ER5,s.EB5], [s.ER2,s.EB5], [s.ER5,s.ER2], # + L2-L5 [s.ER5,s.ER6], [s.ER6,s.ER5], [s.ER6,s.EB5], # + L6 [s.ER2,s.IL2], [s.ER2,s.IF2], [s.ER5,s.IL5], [s.ER5,s.IF5], [s.EB5,s.IL5], [s.EB5,s.IF5]] # + Inh @@ -457,8 +458,8 @@ def createNetwork(): ## Actually make connections - if s.rank==0: print('Making connections (est. time: %i s)...' % (s.performance*s.nconnections/9e2)) - print(' Number of connections on host %i: %i' % (s.rank, s.nconnections)) + if s.rank==0: print(('Making connections (est. time: %i s)...' % (s.performance*s.nconnections/9e2))) + print((' Number of connections on host %i: %i' % (s.rank, s.nconnections))) connstart = time() # See how long connecting the cells takes s.connlist = [] # Create array for storing each of the connections s.stdpconndata = [] # Store data on STDP connections @@ -467,7 +468,7 @@ def createNetwork(): s.precons = [] # Initialize array for presynaptic spike counters s.pstcons = [] # Initialize array for postsynaptic spike counters for con in range(s.nconnections): # Loop over each connection - pregid = s.conndata[0][con] # GID of presynaptic cell + pregid = s.conndata[0][con] # GID of presynaptic cell pstgid = s.conndata[1][con] # Index of postsynaptic cell pstid = s.gidDic[pstgid]# Index of postynaptic cell -- convert from GID to local newcon = s.pc.gid_connect(pregid, s.cells[pstid]) # Create a connection @@ -496,15 +497,15 @@ def createNetwork(): stdpmech.RLantiwt = s.RLrates[s.EorI[pregid],1] # Depression rate stdpmech.tauhebb = stdpmech.tauanti = s.stdpwin # stdp time constant(ms) stdpmech.RLwindhebb = stdpmech.RLwindhebb = s.eligwin # RL eligibility trace window length (ms) - stdpmech.useRLexp = s.useRLexp # RL + stdpmech.useRLexp = s.useRLexp # RL stdpmech.softthresh = s.useRLsoft # RL soft-thresholding else: stdpmech.RLon = 0 # make sure RL is off - + s.nstdpconns = len(s.stdpconndata) # Get number of STDP connections conntime = time()-connstart # See how long it took - if s.usestdp: print(' Number of STDP connections on host %i: %i' % (s.rank, s.nstdpconns)) - if s.rank==0: print(' Done; time = %0.1f s' % conntime) + if s.usestdp: print((' Number of STDP connections on host %i: %i' % (s.rank, s.nstdpconns))) + if s.rank==0: print((' Done; time = %0.1f s' % conntime)) ############################################################################### @@ -519,7 +520,7 @@ def addStimulation(): s.stimconns=[] # Create input connections s.stimtimevecs = [] # Create array for storing time vectors s.stimweightvecs = [] # Create array for holding weight vectors - if s.saveraw: + if s.saveraw: s.stimspikevecs=[] # A list for storing actual cell voltages (WARNING, slow!) s.stimrecorders=[] # And for recording spikes for stim in range(len(s.stimpars)): # Loop over each stimulus type @@ -528,44 +529,44 @@ def addStimulation(): stimvecs = s.makestim(ts.isi, ts.var, ts.width, ts.weight, ts.sta, ts.fin, ts.shape) # Time-probability vectors s.stimstruct.append([ts.name, stimvecs]) # Store for saving later s.stimtimevecs.append(h.Vector().from_python(stimvecs[0])) - + for c in range(s.cellsperhost): - gid = s.cellsperhost*int(s.rank)+c # For deciding E or I + gid = s.cellsperhost*int(s.rank)+c # For deciding E or I seed(s.id32('%d'%(s.randseed+gid))) # Reset random number generator for this cell if ts.fraction>rand(): # Don't do it for every cell necessarily if any(s.cellpops[gid]==ts.pops) and s.xlocs[gid]>=ts.loc[0,0] and s.xlocs[gid]<=ts.loc[0,1] and s.ylocs[gid]>=ts.loc[1,0] and s.ylocs[gid]<=ts.loc[1,1]: - + maxweightincrease = 20 # Otherwise could get infinitely high, infinitely close to the stimulus distancefromstimulus = sqrt(sum((array([s.xlocs[gid], s.ylocs[gid]])-s.modelsize*ts.falloff[0])**2)) fallofffactor = min(maxweightincrease,(ts.falloff[1]/distancefromstimulus)**2) s.stimweightvecs.append(h.Vector().from_python(stimvecs[1]*fallofffactor)) # Scale by the fall-off factor - + stimrand = h.Random() stimrand.MCellRan4() # If everything has the same seed, should happen at the same time stimrand.negexp(1) stimrand.seq(s.id32('%d'%(s.randseed+gid))*1e3) # Set the sequence i.e. seed s.stimrands.append(stimrand) - + stimsource = h.NetStim() # Create a NetStim stimsource.interval = ts.rate**-1*1e3 # Interval between spikes stimsource.number = 1e9 # Number of spikes stimsource.noise = ts.noise # Fractional noise in timing stimsource.noiseFromRandom(stimrand) # Set it to use this random number generator s.stimsources.append(stimsource) # Save this NetStim - + stimconn = h.NetCon(stimsource, s.cells[c]) # Connect this noisy input to a cell for r in range(s.nreceptors): stimconn.weight[r]=0 # Initialize weights to 0, otherwise get memory leaks s.stimweightvecs[-1].play(stimconn._ref_weight[0], s.stimtimevecs[-1]) # Play most-recently-added vectors into weight stimconn.delay=s.mindelay # Specify the delay in ms -- shouldn't make a spot of difference s.stimconns.append(stimconn) # Save this connnection - + if s.saveraw: # and c <=100: stimspikevec = h.Vector() # Initialize vector s.stimspikevecs.append(stimspikevec) # Keep all those vectors stimrecorder = h.NetCon(stimsource, None) stimrecorder.record(stimspikevec) # Record simulation time s.stimrecorders.append(stimrecorder) - print(' Number of stimuli created on host %i: %i' % (s.rank, len(s.stimsources))) + print((' Number of stimuli created on host %i: %i' % (s.rank, len(s.stimsources)))) ############################################################################### @@ -580,7 +581,7 @@ def addBackground(): if s.savebackground: s.backgroundspikevecs=[] # A list for storing actual cell voltages (WARNING, slow!) s.backgroundrecorders=[] # And for recording spikes - for c in range(s.cellsperhost): + for c in range(s.cellsperhost): gid = s.gidVec[c] if s.cellnames[gid] == 'ASC' or s.cellnames[gid] == 'PMd' : # These pops won't receive background stimulations. pass @@ -590,11 +591,11 @@ def addBackground(): backgroundrand.negexp(1) s.backgroundrands.append(backgroundrand) if s.cellnames[gid] == 'EDSC' or s.cellnames[gid] == 'IDSC': - backgroundsource = h.NSLOC() # Create a NSLOC + backgroundsource = h.NSLOC() # Create a NSLOC backgroundsource.interval = s.backgroundrateMin**-1*1e3 # Take inverse of the frequency and then convert from Hz^-1 to ms backgroundsource.noise = 0.3 # Fractional noise in timing elif s.cellnames[gid] == 'EB5': - backgroundsource = h.NSLOC() # Create a NSLOC + backgroundsource = h.NSLOC() # Create a NSLOC backgroundsource.interval = s.backgroundrate**-1*1e3 # Take inverse of the frequency and then convert from Hz^-1 to ms backgroundsource.noise = s.backgroundnoise # Fractional noise in timing else: @@ -606,25 +607,25 @@ def addBackground(): backgroundsource.number = s.backgroundnumber # Number of spikes s.backgroundsources.append(backgroundsource) # Save this NetStim s.backgroundgid.append(gid) # append cell gid associated to this netstim - + backgroundconn = h.NetCon(backgroundsource, s.cells[c]) # Connect this noisy input to a cell for r in range(s.nreceptors): backgroundconn.weight[r]=0 # Initialize weights to 0, otherwise get memory leaks if s.cellnames[gid] == 'EDSC' or s.cellnames[gid] == 'IDSC': backgroundconn.weight[s.backgroundreceptor] = s.backgroundweightExplor # Specify the weight for the EDSC, IDSC and PMd background input - elif s.cellnames[gid] == 'EB5' and s.explorMovs == 2: - backgroundconn.weight[s.backgroundreceptor] = s.backgroundweightExplor # Weight for EB5 input if explor movs via EB5 + elif s.cellnames[gid] == 'EB5' and s.explorMovs == 2: + backgroundconn.weight[s.backgroundreceptor] = s.backgroundweightExplor # Weight for EB5 input if explor movs via EB5 else: backgroundconn.weight[s.backgroundreceptor] = s.backgroundweight[s.EorI[gid]] # Specify the weight -- 1 is NMDA receptor for smoother, more summative activation backgroundconn.delay=2 # Specify the delay in ms -- shouldn't make a spot of difference s.backgroundconns.append(backgroundconn) # Save this connnection - + if s.savebackground: backgroundspikevec = h.Vector() # Initialize vector s.backgroundspikevecs.append(backgroundspikevec) # Keep all those vectors backgroundrecorder = h.NetCon(backgroundsource, None) backgroundrecorder.record(backgroundspikevec) # Record simulation time s.backgroundrecorders.append(backgroundrecorder) - print(' Number created on host %i: %i' % (s.rank, len(s.backgroundsources))) + print((' Number created on host %i: %i' % (s.rank, len(s.backgroundsources)))) s.pc.barrier() @@ -654,7 +655,7 @@ def setupSim(): for c in range(s.cellsperhost): # Loop over each cell and decide which LFP population, if any, it belongs to gid = s.gidVec[c] # Get this cell's GID if s.cellnames[gid] == 'ASC' or s.cellnames[gid] == 'PMd': # 'ER2' won't be fired by background stimulations. - continue + continue for pop in range(s.nlfps): # Loop over each LFP population thispop = s.cellpops[gid] # Population of this cell if sum(s.lfppops[pop]==thispop)>0: # There's a match @@ -663,14 +664,14 @@ def setupSim(): ## Set up raw recording s.rawrecordings = [] # A list for storing actual cell voltages (WARNING, slow!) - if s.saveraw: + if s.saveraw: if s.rank==0: print('\nSetting up raw recording...') s.nquantities = 5 # Number of variables from each cell to record from # Later this part should be modified because NSLOC doesn't have V, u and I. for c in range(s.cellsperhost): gid = s.gidVec[c] # Get this cell's GID if s.cellnames[gid] == 'ASC' or s.cellnames[gid] == 'PMd': # NSLOC doesn't have V, u and I - continue + continue recvecs = [h.Vector() for q in range(s.nquantities)] # Initialize vectors recvecs[0].record(h._ref_t) # Record simulation time recvecs[1].record(s.cells[c]._ref_V) # Record cell voltage @@ -687,7 +688,7 @@ def setupSim(): ## Set up virtual arm if s.useArm != 'None': if s.rank==0: print('\nSetting up virtual arm...') - s.arm = Arm(s.useArm, s.animArm, s.graphsArm) + s.arm = Arm(s.useArm, s.animArm, s.graphsArm) s.arm.targetid = s.targetid s.arm.setup(s)#duration, loopstep, RLinterval, pc, scale, popnumbers, p) @@ -701,18 +702,18 @@ def setupSim(): ''') if s.isOriginal == 0: # With communication program - if s.rank == 0: + if s.rank == 0: #serverManager = s.server.Manager() # isDp in confis.py = 0 s.server.Manager.start() # launch sever process - print "Server process completed and callback function initalized" + print("Server process completed and callback function initalized") e = s.server.Manager.Event() # Queue callback function in the NEURON queue - - # Wait for external spikes for PMd from Plexon + + # Wait for external spikes for PMd from Plexon if s.rank == 0: if s.server.isCommunication == 1: s.server.getServerInfo() # show parameters of the server process - print "[Waiting for spikes; run the client on Windows machine...]" + print("[Waiting for spikes; run the client on Windows machine...]") while s.server.queue.empty(): # only Rank 0 is waiting for spikes in the queue. pass s.pc.barrier() # other workers are waiting here. @@ -726,7 +727,7 @@ def setupSim(): # implement lesion if s.PMdlesion > 0: numLesionedCells = int(round(s.PMdlesion*numrawcells)) - if s.rank==0: print "Lesioning PMd input... removed spike times of %d (%d%%) cells"%(numLesionedCells, s.PMdlesion*100) + if s.rank==0: print("Lesioning PMd input... removed spike times of %d (%d%%) cells"%(numLesionedCells, s.PMdlesion*100)) for itarget in range(len(rawSpikesPMd)): for itrial in range(len(rawSpikesPMd[0])): for icell in range(numLesionedCells): @@ -734,14 +735,14 @@ def setupSim(): # generate spike vectors based on training time, trial duration, and target presentation (eg. alternating trials) if s.duration == 1e3: # if sim duration=1sec, assume its the test trial and select PMd spikes based on targetid - spktPMd = [rawSpikesPMd[s.targetid][s.repeatSingleTrials[s.targetid]][icell]['spkt'][0] for icell in range(numrawcells)] + spktPMd = [rawSpikesPMd[s.targetid][s.repeatSingleTrials[s.targetid]][icell]['spkt'][0] for icell in range(numrawcells)] elif s.repeatSingleTrials[0] > -1: # use single trials for each target during training spktPMd = [] for icell in range(numrawcells): spkt = [] [spkt.extend(list(rawSpikesPMd[itarget][s.repeatSingleTrials[itarget]][icell]['spkt'][0] + (1000*itrial))) \ for itrial,itarget in enumerate(s.trialTargets[:-1])] # replicate spike times over trials - spktPMd.append(spkt) + spktPMd.append(spkt) else: # use all available trials for each target during training pass # play back PMd spikes using VecStims @@ -773,9 +774,9 @@ def runSim(): while round(h.t) < s.duration: run(min(s.duration,h.t+s.loopstep)) # MPI: Get ready to run the simulation (it isn't actually run until pc.runworker() is called I think) if s.server.simMode == 0: - if s.rank==0 and (round(h.t) % s.progupdate)==0: print(' t = %0.1f s (%i%%; time consumed: %0.1f s)' % (h.t/1e3, int(h.t/s.duration*100), (time()-runstart))) + if s.rank==0 and (round(h.t) % s.progupdate)==0: print((' t = %0.1f s (%i%%; time consumed: %0.1f s)' % (h.t/1e3, int(h.t/s.duration*100), (time()-runstart)))) else: - if s.rank==0: print(' t = %0.1f s (%i%%; time consumed: %0.1f s)' % (h.t/1e3, int(h.t/s.duration*100), (time()-runstart))) + if s.rank==0: print((' t = %0.1f s (%i%%; time consumed: %0.1f s)' % (h.t/1e3, int(h.t/s.duration*100), (time()-runstart)))) # Calculate LFP -- WARNING, need to think about how to optimize if s.savelfps: @@ -787,11 +788,11 @@ def runSim(): tmplfps[pop] += s.cells[id].V # Add voltage to LFP estimate if s.verbose: if s.server.Manager.ns.isnan(tmplfps[pop]) or s.server.Manager.ns.isinf(tmplfps[pop]): - print "Nan or inf" + print("Nan or inf") s.hostlfps.append(tmplfps) # Add voltages # Periodic weight saves - if s.usestdp: + if s.usestdp: timesincelastsave = h.t - s.timeoflastsave if timesincelastsave >= s.timebetweensaves: s.timeoflastsave = h.t @@ -799,8 +800,8 @@ def runSim(): for ps in range(s.nstdpconns): if s.stdpmechs[ps].synweight != s.weightchanges[ps][-1][-1]: # Only store connections that changed; [ps] = this connection; [-1] = last entry; [-1] = weight s.weightchanges[ps].append([s.timeoflastsave, s.stdpmechs[ps].synweight]) - - ## Virtual arm + + ## Virtual arm if s.useArm != 'None': armStart = time() s.arm.run(h.t, s) # run virtual arm apparatus (calculate command, move arm, feedback) @@ -824,30 +825,30 @@ def runSim(): #sleep(0.001) #print 'stdp_after: ', stdp.synweight # Synaptic scaling? - + #print(' Arm time = %0.4f s') % (time() - armStart) ## Time adjustment for online mode simulation - if s.PMdinput == 'Plexon' and s.server.simMode == 1: + if s.PMdinput == 'Plexon' and s.server.simMode == 1: # To avoid izhi cell's over shooting when h.t moves forward because sim is slow. - for c in range(s.cellsperhost): + for c in range(s.cellsperhost): gid = s.gidVec[c] if s.cellnames[gid] == 'PMd': # 'PMds don't have t0 variable. continue - s.cells[c].t0 = s.server.newCurrTime.value - h.dt + s.cells[c].t0 = s.server.newCurrTime.value - h.dt dtSave = h.dt # save original dt h.dt = s.server.newCurrTime.value - h.t # new dt active = h.cvode.active() if active != 0: - h.cvode.active(0) + h.cvode.active(0) h.fadvance() # Integrate with new dt if active != 0: - h.cvode.active(1) - h.dt = dtSave # Restore orignal dt - - if s.rank==0: + h.cvode.active(1) + h.dt = dtSave # Restore orignal dt + + if s.rank==0: s.runtime = time()-runstart # See how long it took - print(' Done; run time = %0.1f s; real-time ratio: %0.2f.' % (s.runtime, s.duration/1000/s.runtime)) + print((' Done; run time = %0.1f s; real-time ratio: %0.2f.' % (s.runtime, s.duration/1000/s.runtime))) s.pc.barrier() # Wait for all hosts to get to this point @@ -855,7 +856,7 @@ def runSim(): ### Finalize Simulation (gather data from nodes, etc.) ############################################################################### def finalizeSim(): - + ## Variables to unpack data from all hosts ## Pack data from all hosts @@ -891,7 +892,7 @@ def finalizeSim(): s.totalspikes = 0 # Keep a running tally of the number of spikes s.totalconnections = 0 # Total number of connections s.totalstdpconns = 0 # Total number of stdp connections - if s.saveraw: s.allraw = [] + if s.saveraw: s.allraw = [] for host in range(s.nhosts): # Loop over hosts s.pc.take(host) # Get the last message hostdata = s.pc.upkpyobj() # Unpack them @@ -908,7 +909,7 @@ def finalizeSim(): s.totalspikes = len(s.allspiketimes) # Keep a running tally of the number of spikes s.totalconnections = len(s.allconnections[0]) # Total number of connections s.totalstdpconns = len(s.allstdpconndata) # Total number of STDP connections - + # Record background spike data (cliff: only for one node since takes too long to pack for all and just needed for debugging) if s.savebackground and s.usebackground: @@ -920,7 +921,7 @@ def finalizeSim(): s.allbackgroundspikecells = concatenate((s.allbackgroundspikecells, c+zeros(len(thesespikes)))) # Add this cell's ID to the list s.backgrounddata = transpose(vstack([s.allbackgroundspikecells,s.allbackgroundspiketimes])) else: s.backgrounddata = [] # For saving s no error - + if s.saveraw and s.usestims: s.allstimspikecells=array([]) s.allstimspiketimes=array([]) @@ -932,7 +933,7 @@ def finalizeSim(): else: s.stimspikedata = [] # For saving so no error gathertime = time()-gatherstart # See how long it took - if s.rank==0: print(' Done; gather time = %0.1f s.' % gathertime) + if s.rank==0: print((' Done; gather time = %0.1f s.' % gathertime)) s.pc.barrier() #mindelay = s.pc.allreduce(s.pc.set_maxstep(10), 2) # flag 2 returns minimum value @@ -955,11 +956,11 @@ def finalizeSim(): print('\nAnalyzing...') s.firingrate = float(s.totalspikes)/s.ncells/s.duration*1e3 # Calculate firing rate -- confusing but cool Python trick for iterating over a list s.connspercell = s.totalconnections/float(s.ncells) # Calculate the number of connections per cell - print(' Run time: %0.1f s (%i-s sim; %i scale; %i cells; %i workers)' % (s.runtime, s.duration/1e3, s.scale, s.ncells, s.nhosts)) - print(' Spikes: %i (%0.2f Hz)' % (s.totalspikes, s.firingrate)) - print(' Connections: %i (%i STDP; %0.2f per cell)' % (s.totalconnections, s.totalstdpconns, s.connspercell)) - print(' Mean connection distance: %0.2f um' % mean(s.allconnections[2])) - print(' Mean connection delay: %0.2f ms' % mean(s.allconnections[3])) + print((' Run time: %0.1f s (%i-s sim; %i scale; %i cells; %i workers)' % (s.runtime, s.duration/1e3, s.scale, s.ncells, s.nhosts))) + print((' Spikes: %i (%0.2f Hz)' % (s.totalspikes, s.firingrate))) + print((' Connections: %i (%i STDP; %0.2f per cell)' % (s.totalconnections, s.totalstdpconns, s.connspercell))) + print((' Mean connection distance: %0.2f um' % mean(s.allconnections[2]))) + print((' Mean connection delay: %0.2f ms' % mean(s.allconnections[3]))) ############################################################################### @@ -968,26 +969,26 @@ def finalizeSim(): def saveData(): if s.rank == 0: ## Save to txt file (spikes and conn) - if s.savetxt: + if s.savetxt: filename = '../data/m1ms-spk.txt' fd = open(filename, "w") for c in range(len(s.allspiketimes)): - print >> fd, int(s.allspikecells[c]), s.allspiketimes[c], s.popNamesDic[s.cellnames[int(s.allspikecells[c])]] + print(int(s.allspikecells[c]), s.allspiketimes[c], s.popNamesDic[s.cellnames[int(s.allspikecells[c])]], file=fd) fd.close() - print "[Spikes are stored in", filename, "]" + print("[Spikes are stored in", filename, "]") if s.verbose: filename = 'm1ms-conn.txt' fd = open(filename, "w") for c in range(len(s.allconnections[0])): - print >> fd, int(s.allconnections[0][c]), int(s.allconnections[1][c]), s.allconnections[2][c], s.allconnections[3][c], s.allconnections[4][c] + print(int(s.allconnections[0][c]), int(s.allconnections[1][c]), s.allconnections[2][c], s.allconnections[3][c], s.allconnections[4][c], file=fd) fd.close() - print "[Connections are stored in", filename, "]" + print("[Connections are stored in", filename, "]") ## Save to mat file if s.savemat: savestart = time() # See how long it takes to save - + # Save simulation code filestosave = [] #'main.py', 'shared.py', 'network.py', 'arm.py', 'arminterface.py', 'server.py', 'izhi.py', 'izhi2007.mod', 'stdp.mod', 'nsloc.py', 'nsloc.mod'] # Files to save argv = []; @@ -996,7 +997,7 @@ def saveData(): fobj = open(filestosave[f]) # Open it for reading simcode.append(fobj.readlines()) # Append to list of code to save fobj.close() # Close file object - + # Tidy variables spikedata = vstack([s.allspikecells,s.allspiketimes]).T # Put spike data together connections = vstack([s.allconnections[0],s.allconnections[1]]).T # Put connection data together @@ -1008,13 +1009,13 @@ def saveData(): # Save variables info = {'timestamp':datetime.today().strftime("%d %b %Y %H:%M:%S"), 'runtime':s.runtime, 'popnames':s.popnames, 'popEorI':s.popEorI} # Save date, runtime, and input arguments - + targetPos = s.arm.targetPos handPosAll = s.arm.handPosAll angAll = s.arm.angAll - motorCmdAll = s.arm.motorCmdAll - targetidAll = s.arm.targetidAll - errorAll = s.arm.errorAll + motorCmdAll = s.arm.motorCmdAll + targetidAll = s.arm.targetidAll + errorAll = s.arm.errorAll criticAll = s.arm.criticAll if not hasattr(s, 'phase'): s.phase = '' s.filename = s.outfilestem+'_target_'+str(s.arm.targetid)+s.phase @@ -1022,25 +1023,25 @@ def saveData(): variablestosave = ['targetPos', 'angAll', 'motorCmdAll', 'errorAll'] else: variablestosave = ['info', 'targetPos', 'angAll', 'motorCmdAll', 'errorAll', 'simcode', 'spikedata', 's.cellpops', 's.cellnames', 's.cellclasses', 's.xlocs', 's.ylocs', 's.zlocs', 'connections', 'distances', 'delays', 'weights', 's.EorI'] - - if s.savelfps: - variablestosave.extend(['s.lfptime', 's.lfps']) - if s.usestdp: + + if s.savelfps: + variablestosave.extend(['s.lfptime', 's.lfps']) + if s.usestdp: variablestosave.extend(['stdpdata', 's.allweightchanges']) if s.savebackground: variablestosave.extend(['s.backgrounddata']) - if s.saveraw: + if s.saveraw: variablestosave.extend(['s.stimspikedata', 's.allraw']) if s.usestims: variablestosave.extend(['stimdata']) savecommand = "savemat(s.filename, {" for var in range(len(variablestosave)): savecommand += "'" + variablestosave[var].replace('s.','') + "':" + variablestosave[var] + ", " # Create command out of all the variables savecommand = savecommand[:-2] + "}, oned_as='column')" # Omit final comma-space and complete command - - print('Saving output as %s...' % s.filename) + + print(('Saving output as %s...' % s.filename)) exec(savecommand) # Actually perform the save savetime = time()-savestart # See how long it took to save - print(' Done; time = %0.1f s' % savetime) + print((' Done; time = %0.1f s' % savetime)) ############################################################################### @@ -1050,12 +1051,12 @@ def plotData(): ## Plotting if s.rank == 0: if s.plotraster: # Whether or not to plot - if (s.totalspikes>s.maxspikestoplot): + if (s.totalspikes>s.maxspikestoplot): disp(' Too many spikes (%i vs. %i)' % (s.totalspikes, s.maxspikestoplot)) # Plot raster, but only if not too many spikes - elif s.nhosts>1: - disp(' Plotting raster despite using too many cores (%i)' % s.nhosts) + elif s.nhosts>1: + disp(' Plotting raster despite using too many cores (%i)' % s.nhosts) analysis.plotraster()#;allspiketimes, allspikecells, EorI, ncells, connspercell, backgroundweight, firingrate, duration) - else: + else: print('Plotting raster...') analysis.plotraster()#allspiketimes, allspikecells, EorI, ncells, connspercell, backgroundweight, firingrate, duration) @@ -1081,6 +1082,3 @@ def plotData(): analysis.plot3darch() show(block=False) - - - diff --git a/nsloc.py b/nsloc.py index 6b931ae..7950526 100644 --- a/nsloc.py +++ b/nsloc.py @@ -3,7 +3,7 @@ Usage example: from neuron import h - from nsloc import nsloc + from nsloc import nsloc cell = nsloc(cellid) Version: 2014July23 by salvadordura@gmail.com diff --git a/server.py b/server.py index c1e04eb..7ed90f9 100644 --- a/server.py +++ b/server.py @@ -1,7 +1,7 @@ """ -server.py +server.py -Interface between model and plexon +Interface between model and plexon Developed by giljael Version: 2015feb14 modified by salvadord @@ -17,7 +17,7 @@ import time import os import sys -import Queue +import queue import array import math import numpy as np @@ -35,10 +35,10 @@ lock = Lock() # When there is no data from the plexon system, -# the client will send 'exit' to the server and +# the client will send 'exit' to the server and # terminate itserlf after timeOut * timeWnd time. # ex) 100 * 500 = 50 seconds -timeOut = 500 #10000 +timeOut = 500 #10000 # header ID NODATA = 0 @@ -50,7 +50,7 @@ # channel start CH_START = 1 -# channel end +# channel end CH_END = 96 # unit end UNIT_END = 4 @@ -71,11 +71,11 @@ ##################################################################### # Variables user may modify # ##################################################################### -# port number for sending initial data to the client from the server. +# port number for sending initial data to the client from the server. # When it is changed, initPort in client.m should be changed as well. -initPort = 7869 +initPort = 7869 # port number for receiving spikes from the client -comPort = 9999 +comPort = 9999 # Time window for binning and summing spikes binWnd = 100 # (ms) autorun = 1 # autorun- 0: interactive mode 1: autorun upto @@ -95,7 +95,7 @@ isUdp = 0 # 0: TCP, 1: UDP simMode = 0 # simMode- Offline mode(0), online mode(1) LR = 1000e3 -unsortedSpkFilter = 0 # 0: counts unsorted unit, 1: filter unsorted unit +unsortedSpkFilter = 0 # 0: counts unsorted unit, 1: filter unsorted unit syncSpkFilter = 0 # 0: default, no filtering, 1: filter spikes having identical timestamps, 2: filter spikes in a same time window. isDp = 0 # 0: NSLOC, 1:DP isLwc = 0 # 0: Heavy weight client, 1: light weight client @@ -104,17 +104,17 @@ TIMEOUT = 20 * 0.001 # it should be second syncWnd = 0.001 -# For latency measure +# For latency measure qTimeList = [] deqTimeList = [] if verbose: fdtime = open('SpikeTimeStamp.tsv', 'w') - fdtime.write('CH#\tUnit#\tTS\n') + fdtime.write('CH#\tUnit#\tTS\n') fdtime2 = open('SpikeTimeStamp2.tsv', 'w') - fdtime2.write('CH#\tUnit#\tTS\tSimTS\n') - fdtime3 = open('log.tsv', 'w') - + fdtime2.write('CH#\tUnit#\tTS\tSimTS\n') + fdtime3 = open('log.tsv', 'w') + ##################################################################### # Parameter validation # ##################################################################### @@ -124,7 +124,7 @@ if isCommunication == 0 or simMode == 1: localFileRead = 0 # queue test w/ no communication is useless -if isDp == 1: +if isDp == 1: if simMode == 1: LR = 500 # latency requirement (ms) else: # NSLOC @@ -134,17 +134,17 @@ if simMode == 0: LR = 1000e3 -if unsortedSpkFilter == 0: +if unsortedSpkFilter == 0: UNIT_START = 0 else : UNIT_START = 1 # Print out server's information -def getServerInfo(): +def getServerInfo(): if simMode == 1: - fname = "Online/" + fname = "Online/" else: - fname = "Offline/" + fname = "Offline/" if isUdp == 1: fname += "UDP/" else: @@ -160,18 +160,18 @@ def getServerInfo(): if unsortedSpkFilter == 1: fname += "unsortedSpkFilter(" if isLwc == 1: - fname += "Lwc)/" + fname += "Lwc)/" else: fname += "Hwc)/" if syncSpkFilter: fname += "syncSpkFilter(" if isLwc == 1: - fname += "Lwc)/" + fname += "Lwc)/" else: fname += "Hwc)/" if verbose: - fname += "Verbose/" - print fname + "server is launched!!!" + fname += "Verbose/" + print(fname + "server is launched!!!") @@ -182,12 +182,12 @@ def checkLocalIndexbyKey(dicName, key): except: index = -1 return index - + def feedQueue(qItem, theVerbose = 0): if isDp: # DP model currQTime = qItem[CH_END] * binWnd else: # NSLOC model - spkNum = qItem[SPKSZ - 1] + spkNum = qItem[SPKSZ - 1] if spkNum == 0: # timeout item currQTime = qItem[SPKSZ - 2] * 1000 # sec -> ms else: @@ -200,12 +200,12 @@ def feedQueue(qItem, theVerbose = 0): with lock: currSimTime = currNeuronTime.value if verbose: - print "[feedQueue] currQTime:", currQTime, "currSimTime: ", currSimTime, "LR: ", LR + print("[feedQueue] currQTime:", currQTime, "currSimTime: ", currSimTime, "LR: ", LR) if (currQTime - currSimTime) > LR: # Simulation is slow if qsize == queue.qsize(): # callback doesn't get an item in Q so far queue.get() # discard an item if verbose: - print "[feedQueue] An item is discared in Q" + print("[feedQueue] An item is discared in Q") else: # simulation is fast break @@ -216,15 +216,15 @@ def __init__(self, dir="data"): if isDp == 1: # binned data self.qItem = [0] * (CH_END + 1) self.invl = binWnd - self.vec = h.mua # mua : vector in Hoc + self.vec = h.mua # mua : vector in Hoc else: # spikes self.qItem = [0] * (SPKSZ + 1) # = 3 * SPKNUM + 1 + Serial Number - self.invl = 0.025 + self.invl = 0.025 if pc.id() == 0: self.fih = h.FInitializeHandler(1, self.callback) else: self.fih = h.FInitializeHandler(1, self.callbackSlave) - # for timing measurement + # for timing measurement self.num = 0 self.write = 1 @@ -233,10 +233,10 @@ def __init__(self, dir="data"): if isDp == 0: # read spikes from a file by the server w/o communication self.fname = self.dir + "/spikePMd-6sec.tsv" self.fd = open(self.fname, "r") - print "[Open ", self.fname, " for test]" + print("[Open ", self.fname, " for test]") else: # DP pass - + # callback function for the DP (isDp = 1) and NSLOC (isDp = 0) model def callbackSlave(self): try: @@ -246,11 +246,11 @@ def callbackSlave(self): else: # DP pass # get an item from the queue. Raise exception when it fails - else: # DP and NSLOC with online and offline mode - try: + else: # DP and NSLOC with online and offline mode + try: nextInvl = self.invl - currSimTime = h.t - newCurrTime.value = currSimTime + currSimTime = h.t + newCurrTime.value = currSimTime if isDp == 1: self.qItem = queue.get() else: @@ -258,7 +258,7 @@ def callbackSlave(self): #except Queue.Empty as e: # No item in Q except: # No item in Q if verbose: - print "[callbackSlave] No item in Q" + print("[callbackSlave] No item in Q") else: if isDp == 0 and not n == 0: #NSLOC #self.qItem = vec.to_python() @@ -274,13 +274,13 @@ def callbackSlave(self): for item in deqTimeList: f1.write("%s\n" % item) f1.close() - print "=======NslocPulltime printed!!!" + print("=======NslocPulltime printed!!!") spkNum = int(s.vec[SPKSZ - 1]) # update current neuron time if spkNum == 0: # timeout item currSimTime = self.qItem[SPKSZ - 2] * 1000 + h.t # sec -> ms else: - #currSimTime = self.qItem[(spkNum - 1) * 3 + 2] * 1000 # sec -> ms + #currSimTime = self.qItem[(spkNum - 1) * 3 + 2] * 1000 # sec -> ms currSimTime = s.vec[(spkNum - 1) * 3 + 2] * 1000 # sec -> ms # [0]:CH_ID, [1]:Unit_ID, [2]: Time_stamp for i in range(spkNum): @@ -294,37 +294,37 @@ def callbackSlave(self): #h.s.inncl.o(localIndex).event(timeStamp, 1) s.inncl[localIndex].event(timeStamp, 1) if verbose: - print "[CallbackSlave] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2] * 1000, "t: ", h.t + print("[CallbackSlave] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2] * 1000, "t: ", h.t) #fdtime.write(str(int(self.qItem[3*i])) + '\t' +str(int(self.qItem[3*i+1])) + '\t' +str('{:g}'.format(self.qItem[3*i+2])) + '\n') else: # spike ignored. For dubug if verbose: - print "[CallbackSlave] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2], "t: ", h.t + print("[CallbackSlave] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2], "t: ", h.t) #fdtime2.write(str(int(self.qItem[3*i])) + '\t' +str(int(self.qItem[3*i+1])) + '\t' +str('{:g}'.format(self.qItem[3*i+2])) + '\t' + str(h.t/1000) + '\n') # move t to currSimTime #if currSimTime > h.t and currSimTime < duration: if currSimTime > h.t: if simMode == 1: # online mode - newCurrTime.value = currSimTime + newCurrTime.value = currSimTime else: currSimTime = h.t newCurrTime.value = currSimTime finally: if verbose: - print "[CallbackSlave] currNeuronTime: ", currSimTime + print("[CallbackSlave] currNeuronTime: ", currSimTime) except: if isCommunication == 0: self.fd.close() if verbose: - print "[CallbackSlave] exception occurs:", traceback.print_exc() + print("[CallbackSlave] exception occurs:", traceback.print_exc()) finally: # update current neuron time if isCommunication == 0: - pass + pass else: if verbose: - print('[callback.h.t End(id= %i)] %.3f' %(pc.id(), h.t)) + print(('[callback.h.t End(id= %i)] %.3f' %(pc.id(), h.t))) if newCurrTime.value + nextInvl <= duration: - h.cvode.event(newCurrTime.value + nextInvl, self.callbackSlave) - + h.cvode.event(newCurrTime.value + nextInvl, self.callbackSlave) + # callback function for the DP (isDp = 1) and NSLOC (isDp = 0) model def callback(self): try: @@ -335,68 +335,68 @@ def callback(self): localts = locals.split() if localts[0] == '1': # spikes if unsortedSpkFilter and localts[2] == '0': # filter unsorted spikes - continue + continue else: tt = float(localts[3]) * 1000 if tt > h.t: s.inncl[int(localts[1]) - 1].event(tt, 1) if verbose: - fdtime.write(str(1) + '\t' + str(int(localts[1])) + '\t' +str(int(localts[2])) + '\t' +str('{:g}'.format(float(localts[3]))) + '\n') + fdtime.write(str(1) + '\t' + str(int(localts[1])) + '\t' +str(int(localts[2])) + '\t' +str('{:g}'.format(float(localts[3]))) + '\n') else: # DP pass # get an item from the queue. Raise exception when it fails - else: # DP and NSLOC with online and offline mode - try: + else: # DP and NSLOC with online and offline mode + try: nextInvl = self.invl - currSimTime = h.t + currSimTime = h.t newCurrTime.value = currSimTime if timeMeasure: - timer = time.time() + timer = time.time() if isDp == 1: self.qItem = queue.get() else: self.qItem = queue.get(False) # Q has an item? n = pc.broadcast(s.vec.from_python(self.qItem), 0) # convert python list to hoc vector - except Queue.Empty as e: # No item in Q + except queue.Empty as e: # No item in Q n = pc.broadcast(s.emptyVec, 0) if verbose: - print "[callback] No item in Q" + print("[callback] No item in Q") else: if isDp == 1: - if timeMeasure and self.write: - deqTimeList.append(time.time() - timer) - self.num += 1 - if self.num == iteration: - self.write = 0 - if not os.path.exists(self.dir): - os.mkdir(self.dir) - fname = self.dir + "/DpPullTimefile" - f1 = open(fname, "w+") - for item in deqTimeList: - f1.write("%s\n" % item) - f1.close() - print "=======DpPull time printed!!!" + if timeMeasure and self.write: + deqTimeList.append(time.time() - timer) + self.num += 1 + if self.num == iteration: + self.write = 0 + if not os.path.exists(self.dir): + os.mkdir(self.dir) + fname = self.dir + "/DpPullTimefile" + f1 = open(fname, "w+") + for item in deqTimeList: + f1.write("%s\n" % item) + f1.close() + print("=======DpPull time printed!!!") # for debug #if 0 and row > 0:#verbose: # f_handle = file('spk_rcv.txt', 'a') # for item in self.qItem: # f_handle.write("%d " % item) # f_handle.write("\n") - # f_handle.close() + # f_handle.close() currSimTime = self.qItem[CH_END] * binWnd if self.qItem[0] == -1: # timeout message in Q - print "[callback-DP] timeout message" + print("[callback-DP] timeout message") h.setUpdateDPzero() #h.updateDpWithZero() else: # data in a queue item in Q # Convert the python list to the hoc vector (h.mua) - self.vec = self.vec.from_python(self.qItem) + self.vec = self.vec.from_python(self.qItem) h.updateDpWithMua() else: #NSLOC - if timeMeasure and self.write: - deqTimeList.append(time.time() - timer) + if timeMeasure and self.write: + deqTimeList.append(time.time() - timer) self.num += 1 if self.num == iteration: - self.write = 0 + self.write = 0 if not os.path.exists(self.dir): os.mkdir(self.dir) fname = self.dir + "/NslocPullTimefile" @@ -404,7 +404,7 @@ def callback(self): for item in deqTimeList: f1.write("%s\n" % item) f1.close() - print "=======NslocPulltime printed!!!" + print("=======NslocPulltime printed!!!") spkNum = self.qItem[SPKSZ - 1] # update current neuron time if spkNum == 0: # timeout item @@ -413,46 +413,46 @@ def callback(self): currSimTime = self.qItem[(spkNum - 1) * 3 + 2] * 1000 # sec -> ms # [0]:CH_ID, [1]:Unit_ID, [2]: Time_stamp for i in range(spkNum): - timeStamp = self.qItem[3 * i + 2] * 1000 # second -> ms + timeStamp = self.qItem[3 * i + 2] * 1000 # second -> ms if h.t <= timeStamp and timeStamp <= duration: # queue spikes in the NEURON queue and ignore old spikes # Queue netcon event. CH_START = 1, CH_END = 96 localIndex = checkLocalIndexbyKey(s.innclDic, int(self.qItem[3 * i] - 1)) if not localIndex == -1: - s.inncl[localIndex].event(timeStamp, 1) - if verbose: - print "[Callback] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2] * 1000, "t: ", h.t + s.inncl[localIndex].event(timeStamp, 1) + if verbose: + print("[Callback] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2] * 1000, "t: ", h.t) #fdtime.write(str(int(self.qItem[3*i])) + '\t' +str(int(self.qItem[3*i+1])) + '\t' +str('{:g}'.format(self.qItem[3*i+2])) + '\n') else: # spike ignored. For dubug if verbose: - print "[Callback] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2], "t: ", h.t + print("[Callback] ChID:", self.qItem[3 * i], "UnitID:", self.qItem[3 * i + 1], "Time:", self.qItem[3 * i + 2], "t: ", h.t) #fdtime.write(str(int(self.qItem[3*i])) + '\t' +str(int(self.qItem[3*i+1])) + '\t' +str('{:g}'.format(self.qItem[3*i+2])) + '\t' + str('{:g}'.format(h.t/1000)) + '\n') # move t to currSimTime if currSimTime > h.t:# and currSimTime < duration: - if simMode == 1: # online mode + if simMode == 1: # online mode newCurrTime.value = currSimTime else: - currSimTime = h.t + currSimTime = h.t newCurrTime.value = currSimTime finally: with lock: currNeuronTime.value = currSimTime # (ms) if verbose: - print "[Callback] currNeuronTime: ", currSimTime, nextInvl + print("[Callback] currNeuronTime: ", currSimTime, nextInvl) except: if isCommunication == 0: self.fd.close() if verbose: - print "[Callback] exception occurs:", traceback.print_exc() + print("[Callback] exception occurs:", traceback.print_exc()) finally: # update current neuron time if isCommunication == 0: pass else: if newCurrTime.value + nextInvl <= duration: - h.cvode.event(newCurrTime.value + nextInvl, self.callback) + h.cvode.event(newCurrTime.value + nextInvl, self.callback) if verbose: - print('[callback.h.t End(id= %i)] %.3f' %(pc.id(), h.t)) + print(('[callback.h.t End(id= %i)] %.3f' %(pc.id(), h.t))) fdtime2.write('Invl: ' + str(nextInvl) + '\t' + 'h.t + Invl: ' + str(h.t + nextInvl) + '\n') - fdtime3.write('invl, h.t + nextInvl: ' + str(nextInvl) + ' ' + str(h.t + nextInvl) + '\n') + fdtime3.write('invl, h.t + nextInvl: ' + str(nextInvl) + ' ' + str(h.t + nextInvl) + '\n') # Connect python server to plexon client. # It waits for receiving an "ALIVE" message from the client. @@ -469,49 +469,49 @@ def nrn_py_connectPlxClient(): # receive ALIVE from the client data, addr = Sock.recvfrom(buf) - print data + print(data) dataLen = len(data) dlen = dataLen / 2 if dataLen == 2: - cdata = struct.unpack('H' * dlen, data) + cdata = struct.unpack('H' * dlen, data) if cdata[0] == ALIVE: if verbose: - print "I received ALIVE" + print("I received ALIVE") # send INITDATA to the client a = array.array('H') a.extend([INITDATA, comPort, isUdp, isDp, isLwc, binWnd, timeOut, unsortedSpkFilter, syncSpkFilter, verbose]) Sock.sendto(a.tostring(), addr) if verbose: - print "I sent INITDATA" + print("I sent INITDATA") data, addr = Sock.recvfrom(buf) dataLen = len(data) dlen = dataLen / 2 cdata = struct.unpack('H' * dlen, data) - if cdata[0] == ACK: - print "Connection success!!!" + if cdata[0] == ACK: + print("Connection success!!!") else: - print "ERROR: Connection failure!!!" + print("ERROR: Connection failure!!!") raise - Sock.close() - + Sock.close() + # Create socket and bind to address for spike exchange addr = (host, comPort) if isUdp: Sock = socket(AF_INET, SOCK_DGRAM) - Sock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) - Sock.bind(addr) + Sock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) + Sock.bind(addr) return (Sock, 0, 0) else: Sock = socket(AF_INET, SOCK_STREAM) - Sock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) + Sock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) Sock.bind(addr) Sock.listen(1) conn, addr = Sock.accept() return (Sock, conn, addr) except: - print "[nrn_py_connectPlxClient] exception occurs:", sys.exc_info()[0] + print("[nrn_py_connectPlxClient] exception occurs:", sys.exc_info()[0]) sys.exit(0) - + # Clarify spikes and external events. # length(cdata) = drows * 4 fields : |data type|ch#|unit#|TS| @@ -525,8 +525,8 @@ def nrn_py_clarifyEvt(cdataNp): spkArr = cdataNp[cdataNp[:, 0] == 1] # spikes evtArr = cdataNp[cdataNp[:, 0] == 4] # external events except: - print "[nrn_py_clarifyEvt] exception occurs:", sys.exc_info()[0] - finally: + print("[nrn_py_clarifyEvt] exception occurs:", sys.exc_info()[0]) + finally: return (spkArr, evtArr) # filter system events @@ -537,18 +537,18 @@ def nrn_py_filterSysEvt(evtArr): if row > 0: newEvtArr = evtArr[evtArr[:,1] == 257] if verbose: - print "nrn_py_filterSysEvt returns" + print("nrn_py_filterSysEvt returns") except: - print "[nrn_py_filterSysEvt] exception occurs:", sys.exc_info()[0] - finally: + print("[nrn_py_filterSysEvt] exception occurs:", sys.exc_info()[0]) + finally: return newEvtArr - - -# reorder spikes + + +# reorder spikes # orderedArr = nrn_py_reorder(spkArr) def nrn_py_reorder(spkArr): if not hasattr(nrn_py_reorder, "keptArr"): - nrn_py_reorder.keptArr = np.zeros((0, 4)) # this is a static numpy array + nrn_py_reorder.keptArr = np.zeros((0, 4)) # this is a static numpy array try: row, col = nrn_py_reorder.keptArr.shape orderedArr = np.zeros((0, 4)) @@ -556,7 +556,7 @@ def nrn_py_reorder(spkArr): maxTS = nrn_py_reorder.keptArr[row - 1][3] else: maxTS = 0 - #concatenates + #concatenates newArr = np.concatenate((nrn_py_reorder.keptArr, spkArr), axis = 0) row, col = newArr.shape if row > 0: @@ -565,10 +565,10 @@ def nrn_py_reorder(spkArr): nrn_py_reorder.keptArr = newArr[newArr[:,3] > maxTS] orderedArr = newArr[newArr[:, 3] <= maxTS] except: - print "[nrn_py_reorder] exception occurs:", sys.exc_info()[0] + print("[nrn_py_reorder] exception occurs:", sys.exc_info()[0]) finally: return orderedArr - + # Filter sync spikes # syncSpkFilter = 1: filter spikes having same timestamp # syncSpkFilter = 2: filter spikes in a sync window [syncWnd * x, syncWnd * (x + 1)) @@ -584,18 +584,18 @@ def nrn_py_filterSyncSpk(spkArr): if row > 0: if row == 1: # one element in spkArr newFltdArr = spkArr - else: # more than one in spkArr + else: # more than one in spkArr sync = 0 - baseValue = np.array(spkArr[0], ndmin = 2) # first elem in spkArr + baseValue = np.array(spkArr[0], ndmin = 2) # first elem in spkArr for ii in range(1, row): # 1, 2, 3, ..., m-1 if baseValue[0][3] != spkArr[ii][3]: if sync == 0: - newFltdArr = np.concatenate((newFltdArr, baseValue), axis = 0) + newFltdArr = np.concatenate((newFltdArr, baseValue), axis = 0) else: sync = 0 else: sync = 1 - baseValue = np.array(spkArr[ii], ndmin = 2) + baseValue = np.array(spkArr[ii], ndmin = 2) if sync == 0: newFltdArr = np.concatenate((newFltdArr, baseValue), axis = 0) else: @@ -606,26 +606,26 @@ def nrn_py_filterSyncSpk(spkArr): if row > 0: # sync window size is syncWnd # find the biggest "j" - maxTS = newArr[row - 1][3] + maxTS = newArr[row - 1][3] newJ = nrn_py_filterSyncSpk.jValue + 1 while 1: - if syncWnd * (newJ - 1) <= maxTS and maxTS < (syncWnd * newJ): + if syncWnd * (newJ - 1) <= maxTS and maxTS < (syncWnd * newJ): break else: newJ += 1 nrn_py_filterSyncSpk.FRem = newArr[newArr[:, 3] >= syncWnd * (newJ - 1)] for ii in range(nrn_py_filterSyncSpk.jValue, newJ - 1): - newArr = newArr[newArr[:, 3] >= syncWnd * ii] - temp = newArr[newArr[:, 3] < syncWnd * (ii + 1)] - row, col = temp.shape - if row == 1: - newFltdArr = np.concatenate((newFltdArr, temp), axis = 0) + newArr = newArr[newArr[:, 3] >= syncWnd * ii] + temp = newArr[newArr[:, 3] < syncWnd * (ii + 1)] + row, col = temp.shape + if row == 1: + newFltdArr = np.concatenate((newFltdArr, temp), axis = 0) nrn_py_filterSyncSpk.jValue = newJ - 1 else: - print "Wrong filtering options" + print("Wrong filtering options") except: - print "[nrn_py_filterSyncSpk] exception occurs:", sys.exc_info()[0] - finally: + print("[nrn_py_filterSyncSpk] exception occurs:", sys.exc_info()[0]) + finally: return newFltdArr @@ -636,49 +636,49 @@ def nrn_py_filterUnsortedSpk(spkArr): if row > 0: newFltdArr = spkArr[spkArr[:, 2] > 0] except: - print "[nrn_py_filterUnsortedSpk] exception occurs:", sys.exc_info()[0] + print("[nrn_py_filterUnsortedSpk] exception occurs:", sys.exc_info()[0]) finally: - return newFltdArr + return newFltdArr ## # mua, binningComplete = nrn_py_binSpk(spkArr, binStart, timeoutFlag) # Here mua is a numpy array, so it should be converted to python list. -def nrn_py_binSpk(spkArr, binStart, timeoutFlag): +def nrn_py_binSpk(spkArr, binStart, timeoutFlag): try: # MUA in 1st, and 2nd time interval if not hasattr(nrn_py_binSpk, "binnedMsg"): - nrn_py_binSpk.binnedMsg = np.zeros((CH_END + 1, 1)) # 96 channel + binWnd + nrn_py_binSpk.binnedMsg = np.zeros((CH_END + 1, 1)) # 96 channel + binWnd if not hasattr(nrn_py_binSpk, "dataHave"): nrn_py_binSpk.dataHave = 0 if not hasattr(nrn_py_binSpk, "BRem"): nrn_py_binSpk.BRem = np.zeros((0, 4)) # spikes - + #if unsortedSpkFilter == 0: # UNIT_START = 0 #else: - # UNIT_START = 1 - + # UNIT_START = 1 + newBinnedMsg = np.zeros((CH_END + 1, 1)) binningComplete = 0 binEnd = (binStart + binWnd)/1000.0 binStart = binStart/1000.0 - + # concatenate - newArr = np.concatenate((nrn_py_binSpk.BRem, spkArr), axis = 0) + newArr = np.concatenate((nrn_py_binSpk.BRem, spkArr), axis = 0) nrn_py_binSpk.BRem = np.zeros((0, 4)) row, col = newArr.shape if row > 0: nrn_py_binSpk.BRem = newArr[newArr[:, 3] >= binEnd] newArr = newArr[newArr[:, 3]< binEnd] - row, col = newArr.shape + row, col = newArr.shape if row > 0: - for j in range(row): + for j in range(row): channelID = newArr[j][1] unitID = newArr[j][2] timeStamp = newArr[j][3] #(ms) - #Timestamp is in the current time interval + #Timestamp is in the current time interval if binStart <= timeStamp and timeStamp < binEnd: - if CH_START <= channelID and channelID <= CH_END: + if CH_START <= channelID and channelID <= CH_END: if UNIT_START <= unitID and unitID <= UNIT_END: nrn_py_binSpk.dataHave += 1 # Update spike counts/ch (MUA) @@ -691,33 +691,33 @@ def nrn_py_binSpk(spkArr, binStart, timeoutFlag): nrn_py_binSpk.dataHave = 0 nrn_py_binSpk.binnedMsg = np.zeros((CH_END + 1, 1)) except: - print "[nrn_py_binSpk] exception occurs:", sys.exc_info()[0] - finally: + print("[nrn_py_binSpk] exception occurs:", sys.exc_info()[0]) + finally: return (newBinnedMsg, binningComplete) ## # The LWC server with the Lightweight client which sends raw data to the server # port: port number for UDP server -# binWnd: choose proper time window in ms. Ex) 100, 200, ... +# binWnd: choose proper time window in ms. Ex) 100, 200, ... # verbose: 1 prints information useful for debugging # Received data: <> def nrn_py_interfaceDpLwc(dir): - print "nrn_py_interfaceDpLwc running..." - + print("nrn_py_interfaceDpLwc running...") + # MUA per channel in a binWnd: [spikes in ch0, spikes in ch1, ..., spikes in ch96, timeInterval] mua = [0] * (CH_END + 1) buf = 4096 if verbose: # Remove previous data files if isUdp: - fname = "Udp" + fname = "Udp" else: - fname = "Tcp" + fname = "Tcp" for fileRemove in glob('*LWC*.tsv'): os.unlink(fileRemove) f1 = open(fname + "LWCTS.tsv", "w+") - f1.write('Type\tCH#\tUnit#\tTS\n') + f1.write('Type\tCH#\tUnit#\tTS\n') f2 = open(fname + "LWC2.tsv", "w+") f2.write('IntervalEnd\tSpikes in all CHs\tTotal spikes\n') f3 = open(fname + "LWCMua.tsv", "w+") @@ -729,13 +729,13 @@ def nrn_py_interfaceDpLwc(dir): spikePerTS = 0 # spike counts per timestamp totalRcvBytes = 0 a = time.time() - + # connect to Plexon client - Sock, conn, addr = nrn_py_connectPlxClient() + Sock, conn, addr = nrn_py_connectPlxClient() SN = 1 # serial number lastEndTS = 0.0 # second binStart = 0 - + if timeMeasure: mAddr0 = (ackServerIP, measurePort0) mAddr1 = (ackServerIP, measurePort1) @@ -744,32 +744,32 @@ def nrn_py_interfaceDpLwc(dir): totalMUA = 0 mSock = socket(AF_INET, SOCK_DGRAM) - mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) - mA = array.array('H') - mA.extend([0]) + mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) + mA = array.array('H') + mA.extend([0]) fTimeNum = 1 bTimeList = [] # Rcv messages from the client while 1: if verbose: - getTime1 = time.time() + getTime1 = time.time() if isUdp: - data, addr = Sock.recvfrom(buf) + data, addr = Sock.recvfrom(buf) else: data = conn.recv(buf) if verbose: f3.write("gTime1: " + str(time.time() - getTime1) + "\t") - dataHave = 0 - dataLen = len(data) - if data == 'exit': - print "Client has exited!" + dataHave = 0 + dataLen = len(data) + if data == 'exit': + print("Client has exited!") if verbose: - print "Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt + print("Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt) f1.write('Total spikes: ' + str(totalSpikeCnt) + '\tValid spikes: ' + str(validSpikeCnt) + '\n') - f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') + f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') f1.close() - f2.close() + f2.close() # close socket & generating exit and queue it Sock.close() sys.exit(0) @@ -782,38 +782,38 @@ def nrn_py_interfaceDpLwc(dir): cdataNp = np.zeros((0, 4)) timeoutFlag = 1 if verbose: - print "No data" + print("No data") else: if not isUdp: - while dataLen % 32 != 0: # For avoiding errors because of MSS in TCP + while dataLen % 32 != 0: # For avoiding errors because of MSS in TCP data += conn.recv(buf) - dataLen = len(data) - dlen = dataLen / 8 # double : 8 bytes + dataLen = len(data) + dlen = dataLen / 8 # double : 8 bytes drows = dlen / 4 # 4 fields : |data type|ch#|unit#|TS| cdata = struct.unpack('d' * dlen, data) cdataNp = np.asarray(cdata).reshape(drows, 4) # convert tuple to 2D numpy array - if timeMeasure: - mA[0] = O_RCV + if timeMeasure: + mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr3) #### for filtering - mA[0] = O_SND + mA[0] = O_SND mSock.sendto(mA.tostring(), mAddr1) #### - # filtering + # filtering if syncSpkFilter > 0 and isLwc == 1: cdataNp = nrn_py_filterSyncSpk(cdataNp) if unsortedSpkFilter == 1 and isLwc == 1: cdataNp = nrn_py_filterUnsortedSpk(cdataNp); - if timeMeasure: + if timeMeasure: #### for filtering - mA[0] = O_RCV + mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr1) #### # binning (Lwc) if isDp == 1 and isLwc == 1: # binning in the DP model with the HWC mode if timeMeasure: - mA[0] = O_SND + mA[0] = O_SND mSock.sendto(mA.tostring(), mAddr2) mua, binningComplete = nrn_py_binSpk(cdataNp, binStart, timeoutFlag) if binningComplete == 1: @@ -824,24 +824,24 @@ def nrn_py_interfaceDpLwc(dir): mua[CH_END] = (binStart + binWnd) / binWnd # binning window binStart = binStart + binWnd; if verbose: - print "Client sends a NODATA message!" + print("Client sends a NODATA message!") timeStamp = cdata[1] - totalRcvBytes += NODATA_SIZE - f1.write(str(timeStamp) + '\n') + totalRcvBytes += NODATA_SIZE + f1.write(str(timeStamp) + '\n') f2.write(str(timeStamp) + '\t' + str(spikePerTS) + '\t' + str(totalSpikeCnt) + '\t' + str(NODATA_SIZE) + '\n') qtime = time.time() - if timeMeasure: - mA[0] = O_RCV + if timeMeasure: + mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr2) - #mA[0] = O_END - #mSock.sendto(mA.tostring(), mAddr) + #mA[0] = O_END + #mSock.sendto(mA.tostring(), mAddr) timer = time.time() feedQueue(mua, verbose) # add 0 + timeStamp to the queue - if timeMeasure: # timing measure + if timeMeasure: # timing measure totalMUA +=1 - if totalMUA <= iteration: + if totalMUA <= iteration: qTimeList.append(time.time() - timer) if totalMUA == iteration + 1: if not os.path.exists(dir): @@ -850,29 +850,29 @@ def nrn_py_interfaceDpLwc(dir): f1 = open(fname, "w+") for item in qTimeList: f1.write("%s\n" % item) - f1.close() - print "====LwcDp push" + f1.close() + print("====LwcDp push") if verbose: - a = time.time() - qtime + a = time.time() - qtime if 0: - print "No data" - + print("No data") + ###################################################################### # For DP-HWC mode # The HWC server with the Heavyweight client which sends processed data to the server # Received data: <> def nrn_py_interfaceDpHwc(dir): - print "nrn_py_interfaceDpHwc running..." + print("nrn_py_interfaceDpHwc running...") # MUA per channel in a binWnd: [spikes in ch0, spikes in ch1, ..., spikes in ch96, timeInterval] mua = [0] * (CH_END + 1) buf = 4096 - - if verbose: + + if verbose: # Remove previous data files if isUdp: - fname = "Udp" - else: + fname = "Udp" + else: fname = "Tcp" for fileRemove in glob('*HWC*.tsv'): os.unlink(fileRemove) @@ -881,49 +881,49 @@ def nrn_py_interfaceDpHwc(dir): f2 = open(fname + "HWC2.tsv", "w+") f2.write('IntervalEnd\tSpikes in all CHs\tTotal spikes\n') f3 = open(fname + "HWCSua.tsv", "w+") - f3.write('IntervalEnd\tSUA/Ch\n') + f3.write('IntervalEnd\tSUA/Ch\n') totalSpikeCnt = 0 spikePerTS = 0 # spike counts per timestamp totalRcvBytes = 0 validSpikeCnt = 0 a = time.time() - + # connect to Plexon client - Sock, conn, addr = nrn_py_connectPlxClient() + Sock, conn, addr = nrn_py_connectPlxClient() SN = 1 # serial number lastEndTS = 0.0 # second - + # Rcv messages from the client if timeMeasure: - totalMUA = 0 + totalMUA = 0 mAddr0 = (ackServerIP, measurePort0) mAddr1 = (ackServerIP, measurePort1) mAddr2 = (ackServerIP, measurePort2) mAddr3 = (ackServerIP, measurePort3) mSock = socket(AF_INET, SOCK_DGRAM) - mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) - mA = array.array('H') - mA.extend([0]) - number = 1 + mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) + mA = array.array('H') + mA.extend([0]) + number = 1 while 1: if verbose: - getTime1 = time.time() + getTime1 = time.time() if isUdp: - data, addr = Sock.recvfrom(buf) + data, addr = Sock.recvfrom(buf) else: data = conn.recv(buf) if verbose: f3.write("gTime1: " + str(time.time() - getTime1) + "\t") - dataHave = 0 - dataLen = len(data) - if data == 'exit': - print "Client has exited!" + dataHave = 0 + dataLen = len(data) + if data == 'exit': + print("Client has exited!") if verbose: - print "Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt + print("Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt) f1.write('Total spikes: ' + str(totalSpikeCnt) + '\tValid spikes: ' + str(validSpikeCnt) + '\n') - f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') + f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') f1.close() - f2.close() + f2.close() # close socket & generating exit and queue it Sock.close() sys.exit(0) @@ -932,37 +932,37 @@ def nrn_py_interfaceDpHwc(dir): if verbose: getTime2 = time.time() if dataLen == 4: - #'H': unsigned short (2bytes) - cdata = struct.unpack('H' * 2, data) + #'H': unsigned short (2bytes) + cdata = struct.unpack('H' * 2, data) mua = [-1] * (CH_END + 1) mua[CH_END] = cdata[1] # binning window if verbose: - print "Client sends a NODATA message!" + print("Client sends a NODATA message!") timeStamp = cdata[1] - totalRcvBytes += NODATA_SIZE - f1.write(str(timeStamp) + '\n') + totalRcvBytes += NODATA_SIZE + f1.write(str(timeStamp) + '\n') f2.write(str(timeStamp) + '\t' + str(spikePerTS) + '\t' + str(totalSpikeCnt) + '\t' + str(NODATA_SIZE) + '\n') qtime = time.time() feedQueue(mua, verbose) # add 0 + timeStamp to the queue if verbose: - a = time.time() - qtime + a = time.time() - qtime if 0: - print "No data" + print("No data") else: if not isUdp: # DATA_SIZE = (CH_END + 1) * 8 # 1 for binning window. 8 is double in Matlab - while dataLen < DATA_SIZE: # For avoiding errors because of MSS in TCP + while dataLen < DATA_SIZE: # For avoiding errors because of MSS in TCP data += conn.recv(buf) - dataLen = len(data) - dlen = dataLen / 8 # double : 8 bytes # [ch1|ch2|...|ch96|binning window] + dataLen = len(data) + dlen = dataLen / 8 # double : 8 bytes # [ch1|ch2|...|ch96|binning window] mua = struct.unpack('d' * dlen, data) - if timeMeasure: - mA[0] = O_RCV - mSock.sendto(mA.tostring(), mAddr3) - timer = time.time() + if timeMeasure: + mA[0] = O_RCV + mSock.sendto(mA.tostring(), mAddr3) + timer = time.time() - feedQueue(mua, verbose) + feedQueue(mua, verbose) if timeMeasure: # timing measure totalMUA +=1 @@ -975,61 +975,61 @@ def nrn_py_interfaceDpHwc(dir): f1 = open(fname, "w+") for item in qTimeList: f1.write("%s\n" % item) - f1.close() - print "==========HwcDpPushTimeFile" + f1.close() + print("==========HwcDpPushTimeFile") if verbose: totalRcvBytes += DATA_SIZE - # for verification between client and server - timeStamp = cdata[dataOffset + 1 + CH_END * 6] * binWnd - f3.write(str(timeStamp)) + # for verification between client and server + timeStamp = cdata[dataOffset + 1 + CH_END * 6] * binWnd + f3.write(str(timeStamp)) # print MUA/ch only f1.write(str(timeStamp)) for i in range(CH_END): spikePerTS += mua[i] f1.write('\t' + str(mua[i])) # 6=> ch:unit0:unit1:unit2:unit3:unit4 for j in range(UNIT_END + 1): - f3.write('\t' + str(cdata[dataOffset + 1 + i * 6 + 1 + j])) - f1.write('\n') - f3.write('\n') + f3.write('\t' + str(cdata[dataOffset + 1 + i * 6 + 1 + j])) + f1.write('\n') + f3.write('\n') totalSpikeCnt += spikePerTS f2.write(str(timeStamp) + '\t' + str(spikePerTS) + '\t' + str(totalSpikeCnt) + '\t' + str(DATA_SIZE) + '\n') - spikePerTS = 0 + spikePerTS = 0 ######################################################################################### # The server with the Lightweight client which sends raw data to the server -# It does't do the binning process, and just queues a chunk of spikes (CHID, UNITID, TS). +# It does't do the binning process, and just queues a chunk of spikes (CHID, UNITID, TS). # Received data: <> def nrn_py_interfaceNsloc(dir): - print "[nrn_py_interfaceNsloc running...]" + print("[nrn_py_interfaceNsloc running...]") # Monkey spike: [Channel ID, Unit ID, timestamp] - spk = [0] * (SPKSZ + 1) # = 3 * SPKNUM + 1 + Serial Number + spk = [0] * (SPKSZ + 1) # = 3 * SPKNUM + 1 + Serial Number buf = 4096 - + if verbose: # Remove previous data files if isUdp: - fname = "Udp" + fname = "Udp" else: - fname = "Tcp" + fname = "Tcp" for fileRemove in glob('*NS*.tsv'): os.unlink(fileRemove) f1 = open(fname + "NSTS.tsv", "w+") - f1.write('Type\tCH#\tUnit#\tTS\n') + f1.write('Type\tCH#\tUnit#\tTS\n') f2 = open(fname + "NSValidTS.tsv", "w+") - f2.write('CH#\tUnit#\tTS\n') + f2.write('CH#\tUnit#\tTS\n') f3 = open(fname + "NSValidTS2.tsv", "w+") - f3.write('CH#\tUnit#\tTS\n') + f3.write('CH#\tUnit#\tTS\n') totalSpikeCnt = 0 validSpikeCnt = 0 a = time.time() - - + + # connect to Plexon client - Sock, conn, addr = nrn_py_connectPlxClient() + Sock, conn, addr = nrn_py_connectPlxClient() SN = 1 # serial number lastEndTS = 0.0 # second - + if timeMeasure: totalMUA = 0 mAddr0 = (ackServerIP, measurePort0) @@ -1037,31 +1037,31 @@ def nrn_py_interfaceNsloc(dir): mAddr2 = (ackServerIP, measurePort2) mAddr3 = (ackServerIP, measurePort3) mSock = socket(AF_INET, SOCK_DGRAM) - mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) - mA = array.array('H') - mA.extend([0]) - bTimeList = [] - number = 1 + mSock.setsockopt(SOL_SOCKET, SO_REUSEADDR, 1) + mA = array.array('H') + mA.extend([0]) + bTimeList = [] + number = 1 # Rcv messages from the client while 1: if verbose: - getTime1 = time.time() + getTime1 = time.time() if isUdp: - data, addr = Sock.recvfrom(buf) + data, addr = Sock.recvfrom(buf) else: data = conn.recv(buf) if verbose: f3.write("gTime1: " + str(time.time() - getTime1) + "\t") - dataHave = 0 - dataLen = len(data) - if data == 'exit': - print "Client has exited!" + dataHave = 0 + dataLen = len(data) + if data == 'exit': + print("Client has exited!") if verbose: - print "Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt + print("Total spikes: ", totalSpikeCnt, "Valid spikes: ", validSpikeCnt) f1.write('Total spikes: ' + str(totalSpikeCnt) + '\tValid spikes: ' + str(validSpikeCnt) + '\n') - f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') + f2.write('Valid spikes: ' + str(validSpikeCnt) + '\n') f1.close() - f2.close() + f2.close() # close socket & generating exit and queue it Sock.close() sys.exit(0) @@ -1074,24 +1074,24 @@ def nrn_py_interfaceNsloc(dir): cdataNp = np.zeros((0, 4)) timeoutFlag = 1 if 0: - print "No data" + print("No data") else: if not isUdp: - while dataLen % 32 != 0: # For avoiding errors because of MSS in TCP + while dataLen % 32 != 0: # For avoiding errors because of MSS in TCP data += conn.recv(buf) - dataLen = len(data) - dlen = dataLen / 8 # double : 8 bytes + dataLen = len(data) + dlen = dataLen / 8 # double : 8 bytes drows = dlen / 4 # 4 fields : |data type|ch#|unit#|TS| cdata = struct.unpack('d' * dlen, data) cdataNp = np.asarray(cdata).reshape(drows, 4) # convert tuple to 2D numpy array - if timeMeasure: - mA[0] = O_RCV + if timeMeasure: + mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr3) - #number += 1 + #number += 1 if isLwc == 1: if timeMeasure: #### ack_rcv1 for filtering - mA[0] = O_SND + mA[0] = O_SND mSock.sendto(mA.tostring(), mAddr1) if syncSpkFilter > 0: @@ -1099,22 +1099,22 @@ def nrn_py_interfaceNsloc(dir): if unsortedSpkFilter == 1: cdataNp = nrn_py_filterUnsortedSpk(cdataNp); - if timeMeasure: + if timeMeasure: #### ack_rcv1 for filtering mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr1) - - drows, col = cdataNp.shape - if timeMeasure: + + drows, col = cdataNp.shape + if timeMeasure: #mA[0] = O_FLT - #mSock.sendto(mA.tostring(), mAddr) - #### ack_rcv2 for buffering - mA[0] = O_SND + #mSock.sendto(mA.tostring(), mAddr) + #### ack_rcv2 for buffering + mA[0] = O_SND mSock.sendto(mA.tostring(), mAddr2) # "NODATA" message if drows == 0 and timeoutFlag == 1: - if verbose: - print "Client sends a NODATA message due to timeout!" + if verbose: + print("Client sends a NODATA message due to timeout!") if verbose: getTime2 = time.time() # Generating timeout queue item (Itimeout) @@ -1123,33 +1123,33 @@ def nrn_py_interfaceNsloc(dir): spk[SPKSZ - 1] = 0 # number of spikes spk[SPKSZ] = SN # serial number from 1 lastEndTS += TIMEOUT - SN += 1 + SN += 1 if verbose: f3.write("loopTime: " + str(time.time()- a) + "\t") qtime = time.time() feedQueue(spk, verbose) #queue.put(spk) if verbose: - f3.write("qTime: " + str(time.time()- qtime) + "\n") - a = time.time() - #f1.write("========\n") + f3.write("qTime: " + str(time.time()- qtime) + "\n") + a = time.time() + #f1.write("========\n") else: if verbose and drows > 0:#verbose: f_handle = file('spk_rcv.txt', 'a') np.savetxt(f_handle, cdataNp, fmt= '%4d %4d %4d %10.4f') f_handle.close() - #f3.write("gTime2: " + str(time.time() - getTime2) + "\n") - #a = time.time() # for measureing loop time + #f3.write("gTime2: " + str(time.time() - getTime2) + "\n") + #a = time.time() # for measureing loop time for j in range(drows): - if verbose: + if verbose: totalSpikeCnt += 1 - #print int(cdataNp[j][1]), int(cdataNp[j][2]), cdata[j][3] - #f1.write(str(int(cdataNp[j][0])) + '\t' + str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') - f1.write(str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') - if cdataNp[j][0] == 1: # neuron types + #print int(cdataNp[j][1]), int(cdataNp[j][2]), cdata[j][3] + #f1.write(str(int(cdataNp[j][0])) + '\t' + str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') + f1.write(str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') + if cdataNp[j][0] == 1: # neuron types channelID = int(cdataNp[j][1]) unitID = int(cdataNp[j][2]) timeStamp = cdataNp[j][3] # (second) - if CH_START <= channelID and channelID <= CH_END: + if CH_START <= channelID and channelID <= CH_END: if UNIT_START <= unitID and unitID <= UNIT_END: spk[dataHave * 3] = channelID spk[dataHave * 3 + 1] = unitID @@ -1159,12 +1159,12 @@ def nrn_py_interfaceNsloc(dir): validSpikeCnt += 1 #print "ChID:", spk[0], "UnitID:", spk[1], "Time:", spk[2] f2.write(str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') #'\t' + str(dataHave) + '\t' + str(j) + '\t' + str(drows -1) + '\n') - f3.write(str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') #'\t' + str(dataHave) + '\t' + str(j) + '\t' + str(drows -1) + '\n') - if dataHave == SPKNUM or j == drows - 1: + f3.write(str(int(cdataNp[j][1])) + '\t' +str(int(cdataNp[j][2])) + '\t' +str('{:g}'.format(cdataNp[j][3])) + '\n') #'\t' + str(dataHave) + '\t' + str(j) + '\t' + str(drows -1) + '\n') + if dataHave == SPKNUM or j == drows - 1: spk[SPKSZ - 1] = dataHave lastEndTS = spk[(dataHave - 1) * 3 + 2] spk[SPKSZ] = SN # serial number from 1 - SN += 1 + SN += 1 if verbose: f3.write("loopTime: " + str(time.time()- a) + "\t") qtime = time.time() @@ -1174,18 +1174,18 @@ def nrn_py_interfaceNsloc(dir): mA[0] = O_RCV mSock.sendto(mA.tostring(), mAddr2) #### ack_rcv2 for buffering - mA[0] = O_SND + mA[0] = O_SND mSock.sendto(mA.tostring(), mAddr2) - #mA[0] = O_END - #mSock.sendto(mA.tostring(), mAddr) - - timer = time.time() - + #mA[0] = O_END + #mSock.sendto(mA.tostring(), mAddr) + + timer = time.time() + feedQueue(spk, verbose) #queue.put(spk) - dataHave = 0 + dataHave = 0 if timeMeasure: # timing measure - totalMUA +=1 + totalMUA +=1 if totalMUA <= iteration: qTimeList.append(time.time() - timer) if totalMUA == iteration + 1: @@ -1195,21 +1195,21 @@ def nrn_py_interfaceNsloc(dir): f1 = open(fname, "w+") for item in qTimeList: f1.write("%s\n" % item) - f1.close() - print "Nsloc push" - btime = time.time() + f1.close() + print("Nsloc push") + btime = time.time() if verbose: - f3.write("qTime: " + str(time.time()- qtime) + "\n") - a = time.time() - #f1.write("========\n") - + f3.write("qTime: " + str(time.time()- qtime) + "\n") + a = time.time() + #f1.write("========\n") + # For benchmark of no communication case with DP or NSLOC model def serverNoComm(dir): # MUA per channel in a binWnd: [spikes in ch0, spikes in ch1, ..., spikes in ch96, timeInterval] - print "[serverNoComm is running...]" + print("[serverNoComm is running...]") try: - if isDp == 1: # inputs through DP cells + if isDp == 1: # inputs through DP cells mua = [0] * (CH_END + 1) fname = dir + "/StoredMua.tsv" fd = open(fname, 'r') @@ -1222,9 +1222,9 @@ def serverNoComm(dir): feedQueue(mua, verbose) time.sleep(0.1) # emulating 100 ms delay #print "%s" % (line) - else: # inputs through NSLOC units + else: # inputs through NSLOC units # Monkey spike: [Channel ID, Unit ID, timestamp] - spk = [0] * (SPKSZ + 1) # = 3 * SPKNUM + 1 + Serial Number + spk = [0] * (SPKSZ + 1) # = 3 * SPKNUM + 1 + Serial Number dataHave = 0 SN = 1 fname = dir + "/spikePMd-6sec.tsv" @@ -1240,27 +1240,27 @@ def serverNoComm(dir): # continue #else: if True: - if CH_START <= channelID and channelID <= CH_END: + if CH_START <= channelID and channelID <= CH_END: if UNIT_START <= unitID and unitID <= UNIT_END: spk[dataHave * 3] = channelID spk[dataHave * 3 + 1] = unitID spk[dataHave * 3 + 2] = timeStamp dataHave += 1 - if dataHave == SPKNUM: + if dataHave == SPKNUM: spk[SPKSZ - 1] = dataHave spk[SPKSZ] = SN # serial number from 1 - SN += 1 + SN += 1 feedQueue(spk, verbose) #queue.put(spk) - dataHave = 0 + dataHave = 0 if dataHave > 0: spk[SPKSZ - 1] = dataHave spk[SPKSZ] = SN # serial number from 1 feedQueue(spk, verbose) #queue.put(spk) - dataHave = 0 + dataHave = 0 except: - print "[serverNoComm] exception occurs:", sys.exc_info()[0] + print("[serverNoComm] exception occurs:", sys.exc_info()[0]) finally: - print "[ServerNoComm is terminated!!!]" + print("[ServerNoComm is terminated!!!]") fd.close() sys.exit(0) @@ -1283,10 +1283,10 @@ def start(self): self.workers.append(Process(target=serverNoComm, args = (self.dir, ))) for w in self.workers: - w.start() - + w.start() + def stop(self): for w in self.workers: - w.terminate() - print "[Server process is terminated!!!]" + w.terminate() + print("[Server process is terminated!!!]") diff --git a/shared.py b/shared.py index 4fcac1f..8597793 100644 --- a/shared.py +++ b/shared.py @@ -1,5 +1,5 @@ """ -globals.py +globals.py list of model objects (paramaters and variables) to be shared across modules Can modified manually or via arguments from main.py @@ -19,7 +19,7 @@ from time import time from math import radians import hashlib -def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xffffffff # for random seeds (bitwise AND to retain only lower 32 bits) +def id32(obj): return int(hashlib.md5(obj.encode("utf-8")).hexdigest()[0:8],16)# hash(obj) & 0xffffffff # for random seeds (bitwise AND to retain only lower 32 bits) ## MPI @@ -27,7 +27,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf nhosts = int(pc.nhost()) # Find number of hosts rank = int(pc.id()) # rank 0 will be the master -if rank==0: +if rank==0: pc.gid_clear() print('\nSetting parameters...') @@ -41,7 +41,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf popnames = ['PMd', 'ASC', 'EDSC', 'IDSC', 'ER2', 'IF2', 'IL2', 'ER5', 'EB5', 'IF5', 'IL5', 'ER6', 'IF6', 'IL6'] popclasses = [-1, -1, 1, 4, 1, 5, 4, 1, 1, 5, 4, 1, 5, 4] # Izhikevich population type popEorI = [ 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 1] # Whether it's excitatory or inhibitory -popratios = [96, 64, 64, 64, 150, 25, 25, 168, 72, 40, 40, 192, 32, 32] # Cell population numbers +popratios = [96, 64, 64, 64, 150, 25, 25, 168, 72, 40, 40, 192, 32, 32] # Cell population numbers popyfrac = [[-1,-1], [-1,-1], [-1,-1], [-1,-1], [0.1,0.31], [0.1,0.31], [0.1,0.31], [0.31,0.52], [0.52,0.77], [0.31,0.77], [0.31,0.77], [0.77,1.0], [0.77,1.0], [0.77,1.0]] # data from Weiler et. al 2008 (updated from Ben's excel sheet) receptornames = ['AMPA', 'NMDA', 'GABAA', 'GABAB', 'opsin'] # Names of the different receptors @@ -52,17 +52,17 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf popnumbers = [] PMdinput = 'spikes' # 'Plexon', 'spikes', 'SSM', 'targetSplit' - -# Define params for each cell: cellpops, cellnames, cellclasses, EorI + +# Define params for each cell: cellpops, cellnames, cellclasses, EorI cellpops = [] # Store list of populations for each cell -- e.g. 1=ER2 vs. 2=IF2 cellnames = [] # Store list of names for each cell -- e.g. 'ER2' vs. 'IF2' cellclasses = [] # Store list of classes types for each cell -- e.g. pyramidal vs. interneuron EorI = [] # Store list of excitatory/inhibitory for each cell popnumbers = scale*array(popratios) # Number of neurons in each population -if PMdinput == 'Plexon' and 'PMd' in popnames: +if PMdinput == 'Plexon' and 'PMd' in popnames: popratios[popnames.index('PMd')] = server.numPMd popnumbers[popnames.index('PMd')] = server.numPMd # Number of PMds is fixed. -ncells = int(sum(popnumbers))# Calculate the total number of cells +ncells = int(sum(popnumbers))# Calculate the total number of cells for c in range(len(popnames)): start = sum(popnumbers[:c]) popGidStart.append(start) @@ -73,23 +73,23 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf cellnames.append(popnames[pop]) # Append this cell's population type to the list cellclasses.append(popclasses[pop]) # Just use integers to index them EorI.append(popEorI[pop]) # Append this cell's excitatory/inhibitory state to the list -cellpops = array(cellpops) +cellpops = array(cellpops) cellnames = array(cellnames) EorI = array(EorI) # Assign numbers to each of the different variables so they can be used in the other functions # initializes the following variables: -# ASC, EDSC, ER2, IF2, IL2, ER5, EB5, IF5, IL5, ER6, IF6, IL6, AMPA, NMDA, GABAA, GABAB, opsin, Epops, Ipops, allpops +# ASC, EDSC, ER2, IF2, IL2, ER5, EB5, IF5, IL5, ER6, IF6, IL6, AMPA, NMDA, GABAA, GABAB, opsin, Epops, Ipops, allpops for i,name in enumerate(popnames): exec(name+'=i') # Set population names to the integers for i,name in enumerate(receptornames): exec(name+'=i') # Set population names to the integers -allpops = array(range(npops)) # Create an array with all the population numbers +allpops = array(list(range(npops))) # Create an array with all the population numbers Epops = allpops[array(popEorI)==0] # Pick out numbers corresponding to excitatory populations Ipops = allpops[array(popEorI)==1] # Pick out numbers corresponding to inhibitory populations # for creating natural and artificial stimuli (need to import after init of population and receptor indices) -from stimuli import touch, stimmod, makestim +from stimuli import touch, stimmod, makestim PMdconnprob = 2.0 @@ -105,14 +105,14 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf connprobs[ER2,IL2]=0.51 connprobs[ER2,IF2]=0.43 connprobs[EB5,ER2]=0 # none by wiring matrix in (Weiler et al., 2008) -connprobs[EB5,EB5]=0.04 * 4 # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 -connprobs[EB5,ER5]=0 # set using (Kiritani et al., 2012) Fig. 6D, Table 1 +connprobs[EB5,EB5]=0.04 * 4 # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 +connprobs[EB5,ER5]=0 # set using (Kiritani et al., 2012) Fig. 6D, Table 1 connprobs[EB5,ER6]=0 # none by suggestion of Ben and Gordon over phone connprobs[EB5,IL5]=0 # ruled out by (Apicella et al., 2012) Fig. 7 connprobs[EB5,IF5]=0.43 # CK: was 0.43 but perhaps the cause of the blow-up? connprobs[ER5,ER2]=0.2 # weak by wiring matrix in (Weiler et al., 2008) connprobs[ER5,EB5]=0.21 * 4 # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 -connprobs[ER5,ER5]=0.11 * 4 # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 +connprobs[ER5,ER5]=0.11 * 4 # set using (Kiritani et al., 2012) Fig. 6D, Table 1, value x 4 connprobs[ER5,ER6]=0.2 # weak by wiring matrix in (Weiler et al., 2008) connprobs[ER5,IL5]=0 # ruled out by (Apicella et al., 2012) Fig. 7 connprobs[ER5,IF5]=0.43 @@ -124,7 +124,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf connprobs[ER6,IF6]=0.43 connprobs[IL2,ER2]=0.35 connprobs[IL2,IL2]=0.09 -connprobs[IL2,IF2]=0.53 +connprobs[IL2,IF2]=0.53 connprobs[IF2,ER2]=0.44 connprobs[IF2,IL2]=0.34 connprobs[IF2,IF2]=0.62 @@ -142,9 +142,9 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf connprobs[IF6,ER6]=0.44 connprobs[IF6,IL6]=0.34 connprobs[IF6,IF6]=0.62 -connprobs[ASC,ER2]=0.6 +connprobs[ASC,ER2]=0.6 connprobs[EB5,EDSC]=1.0 #0.6 -connprobs[EB5,IDSC]=0.0 # hard-wire so receives same input as EB5->EDSC +connprobs[EB5,IDSC]=0.0 # hard-wire so receives same input as EB5->EDSC connprobs[IDSC,EDSC]=0.0 # hard-wire so projects to antagonist muscle subpopulation connprobs[PMd,ER5]=PMdconnprob @@ -199,7 +199,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf connweights[ASC,ER2,AMPA]=1.0 connweights[EB5,EDSC,AMPA]=2.0 connweights[EB5,IDSC,AMPA]=0.5 -connweights[IDSC,EDSC,GABAA]=1.0 +connweights[IDSC,EDSC,GABAA]=1.0 connweights[PMd,ER5,AMPA]=PMdconnweight @@ -231,7 +231,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf verbose = 0 # Whether to write nothing (0), diagnostic information on events (1), or everything (2) a file directly from izhi.mod filename = 'm1ms' # Set file output name plotraster = False # Whether or not to plot a raster -plotpeth = False # plot perievent time histogram +plotpeth = False # plot perievent time histogram plotpsd = False # plot power spectral density maxspikestoplot = 3e8 # Maximum number of spikes to plot plotconn = False # whether to plot conn matrix @@ -273,7 +273,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf stdpwin = 10 # length of stdp window (ms) (scholarpedia=10; Frem13=20(+),40(-)) eligwin = 50 # length of RL eligibility window (ms) (Frem13=500ms) useRLexp = 0 # Use binary or exp decaying eligibility trace -useRLsoft = 1 # Use soft thresholding +useRLsoft = 1 # Use soft thresholding maxweight = 8 # Maximum synaptic weight timebetweensaves = 5*1e3 # How many ms between saving weights(can't be smaller than loopstep) timeoflastsave = -inf # Never saved @@ -295,9 +295,9 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf useArm = 'dummyArm' # what type of arm to use: 'randomOutput', 'dummyArm' (simple python arm), 'musculoskeletal' (C++ full arm model) animArm = False # shows arm animation graphsArm = False # shows graphs (arm trajectory etc) when finisheds -targetid = 1 # initial target +targetid = 1 # initial target minRLerror = 0.002 # minimum error change for RL (m) -armLen = [0.4634 - 0.173, 0.7169 - 0.4634] # elbow - shoulder from MSM;radioulnar - elbow from MSM; +armLen = [0.4634 - 0.173, 0.7169 - 0.4634] # elbow - shoulder from MSM;radioulnar - elbow from MSM; nMuscles = 4 # number of muscles startAng = [0.62,1.53] # starting shoulder and elbow angles (rad) = natural rest position targetDist = 0.15 # target distance from center (15 cm) @@ -305,16 +305,16 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf initArmMovement = 250 # time after which to start moving arm (adds initial delay to avoid using initial burst of activity due to background noise init) motorCmdStartCell = popGidStart[EDSC] # start cell for motor command motorCmdEndCell = popGidStart[EDSC] + popnumbers[EDSC] # end cell for motor command -if useArm == 'dummyArm': +if useArm == 'dummyArm': cmdmaxrate = scale*20.0 # maximum spikes for motor command (normalizing value) cmdmaxrateTest = scale*20.0 # maximum spikes for motor command (normalizing value) -elif useArm == 'musculoskeletal': +elif useArm == 'musculoskeletal': cmdmaxrate = 8*scale*20.0 # maximum spikes for motor command (normalizing value 0 to 1) cmdmaxrateTest = 8*scale*20.0 # maximum spikes for motor command (normalizing value) cmdtimewin = 100 # spike time window for motor command (ms) antagInh = 0 # antagonist muscle inhibition # proprioceptive encoding -pStart = popGidStart[ASC] +pStart = popGidStart[ASC] numPcells = popnumbers[ASC] # number of proprioceptive (P) cells to encode shoulder and elbow angles minPval = radians(-30) # min angle to encode maxPval = radians(135) # max angle to encode @@ -334,7 +334,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf timeoflastreset = 0 # time when arm was last reseted (start from center etc. to simulate trials) -## PMd inputs +## PMd inputs if PMdinput == 'Plexon': vec = h.Vector() # temporary Neuron vectors emptyVec = h.Vector() @@ -345,7 +345,7 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf targetPMdInputs = [] # PMd units active for each trial maxPMdRate = 100.0 minPMdRate = 0.01 - PMdNoiseRatio = 0.1 + PMdNoiseRatio = 0.1 if PMdinput == 'spikes': pmdnc = [] pmdncDic = {} # dictionary to relate global and local PMd netcons @@ -373,7 +373,4 @@ def id32(obj): return int(hashlib.md5(obj).hexdigest()[0:8],16)# hash(obj) & 0xf benchstart = time() for i in range(int(1.36e6)): tmp=0 # Number selected to take 0.1 s on my machine performance = 1/(10*(time() - benchstart))*100 - print(' Running at %0.0f%% default speed (%0.0f%% total)' % (performance, performance*nhosts)) - - - + print((' Running at %0.0f%% default speed (%0.0f%% total)' % (performance, performance*nhosts))) diff --git a/stimuli.py b/stimuli.py index b31aece..ad60d01 100644 --- a/stimuli.py +++ b/stimuli.py @@ -50,7 +50,7 @@ class opto: # Allow defaults to be overriden by name, e.g. instead of touch(), stimmod(touch,sta=1,fin=3) def stimmod(stimclass, name=None, receptor=None, isi=None, var=None, width=None, weight=None, sta=None, fin=None, shape=None, pops=None, fraction=None, rate=None, noise=None, loc=None, falloff=None): - stiminstance = stimclass() + stiminstance = stimclass() stimproperties = ['name', 'receptor', 'isi', 'var', 'width', 'weight', 'sta', 'fin', 'shape', 'pops', 'fraction', 'rate', 'noise', 'loc', 'falloff']; for prop in stimproperties: thisproperty = eval(prop) @@ -62,7 +62,7 @@ def stimmod(stimclass, name=None, receptor=None, isi=None, var=None, width=None, ## Define stimulus-making code def makestim(isi=1, variation=0, width=0.05, weight=10, start=0, finish=1, stimshape='gaussian'): from pylab import r_, convolve, shape - + # Create event times timeres = 0.005 # Time resolution = 5 ms = 200 Hz pulselength = 10 # Length of pulse in units of width @@ -73,19 +73,19 @@ def makestim(isi=1, variation=0, width=0.05, weight=10, start=0, finish=1, stims while currenttime=0 and currenttime Date: Thu, 7 Jul 2022 13:45:33 +0200 Subject: [PATCH 2/3] update model for python3 and modeldb-ci support --- network.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/network.py b/network.py index 9eaa711..1e47fe1 100644 --- a/network.py +++ b/network.py @@ -36,7 +36,6 @@ import analysis from arm import Arm # Class with arm methods and variables -import numpy as np ############################################################################### ### Sequences of commands to run full model @@ -353,7 +352,7 @@ def createNetwork(): nCells = s.motorCmdEndCell - s.motorCmdStartCell s.motorCmdCellRange = [] for i in range(s.nMuscles): - s.motorCmdCellRange.append(list(np.arange(s.motorCmdStartCell + (nCells/s.nMuscles)*i, s.motorCmdStartCell + (nCells/s.nMuscles)*i + (nCells/s.nMuscles),dtype=int))) # cells used to for shoulder motor command + s.motorCmdCellRange.append(list(range(s.motorCmdStartCell + int(nCells/s.nMuscles)*i, s.motorCmdStartCell + int(nCells/s.nMuscles)*i + int(nCells/s.nMuscles)))) # cells used to for shoulder motor command ## Calculate distances and probabilities @@ -867,7 +866,7 @@ def finalizeSim(): hostspikecells=array([]) hostspiketimes=array([]) for c in range(len(s.hostspikevecs)): # fails when saving raw - thesespikes = array(s.hostspikevecs[c]) # Convert spike times to an array + thesespikes = array([s.hostspikevecs[c].x[i] for i in range(s.hostspikevecs[c].size())]) # Convert spike times to an array nthesespikes = len(thesespikes) # Find out how many of spikes there were for this cell hostspiketimes = concatenate((hostspiketimes, thesespikes)) # Add spikes from this cell to the list #hostspikecells = concatenate((hostspikecells, (c+host*s.cellsperhost)*ones(nthesespikes))) # Add this cell's ID to the list From c7eeafb3d3d1b1b67681ba3b4385cf869e18929e Mon Sep 17 00:00:00 2001 From: Christos Kotsalos Date: Thu, 7 Jul 2022 14:11:15 +0200 Subject: [PATCH 3/3] update model for python3 and modeldb-ci support --- main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main.py b/main.py index c4d6d5f..80930f3 100644 --- a/main.py +++ b/main.py @@ -1,5 +1,5 @@ import matplotlib -# matplotlib.use('Agg') # needed for hpc batch sims +matplotlib.use('Agg') # needed for hpc batch sims import sys from numpy import mean, zeros