# analysis.py -- Python script for running the simulation interactively
#
# Usage:
#   ipython --pylab -i sim.py
# 
# Last update: 13/09/10 (salvadord)

#
# Global parameters
#

# Load the sys and os packages.
import sys, os
from neuron import h # Load the interpreter object.
from numpy import *
from pylab import *
#from pylab import figure, array, size, zeros, subplot, ion, show, pause, hold, sum
from os import system,listdir
from copy import deepcopy
import scipy.io
import numpy as np
import pickle
from bicolormap import bicolormap
import time
import analyse_funcs
    
# for ipython debugging
#from IPython.core.debugger import Tracer; debug_here = Tracer()


# show graphs with batch simulations error results 
def batchErrors(simdatadir, param1Arg=None):
    # Set up the simulation data directory.
    #simdatadir = 'data/13sep06_sim1'

    invert_axis = 1;
    
    # seed values
    wseedvals =[120456, 398115]#, 534031, 796321, 895199]
    iseedvals = [1235, 2837]#, 3955, 4506, 6789]
    
    # parameter values
    #arange(10,110,10)# #[0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15]#[10,20,30,40,50,60,70,80,90,100]#[50, 100,250,500,750,1000]#arange(20,180,20)#[25, 50,75, 100, 125, 150, 175, 200, 225, 250]#
    if param1Arg==None:
        param1_range =arange(0.8,1.24,0.04)#arange(1,9)# [0,1]#arange(50,550,50)
    else:
        param1_range=param1Arg
    param2_range = [0,1,2,3,4,5,6,7]
    
    # number of errror values = 6 = (ang vs xy) x (10%, 50%, 90%)
    err_vals = 6;
    
    # create array for all results (2 errors, 10 train times, 5 wseeds, 5 inseed, 4 error values)
    error_all = zeros((len(param1_range),len(param2_range),len(wseedvals),len(iseedvals),err_vals))
    
    # create array for results avg'd over seeds (2 errors, 10 train times, 4 error values)
    error_avg = zeros((len(param1_range),len(param2_range),err_vals))
    error_std = zeros((len(param1_range),len(param2_range),err_vals))
    error_p1_avg = zeros((len(param1_range),err_vals))
    
    # Loop over param1 values
    iparam1 = -1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        iparam2 = -1
        skipValue = 0
        # Loop over param2 values
        for param2 in param2_range:
            iparam2 = iparam2 + 1
            iwseed = -1
            # Loop over wiring seed...
            for wseed in wseedvals:
                iwseed = iwseed + 1
                iiseed = -1
                # Loop over input seed...
                for iseed in iseedvals:
                    iiseed = iiseed + 1
                # set filename
                outfilestem = '%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                if os.path.isfile(outfilestem):
                    print "loading "+str(outfilestem)
                else:
                    outfilestem = '%s/p1-%.2f_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                    if os.path.isfile(outfilestem):
                        print "loading "+str(outfilestem)
                    else: 
                        outfilestem = '%s/p1-%.3f_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                        if os.path.isfile(outfilestem):
                            print "loading "+str(outfilestem)
                        else:
                            skipValue=1

                # get errors from nqs
                if skipValue==0:
                    tmp = h.calcErrorFromNQ(outfilestem).to_python()
                    tmp =array(tmp)
                    error_all[iparam1, iparam2, iwseed, iiseed, :] = tmp
                #print tmp
            # calculate avg error over seeds
            if skipValue==0:
                for i in arange(err_vals):
                    error_avg[iparam1, iparam2, i] = mean(error_all[iparam1, iparam2, :, :, i])
                    error_std[iparam1, iparam2, i]= std(error_all[iparam1, iparam2, :, :, i])
    
        # calculate average over all targets (param2) for each value of param1
        if skipValue==0:
            for i in arange(err_vals):
                error_p1_avg[iparam1, i] = mean(error_avg[iparam1,:,i])
                
    # plot results    
    figsize=[800,500]
    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    ax1 = fig.add_subplot(211)
    ax2 = fig.add_subplot(212)
    param_xaxis = 'param2'
    #debug_here()
    
    if invert_axis:
        param_xaxis='param1'
        tmp=param1_range
        param1_range = param2_range
        param2_range = tmp
        error_avg = swapaxes(error_avg, 0, 1)
        
    xaxis=array(param2_range)
    colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
    
    # plot cartesian error
    iparam1=-1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        #ax1.errorbar(xaxis, error_avg[iparam1,:,0], yerr=error_std[0,:,0], fmt='s', color = colorlist[iparam1%10],)
        #ax1.errorbar(xaxis, error_avg[iparam1,:,1], yerr=error_std[0,:,1], fmt='o', color = colorlist[iparam1%10], )
        #ax1.errorbar(xaxis, error_avg[iparam1,:,2], yerr=error_std[0,:,2], fmt='x', color = colorlist[iparam1%10], )
        #ax1.plot(xaxis, error_avg[iparam1,:,0], color = colorlist[iparam1%10], label=str(param1)+', last 10%')
        ax1.plot(xaxis, error_avg[iparam1,:,0], linestyle = "--",color = colorlist[iparam1%10],label=str(param1)+', last 20%')
        ax1.plot(xaxis, error_avg[iparam1,:,2], linestyle = "-", color = colorlist[iparam1%10], label=str(param1)+', last 90%')
    ax1.plot(xaxis, error_p1_avg[:,0], linestyle = "--",color = 'black',linewidth=4, label='AVG - 20%')
    ax1.plot(xaxis, error_p1_avg[:,2], linestyle = "-", color = 'black', linewidth=4, label='AVG - 90%')
    
    ax1.set_ylabel('cartesian error (m)')
    ax1.set_xlabel(param_xaxis)
    ax1.set_title('Cartesian error as a func of param1 and param2')
    #ax1.legend(loc='upper center', bbox_to_anchor=(1, 0.2),  borderaxespad=0., prop={'size':12})
    ax1.grid(True)
    
    # plot angular error
    iparam1=-1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        #ax2.errorbar(xaxis, error_avg[iparam1,:,3], yerr=error_std[0,:,3], fmt='s', color = colorlist[iparam1%10],)
        #ax2.errorbar(xaxis, error_avg[iparam1,:,4], yerr=error_std[0,:,4], fmt='o', color = colorlist[iparam1%10],)
        #ax2.errorbar(xaxis, error_avg[iparam1,:,5], yerr=error_std[0,:,5], fmt='x', color = colorlist[iparam1%10], )
        #ax2.plot(xaxis, error_avg[iparam1,:,3], color = colorlist[iparam1%10], label=str(param1)+', last 10%')
        ax2.plot(xaxis, error_avg[iparam1,:,3], linestyle = "--", color = colorlist[iparam1%10], label=str(param1)+', last 20%')
        ax2.plot(xaxis, error_avg[iparam1,:,5], linestyle = "-",color = colorlist[iparam1%10], label=str(param1)+', last 90%')
    ax2.plot(xaxis, error_p1_avg[:,3], linestyle = "--",color = 'black',linewidth=4, label='AVG - last 20%')
    ax2.plot(xaxis, error_p1_avg[:,5], linestyle = "-", color = 'black', linewidth=4, label=' AVG - last 90%')
        
    ax2.set_ylabel('angular error (rad)')
    ax2.set_xlabel(param_xaxis)
    ax2.set_title('Angular error as a func of param1 and param2')
    ax2.legend(loc='upper right', bbox_to_anchor=(1.1, 2.0),  borderaxespad=0., prop={'size':12})
    ax2.grid(True)

    fig.tight_layout()    
    show()

# show the receptive field of a neuron or group of neurons (targetCell); 
# nSyn = number of synaptic connections to take into account; eg. if 1: ES->EM; if 2, DP->ES->EM
# animate = include synaptic changes over time; make graph

# show error graphs for 8 sims of same day
def batchErrors8(simdatadirRoot):
    numSims = 8

    param1_range =  []#[None] * numSims
    param1_range.append([10,20,30,40,50])
    param1_range.append([10,20,30,40,50,60,70,80,90,100])
    param1_range.append(arange(0.8,1.24,0.04))#([0.012, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15])
    param1_range.append([10,20,30,40,50,60,70,80,90,100])
    param1_range.append([25, 50,75, 100, 125, 150, 175, 200, 225, 250])
    param1_range.append([20,40,60,80,100,120,140,160])
    param1_range.append([200,225,250,275,300,325,350])
    param1_range.append([50, 100,250,500,750,1000])
    
    for i in range(numSims):
        batchErrors(simdatadirRoot+'_sim'+str(i+1), param1_range[i])
    
# show the RF of targetCells
# if nSyn2: use polysynaptic chains of length =2
# if animate: show RF change over time, overimposed with spikes
def receptiveFieldAndSpikes(targetCells, nSyn2, animate, dosave):    
    #dosave=1<
    if dosave:
        moviedir='data/movies/'
        moviename='13oct15b.mpg'
        maxmovieframes = 500 # Maximum number of movie frames
    view3d=0
    showSpikes=1
    
    print "Calculating receptive field of cell group starting at id %d " % targetCells[0]
    targetCells = array(targetCells)
    popnames=['P','ES','IS','ILS','EM','IM','IML'];
    popinds=[ 2, 41, 42, 43, 44, 45, 46];
    popshape=['o',  '^',  'h', '*',  '^',  'h',  '*', '8'];
    
    # convert connectivity to python arrays
    conpreid=array(h.col.connsnq.getcol("id1"), 'i')
    conpostid=array(h.col.connsnq.getcol("id2"), 'i')
    #condelay=array(h.col.connsnq.getcol("del"))
    #condistance=array(h.col.connsnq.getcol("dist"))
    conweight1=array(h.col.connsnq.getcol("wt1"))
    conweight2=array(h.col.connsnq.getcol("wt2"))
    
    order = lexsort((conpreid, conpostid)) # order by post and then pre
    conpreid = conpreid[order]
    conpostid = conpostid[order]
    #condelay = condelay[order]
    #condistance = condistance[order]
    conweight1 = conweight1[order]
    conweight2 = conweight2[order]
    
    maxW = max(conweight1.max(), conweight2.max())
    conweight1 = conweight1/maxW    # normalize so can multiply together
    conweight2 = conweight2/maxW # 
    
    conweight1Original = deepcopy(conweight1) # make copy of weights at t=0
    
    # convert locations to python arrays
    h('objref cellList')
    h('cellList=col.ce')
    h('objref cellListDP') 
    h('cellListDP = cedp')
    n=h.cellList.count() # Number of cells
    ndp=h.cellListDP.count() # number of DP cells
    #cellLocations = zeros((n+ndp,5)) # number of cells and attributes
    cellids =zeros(n+ndp)
    celltypes = zeros(n+ndp)
    xlocs =zeros(n+ndp)
    ylocs =zeros(n+ndp)
    zlocs =zeros(n+ndp)
    
    for i in range(int(n)): # Loop over each cell
        cellids[i]=i # Cell ID
        celltypes[i]=h.cellList.o(i).type() # Cell population
        xlocs[i]=h.cellList.o(i).xloc # X position
        ylocs[i]=h.cellList.o(i).yloc # Y position
        zlocs[i]=h.cellList.o(i).zloc # Z position

    for i in range(int(ndp)): # Loop over each cell
        cellids[n+i]=n+i # Cell ID
        celltypes[n+i]=h.cellListDP.o(i).type # Cell population
        xlocs[n+i]=h.cellListDP.o(i).xloc # X position
        ylocs[n+i]=h.cellListDP.o(i).yloc # Y position
        zlocs[n+i]=h.cellListDP.o(i).zloc # Z position
    
    # arrange neurons in spatial grid and according to muscle group
    spatialGrid = 1
    if spatialGrid:
        x = [0, 0, 0, 0, 0] # initial starting positions
        y = [0, 0, 0, 0, 0]        
        xstep = 1 # step increase
        ystep = 1.5
        xmax = 10
        for i in range(len(xlocs)):
            if celltypes[i] == popinds[1]:    #ES
                xmax = 24
                xoffset = [0, 0, 0, 0, 0]
                yoffset = [0, 0, 0, 0, 14]
                zlocs[i] = 0 # fix to same muscle group
            elif celltypes[i] == popinds[2]: #IS
                zlocs[i] = 4 # do not divide in muscle groups
            elif celltypes[i] == popinds[3]:    #ISL
                zlocs[i] = 4 # do not divide in muscle groups
            elif celltypes[i] == popinds[4]:    #EM
                xoffset = [0, 14, 0, 14, 5] #offset for different muscle populations + inhibitory cells
                yoffset = [0, 0, 14, 14, 7]
                xmax = 12
            elif celltypes[i] == popinds[5]: #IM
                zlocs[i] = 4 # do not divide in muscle groups
                xmax = 16
            elif celltypes[i] == popinds[6]:    #IML
                zlocs[i] = 4 # do not divide in muscle groups
                xmax = 16
            elif celltypes[i] == popinds[0]: #P
                xmax = 12
                xoffset = [0, 14, 0, 14, 3] #offset for different muscle populations + inhibitory cells
                yoffset = [0, 0, 12, 12, 9]    
            
            zloc = int(zlocs[i])
            x[zloc] = x[zloc] + xstep # increase step

            if x[zloc] > xmax: # if reached x limit
                y[zloc] = y[zloc] + ystep
                x[zloc] = 1
            
            xlocs[i] = x[zloc] + xoffset[zloc] # set positions
            ylocs[i] = y[zloc] + yoffset[zloc]
            
            if (i==(192+44+20-1) or (i==(2*(192+44+20)-1))): # when change subplot, restart xy positions
                x = [0, 0, 0, 0, 0] # initial starting positions
                y = [0, 0, 0, 0, 0]
                
    ########################
    # animation data and params
    if animate:
        # convert synaptic changes over time (since connsnq doesn't change with t, can add)
        h('objref synChanges') 
        #h('nqsy.tog()')
        h('nqsy.select("t",">=", 0)')
        h('synChanges = nqsy')
        synpreid=array(h.synChanges.getcol("id1"), 'i')
        synpostid=array(h.synChanges.getcol("id2"), 'i')
        synweight=array(h.synChanges.getcol("wg"))
        syntime=array(h.synChanges.getcol("t"))

        order = lexsort((synpreid, synpostid)) # order syn by post and then pre
        synpreid = synpreid[order]
        synpostid = synpostid[order]
        synweight=synweight[order]        
        syntime=syntime[order]    
        
        # convert spiking data to python arrays
        h('objref spikes')
        h('snq.select("t",">=", 0)')
        h('spikes = snq')
        spikeId=array(h.spikes.getcol("id"))
        spikeType=array(h.spikes.getcol("type"))
        spikeTime=array(h.spikes.getcol("t"))
        spikeMid=array(h.spikes.getcol("mid"))
        
        weightInc = h.plastEEinc / maxW # set weight increases and normalize
    
    def makeplot():
        # visualization options
        figsize = [720,800] # Figure size in pixels
        targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
        normalColor = 'white'#'lightyellow'
        targetSizeFactor = 2
        weightsCmap = 'YlOrRd' #'jet'#'hot'#'autumn'
        cellSize = 50;
        cellBorder = 1;
        
        # create subplots for M, S and P populations
        ion()
        fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
        fig.subplots_adjust(left=0.02) # Less space on left
        fig.subplots_adjust(right=0.93) # Less space on right
        fig.subplots_adjust(bottom=0.08) # Less space on bottom
        #fig.subplots_adjust(wspace=0.25) # More space between
        #fig.subplots_adjust(hspace=0.30) # More space between
        if view3d: 
            proj='3d'
        else:
            proj='rectilinear'
            
        rasterMsubplot = subplot(311, projection=proj) # create subplots
        rasterSsubplot = subplot(312, projection=proj)
        rasterPsubplot = subplot(313, projection=proj)
        
        rasterMsubplot.set_title('Motor population', fontsize=12) # titles
        rasterSsubplot.set_title('Somatosensory population', fontsize=12)
        rasterPsubplot.set_title('Proprioceptive population', fontsize=12)
        
        rasterMsubplot.text(-0.7, 1, "sho\next", fontsize=10, fontweight='bold') # muscle labels M
        rasterMsubplot.text(26.6, 1, "sho\nflex", fontsize=10, fontweight='bold')
        rasterMsubplot.text(-0.7, 15, "elb\next", fontsize=10, fontweight='bold')
        rasterMsubplot.text(26.6, 15, "elb\nflex", fontsize=10, fontweight='bold')
        
        rasterPsubplot.text(-0.7, 1, "sho\next", fontsize=10, fontweight='bold') # muscle labels P
        rasterPsubplot.text(26.6, 1, "sho\nflex", fontsize=10, fontweight='bold')
        rasterPsubplot.text(-0.7, 13, "elb\next", fontsize=10, fontweight='bold')
        rasterPsubplot.text(26.6, 13, "elb\nflex", fontsize=10, fontweight='bold')
         
        setp(rasterMsubplot.get_xticklabels(), visible=False) # hide x and y ticks
        setp(rasterSsubplot.get_xticklabels(), visible=False)
        setp(rasterPsubplot.get_xticklabels(), visible=False)
        setp(rasterMsubplot.get_yticklabels(), visible=False)
        setp(rasterSsubplot.get_yticklabels(), visible=False)
        setp(rasterPsubplot.get_yticklabels(), visible=False)
        
        border = 2
        rasterMsubplot.set_xlim([xlocs.min()-border, xlocs.max()+border]) # set x-y lims
        rasterMsubplot.set_ylim([ylocs.min()-border, ylocs.max()+border])
        rasterSsubplot.set_xlim([xlocs.min()-border, xlocs.max()+border])
        rasterSsubplot.set_ylim([ylocs.min()-border, ylocs.max()+border])
        rasterPsubplot.set_xlim([xlocs.min()-border, xlocs.max()+border])
        rasterPsubplot.set_ylim([ylocs.min()-border, ylocs.max()+border])
        
        ######################
        # plot cells in gray (static)
        cellsScatter = [None]*len(popinds);
        for p in range(len(popinds)):
            #cells=where(cellLocations[:,1]==popinds[p])
            cells=(celltypes==popinds[p])
            if view3d:
                if (p == 0):
                    rasterPsubplot.scatter(xlocs[cells], ylocs[cells], zlocs[cells], c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize) 
                elif (p>=1 and p<=3):
                    rasterSsubplot.scatter(xlocs[cells], ylocs[cells], zlocs[cells], c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize) 
                elif (p>=4 and p<=6):
                    rasterMsubplot.scatter(xlocs[cells], ylocs[cells], zlocs[cells], c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize) 
            else:
                if (p == 0):
                    cellsScatter[p] = rasterPsubplot.scatter(xlocs[cells], ylocs[cells],  c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize, label=popnames[p]) 
                elif (p>=1 and p<=3):
                    cellsScatter[p] = rasterSsubplot.scatter(xlocs[cells], ylocs[cells],  c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize, label=popnames[p]) 
                elif (p>=4 and p<=6):
                    cellsScatter[p] = rasterMsubplot.scatter(xlocs[cells], ylocs[cells],  c=normalColor, marker=popshape[p], linewidth=cellBorder, s=cellSize, label=popnames[p]) 
            
        #############################
        # show target cells in different color
        targetLabel = 'target cell'
        for i in range(len(targetCells)):
            for p in range(len(popinds)):
                # find target cells for each layer
                #targetLayerCells = intersect1d(targetCells, where(cellLocations[:,1]==popinds[p])[0])
                targetLayerCells = (cellids==targetCells[i]) * (celltypes==popinds[p]) # obtain indices for population p
                if targetLayerCells.any(): # it at least one cell was found
                    if (p == 0):
                        rasterPsubplot.scatter(xlocs[targetLayerCells], ylocs[targetLayerCells], c=targetColor, marker=popshape[7], linewidth=cellBorder*2, s=cellSize*targetSizeFactor) 
                    elif (p>=1 and p<=3):
                        rasterSsubplot.scatter(xlocs[targetLayerCells], ylocs[targetLayerCells], c=targetColor, marker=popshape[7],  linewidth=cellBorder*2, s=cellSize*targetSizeFactor) 
                    elif (p>=4 and p<=6):
                        if i == 0: targetPlot =rasterMsubplot.scatter(xlocs[targetLayerCells], ylocs[targetLayerCells], c=targetColor, marker=popshape[7], linewidth=cellBorder*2, s=cellSize*targetSizeFactor, label = targetLabel) 
                        else: rasterMsubplot.scatter(xlocs[targetLayerCells], ylocs[targetLayerCells], c=targetColor, marker=popshape[7], linewidth=cellBorder*2, s=cellSize*targetSizeFactor) 
                    
        legendLabels=['E', 'I', 'ILS', 'P', 'target']
        rasterSsubplot.legend([cellsScatter[1], cellsScatter[2],  cellsScatter[3], cellsScatter[0], targetPlot], legendLabels, columnspacing=0,markerscale=0.8, loc='upper center', prop={'size':9}, scatterpoints=1, bbox_to_anchor=(0.94, 1.05),  borderpad=0.07)#, borderpad=0.05, labelspacing=0.1)
        
        ######################
        # show afferent cells weights (RF)
        rf = [None]*len(popinds);
        
        def calculateRF(mode,rf):
            onlyAMPA = 1
            adaptScale = 0
            scaleMax =[0.5,1,1] 
            # find all cells projecting to the target cell for each layer
            afferentCells = array([], 'i')
            weightColors = []
            for i in range(len(targetCells)):
                targetAfferents = array(conpreid[where(conpostid==targetCells[i])], 'i') # get afferent of next target cell
                if targetAfferents.any(): # if any 
                    for iAfferent in targetAfferents: # for each new afferent 
                        if (iAfferent in afferentCells): # check if already included in list
                            if onlyAMPA: weightColors2 = (conweight1[iAfferent])
                            else: weightColors2 = ((conweight1[iAfferent]) +(conweight2[iAfferent]))
                            weightColors[afferentCells==iAfferent] +=  weightColors2
                        else: # if not included in list
                            afferentCells=append(afferentCells, iAfferent) # add to list of afferent cells
                            if onlyAMPA: weightColors = append(weightColors, conweight1[targetAfferents]); # calculate weigths
                            else: weightColors = append(weightColors, conweight1[afferentCells]-conweight2[afferentCells]); # calculate weigths
                    
                    if nSyn2: # 2 synaptic connections
                        for targetCell2 in targetAfferents:  # loop through all afferent cells (which now become target cells)
                            targetAfferents2 = array(conpreid[where(conpostid==targetCell2)], 'i') # get new set of afferent cells
                            if targetAfferents2.any():
                                for iAfferent in targetAfferents2: # loop through all afferent cells
                                    if (iAfferent in afferentCells): # if already included in list
                                        if onlyAMPA: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                        else: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                        weightColors[afferentCells==iAfferent] +=  weightColors2
                                    else:
                                        afferentCells=append(afferentCells, iAfferent) # add to list of afferent cells
                                        if onlyAMPA: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                        else: weightColors2 = (conweight1[targetCell2]+conweight2[targetCell2]) * (conweight1[iAfferent]+conweight2[iAfferent])
                                        weightColors=append(weightColors, weightColors2)
                                            
            weightColors = weightColors/weightColors.max() # normalize weights to 1
            #print afferentCells
                    
            for p in range(len(popinds)):
                # find all cells projecting to the target cell for each layer
                afferentCellsL = intersect1d(afferentCells, where(celltypes==popinds[p])[0]) # find indices for population p
                #print afferentCellsL
                if afferentCellsL.any():
                    weightColorsL = weightColors[in1d(afferentCells, afferentCellsL)]; # obtain weight colors for this population
                    #print weightColorsL
                    if mode == "new":
                        if (p == 0):
                            rf[p] = rasterPsubplot.scatter(xlocs[afferentCellsL], ylocs[afferentCellsL], cmap=weightsCmap, c=weightColorsL, vmin=0, vmax=scaleMax[0], marker=popshape[p], linewidth=cellBorder, s=cellSize) 
                            if adaptScale and ('c1' in locals() or 'c1' in globals()):
                                rf[p].set_clim([weightColorsL.min(), weightColorsL.max()])
                                c1.set_clim([weightColorsL.min(), weightColorsL.max()])
                                c1.draw_all()
                        elif (p>=1 and p<=3):
                            rf[p] = rasterSsubplot.scatter(xlocs[afferentCellsL], ylocs[afferentCellsL], cmap=weightsCmap,  c=weightColorsL, vmin=0, vmax=scaleMax[1], marker=popshape[p], linewidth=cellBorder, s=cellSize) 
                            if adaptScale and ('c2' in locals() or 'c2' in globals()):
                                rf[p].set_clim([weightColorsL.min(), weightColorsL.max()])
                                c2.set_clim([weightColorsL.min(), weightColorsL.max()])
                                c2.draw_all()
                        elif (p>=4 and p<=6):
                            rf[p] = rasterMsubplot.scatter(xlocs[afferentCellsL], ylocs[afferentCellsL], cmap=weightsCmap,  c=weightColorsL,vmin=0, vmax=scaleMax[0], marker=popshape[p], linewidth=cellBorder, s=cellSize) 
                            if adaptScale and ('c3' in locals() or 'c1' in globals()):
                                rf[p].set_clim([weightColorsL.min(), weightColorsL.max()])
                                c3.set_clim([weightColorsL.min(), weightColorsL.max()])
                                c3.draw_all()
                    elif mode == "update":
                            rf[p].set_color(weightColorsL)
                            rf[p].set_cmap(weightsCmap)
                            
                    #print afferentCells
            
        calculateRF("new",rf)
        tight_layout()
        
        cfraction = 0.09
        cpad = 0.01
        cshrink = 0.75
        cfontsize=9
        if rf[0] != None: 
            c1=colorbar(rf[0], ax=rasterPsubplot, fraction=cfraction, pad=cpad, shrink=cshrink)
            c1.ax.tick_params(labelsize=cfontsize) 
        c2=colorbar(rf[1], ax=rasterSsubplot, fraction=cfraction, pad=cpad, shrink=cshrink)
        c2.set_label('normalized weight to target cell', fontsize=cfontsize+1)
        c2.ax.tick_params(labelsize=cfontsize) 
        c3=colorbar(rf[4], ax=rasterMsubplot,  fraction=cfraction, pad=cpad, shrink=cshrink)
        c3.ax.tick_params(labelsize=cfontsize) 
        
        ########################################
        # ANIMATE SPIKES AND SYNAPTIC CHANGES
        if animate:
            # spike parameters
            raster = [None]*len(popinds);
            spikeColor = 'chartreuse'#'DarkMagenta'#'lime'#'purple'
            spikeSizeFactor = 1.5

            maxframes = 1000 # Maximum number of non-movie frames
            binsize = 50 # Bin size in ms
            ncells = len(xlocs)
            xmax = round(xlocs.max())
            ymax = round(xlocs.max())
            tmax = spikeTime.max()
            totalspikes = len(spikeTime)

            nbins = int(tmax/binsize)+1 # Calculate number of bins; +1 for edge effects
            spikecounts=zeros((ncells,nbins)) # Number of spike counts
            timebins = array(range(nbins))*binsize/1e3 # Calculate time bins
            spikebins = array(spikeTime/binsize,dtype=int) # Convert times to bins

            for s in range(totalspikes): # Loop over each spike
                spikecounts[spikeId[s],spikebins[s]] += 1 # Increment this bin
            spikingrate = sum(spikecounts,axis=0)*1e3/binsize/ncells # Get the total spike counts over time
            
            maxspikes = 5*sum(spikecounts,axis=0).mean() # Find out the maximum number of spikes in any given timestep
            maxiters = int(min(maxframes,maxmovieframes,len(timebins))) if dosave else int(min(maxframes,len(timebins))) # See how many iterations to go for

            # plot spikes and syn changes for every frame
            for b in range(maxiters): # for each bin
                
                rasterMsubplot.set_title('Motor population, t = %0.3f s' % timebins[b],  fontsize=12)
                rasterSsubplot.set_title('Somatosensory population, t = %0.3f s' % timebins[b], fontsize=12)
                rasterPsubplot.set_title('Proprioceptive population = %0.3f s' % timebins[b], fontsize=12)
                
                #############################
                # show spikes
                if showSpikes:
                    fired = array(spikecounts[:,b]>0).flatten() # Find which neurons fired in this timestep
                    print('  bin = %i of %i; %i spikes' % (b, maxiters, sum(fired)))
                    numfired = sum(fired) # Total number of cells that fired
                    
                    for p in range(len(popinds)):
                        firedL = fired * (celltypes==popinds[p]) # calculate cells that fired for population p
                        #print firedL
                        if (p == 0):
                            raster[p] = rasterPsubplot.scatter(xlocs[firedL], ylocs[firedL],  c=spikeColor, marker=popshape[p], linewidth=cellBorder-1, s=cellSize*spikeSizeFactor) 
                        elif (p>=1 and p<=3):
                            raster[p] = rasterSsubplot.scatter(xlocs[firedL], ylocs[firedL],  c=spikeColor, marker=popshape[p], linewidth=cellBorder-1, s=cellSize*spikeSizeFactor) 
                        elif (p>=4 and p<=6):
                            raster[p] = rasterMsubplot.scatter(xlocs[firedL], ylocs[firedL],  c=spikeColor, marker=popshape[p], linewidth=cellBorder-1, s=cellSize*spikeSizeFactor) 
                    
                    show()
                else:
                    print('  bin = %i of %i' % (b, maxiters))
                    
                ############################
                # show changes in synaptic weights
                synBin = where(syntime == b*binsize)[0]
                if (size(synBin)>0 and size(synBin)<size(conweight1)):
                    # update connection matrix based on syn wieght gains
                    conweight1[range(len(synBin))] = conweight1Original[range(len(synBin))] * synweight[synBin]
                    
                    for p in range(len(popinds)): 
                        if rf[p] != None: rf[p].remove() # remove previous scatter    
                    calculateRF("new",rf) # plot new RF

                if dosave: fig.savefig('tmpmovie%04i.png' % b) # Save PNG files for movie
                
                if dosave!=1: pause(1e-10) 
                if showSpikes: 
                    for p in range(len(popinds)): raster[p].remove() # remove previous scatter
                if dosave!=1: pause(1e-10)
    
    makeplot()

    if dosave:
        print('Making movie...')
        system('mkdir -p ' + moviedir)
        system("mencoder 'mf://tmpmovie*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" % (moviedir+moviename))
        print('  Cleaning up...')
        system('rm tmpmovie*.png')

# save potentiated/depressed synapses at each time step (for video)
def saveRLsyn(t1,t2):
    # convert synaptic changes over time (since connsnq doesn't change with t, can add)
        h('objref synChanges') 
        #h('nqsy.tog()')
        #h('nqsy.select("t",">", 0)')
        h('nqsy.select("t","[]", %d, %d)'%(t1,t2))
        h('synChanges = nqsy')
        synpreid=array(h.synChanges.getcol("id1"), 'i')
        synpostid=array(h.synChanges.getcol("id2"), 'i')
        synweight=array(h.synChanges.getcol("wg"))
        syntime=array(h.synChanges.getcol("t"))
        
        savetime=h.syDT
        prevWeight=synweight[syntime==t1]
        if t1<50:
            prevWeight=prevWeight[:len(prevWeight)/2]
            t1=2*savetime
        else: 
            t1=t1+savetime
        synChanged=[0,0,0,0]
        bins=0
        
        for t in arange(t1, syntime.max(), savetime):
            print t
            bins=bins+1
            # get data for this time
            currsynpreid = synpreid[syntime==t] 
            currsynpostid = synpostid[syntime==t]
            currsyntime = syntime[syntime==t]
            
            # find and save connections that have changed
            currWeight = synweight[syntime==t]
            diffWeight = currWeight-prevWeight
            changedWeight = where(diffWeight != 0)[0]
            # CHECK IF CAN CHANGE LINE BELOW TO AVOID FOR!
            #for i in range(len(changedWeight)):
                #synChanged = vstack((synChanged, [currsyntime[changedWeight[i]], currsynpreid[changedWeight[i]], currsynpostid[changedWeight[i]], diffWeight[changedWeight[i]]]))
            synChanged = vstack((synChanged, array([currsyntime[changedWeight], currsynpreid[changedWeight], currsynpostid[changedWeight], diffWeight[changedWeight]]).T))
            
            prevWeight=currWeight
        
        scipy.io.savemat('arch3d/synChanged.mat', mdict={'synChanged': synChanged})

# Calculate average weight of all incoming connections (2nd order) to each EM muscle subpopulation over time
# Useful to analyse the effect that RL has on the weights as a function of  - WORK IN PROGRESS--
def synWeightsEM():
    
    print "Calculating average weight to each EM muscle subpopulation"
    
    EMstart=256 # starting cell index of EM population
    EMlength=192 # number of EM cells
    EMsubpopo = [None]*4
    for i in range(4): EMsubpop[i] = [arange(EMstart+i,EMstart+EMlength,4)] # create array with list of ids for each EM cell subpopulation  
    popnames=['P','ES','IS','ILS','EM','IM','IML'];
    popinds=[ 2, 41, 42, 43, 44, 45, 46];
    popshape=['o',  '^',  'h', '*',  '^',  'h',  '*', '8'];
    
    # convert connectivity to python arrays
    conpreid=array(h.col.connsnq.getcol("id1"), 'i')
    conpostid=array(h.col.connsnq.getcol("id2"), 'i')
    #condelay=array(h.col.connsnq.getcol("del"))
    #condistance=array(h.col.connsnq.getcol("dist"))
    conweight1=array(h.col.connsnq.getcol("wt1"))
    conweight2=array(h.col.connsnq.getcol("wt2"))
    
    order = lexsort((conpreid, conpostid)) # order by post and then pre
    conpreid = conpreid[order]
    conpostid = conpostid[order]
    #condelay = condelay[order]
    #condistance = condistance[order]
    conweight1 = conweight1[order]
    conweight2 = conweight2[order]
    
    maxW = max(conweight1.max(), conweight2.max())
    conweight1 = conweight1/maxW    # normalize so can multiply together
    conweight2 = conweight2/maxW 
    
    # convert synaptic changes over time (since connsnq doesn't change with t, can add)
    h('objref synChanges') 
    #h('nqsy.tog()')
    h('nqsy.select("t",">=", 0)')
    h('synChanges = nqsy')
    synpreid=array(h.synChanges.getcol("id1"), 'i')
    synpostid=array(h.synChanges.getcol("id2"), 'i')
    synweight=array(h.synChanges.getcol("wg"))
    syntime=array(h.synChanges.getcol("t"))
    
    order = lexsort((synpreid, synpostid)) # order syn by post and then pre
    synpreid = synpreid[order]
    synpostid = synpostid[order]
    synweight=synweight[order]        
    syntime=syntime[order]    
    
    # define t
    tMax = syntime.max()
    tstep = 50
    T = arange(0, tMax, tstep)
    
    # create array to store average weight change for each subpopulation
    weightSubpops = array((T.len(), 4))        
    
    def showRF(targetCells):
        onlyAMPA = 1
        # find all cells projecting to the target cell for each layer
        afferentCells = array([], 'i')
        weightColors = []
        for i in range(len(targetCells)):
            targetAfferents = array(conpreid[where(conpostid==targetCells[i])], 'i') # get afferent of next target cell
            if targetAfferents.any(): # if any 
                for iAfferent in targetAfferents: # for each new afferent 
                    if (iAfferent in afferentCells): # check if already included in list
                        if onlyAMPA: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                        else: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                        weightColors[afferentCells==iAfferent] +=  weightColors2
                    else: # if not included in list
                        afferentCells=append(afferentCells, iAfferent) # add to list of afferent cells
                        if onlyAMPA: weightColors = append(weightColors, conweight1[targetAfferents]); # calculate weigths
                        else: weightColors = append(weightColors, conweight1[afferentCells]-conweight2[afferentCells]); # calculate weigths
                
                if nSyn2: # 2 synaptic connections
                    for targetCell2 in targetAfferents:  # loop through all afferent cells (which now become target cells)
                        targetAfferents2 = array(conpreid[where(conpostid==targetCell2)], 'i') # get new set of afferent cells
                        if targetAfferents2.any():
                            for iAfferent in targetAfferents2: # loop through all afferent cells
                                if (iAfferent in afferentCells): # if already included in list
                                    if onlyAMPA: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                    else: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                    weightColors[afferentCells==iAfferent] +=  weightColors2
                                else:
                                    afferentCells=append(afferentCells, iAfferent) # add to list of afferent cells
                                    if onlyAMPA: weightColors2 = (conweight1[targetCell2]) * (conweight1[iAfferent])
                                    else: weightColors2 = (conweight1[targetCell2]+conweight2[targetCell2]) * (conweight1[iAfferent]+conweight2[iAfferent])
                                    weightColors=append(weightColors, weightColors2)
                
                #weightColors = weightColors/weightColors.max() # normalize weights to 1
                weight = weightColors.sum()
                return weight
    

    for b in T:
        # calculate sum of weight for each subpopulation
        
        for i in range(4):
            targetCells = EMsubpop[i]
            weight=calculateRF(targetCells) # plot new RF
            weightSubpops[b, i] = weight
            
            # update weights 
            synBin = where(syntime == b )[0]
            if (size(synBin)>0 and size(synBin)<size(conweight1)):
                # update connection matrix based on syn wieght gains
                conweight1[range(len(synBin))] = conweight1[range(len(synBin))] * synweight[synBin]
                                            
    figure()
    plot(T, weightSubpops)
    show()

# plot raster of saved sim
def raster(nqsdir, nqsparams, tstop=4):
    wseedvals =[120456, 398115, 534031, 796321, 895199] # seed values for filename
    iseedvals = [1235, 2837, 3955, 4506, 6789]
    # load connectivity file
    outfilestem = '"%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-spk.nqs"' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])
    filename = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-spk.nqs' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])

    if os.path.isfile(filename):
        print "loading "+outfilestem
        # get errors from nqs
        h('objref nqaload')
        h('nqaload = new NQS(%s)'%outfilestem)
        h('snq[0]=nqaload')
        h.tstop=tstop;
        h('{gg() drit(0, tstop) grlines()}')
    else:
        print "file not found: "+outfilestem

# plot motor spike rate for each muscls subpopulation, from saved sim
def plotvEM(nqsdir, nqsparams,xmin=4,xmax=400, ymax=60):
    wseedvals =[120456, 398115, 534031, 796321, 895199] # seed values for filename
    iseedvals = [1235, 2837, 3955, 4506, 6789]
    # load connectivity file
    outfilestem = '"%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-nqaupd.nqs"' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])
    filename = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-nqaupd.nqs' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])

    if os.path.isfile(filename):
        print "loading "+outfilestem
        # get errors from nqs
        h('objref nqaload')
        h('nqaload = new NQS(%s)'%outfilestem)
        shext=array(h.nqaload.getcol("shext"))
        shflex=array(h.nqaload.getcol("shflex"))
        elext=array(h.nqaload.getcol("elext"))
        elflex=array(h.nqaload.getcol("elflex"))

        figsize=[800,500]
        fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
        colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
        #plot(shext,color = colorlist[3],label='shext')
        #plot(shflex,color = colorlist[5],label='shflex')
        #plot(elext,color = colorlist[6],label='elext')
        #plot(elflex,color = colorlist[9],label='elflex')
        x=arange(xmin,xmax)
        width=1
        bar(x,shext[xmin:xmax],width, color = colorlist[2],label='shext')
        bar(x,shflex[xmin:xmax],width, bottom=shext[xmin:xmax],color = colorlist[4],label='shflex')
        bar(x,elext[xmin:xmax],width, bottom=shflex[xmin:xmax]+shext[xmin:xmax],color = colorlist[6],label='elext')
        bar(x,elflex[xmin:xmax],width,bottom=elext[xmin:xmax]+shflex[xmin:xmax]+shext[xmin:xmax], color = colorlist[9],label='elflex')

        #xlim([xmin, xmax])
        ylim([0, ymax])
        legend(loc='upper right')

    else:
        print "file not found: "+outfilestem

# Show graphs of weight vs. time for each cell subpopulation (subdivided into muscle groups eg. Pshext->ES, Pshflex->ES)
# One subplot for each layer/pop - 1) Pshflex->,Pshext->, ... 2) ES-> 3) IS, ISL, 4) EMshflex->, EMshext->, 5) IM, IML 
def popWeights(loadnqs=0, nqsfile = "", nqsdir="", nqsparams=[0,0,0,0,0,0,0], savefig=0, savenpy=0, loadnpy=0, animate=0, saveanim=0, tinterval=0):
    
    def index2d(myList, v):
        for i, x in enumerate(myList):
            if v in x:
                return (i)
        
    def calculateWeightMatrix():
        print "calculating final weight matrix..."
        #weightMatrix = ndarray((sum(popslen), sum(popslen)))        
        w, h = len(p), len(p)
        weightMatrix = [[0] * w for i in range(h)]
        onlyweight1 = 1
        for i in range(len(synBin)):
            if (conweight1[i]>0):# || conweight2[i]>0):
                if onlyweight1: weightMatrix[index2d(subpops,conpreid[i])][index2d(subpops,conpostid[i])] += conweight1[i]
                else: weightMatrix[index2d(subpops,conpreid[i])][index2d(subpops,conpostid[i])]  += conweight1[i] + conweight2[i]
        weightMatrix=array(weightMatrix)
        for i in range(len(subpops)):
            for j in range(len(subpops)):
                weightMatrix[i][j]=weightMatrix[i][j]/(len(subpops[i])*len(subpops[j]))
        return weightMatrix
        
    # alternative method under construction
    def calculateWeightMatrix2():
        #weightMatrix = ndarray((sum(popslen), sum(popslen)))        
        w, h = len(p), len(p)
        weightMatrix = [[0] * w for i in range(h)]
        onlyAMPA = 1
        for pre in range(len(subpops)):
            for post in range(len(subpops)):
                indicespre=[]
                indicespost=[]
                for i in conpreid: indicespre.append(i in subpops[pre])
                for i in conpostid: indicespost.append(i in subpops[post])
                indices = indicespre * indicespost
                if (conweight1[i]>0):# || conweight2[i]>0):
                    if onlyAMPA: weightMatrix[pre][post] += conweight1[indices].sum()
                    else: weightMatrix[pre][post]  += conweight1[indices].sum() + conweight2[indices].sum()
        return weightMatrix

    # visualization options
    figsize = [1000,800] # Figure size in pixels
    weightsCmap = 'hot_r'# 'jet'#'YlOrRd' #'jet'#'hot'#'autumn'

    # create figures
    ion()
    if animate: # only create figure 1 (evolution over time) if aniamte=1
        fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
        #fig.subplots_adjust(left=0.02) # Less space on left
        #fig.subplots_adjust(right=0.93) # Less space on right
        fig.subplots_adjust(bottom=0.08) # Less space on bottom    
        subplot1 = subplot(221) # create subplots
        subplot2 = subplot(222)
        subplot3 = subplot(223)
        subplot4 = subplot(224)

    #border = 2
    #weightplot.set_xlim([xlocs.min()-border, xlocs.max()+border]) # set x-y lims
    
    # cell populations parameters
    popslen = [192, 44, 20, 192, 44, 20,192]
    #popnames=['ES','IS','ILS','EM','IM','IML', 'P'];
    popinds=[41, 42, 43, 44, 45, 46, 2];
    popshape=['o',  '^',  'h', '*',  '^',  'h',  '*', '8'];
    
    # cell subpopulations (incuding muscle groups)
    subpops= [None]*13
    p = {'Pse':0, 'Psf':1,'Pee':2,'P':3,'ES':4,'IS':5,'ILS':6,'EMse':7,'EMsf':8,'EMee':9,'EMef':10,'IM':11,'IML':12}
    p2 = {0:'Pse', 1:'Psf',2:'Pee',3:'Pef',4:'ES',5:'IS',6:'ILS',7:'EMse',8:'EMsf',9:'EMee',10:'EMef',11:'IM',12:'IML'}
    
    cellslist=[] # list of cells with muscle groups together
    popstart=0
    index=0
    for i in range(len(popslen)):
        popend = popstart+popslen[i]
        if i==3 or i==6: # for P and EM divide into 4 groups 
            for j in range(4): 
                subpops[index] = list(arange(popstart+j,popend,4)) # create array with list of ids for each cell subpopulation  
                index+=1
        else:
            subpops[index]=list(arange(popstart,popend))
            index+=1
        popstart+=popslen[i] # increase popstart
    
    # reorder to have P first
    neworder=[9,10,11,12,0,1,2,3,4,5,6,7,8]
    subpops = [ subpops[i] for i in neworder]

    # load data from neuron variables
    if (not loadnpy):
        wseedvals =[120456, 398115, 534031, 796321, 895199] # seed values for filename
        iseedvals = [1235, 2837, 3955, 4506, 6789]
        if loadnqs: # load data from nqs saved file
            print "loading data from nqs files..."
            # load connectivity file
            if (nqsfile !=""):
                outfilestem = '"%s-con.nqs"' % (nqsfile)
                filename = '%s-con.nqs' % (nqsfile) 
            else:
                outfilestem = '"%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-con.nqs"' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])
                filename = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-con.nqs' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])

            if os.path.isfile(filename):
                print "loading "+outfilestem
                # get errors from nqs
                h('objref nqaload')
                h('nqaload = new NQS(%s)'%outfilestem)
                # convert connectivity to python arrays
                conpreid=array(h.nqaload.getcol("id1"), 'i')
                conpostid=array(h.nqaload.getcol("id2"), 'i')
                #condelay=array(h.col.connsnq.getcol("del"))
                #condistance=array(h.col.connsnq.getcol("dist"))
                conweight1=array(h.nqaload.getcol("wt1"))
                conweight2=array(h.nqaload.getcol("wt2"))

            # load weights file
            if (nqsfile != ""):
                outfilestem = '"%s-syn.nqs"' % (nqsfile)
                filename = '%s-syn.nqs' % (nqsfile) 
            else:
                outfilestem = '"%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-syn.nqs"' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])
                filename = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-syn.nqs' % (nqsdir, nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], iseedvals[nqsparams[5]], wseedvals[nqsparams[6]])
    
            if os.path.isfile(filename):
                print "loadding "+outfilestem
                # get errors from nqs
                h('objref nqaload2')
                h('nqaload2 = new NQS(%s)'%outfilestem)
                # convert connectivity to python arrays
                #h('objref synChanges') 
                #h('nqaload2.select("t",">=", 0)')
                #h('synChanges = nqaload2')
                synpreid=array(h.nqaload2.getcol("id1"), 'i')
                synpostid=array(h.nqaload2.getcol("id2"), 'i')
                synweight=array(h.nqaload2.getcol("wg"))
                syntime=array(h.nqaload2.getcol("t"))


            else: # load from current simulation environment
                print "loading data from current sim..."
                # convert connectivity to python arrays
                conpreid=array(h.col.connsnq.getcol("id1"), 'i')
                conpostid=array(h.col.connsnq.getcol("id2"), 'i')
                #condelay=array(h.col.connsnq.getcol("del"))
                #condistance=array(h.col.connsnq.getcol("dist"))
                conweight1=array(h.col.connsnq.getcol("wt1"))
                conweight2=array(h.col.connsnq.getcol("wt2"))
    
                h('objref synChanges') 
                h('nqsy.select("t",">=", 0)')
                h('synChanges = nqsy')
                synpreid=array(h.synChanges.getcol("id1"), 'i')
                synpostid=array(h.synChanges.getcol("id2"), 'i')
                synweight=array(h.synChanges.getcol("wg"))
                syntime=array(h.synChanges.getcol("t"))

            if 'conpreid' in locals():
                order = lexsort((conpostid, conpreid)) # order by post and then pre
                #condelay = condelay[order]
                #condistance = condistance[order]
                conpreid = conpreid[order]
                conpostid = conpostid[order]
                conweight1 = conweight1[order]
                conweight2 = conweight2[order]
                
                # normalize conweight1 and conweight2 separately
                #maxW = max(conweight1.max(), conweight2.max())
                #conweight1 = conweight1/maxW    # normalize so can multiply together
                #conweight2 = conweight2/maxW 
                
                #conweight1Original = deepcopy(conweight1) # make copy of weights at t=0
                
                # sum conweight1+conweight2 and normalize
                conweight = conweight1+conweight2
                conweightOriginal = conweight / max(conweight)
                
                # convert synaptic changes over time (since connsnq doesn't change with t, can add)
                #synpreid=synpreid[0:len(conpreid)]
                order = lexsort((synpostid, synpreid, syntime)) # order syn by post and then pre
                synpreid = synpreid[order]
                synpostid = synpostid[order]
                synweight=synweight[order]        
                syntime=syntime[order]    
            else: 
                print "Coulndn't load file: "+outfilestem
        

    fig2 = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    fontsiz=14
    matrixsubplot = subplot(111)
    matrixsubplot.set_xlabel('postsynaptic population',fontsize=fontsiz)
    matrixsubplot.set_ylabel('presynaptic population',fontsize=fontsiz)

    # plot static weight matrix - use last set of weights
    if (not animate): 
        if loadnpy:
            weightMatrix=np.load('weightMatrix.npy')
        else:
            # load last set of weights
            synBin = where(syntime == syntime.max())[0]
            if (size(synBin)>0 and size(synBin)<size(conweight1)):
                # update connection matrix based on syn wieght gains
                #conweight1[range(len(synBin))] = (conweight1Original[range(len(synBin))] * synweight[synBin]) - conweight1Original[range(len(synBin))] # absolute weight increase
                conweight1[range(len(synBin))] = synweight[synBin] # relative increase 
                #conweight1[range(len(synBin))] = conweight1Original[range(len(synBin))] * synweight[synBin] # final weight
                #conweight1 = conweightOriginal # original weights (includes P population)
                #synBin = zeros(len(conweight1)) # TEMPORARY LINE to include P population weights - remove after!

            # calculate sum of weight for each subpopulation
            weightMatrix=calculateWeightMatrix()
            weightMatrix = weightMatrix/weightMatrix.max()
            EorI = [1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1] # define excitatory vs inhibitory pops
            for (x,y),value in ndenumerate(weightMatrix): # make weight negative if inhib
                if (EorI[x] == -1):
                    weightMatrix[x,y] = -1 * value 
        

        
        # plot 2dmap
        #weightMatrix[6,6] = weightMatrix[5,5]
        im=matrixsubplot.pcolor(array(weightMatrix), cmap=bicolormap(gap=0.05, mingreen=-0.6,redbluemix=0.0,epsilon=0.005), edgecolors='k', linewidths=1, vmin=-1, vmax=1)
        xticks([float(i)+0.5 for i in p.values()], p.keys())
        yticks([float(i)+0.5 for i in p.values()], p.keys())
            
        matrixsubplot.axis([3, len(p), 3,len(p)])  # set the limits of the plot to the limits of the data
        #matrixsuplot.tick_params(labelsize=8)
        
        c=colorbar(im, ax=matrixsubplot, label='normalized effective connectivity',fontsize=fontsiz)#, fraction=0.09,pad=0.01)     
        try:
            filename
        except:    
            fig2.canvas.set_window_title('Data loaded from saved array')
        else:
            fig2.canvas.set_window_title(filename)
        show()
        #fig.tight_layout()
        if savenpy:
            np.save("weightMatrix.npy", weightMatrix)
        if savefig:
            fig2.savefig("connectivity.pdf",format='pdf')            
            #fig2.savefig('gif/%s_p1-%d_p2-%d_p3-%d_p4-%d_p5-%d_i-%d_w-%d_wmat.png' % (nqsdir[5:], nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], nqsparams[5], nqsparams[6])) # Save PNG file
        
    if animate:
        # define t
        if 'syntime' in globals():
            tMax = syntime.max()
        else:
            tMax = 35000
        tstep = tinterval
        T = arange(0, tMax+tstep, tstep)
        print T
        if load: 
            weightMatrixT=np.load('weightMatrixT.npy')
        else:
            weightMatrixT = [None] * len(T)
        t=0
        for b in T:
            if load:
                weightMatrix=weightMatrixT[t]
            else:
                # calculate sum of weight for each subpopulation
                weightMatrix=calculateWeightMatrix() # plot new RF
                weightMatrixT[t]=weightMatrix
            
            t+=1    
            if t==1: #fixed max value 
                plotvmax=1.2*weightMatrix.max()
            
            #plot 2dmap  , save and wait ; make video
            im=matrixsubplot.pcolor(weightMatrix, cmap=weightsCmap, vmin=0, vmax=plotvmax)
            xticks(p.values(), p.keys())
            yticks(p.values(), p.keys())
            matrixsubplot.axis([0, len(p), 0,len(p)]) 
            if t==1:
                c=colorbar(im, ax=matrixsubplot)
            else:
                c.set_clim([0,plotvmax])
                c.draw_all()
            show()
            #fig.tight_layout()
            if dosave: fig.savefig('tmpmovie%04i.png' % b) # Save PNG files for movie
            if dosave!=1: pause(1)#(1e-10) 
            # update weights 
            if (not load):
                synBin = where(syntime == b )[0]
                if (size(synBin)>0 and size(synBin)<size(conweight1)):
                    # update connection matrix based on syn wieght gains
                    conweight1[range(len(synBin))] = conweight1Original[range(len(synBin))] * synweight[synBin]
                    #conweight1[range(len(synBin))] = synweight[synBin]
        
        # create video from images            
        if saveanim:
            moviedir='data/movies/'
            moviename='14jan28_weights.mpg'
            print('Making movie...')
            system('mkdir -p ' + moviedir)
            system("mencoder 'mf://tmpmovie*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" % (moviedir+moviename))
            print('  Cleaning up...')
            system('rm tmpmovie*.png')

        if savenpy:
            np.save("weightMatrixT.npy", weightMatrixT)
    
        # show plots with evolution of weights over time
        if (not load): weightMatrixT=array(weightMatrixT);
        colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
        legendx=0.3
        legendy=1.0
        legendfontsize=8
        xlabel = 'time (ms)'
        ylabel = 'summed weights'
        #p = {'Pse':0, 'Psf':1,'Pee':2,'Pef':3,'ES':4,'IS':5,'ILS':6,'EMse':7,'EMsf':8,'EMee':9,'EMef':10,'IM':11,'IML':12}
        
        #subplot 1: Px4 (0:3) -> ES (4); IS (5)->ES (4); IM (11)->EM x4 (7:10) total: 9
        subplot1.set_title('P->ES, IS->ES, IM->EM', fontsize=12) # titles
        icolor=0
        for pre in arange(0,4):
            for post in arange(4,5):
                subplot1.plot(T,weightMatrixT[:, pre, post], colorlist[icolor], label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
        for pre in arange(5,6):
            for post in arange(4,5):
                subplot1.plot(T,weightMatrixT[:, pre, post], colorlist[icolor], label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
        for pre in arange(11,12):
            for post in arange(7,11):
                subplot1.plot(T,weightMatrixT[:, pre, post], colorlist[icolor], label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
        subplot1.legend(loc='upper right', bbox_to_anchor=(legendx, legendy),  borderaxespad=0., prop={'size':legendfontsize})
        subplot1.set_ylabel(ylabel)
        subplot1.set_xlabel(xlabel)
        subplot1.grid(True)

        # ES (4)-> ES, EMx4, IS, ISL (4:10) total: 8
        subplot2.set_title('ES', fontsize=12)
        icolor=0
        for pre in arange(4,5):
            for post in arange(4,11):
                subplot2.plot(T,weightMatrixT[:, pre, post], colorlist[icolor], label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
        subplot2.legend(loc='upper right', bbox_to_anchor=(legendx, legendy),  borderaxespad=0., prop={'size':legendfontsize})
        subplot2.set_ylabel(ylabel)
        subplot2.set_xlabel(xlabel)    
        subplot2.grid(True)

        # EMx4 (7:10)-> EMx4 (7:10) total:16 (use dotted vs line))
        subplot3.set_title('EM recurrent', fontsize=12)
        icolor=0
        linesty='-'
        for pre in arange(7,11):
            for post in arange(7,11):
                subplot3.plot(T,weightMatrixT[:, pre, post], colorlist[icolor],  linestyle = linesty, label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
                if icolor>7:
                    icolor=0
                    linesty='--'
        subplot3.legend(loc='upper right', bbox_to_anchor=(legendx, legendy),  borderaxespad=0., prop={'size':legendfontsize})
        subplot3.set_ylabel(ylabel)
        subplot3.set_xlabel(xlabel)
        subplot3.grid(True)

        # EMx4 (7:10)->  ES, IM, IML (3, 11:12) total:12
        subplot4.set_title('EM', fontsize=12)
        icolor=0
        linesty='-'
        for pre in arange(7,11):
            for post in (4,11,12):
                subplot4.plot(T,weightMatrixT[:, pre, post], colorlist[icolor], linestyle = linesty, label=str(p2[pre])+'->'+str(p2[post]))
                icolor+=1
                if icolor>5:
                    icolor=0
                    linesty='--'
        subplot4.legend(loc='upper right', bbox_to_anchor=(legendx, legendy),  borderaxespad=0., prop={'size':legendfontsize})
        subplot4.set_ylabel(ylabel)
        subplot4.set_xlabel(xlabel)
        subplot4.grid(True)
        show()

# call popWeight to plot the weights of the 4 different targets in a specific sim (with specific params)
def popWeights4Targ(nqsdir, nqsparams=[0,0,0,0], savefig=0):
    for i in range(4):
        popWeights(1, nqsdir, [nqsparams[0],nqsparams[1],nqsparams[2],nqsparams[3],i, 0,0], savefig)


# Compare population weights for each target after training
def popMuscles(nqsdir="", plotfig = 0, savefig = 0, save2Matlab = 1):
    
    def index2d(myList, v):
        for i, x in enumerate(myList):
            if v in x:
                return (i)
        
    def calculateWeightMatrix():
        print "calculating final weight matrix..."
        #weightMatrix = ndarray((sum(popslen), sum(popslen)))        
        w, h = len(p), len(p)
        weightMatrix = [[0] * w for i in range(h)]
        onlyweight1 = 1
        for i in range(len(synBin)):
            if (conweight1[i]>0):# || conweight2[i]>0):
                if onlyweight1: weightMatrix[index2d(subpops,conpreid[i])][index2d(subpops,conpostid[i])] += conweight1[i]
                else: weightMatrix[index2d(subpops,conpreid[i])][index2d(subpops,conpostid[i])]  += conweight1[i] + conweight2[i]
        weightMatrix=array(weightMatrix)
        for i in range(len(subpops)):
            for j in range(len(subpops)):
                weightMatrix[i][j]=weightMatrix[i][j]/(len(subpops[i])*len(subpops[j]))
        return weightMatrix
        
    # visualization options
    figsize = [1000,800] # Figure size in pixels
    weightsCmap = 'hot_r'# 'jet'#'YlOrRd' #'jet'#'hot'#'autumn'

    # create figures
    ion()

    # cell populations parameters
    popslen = [192, 44, 20, 192, 44, 20,192]
    #popnames=['ES','IS','ILS','EM','IM','IML', 'P'];
    popinds=[41, 42, 43, 44, 45, 46, 2];
    popshape=['o',  '^',  'h', '*',  '^',  'h',  '*', '8'];
    
    # cell subpopulations (incuding muscle groups)
    subpops= [None]*13
    p = {'Pse':0, 'Psf':1,'Pee':2,'P':3,'ES':4,'IS':5,'ILS':6,'EMse':7,'EMsf':8,'EMee':9,'EMef':10,'IM':11,'IML':12}
    p2 = {0:'Pse', 1:'Psf',2:'Pee',3:'Pef',4:'ES',5:'IS',6:'ILS',7:'EMse',8:'EMsf',9:'EMee',10:'EMef',11:'IM',12:'IML'}
    
    cellslist=[] # list of cells with muscle groups together
    popstart=0
    index=0
    for i in range(len(popslen)):
        popend = popstart+popslen[i]
        if i==3 or i==6: # for P and EM divide into 4 groups 
            for j in range(4): 
                subpops[index] = list(arange(popstart+j,popend,4)) # create array with list of ids for each cell subpopulation  
                index+=1
        else:
            subpops[index]=list(arange(popstart,popend))
            index+=1
        popstart+=popslen[i] # increase popstart
    
    # reorder to have P first
    neworder=[9,10,11,12,0,1,2,3,4,5,6,7,8]
    subpops = [ subpops[i] for i in neworder]

    # load data from neuron variables
    wseedvals =[120456, 398115, 534031, 796321, 895199] # seed values for filename
    iseedvals = [1235, 1235+(2*17),2837, 3955, 4506, 6789]

    print "loading data from nqs files..."
    # load connectivity file
    targets=[0,1,2,3]
    loaded = 0
    for itarget in targets:
        iseedval = iseedvals[0]
        for wseedval in wseedvals:
            loaded = 0
            outfilestem = '"%s/target-%d_i-%d_w-%d_train-con.nqs"' % (nqsdir, itarget, iseedval, wseedval)
            filename = '%s/target-%d_i-%d_w-%d_train-con.nqs' % (nqsdir,itarget, iseedval, wseedval)

            if os.path.isfile(filename):
                print "loading "+outfilestem
                # get errors from nqs
                h('objref nqaload')
                h('nqaload = new NQS(%s)'%outfilestem)
                # convert connectivity to python arrays
                conpreid=array(h.nqaload.getcol("id1"), 'i')
                conpostid=array(h.nqaload.getcol("id2"), 'i')
                #condelay=array(h.col.connsnq.getcol("del"))
                #condistance=array(h.col.connsnq.getcol("dist"))
                conweight1=array(h.nqaload.getcol("wt1"))
                conweight2=array(h.nqaload.getcol("wt2"))
                loaded += 1

            # load weights file

            outfilestem = '"%s/target-%d_i-%d_w-%d_train-syn.nqs"' % (nqsdir, itarget, iseedval, wseedval)
            filename = '%s/target-%d_i-%d_w-%d_train-syn.nqs' % (nqsdir, itarget, iseedval, wseedval)

            if os.path.isfile(filename):
                print "loading "+outfilestem
                # get errors from nqs
                h('objref nqaload2')
                h('nqaload2 = new NQS(%s)'%outfilestem)
                # convert connectivity to python arrays
                #h('objref synChanges') 
                #h('nqaload2.select("t",">=", 0)')
                #h('synChanges = nqaload2')
                synpreid=array(h.nqaload2.getcol("id1"), 'i')
                synpostid=array(h.nqaload2.getcol("id2"), 'i')
                synweight=array(h.nqaload2.getcol("wg"))
                syntime=array(h.nqaload2.getcol("t"))
                loaded += 1


            if loaded == 2: #'conpreid' in locals():
                order = lexsort((conpostid, conpreid)) # order by post and then pre
                #condelay = condelay[order]
                #condistance = condistance[order]
                conpreid = conpreid[order]
                conpostid = conpostid[order]
                conweight1 = conweight1[order]
                conweight2 = conweight2[order]
                
                # normalize conweight1 and conweight2 separately
                #maxW = max(conweight1.max(), conweight2.max())
                #conweight1 = conweight1/maxW    # normalize so can multiply together
                #conweight2 = conweight2/maxW 
                
                #conweight1Original = deepcopy(conweight1) # make copy of weights at t=0
                
                # sum conweight1+conweight2 and normalize
                conweight = conweight1+conweight2
                conweightOriginal = conweight / max(conweight)
                
                # convert synaptic changes over time (since connsnq doesn't change with t, can add)
                #synpreid=synpreid[0:len(conpreid)]
                order = lexsort((synpostid, synpreid, syntime)) # order syn by post and then pre
                synpreid = synpreid[order]
                synpostid = synpostid[order]
                synweight=synweight[order]        
                syntime=syntime[order]    

        
                if plotfig:
                    fig2 = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
                    fontsiz=14
                    matrixsubplot = subplot(111)
                    matrixsubplot.set_xlabel('postsynaptic population',fontsize=fontsiz)
                    matrixsubplot.set_ylabel('presynaptic population',fontsize=fontsiz)

                # plot static weight matrix - use last set of weights
                # load last set of weights
                synBin = where(syntime == syntime.max())[0]
                if (size(synBin)>0 and size(synBin)>size(conweight1)):
                    synBin = synBin[0:(len(synBin)/2)-1]
                #synBin = where(syntime == 4000)[0]
                if (size(synBin)>0 and size(synBin)<size(conweight1)):
                    # update connection matrix based on syn wieght gains
                    conweight1[range(len(synBin))] = (conweightOriginal[range(len(synBin))] * synweight[synBin]) - conweightOriginal[range(len(synBin))] # absolute weight increase
                    #conweight1[range(len(synBin))] = synweight[synBin] # relative increase 
                    #conweight1[range(len(synBin))] = conweight1Original[range(len(synBin))] * synweight[synBin] # final weight
                    #conweight1 = conweightOriginal # original weights (includes P population)
                    #synBin = zeros(len(conweight1)) # TEMPORARY LINE to include P population weights - remove after!

                    # calculate sum of weight for each subpopulation
                    weightMatrix=calculateWeightMatrix()
                    #weightMatrix = weightMatrix/weightMatrix.max()
                    EorI = [1, 1, 1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1] # define excitatory vs inhibitory pops
                    for (x,y),value in ndenumerate(weightMatrix): # make weight negative if inhib
                        if (EorI[x] == -1):
                            weightMatrix[x,y] = -1 * value 

                    # plot 2dmap
                    if plotfig:
                        im=matrixsubplot.pcolor(array(weightMatrix), cmap=bicolormap(gap=0.05, mingreen=-0.6,redbluemix=0.0,epsilon=0.005), edgecolors='k', linewidths=1)#, vmin=-1, vmax=1)
                        xticks([float(i)+0.5 for i in p.values()], p.keys())
                        yticks([float(i)+0.5 for i in p.values()], p.keys())
                            
                        matrixsubplot.axis([3, len(p), 3,len(p)])  # set the limits of the plot to the limits of the data
                        #matrixsuplot.tick_params(labelsize=8)
                        
                        c=colorbar(im, ax=matrixsubplot, label='normalized effective connectivity')#,fontsize=fontsiz)#, fraction=0.09,pad=0.01)     
                        try:
                            filename
                        except:    
                            fig2.canvas.set_window_title('Data loaded from saved array')
                        else:
                            fig2.canvas.set_window_title(filename)
                        show()
                    #fig.tight_layout()
                    
                    if savefig:
                        fig2.savefig("connectivity.pdf",format='pdf')            
                        #fig2.savefig('gif/%s_p1-%d_p2-%d_p3-%d_p4-%d_p5-%d_i-%d_w-%d_wmat.png' % (nqsdir[5:], nqsparams[0], nqsparams[1], nqsparams[2], nqsparams[3], nqsparams[4], nqsparams[5], nqsparams[6])) # Save PNG file

                    if save2Matlab:
                        print "saving matlab file with weights..."
                        scipy.io.savemat(('%s/target-%d_i-%d_w-%d_weights.mat'%(nqsdir, itarget, iseedval, wseedval)), \
                                        mdict={'weightMatrix': weightMatrix})
            else: 
                print "Coulndn't load file: "+outfilestem
                              


# Show 2D pcolor of 704x704 cells with connectivity/weights - over time?
def cellWeightsMatrix(animate, dosave, tinterval):
    # visualization options
    figsize = [800,800] # Figure size in pixels
    weightsCmap = 'jet'#'YlOrRd' #'jet'#'hot'#'autumn'

    # create figure
    ion()
    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    #fig.subplots_adjust(left=0.02) # Less space on left
    #fig.subplots_adjust(right=0.93) # Less space on right
    fig.subplots_adjust(bottom=0.08) # Less space on bottom    
    weightplot = subplot(111) # create subplot        
    weightplot.set_title('Weight Matrix', fontsize=12) # titles
    weightplot.set_xlabel('post (Pse,Psf,Pee,Pef,ES,IS,ILS,EMse,EMsf,EMee,EMef,IM,IML)')
    weightplot.set_ylabel('pre (Pse,Psf,Pee,Pef,ES,IS,ILS,EMse,EMsf,EMee,EMef,IM,IML)')
    
    #setp(weightplot.get_xticklabels(), visible=False) # hide x and y ticks
    #setp(weightplot.get_yticklabels(), visible=False)
    
    #border = 2
    #weightplot.set_xlim([xlocs.min()-border, xlocs.max()+border]) # set x-y lims
    
    # cell populations parameters
    popslen = [192, 44, 20, 192, 44, 20,192]
    #popnames=['P','ES','IS','ILS','EM','IM','IML'];
    #popinds=[ 2, 41, 42, 43, 44, 45, 46];
    popshape=['o',  '^',  'h', '*',  '^',  'h',  '*', '8'];
    
    # cell subpopulations (incuding muscle groups)
    subpops= [None]*13
    cellslist=[] # list of cells with muscle groups together
    popstart=0
    index=0
    for i in range(len(popslen)):
        popend = popstart+popslen[i]
        if i==3 or i==6: # for P and EM divide into 4 groups 
            for j in range(4): 
                subpops[index] = list(arange(popstart+j,popend,4)) # create array with list of ids for each EM cell subpopulation  
                cellslist=cellslist+subpops[index]
                index+=1
        else:
            subpops[index]=list(arange(popstart,popend))
            cellslist=cellslist+subpops[index]
            index+=1
        # increase popstart
        popstart+=popslen[i]
    #move P to the beginning
    cellslist=cellslist[-192:]+cellslist[:-192]


    # convert connectivity to python arrays
    conpreid=array(h.col.connsnq.getcol("id1"), 'i')
    conpostid=array(h.col.connsnq.getcol("id2"), 'i')
    #condelay=array(h.col.connsnq.getcol("del"))
    #condistance=array(h.col.connsnq.getcol("dist"))
    conweight1=array(h.col.connsnq.getcol("wt1"))
    conweight2=array(h.col.connsnq.getcol("wt2"))
    
    order = lexsort((conpreid, conpostid)) # order by post and then pre
    conpreid = conpreid[order]
    conpostid = conpostid[order]
    #condelay = condelay[order]
    #condistance = condistance[order]
    conweight1 = conweight1[order]
    conweight2 = conweight2[order]
    
    maxW = max(conweight1.max(), conweight2.max())
    conweight1 = conweight1/maxW    # normalize so can multiply together
    conweight2 = conweight2/maxW 
    
    conweight1Original = deepcopy(conweight1) # make copy of weights at t=0
    
    # convert synaptic changes over time (since connsnq doesn't change with t, can add)
    h('objref synChanges') 
    #h('nqsy.tog()')
    h('nqsy.select("t",">=", 0)')
    h('synChanges = nqsy')
    synpreid=array(h.synChanges.getcol("id1"), 'i')
    synpostid=array(h.synChanges.getcol("id2"), 'i')
    synweight=array(h.synChanges.getcol("wg"))
    syntime=array(h.synChanges.getcol("t"))
    
    order = lexsort((synpreid, synpostid)) # order syn by post and then pre
    synpreid = synpreid[order]
    synpostid = synpostid[order]
    synweight=synweight[order]        
    syntime=syntime[order]    
    
        
    def calculateWeightMatrix():
        #weightMatrix = ndarray((sum(popslen), sum(popslen)))        
        w, h = sum(popslen), sum(popslen)
        weightMatrix = [[0] * w for i in range(h)]
        onlyAMPA = 0
        for i in range(len(conweight1)):
            if onlyAMPA: weightMatrix[conpreid[i]][conpostid[i]] = conweight1[i]
            else: weightMatrix[conpreid[i]][conpostid[i]] = conweight1[i] + conweight2[i]
        #weightMatrix = [ mylist[i] for i in myorder] # reorder
        weightMatrix2=deepcopy(weightMatrix)
        for i,sublist in enumerate(weightMatrix):
            weightMatrix2[i]=[sublist[j] for j in cellslist] #reorder rows
        
        weightMatrix2=[weightMatrix2[j] for j in cellslist]  # reorder columns
            
        weightMatrix2=array(weightMatrix2)            
        return weightMatrix2
    
    # calculate sum of weight for each subpopulation
    weightMatrix=calculateWeightMatrix() # plot new RF
    
    # plot 2dmap
    im=weightplot.imshow(weightMatrix, cmap=weightsCmap)
    #im=imshow(weightMatrix, cmap=weightsCmap)
    weightplot.axis([0, sum(popslen), 0,sum(popslen)])  # set the limits of the plot to the limits of the data
    c=colorbar(im, ax=weightplot)#, fraction=0.09,pad=0.01)     
    show()
    #fig.tight_layout()
    
    if animate:
        # define t
        tMax = syntime.max()
        tstep = tinterval
        T = arange(0, tMax, tstep)
        for b in T:
            # calculate sum of weight for each subpopulation
            weightMatrix=calculateWeightMatrix() # plot new RF
            
            #plot 2dmap  , save and wait ; make video
            im=weightplot.imshow(weightMatrix, cmap=weightsCmap, vmin=0, vmax=1.5)
            #im=imshow(weightMatrix, cmap=weightsCmap, vmin=0, vmax=1.5)
            #plt.imshow(z)
            #plt.pcolor(Z)
        
            show()
            #fig.tight_layout()
            if dosave: fig.savefig('tmpmovie%04i.png' % b) # Save PNG files for movie
            if dosave!=1: pause(1e-10) 
            
            # update weights 
            synBin = where(syntime == b )[0]
            if (size(synBin)>0 and size(synBin)<size(conweight1)):
                # update connection matrix based on syn wieght gains
                conweight1[range(len(synBin))] = conweight1Original[range(len(synBin))] * synweight[synBin]
                
    if dosave:
        moviedir='data/movies/'
        moviename='14jan28_weights.mpg'
        print('Making movie...')
        system('mkdir -p ' + moviedir)
        system("mencoder 'mf://tmpmovie*.png' -mf type=png:fps=10 -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o %s" % (moviedir+moviename))
        print('  Cleaning up...')
        system('rm tmpmovie*.png')
                                            
# calculate the position of center-out targets in joe's BMI experiments to use in sims                    
def calculateTargetsCenterOut():
    monkeyArmLen=[0.0839, 0.2038]
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634] # added average hand size (included in monkey dimensions)
    factor = 2*(varmLen[0]+varmLen[1])/(monkeyArmLen[0]+monkeyArmLen[1])
    monkeyCenter = [0.050536, 0.095114]
    varmAngleCenter = [0.62, 1.53] # 0.53, 1.41
    varmCenter = angles2pos(varmAngleCenter[0], varmAngleCenter[1], varmLen[0], varmLen[1])
    targetDist = 0.04
    targetDist45deg = 0.0282

    print "\nvarmCenter:"
    print varmCenter
    
    targets = zeros((8,2))
    targets_pos = zeros((8,2))

    targets_pos[0] = varmCenter[0]+targetDist*factor, varmCenter[1]+0 #right
    targets_pos[1] = varmCenter[0]-targetDist*factor, varmCenter[1]+0 #left
    targets_pos[2] = varmCenter[0], varmCenter[1]+targetDist*factor #top
    targets_pos[3] = varmCenter[0], varmCenter[1]-targetDist*factor #bottom
    targets_pos[4] = varmCenter[0]+targetDist45deg*factor, varmCenter[1]+targetDist45deg*factor #right-top
    targets_pos[5] = varmCenter[0]-targetDist45deg*factor, varmCenter[1]+targetDist45deg*factor #left-top
    targets_pos[6] = varmCenter[0]+targetDist45deg*factor, varmCenter[1]-targetDist45deg*factor #left-bottom
    targets_pos[7] = varmCenter[0]-targetDist45deg*factor, varmCenter[1]-targetDist45deg*factor #rightbottom

    print "\ntarget positions:"
    print targets_pos

    print "\ntarget displacements:"
    print targets_pos - varmCenter

    targets[0] = pos2angles(varmCenter[0]+targetDist*factor, varmCenter[1]+0, varmLen[0], varmLen[1]) #right
    targets[1] = pos2angles(varmCenter[0]-targetDist*factor, varmCenter[1]+0, varmLen[0], varmLen[1])#left
    targets[2] = pos2angles(varmCenter[0], varmCenter[1]+targetDist*factor, varmLen[0], varmLen[1]) #top
    targets[3] = pos2angles(varmCenter[0], varmCenter[1]-targetDist*factor, varmLen[0], varmLen[1]) #bottom
    targets[4] = pos2angles(varmCenter[0]+targetDist45deg*factor, varmCenter[1]+targetDist45deg*factor, varmLen[0], varmLen[1]) #right-top
    targets[5] = pos2angles(varmCenter[0]-targetDist45deg*factor, varmCenter[1]+targetDist45deg*factor, varmLen[0], varmLen[1]) #left-top
    targets[6] = pos2angles(varmCenter[0]+targetDist45deg*factor, varmCenter[1]-targetDist45deg*factor, varmLen[0], varmLen[1]) #left-bottom
    targets[7] = pos2angles(varmCenter[0]-targetDist45deg*factor, varmCenter[1]-targetDist45deg*factor, varmLen[0], varmLen[1]) #rightbottom
    
    print "\ntarget angular displacements:"
    print targets - varmAngleCenter

    print "\nabsolute target angular displacements:"
    print sum(abs(targets - varmAngleCenter),1)
    
    print "\ntarget angles:"
    print targets

    return targets    
    
# convert cartesian position to joint angles
def pos2angles(x,y,l1,l2):
    
    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; 
    
    return [shang,elang]

# convert joint angles to cartesian position
def angles2pos(shang,elang,l1,l2):
    armAng=zeros(2)
    armLen=zeros(2)
    armAng[0]=shang
    armAng[1]=elang
    armLen[0]=l1
    armLen[1]=l2
    
    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])
    
    return wristPosx,wristPosy    

# calculates deviation of trajectory from ideal (straight line) traj to target
def trajDeviation(x,y,center,targetInd, targetXdist):
    # subtract center position from trajectory
    x=x-center[0]
    y=y-center[1]

    # set angle of rotation as a function of target
    targetAngles = [0, - np.pi, -np.pi/2, np.pi/2, -np.pi/4, -3*np.pi/4, np.pi/4, 3*np.pi/4] 
    rotAng = targetAngles[targetInd]
    
    # rotate trajectory to position over positive x-axis (=target 0, center-right)
    xrot = x*cos(rotAng) - y*sin(rotAng)
    yrot = x*sin(rotAng) + y*cos(rotAng)

    # include only trajectory points within trajDist*target x pos
    # enables comparison of different targets irrespective of time to target
    trajDist=1.1
    xlim=where(xrot<trajDist*(targetXdist-center[0]))[0]

    # add penalty for displacements going in wrong direction
    xneg=where(xrot<0)
    yrot[xneg]=yrot[xneg]+abs(xrot[xneg])

    # calculate deviation as the sum of the absolute y displacements
    #ydisp = sum(abs(yrot)) #sum
    ydisp = mean(abs(yrot[xlim])) #mean
    return ydisp

# calculate the average (replaced by min!) distance between the trajectory and the target during trange
def trajAccuracy(x,y,targetPos,varmLen, trange, method):
    #method = 'min'#

    if method == 'sum':
        dist=0
    elif method == 'min':
        dist=1
    for i in trange:
        if method == 'sum':
            dist += sqrt((x[i]-targetPos[0])**2 + (y[i]-targetPos[1])**2) # sum
        elif method == 'min':
            dist = min(dist, sqrt((x[i]-targetPos[0])**2 + (y[i]-targetPos[1])**2)) # min
    #print "min dist:"+str(dist)
    if method == 'sum':
        dist = dist/len(trange) # average 
    return dist

# plot center out reaching trajectory for single target
def plotCenterOutTrajSingleParam(simdatadir, trajStartArg, trajStopArg, niseeds=100):

    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    elbowCenter = [varmLen[0] * cos(0.62), varmLen[0] * sin(0.62)]
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456, 398115, 534031, 796321, 895199]
    iseedvals = arange(1235,  1235+(17*niseeds), 17)
    
    # parameter values
    param1_range = arange(1,2,1)#  arange(0.8,1.24,0.04)#arange(200,375,25)#arange(20,180,20)##arange(10,110,10)# [0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2]##arange(0.8,1.24,0.04)# arange(1,9)#arange(50,550,50)# arange(10,110,10)## #[0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15]#[10,20,30,40,50,60,70,80,90,100]#[50, 100,250,500,750,1000]#arange(20,180,20)##
    param2_range = [0,1,2,3] # [0,1,2,3,4,5,6,7]

    # Visualization options
    figsize =[800,800]# [1300,500] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2*2 #??
    targetSize = 30 *2
    targetMarker = 'x'
    colorlist = ['black','blue', 'red', 'green', 'brown', 'cyan','darkgrey', 'magenta', 'orange','yellow']
    trajLine = 1
    trajSize = 1
    trajMarker = '.'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)
        
    # create subplots for each of the param values tested
    ion()
    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    fig.subplots_adjust(left=0.02) # Less space on left
    fig.subplots_adjust(right=0.93) # Less space on right
    fig.subplots_adjust(bottom=0.08) # Less space on bottom
    subplotsx = 1#5
    subplotsy = 1#2
    maxSubplots = subplotsx*subplotsy
    param1Subplots =  [None] * min(len(param1_range), maxSubplots)
    border = 0.1
    for i in range(min(len(param1_range), maxSubplots)):
        param1Subplots[i] = subplot(subplotsy, subplotsx, i+1)
        #setp(param1Subplots[i].get_xticklabels(), visible=False) # hide x and y ticks
        #setp(param1Subplots[i].get_yticklabels(), visible=False)
        param1Subplots[i].set_xlim([targetsPos[1][0]-border*1.5, targetsPos[0][0]+border]) # set x-y lims
        param1Subplots[i].set_ylim([targetsPos[3][1]-border*2.5, targetsPos[2][1]+border])
    
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    iparam1=0
    param1Subplots[iparam1].set_title('Trajectories for multiple input seeds')
    # Loop over param2 values
    iparam2=-1
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        iwseed = -1
        # plot target location
        param1Subplots[iparam1].scatter(targetsPos[iparam2][0], targetsPos[iparam2][1],  c=colorlist[iparam2], marker=targetMarker, linewidth=targetLine, s=targetSize)
        # Loop over wiring seed...
        for wseed in wseedvals:
            iwseed = iwseed + 1
            iiseed = -1
            # Loop over input seed...
            for iseed in iseedvals:
                iiseed = iiseed + 1
                # set filename
                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                #outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                outfilestem = '"%s/target-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param2, iseed, wseed)
                filename = '%s/target-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param2, iseed, wseed)
                if os.path.isfile(filename):
                    print outfilestem
                    # get errors from nqs
                    h('nqaload = new NQS(%s)'%outfilestem)
                    t=array(h.nqaload.getcol("t"))
                    x=array(h.nqaload.getcol("x"))
                    y=array(h.nqaload.getcol("y"))
                    shAng = array(h.nqaload.getcol("ang0"))
                    #print shAng
                    elAng = array(h.nqaload.getcol("ang1"))
                    #print elAng
                    param1Subplots[iparam1].scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam2], marker=trajMarker, linewidth=trajLine, s=trajSize) 
                
    # draw arm
    plot([0,elbowCenter[0]], [0, elbowCenter[1]],'k-',lw=4) #elbow
    plot([elbowCenter[0], varmCenter[0]], [elbowCenter[1], varmCenter[1]],'k-',lw=4)
    tight_layout()
    fig.savefig('gif/%s_traj_%d.png' % (simdatadir[5:], trajStop))

# plot center out reaching trajectories
def plotCenterOutTraj(simdatadir, trajStartArg, trajStopArg, param1Arg=None):

    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456, 398115]#, 534031, 796321, 895199]
    iseedvals = [1235, 2837]#, 3955, 4506, 6789]
    
    # parameter values
    if param1Arg==None:
        param1_range =  arange(0.8,1.24,0.04)#arange(200,375,25)#arange(20,180,20)##arange(10,110,10)# [0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2]##arange(0.8,1.24,0.04)# arange(1,9)#arange(50,550,50)# arange(10,110,10)## #[0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15]#[10,20,30,40,50,60,70,80,90,100]#[50, 100,250,500,750,1000]#arange(20,180,20)##
    else:
        param1_range = param1Arg
    param2_range = [0,1,2,3,4,5,6,7]

    # Visualization options
    figsize =[1300,500] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2 #
    targetSize = 30
    targetMarker = 'x'
    colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
    trajLine = 1
    trajSize = 1
    trajMarker = '.'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)
        
    # create subplots for each of the param values tested
    ion()
    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    fig.subplots_adjust(left=0.02) # Less space on left
    fig.subplots_adjust(right=0.93) # Less space on right
    fig.subplots_adjust(bottom=0.08) # Less space on bottom
    subplotsx = 5
    subplotsy = 2
    maxSubplots = subplotsx*subplotsy
    param1Subplots =  [None] * min(len(param1_range), maxSubplots)
    border = 0.1
    for i in range(min(len(param1_range), maxSubplots)):
        param1Subplots[i] = subplot(subplotsy, subplotsx, i+1)
        setp(param1Subplots[i].get_xticklabels(), visible=False) # hide x and y ticks
        setp(param1Subplots[i].get_yticklabels(), visible=False)
        param1Subplots[i].set_xlim([targetsPos[1][0]-border, targetsPos[0][0]+border]) # set x-y lims
        param1Subplots[i].set_ylim([targetsPos[3][1]-border, targetsPos[2][1]+border])
    
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    iparam1 = -1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        iparam2 = -1
        param1Subplots[iparam1].set_title('param1= '+str(param1))
        # Loop over param2 values
        for param2 in param2_range:
            iparam2 = iparam2 + 1
            iwseed = -1
            # plot target location
            param1Subplots[iparam1].scatter(targetsPos[iparam2][0], targetsPos[iparam2][1],  c=colorlist[iparam2], marker=targetMarker, linewidth=targetLine, s=targetSize)
            # Loop over wiring seed...
            for wseed in wseedvals:
                iwseed = iwseed + 1
                iiseed = -1
                # Loop over input seed...
                for iseed in iseedvals:
                    iiseed = iiseed + 1
                    skipValue=0
                    # set filename and check if file exists
                    outfilestem = '%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                    if os.path.isfile(outfilestem):
                        print "loading "+str(outfilestem)
                    else:
                        outfilestem = '%s/p1-%.2f_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                        if os.path.isfile(outfilestem):
                            print "loading "+str(outfilestem)
                        else: 
                            outfilestem = '%s/p1-%.3f_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                            if os.path.isfile(outfilestem):
                                print "loading "+str(outfilestem)
                            else:
                                skipValue=1

                            
                    if skipValue==0:
                        # get errors from nqs
                        h('nqaload = new NQS("%s")'%outfilestem)
                        print("extracting trajectory data...")
                        t=array(h.nqaload.getcol("t"))
                        x=array(h.nqaload.getcol("x"))
                        y=array(h.nqaload.getcol("y"))
                        #print t,x,y
                    
                        param1Subplots[iparam1].scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam2], marker=trajMarker, linewidth=trajLine, s=trajSize) 
                
    tight_layout()

# plot center out reaching trajectories for 8 sims with same date
def plotCenterOutTraj8(simdatadirRoot, trajStartArg, trajStopArg):    
# call plotCenterOutTraj for each sim in same date
    
    numSims = 8
    
    param1_range =  []#[None] * numSims
    param1_range.append([10,20,30,40,50])
    param1_range.append([10,20,30,40,50,60,70,80,90,100])
    param1_range.append(arange(0.8,1.24,0.04))#([0.012, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15])
    param1_range.append([10,20,30,40,50,60,70,80,90,100])
    param1_range.append([25, 50,75, 100, 125, 150, 175, 200, 225, 250])
    param1_range.append([20,40,60,80,100,120,140,160])
    param1_range.append([200,225,250,275,300,325,350])
    param1_range.append([50, 100,250,500,750,1000])
    
    for i in range(numSims):
        plotCenterOutTraj(simdatadirRoot+'_sim'+str(i+1), trajStartArg, trajStopArg, param1_range[i])

# plot center out reaching trajectories for 4param batch sims
def plotCenterOutTraj4param(simdatadir, trajStartArg, trajStopArg):
    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    #targetsAng=[[deg2rad(90),deg2rad(70)],[deg2rad(90),deg2rad(100)],[deg2rad(100),deg2rad(80)],[deg2rad(100),deg2rad(70)],[deg2rad(95),deg2rad(95)]]

    #targetsAng=[[deg2rad(90),deg2rad(100)],[deg2rad(100),deg2rad(80)]] # only targets 14 and 15

    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    
    wseedvals =[120456,398115, 534031, 796321, 895199]
    iseedvals = [1235, 2837, 3955, 4506, 6789]

    # set param range by hand
    param1_range = arange(200,400,50) # EMNoiseRate train
    param2_range= arange(0, 8, 2) # EMNoiseMuscleGroup
    param3_range = arange(500,2000,500) # EMMuscleNoiseChangeDT
    param4_range = arange(8,40,8) # exploreTot
    param5_range = arange(4) # only four targets
    wseed_range = arange(2)
    inseed_range = arange(2)

    # read from param file
    with open('%s/params'% (simdatadir)) as f:
        param1_range, param2_range, param3_range, param4_range, param5_range, wseed_range, inseed_range = pickle.load(f)

    # Visualization options
    figsize =[1000,700] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2 #
    targetSize = 30
    targetMarker = 'x'
    #colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
    colorlist = ['blue', 'green', 'red', 'magenta', 'darkgrey']
    trajLine = 1
    trajSize = 1
    trajMarker = '.'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)
        
    # create subplots for each of the param values tested
    ion()
    for param1 in param1_range:
        for param2 in param2_range:
            fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
            fig.canvas.set_window_title('param1 = '+str(param1)+', param2 = '+str(param2))
            #fig.suptitle('param1 = '+str(param1)+', param2 = '+str(param2), fontsize=8)
            fig.subplots_adjust(left=0.02) # Less space on left
            fig.subplots_adjust(right=0.93) # Less space on right
            fig.subplots_adjust(bottom=0.08) # Less space on bottom
            subplotsy = len(param3_range)
            subplotsx = len(param4_range)
            maxSubplots = subplotsx*subplotsy
            paramSubplots =  [None] * min(len(param3_range)*len(param4_range), maxSubplots)
            border = 0.1
            for i in range(len(param3_range)*len(param4_range)):
                paramSubplots[i] = subplot(subplotsy, subplotsx, i+1)
                setp(paramSubplots[i].get_xticklabels(), visible=False) # hide x and y ticks
                setp(paramSubplots[i].get_yticklabels(), visible=False)
                paramSubplots[i].set_xlim([targetsPos[1][0]-border, targetsPos[0][0]+border]) # set x-y lims
                paramSubplots[i].set_ylim([targetsPos[3][1]-border, targetsPos[2][1]+border])
                
            # plot data for each of the 4 starting conditions in same color; different color for each target
            # Loop over param3 values
            h('objref nqaload')
            iparam3 = -1
            for param3 in param3_range:
                iparam3 = iparam3 + 1
                iparam4 = -1
                # Loop over param3 values
                for param4 in param4_range:
                    iparam4 = iparam4 + 1
                    iparam5 = -1
                    paramSubplots[iparam3*len(param4_range)+iparam4].set_title('param3 = '+str(param3)+', param4 = '+str(param4),fontsize=8)
                    # Loop over param5 values (target)
                    for param5 in param5_range:
                        iparam5 = iparam5 + 1
                        iwseed = -1
                        # plot target location
                        
                        paramSubplots[iparam3*len(param4_range)+iparam4].scatter(targetsPos[iparam5][0], targetsPos[iparam5][1],  c=colorlist[iparam5], marker=targetMarker, linewidth=targetLine, s=targetSize)
                        # Loop over wiring seed...
                        for wseed in array(wseedvals)[wseed_range]:
                            iwseed = iwseed + 1
                            iiseed = -1
                            # Loop over input seed...
                            for iseed in array(iseedvals)[inseed_range]:
                                iiseed = iiseed + 1
                                skipValue=0
                                # set filename and check if file exists
                                outfilestem = '%s/p1-%d_p2-%d_p3-%d_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                if os.path.isfile(outfilestem):
                                    print "loading "+str(outfilestem)
                                else:
                                    outfilestem = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                    if os.path.isfile(outfilestem):
                                        print "loading "+str(outfilestem)
                                    else: 
                                        outfilestem = '%s/p1-%.1f_p2-%d_p3-%.3f_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                        if os.path.isfile(outfilestem):
                                            print "loading "+str(outfilestem)
                                        else:
                                            skipValue=1
                                    
                                if skipValue==0:
                                    # get errors from nqs
                                    h('nqaload = new NQS("%s")'%outfilestem)
                                    print("extracting trajectory data...")
                                    t=array(h.nqaload.getcol("t"))
                                    x=array(h.nqaload.getcol("x"))
                                    y=array(h.nqaload.getcol("y"))
                                    #print t,x,y
                            
                                    paramSubplots[iparam3*len(param4_range)+iparam4].scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam5], marker=trajMarker, linewidth=trajLine, s=trajSize) 
                
            # set tight layout and save figure
            tight_layout()
            fig.savefig('gif/%s_p1-%.2f_p2-%.2f_traj.png' % (simdatadir[5:], param1, param2)) # Save PNG files for movie
    close('all')        
        

    
# save spiking data and arm trajectory to matlab file


# plot center out reaching trajectories
def plotCenterOutAng4param(simdatadir, trajStartArg, trajStopArg):
    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    targetsAng=[[deg2rad(90),deg2rad(70)],[deg2rad(90),deg2rad(100)],[deg2rad(100),deg2rad(80)],[deg2rad(100),deg2rad(70)],[deg2rad(95),deg2rad(95)]]

    targetsAng=[[deg2rad(90),deg2rad(100)],[deg2rad(100),deg2rad(80)]] # only targets 14 and 15

    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    
    wseedvals =[120456,398115, 534031, 796321, 895199]
    iseedvals = [1235, 2837, 3955, 4506, 6789]

    # set param range by hand
    param1_range = arange(200,400,50) # EMNoiseRate train
    param2_range= arange(0, 8, 2) # EMNoiseMuscleGroup
    param3_range = arange(500,2000,500) # EMMuscleNoiseChangeDT
    param4_range = arange(8,40,8) # exploreTot
    param5_range = arange(4) # only four targets
    wseed_range = arange(2)
    inseed_range = arange(2)

    # read from param file
    with open('%s/params'% (simdatadir)) as f:
        param1_range, param2_range, param3_range, param4_range, param5_range, wseed_range, inseed_range = pickle.load(f)

    param5_range=arange(14,15,1)
    
    # Visualization options
    figsize =[1000,700] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2 #
    targetSize = 30
    targetMarker = 'x'
    #colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
    colorlist = ['blue', 'green', 'red', 'magenta', 'darkgrey']
    trajLine = 1
    trajSize = 1
    trajSize2 = 1
    trajMarker = '.'
    trajMarker2 = '*'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)
        
    # create subplots for each of the param values tested
    ion()
    for param1 in param1_range:
        for param2 in param2_range:
            fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
            fig.canvas.set_window_title('param1 = '+str(param1)+', param2 = '+str(param2))
            #fig.suptitle('param1 = '+str(param1)+', param2 = '+str(param2), fontsize=8)
            fig.subplots_adjust(left=0.02) # Less space on left
            fig.subplots_adjust(right=0.93) # Less space on right
            fig.subplots_adjust(bottom=0.08) # Less space on bottom
            subplotsy = len(param3_range)
            subplotsx = len(param4_range)
            maxSubplots = subplotsx*subplotsy
            paramSubplots =  [None] * min(len(param3_range)*len(param4_range), maxSubplots)
            border = 0.1
            for i in range(len(param3_range)*len(param4_range)):
                paramSubplots[i] = subplot(subplotsy, subplotsx, i+1)
                setp(paramSubplots[i].get_xticklabels(), visible=False) # hide x and y ticks
                setp(paramSubplots[i].get_yticklabels(), visible=False)
                #paramSubplots[i].set_xlim([targetsPos[1][0]-border, targetsPos[0][0]+border]) # set x-y lims
                #paramSubplots[i].set_ylim([targetsPos[3][1]-border, targetsPos[2][1]+border])
                
            # plot data for each of the 4 starting conditions in same color; different color for each target
            # Loop over param3 values
            h('objref nqaload')
            iparam3 = -1
            for param3 in param3_range:
                iparam3 = iparam3 + 1
                iparam4 = -1
                # Loop over param3 values
                for param4 in param4_range:
                    iparam4 = iparam4 + 1
                    iparam5 = -1
                    paramSubplots[iparam3*len(param4_range)+iparam4].set_title('param3 = '+str(param3)+', param4 = '+str(param4),fontsize=8)
                    # Loop over param5 values (target)
                    for param5 in param5_range:
                        iparam5 = iparam5 + 1
                        iwseed = -1
                        # plot target location
                        
                        #paramSubplots[iparam3*len(param4_range)+iparam4].scatter(targetsPos[iparam5][0], targetsPos[iparam5][1],  c=colorlist[iparam5], marker=targetMarker, linewidth=targetLine, s=targetSize)
                        # Loop over wiring seed...
                        for wseed in array(wseedvals)[wseed_range]:
                            iwseed = iwseed + 1
                            iiseed = -1
                            # Loop over input seed...
                            for iseed in array(iseedvals)[inseed_range]:
                                iiseed = iiseed + 1
                                skipValue=0
                                # set filename and check if file exists
                                outfilestem = '%s/p1-%d_p2-%d_p3-%d_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                if os.path.isfile(outfilestem):
                                    print "loading "+str(outfilestem)
                                else:
                                    outfilestem = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                    if os.path.isfile(outfilestem):
                                        print "loading "+str(outfilestem)
                                    else: 
                                        outfilestem = '%s/p1-%.1f_p2-%d_p3-%.3f_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                        if os.path.isfile(outfilestem):
                                            print "loading "+str(outfilestem)
                                        else:
                                            skipValue=1
                                    
                                if skipValue==0:
                                    # get errors from nqs
                                    h('nqaload = new NQS("%s")'%outfilestem)
                                    print("extracting trajectory data...")
                                    t=array(h.nqaload.getcol("t"))
                                    x=array(h.nqaload.getcol("x"))
                                    y=array(h.nqaload.getcol("y"))
                                    shang=array(h.nqaload.getcol("ang0"))
                                    elang=array(h.nqaload.getcol("ang1"))
                                
                            
                                    #paramSubplots[iparam3*len(param4_range)+iparam4].scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam5], marker=trajMarker, linewidth=trajLine, s=trajSize)  
                                    # plot shoulder ang
                                    paramSubplots[iparam3*len(param4_range)+iparam4].scatter(t[trajStart:trajStop], shang[trajStart:trajStop],  edgecolor=colorlist[(iwseed*2)+iiseed], marker=trajMarker, linewidth=trajLine, s=trajSize)       
                                    # plot elbow ang
                                    paramSubplots[iparam3*len(param4_range)+iparam4].scatter(t[trajStart:trajStop], elang[trajStart:trajStop],  edgecolor=colorlist[(iwseed*2)+iiseed], marker=trajMarker2, linewidth=trajLine, s=trajSize2)   
                                    # plot shoulder+elbow targets
                                    #paramSubplots[iparam3*len(param4_range)+iparam4].plot([t[trajStart], t[trajStop]], [targetsAng[0][0], targetsAng[0][0]], color='black', linestyle='..')
                                    #paramSubplots[iparam3*len(param4_range)+iparam4].plot([t[trajStart], t[trajStop]], [targetsAng[0][1], targetsAng[0][1]], color='darkgrey', linestyle='..')

                                    paramSubplots[iparam3*len(param4_range)+iparam4].scatter(t[trajStart:trajStop], [targetsAng[0][0]]*len(t[trajStart:trajStop]), edgecolor='black', marker=trajMarker, linewidth=trajLine, s=trajSize)
                                    paramSubplots[iparam3*len(param4_range)+iparam4].scatter(t[trajStart:trajStop], [targetsAng[0][1]]*len(t[trajStart:trajStop]), edgecolor='darkgrey', marker=trajMarker, linewidth=trajLine, s=trajSize)



                
            # set tight layout and save figure
            tight_layout()
            fig.savefig('gif/%s_p1-%.2f_p2-%.2f_traj.png' % (simdatadir[5:], param1, param2)) # Save PNG files for movie
    close('all')        
        
# save spiking data and arm trajectory to matlab file


# plot center out reaching trajectories

def errorCenterOutTraj4param(simdatadir, trajStartArg, trajStopArg, accuracyFactor=0, deviationFactor=0):
    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    
    wseedvals =[120456,398115, 534031, 796321, 895199]
    iseedvals = [1235, 2837, 3955, 4506, 6789]
    
    param1_range = arange(200,400,50) # EMNoiseRate train
    param2_range= arange(0, 8, 2) # EMNoiseMuscleGroup
    param3_range = arange(500,2000,500) # EMMuscleNoiseChangeDT
    param4_range = arange(8,40,8) # exploreTot
    param5_range = arange(4) # only four targets
    wseed_range = arange(2)
    inseed_range = arange(2)

    with open('%s/params'% (simdatadir)) as f:
        param1_range, param2_range, param3_range, param4_range, param5_range, wseed_range, inseed_range = pickle.load(f)

    # Visualization options
    figsize =[700,700] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2 #
    targetSize = 30
    targetMarker = 'x'
    #colorlist = ['black','darkgrey','blue', 'cyan', 'green', 'brown', 'red', 'magenta', 'orange','yellow']
    colorlist = ['blue', 'green', 'red', 'magenta']
    trajLine = 1
    trajSize = 1
    trajMarker = '.'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)

    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    #fig.canvas.set_window_title('param1 = '+str(param1)+', param2 = '+str(param2))
    #fig.suptitle('param1 = '+str(param1)+', param2 = '+str(param2), fontsize=8)
    fig.subplots_adjust(left=0.02) # Less space on left
    fig.subplots_adjust(right=0.93) # Less space on right
    fig.subplots_adjust(bottom=0.08) # Less space on bottom
    subplotsy = len(param1_range)
    subplotsx = len(param2_range)
    maxSubplots = subplotsx*subplotsy
    paramSubplots =  [None] * maxSubplots #min(len(param2_range)*len(param2_range), maxSubplots)
    border = 0.1
    for i in range(len(param1_range)*len(param2_range)):
        paramSubplots[i] = subplot(subplotsy, subplotsx, i+1)
        
    # create subplots for each of the param values tested
    errorAvg = zeros((len(param1_range),len(param2_range),len(param3_range),len(param4_range)))
    ion()
    iparam1 =-1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        iparam2 = -1
        for param2 in param2_range:
            iparam2 = iparam2 + 1
            # plot data for each of the 4 starting conditions in same color; different color for each target
            # Loop over param3 values
            
            accuracy = zeros((len(param3_range),len(param4_range),len(param5_range),len(wseed_range),len(inseed_range)))
            deviation = zeros((len(param3_range),len(param4_range),len(param5_range),len(wseed_range),len(inseed_range)))

            accuracy_mean = zeros((len(param3_range),len(param4_range),len(param5_range)))
            deviation_mean = zeros((len(param3_range),len(param4_range),len(param5_range)))

            accuracy_max = zeros((len(param3_range),len(param4_range)))
            deviation_max = zeros((len(param3_range),len(param4_range)))

            
            h('objref nqaload')
            iparam3 = -1
            for param3 in param3_range:
                iparam3 = iparam3 + 1
                iparam4 = -1
                paramSubplots[iparam1*len(param2_range)+iparam2].set_title('param1 = '+str(param1)+', param2 = '+str(param2),fontsize=8)
                # Loop over param3 values
                for param4 in param4_range:
                    iparam4 = iparam4 + 1
                    iparam5 = -1
                    # Loop over param5 values (target)
                    for param5 in param5_range:
                        iparam5 = iparam5 + 1
                        iwseed = -1
                        # Loop over wiring seed...
                        for wseed in array(wseedvals)[wseed_range]:
                            iwseed = iwseed + 1
                            iiseed = -1
                            # Loop over input seed...
                            for iseed in array(iseedvals)[inseed_range]:
                                iiseed = iiseed + 1
                                skipValue=0
                                # set filename and check if file exists
                                outfilestem = '%s/p1-%d_p2-%d_p3-%d_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                if os.path.isfile(outfilestem):
                                    print "loading "+str(outfilestem)
                                else:
                                    outfilestem = '%s/p1-%.2f_p2-%.2f_p3-%.2f_p4-%.2f_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                    if os.path.isfile(outfilestem):
                                        print "loading "+str(outfilestem)
                                    else: 
                                        outfilestem = '%s/p1-%.1f_p2-%d_p3-%.3f_p4-%d_p5-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4, param5, iseed, wseed)
                                        if os.path.isfile(outfilestem):
                                            print "loading "+str(outfilestem)
                                        else:
                                            skipValue=1
                                    
                                if skipValue==0:
                                    # get errors from nqs
                                    h('nqaload = new NQS("%s")'%outfilestem)
                                    print("extracting trajectory data...")
                                    t=array(h.nqaload.getcol("t"))
                                    x=array(h.nqaload.getcol("x"))
                                    y=array(h.nqaload.getcol("y"))
                                    #print t,x,y

                                    accuracy[iparam3,iparam4, iparam5, iwseed, iiseed] = trajAccuracy(x[trajStart:trajStop], y[trajStart:trajStop], targetsPos[param5], varmLen, [int(0.1*(trajStop-trajStart-1)),trajStop-trajStart-1])
                                    deviation[iparam3,iparam4, iparam5, iwseed, iiseed] = trajDeviation(x[trajStart:trajStop], y[trajStart:trajStop], varmCenter, param5, targetsPos[0][0]) 

                        # calculate mean over seeds
                        accuracy_mean[iparam3,iparam4, iparam5]= mean(accuracy[iparam3,iparam4,iparam5,:,:])
                        deviation_mean[iparam3,iparam4, iparam5]= mean(deviation[iparam3,iparam4,iparam5,:,:])

                    # calculate max over targets
                    accuracy_max[iparam3,iparam4]= max(accuracy_mean[iparam3,iparam4,:])
                    deviation_max[iparam3,iparam4]= max(deviation_mean[iparam3,iparam4,:])

                    # calcualate error as max of acc and dev for each param3 and param4
                    #error = accuracyFactor*accuracy+deviationFactor*deviation

                    # store 
                    #errorAvg[iparam1, iparam2, iparam3,iparam4] = mean(error[iparam3,iparam4,:,:,:])
                    errorAvg[iparam1, iparam2, iparam3,iparam4] = max(accuracy_max[iparam3,iparam4], deviation_max[iparam3,iparam4])
                                                        
    #print errorAvg

    iparam1 =-1
    for param1 in param1_range:
        iparam1 = iparam1 + 1
        iparam2 = -1
        for param2 in param2_range:
            iparam2 = iparam2 + 1

            # plot error results
            #imsubplot=paramSubplots[iparam1*len(param2_range)+iparam2].imshow(errorAvg, vmin=0, vmax=0.04, origin='upper', interpolation='none')
            imsubplot=paramSubplots[iparam1*len(param2_range)+iparam2].pcolor(errorAvg[iparam1,iparam2], vmin=nanmin(errorAvg), vmax=nanmax(errorAvg))
            paramSubplots[iparam1*len(param2_range)+iparam2].invert_yaxis()

            paramSubplots[iparam1*len(param2_range)+iparam2].set_xticklabels(param4_range)
            paramSubplots[iparam1*len(param2_range)+iparam2].set_xticks(arange(0.5, len(param4_range), 1))
            paramSubplots[iparam1*len(param2_range)+iparam2].set_yticklabels(param3_range)
            paramSubplots[iparam1*len(param2_range)+iparam2].set_yticks(arange(0.5, len(param3_range), 1))
            paramSubplots[iparam1*len(param2_range)+iparam2].tick_params(labelsize=8)


            #scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam5], marker=trajMarker, linewidth=trajLine, s=trajSize) 
                
    # set tight layout and save figure
    cfraction = 1
    cpad = 1
    cshrink = 1
    cfontsize = 8
    #c=colorbar(paramSubplots[iparam1*len(param2_range)+iparam2, ax=rasterPsubplot, fraction=cfraction, pad=cpad, shrink=cshrink)
    tight_layout()
    c=colorbar(imsubplot,ax=paramSubplots[iparam1*len(param2_range)+iparam2])
    c.ax.tick_params(labelsize=cfontsize) 
    fig.savefig('gif/%s_error.png' % (simdatadir[5:])) # Save PNG files for movie
    
# evaluation function over 4 tragets -- note: not adapted to multiple seeds
def errorCenterOutTraj(simdatafile, trajStart=2, trajStop=200, accFactor=1, devFactor=0):
    
    simdatafiles = []
    simdatafiles.append(simdatafile.replace("target_3","target_0"))
    simdatafiles.append(simdatafile.replace("target_3","target_1"))
    simdatafiles.append(simdatafile.replace("target_3","target_2"))
    simdatafiles.append(simdatafile)

    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    
    accuracy = zeros((len(simdatafiles)))
    deviation = zeros((len(simdatafiles)))

    # create subplots for each of the param values tested
    ifile =-1
    for filename in simdatafiles:
        ifile = ifile + 1
            
        h('objref nqaload')
        full_filename = '%s_test-nqa.nqs' % (filename)
        skipValue=0
        if os.path.isfile(full_filename):
            print "loading "+str(full_filename)
        else:
            skipValue=1
            
        if skipValue==0:
            # get errors from nqs
            h('nqaload = new NQS("%s")'%full_filename)
            print("extracting trajectory data...")
            t=array(h.nqaload.getcol("t"))
            x=array(h.nqaload.getcol("x"))
            y=array(h.nqaload.getcol("y"))
            #print t,x,y

            accuracy[ifile] = trajAccuracy(x[trajStart:trajStop], y[trajStart:trajStop], targetsPos[ifile], varmLen, [int(0.01*(trajStop-trajStart-1)),trajStop-trajStart-1], 'sum')
            deviation[ifile] = trajDeviation(x[trajStart:trajStop], y[trajStart:trajStop], varmCenter, ifile, targetsPos[0][0]) 

    # calculate max over targets
    print "accuracy:"
    print accuracy
    print "deviation:"
    print deviation
    accuracy_max= accFactor*max(accuracy)
    deviation_max= devFactor*max(deviation)

    error = max(accuracy_max, deviation_max)
    #saveFile=simdatafile[0:simdatafile.find("_iter")]
    saveFile=simdatafile[0:simdatafile.find("_target")]
    with open('%s_errortmp'% (saveFile), 'w') as f:
        pickle.dump(error, f)

    return error

def errorCenterOutTrajTmp(simdatafile, trajStart=2, trajStop=200):
    simdatafiles = []
    simdatafiles.append(simdatafile.replace("target_3","target_0"))
    simdatafiles.append(simdatafile.replace("target_3","target_1"))
    simdatafiles.append(simdatafile.replace("target_3","target_2"))
    simdatafiles.append(simdatafile)

    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    
    accuracy = zeros((len(simdatafiles)))
    deviation = zeros((len(simdatafiles)))

    # create subplots for each of the param values tested
    ifile =-1
    for filename in simdatafiles:
        ifile = ifile + 1
            
        h('objref nqaload')
        full_filename = '%s_test-nqa.nqs' % (filename)
        skipValue=0
        if os.path.isfile(full_filename):
            print "loading "+str(full_filename)
        else:
            skipValue=1
            
        if skipValue==0:
            # get errors from nqs
            h('nqaload = new NQS("%s")'%full_filename)
            print("extracting trajectory data...")
            t=array(h.nqaload.getcol("t"))
            x=array(h.nqaload.getcol("x"))
            y=array(h.nqaload.getcol("y"))
            #print t,x,y

            accuracy[ifile] = trajAccuracy(x[trajStart:trajStop], y[trajStart:trajStop], targetsPos[ifile], varmLen, [int(0.1*(trajStop-trajStart-1)),trajStop-trajStart-1])
            deviation[ifile] = trajDeviation(x[trajStart:trajStop], y[trajStart:trajStop], varmCenter, ifile, targetsPos[0][0]) 

    # calculate max over targets
    print "accuracy:"
    print accuracy
    print "deviation:"
    print deviation
    accuracy_max= max(accuracy)
    deviation_max= max(deviation)

    error = max(accuracy_max, deviation_max)
    saveFile=simdatafile[0:simdatafile.find("_iter")]
    #with open('%s_errortmp'% (saveFile), 'w') as f:
    #    pickle.dump(error, f)

    return error
    
    
# save spiking data and arm trajectory to matlab file
def save2Matlab(simdatadir, param1Arg, t1Arg, t2Arg):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    iseedvals =     [1235, 2837, 3955, 4506, 6789, 1236, 2838, 3956, 4507, 6790, 1237, 2839, 3957, 4508, 6791, 1238, 2840, 3958, 4509, 6792]# [1235, 2837]#, 3955, 4506, 6789]
    
    # parameter values
    param2_range = [0,1]#[0,1,2,3,4,5,6,7]

            
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    param1=param1Arg
    iparam2 = -1
    # Loop over param2 values
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        iwseed = -1
        # Loop over wiring seed...
        for wseed in wseedvals:
            iwseed = iwseed + 1
            iiseed = -1
            # Loop over input seed...
            for iseed in iseedvals:
                iiseed = iiseed + 1
                # set filename
                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                #outfilestem = '"%s/target-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param2, iseed, wseed)
                outfilestem_format = '%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param1, param2, iseed, wseed)
                
                # get fields from nqs
                if os.path.isfile(outfilestem_format):
                    try:
                        h('nqaload = new NQS(%s)'%outfilestem)
                        h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        t=array(h.nqaload.getcol("t"))
                        xHand = array(h.nqaload.getcol("x"))
                        yHand = array(h.nqaload.getcol("y"))
                        xElbow = array(h.nqaload.getcol("ex"))
                        yElbow = array(h.nqaload.getcol("ey"))
                        shAng = array(h.nqaload.getcol("ang0"))
                        elAng = array(h.nqaload.getcol("ang1"))
                        shExtMusLength = array(h.nqaload.getcol("ML0"))
                        shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                        elExtMusLength = array(h.nqaload.getcol("ML2"))
                        elFlexMusLength = array(h.nqaload.getcol("ML3"))
                        
                        #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                        outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)

                        h('spikesload = new NQS(%s)' % outfilestem)
                        h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        tspike=array(h.spikesload.getcol("t"))
                        cellid=array(h.spikesload.getcol("id"))
                        celltype=array(h.spikesload.getcol("type"))
                        muscleid=array(h.spikesload.getcol("mid"))
                        
                        #save to matlab
                        #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                        #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                        scipy.io.savemat(('%s_target%d_i%d_w%d_spk.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                        mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    
                                    
                        scipy.io.savemat(('%s_target%d_i%d_w%d_arm.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                        mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    
                    except:
                        pass

# save spiking data and arm trajectory to matlab file
def save2Matlab_iseeds(simdatadir, t1Arg=20, t2Arg=1000, simtype='test'):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    #wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    #iseedvals =     [1235, 2837, 3955, 4506, 6789, 1236, 2838, 3956, 4507, 6790, 1237, 2839, 3957, 4508, 6791, 1238, 2840, 3958, 4509, 6792]# [1235, 2837]#, 3955, 4506, 6789]
    wseedvals =[120456, 398115, 534031, 796321, 895199]
    iseedvals = arange(1235,  1235+(17*100), 17)
    
    # parameter values
    param2_range = [0,1,2,3]#,4,5,6,7]

            
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    iparam2 = -1
    # Loop over param2 values
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        iwseed = -1
        # Loop over wiring seed...
        for wseed in wseedvals:
            iwseed = iwseed + 1
            iiseed = -1
            # Loop over input seed...
            for iseed in iseedvals:
                iiseed = iiseed + 1
                # set filename
                outfilestem = '"%s/target-%d_i-%d_w-%d_%s-nqa.nqs"' % (simdatadir, param2, iseed, wseed, simtype)
                outfilestem_format = '%s/target-%d_i-%d_w-%d_%s-nqa.nqs' % (simdatadir, param2, iseed, wseed, simtype)

                # get fields from nqs
                if os.path.isfile(outfilestem_format):
                    try:
                        # get fields from nqs
                        h('nqaload = new NQS(%s)'%outfilestem)
                        h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        t=array(h.nqaload.getcol("t"))
                        xHand = array(h.nqaload.getcol("x"))
                        yHand = array(h.nqaload.getcol("y"))
                        xElbow = array(h.nqaload.getcol("ex"))
                        yElbow = array(h.nqaload.getcol("ey"))
                        shAng = array(h.nqaload.getcol("ang0"))
                        elAng = array(h.nqaload.getcol("ang1"))
                        shExtMusLength = array(h.nqaload.getcol("ML0"))
                        shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                        elExtMusLength = array(h.nqaload.getcol("ML2"))
                        elFlexMusLength = array(h.nqaload.getcol("ML3"))
                        
                        outfilestem = '"%s/target-%d_i-%d_w-%d_%s-spk.nqs"' % (simdatadir, param2, iseed, wseed, simtype)
                                              
                        h('spikesload = new NQS(%s)' % outfilestem)
                        h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        tspike=array(h.spikesload.getcol("t"))
                        cellid=array(h.spikesload.getcol("id"))
                        celltype=array(h.spikesload.getcol("type"))
                        muscleid=array(h.spikesload.getcol("mid"))
                        
                        #save to matlab
                       
                        scipy.io.savemat(('%s/target%d_i%d_w%d_spk.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                        mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    
                                    
                        scipy.io.savemat(('%s/target%d_i%d_w%d_arm.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                        mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    
                    except:
                        pass


# save spiking data and arm trajectory to matlab file
def save2Matlab_traintimes(simdatadir, t1Arg, t2Arg):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    #wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    #iseedvals =     [1235, 2837, 3955, 4506, 6789, 1236, 2838, 3956, 4507, 6790, 1237, 2839, 3957, 4508, 6791, 1238, 2840, 3958, 4509, 6792]# [1235, 2837]#, 3955, 4506, 6789]
    wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    iseedvals = arange(1235,  1235+(17*100), 17)
    
    # parameter values
    param2_range = [0,1,2,3]#[0,1,2,3,4,5,6,7]
    tphase1_range = [0, 2, 4, 6, 8, 10]  # range of training epochs for phase 1
    tphase2_range = [0, 10, 20, 30, 40, 50, 60] # range of training epochs for phase 2
            
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    iparam2 = -1
    # Loop over param2 values
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        for tphase1 in tphase1_range:
            # loop over epochs for training phase 2
            for tphase2 in tphase2_range:
                iwseed = -1
                # Loop over wiring seed...
                for wseed in wseedvals:
                    iwseed = iwseed + 1
                    iiseed = -1
                    # Loop over input seed...
                    for iseed in iseedvals:
                        iiseed = iiseed + 1
                        # set filename
                        outfilestem = '"%s/target-%d_tp1-%d_tp2-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param2, tphase1, tphase2, iseed, wseed)
                        outfilestem_format = '%s/target-%d_tp1-%d_tp2-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param2, tphase1, tphase2, iseed, wseed)
                        
                        print
                        # get fields from nqs
                        if os.path.isfile(outfilestem_format):
                            try:
                                # get fields from nqs
                                h('nqaload = new NQS(%s)'%outfilestem)
                                h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                t=array(h.nqaload.getcol("t"))
                                xHand = array(h.nqaload.getcol("x"))
                                yHand = array(h.nqaload.getcol("y"))
                                xElbow = array(h.nqaload.getcol("ex"))
                                yElbow = array(h.nqaload.getcol("ey"))
                                shAng = array(h.nqaload.getcol("ang0"))
                                elAng = array(h.nqaload.getcol("ang1"))
                                shExtMusLength = array(h.nqaload.getcol("ML0"))
                                shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                                elExtMusLength = array(h.nqaload.getcol("ML2"))
                                elFlexMusLength = array(h.nqaload.getcol("ML3"))
                                
                                outfilestem = '"%s/target-%d_tp1-%d_tp2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param2, tphase1, tphase2, iseed, wseed)

                                h('spikesload = new NQS(%s)' % outfilestem)
                                h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                tspike=array(h.spikesload.getcol("t"))
                                cellid=array(h.spikesload.getcol("id"))
                                celltype=array(h.spikesload.getcol("type"))
                                muscleid=array(h.spikesload.getcol("mid"))
                                
                                #save to matlab
                               
                                scipy.io.savemat(('%s/target%d_tp1%d_tp2%d_i%d_w%d_spk.mat'%(simdatadir, iparam2, tphase1, tphase2, iiseed, iwseed)), \
                                mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    
                                            
                                scipy.io.savemat(('%s/target%d_tp1%d_tp2%d_i%d_w%d_arm.mat'%(simdatadir, iparam2, tphase1, tphase2, iiseed, iwseed)), \
                                mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    
                            except:
                                pass



                    
def save2MatlabTrain(simdatadir, param1Arg, t1Arg, t2Arg):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    iseedvals = [1235]#, 2837]#, 3955, 4506, 6789]
    
    # parameter values
    param2_range = [0,1,2,3,4,5,6,7]

            
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    param1=param1Arg
    iparam2 = -1
    # Loop over param2 values
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        iwseed = -1
        # Loop over wiring seed...
        for wseed in wseedvals:
            iwseed = iwseed + 1
            iiseed = -1
            # Loop over input seed...
            for iseed in iseedvals:
                iiseed = iiseed + 1
                # set filename
                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)

                # get fields from nqs
                h('nqaload = new NQS(%s)'%outfilestem)
                h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                t=array(h.nqaload.getcol("t"))
                xHand = array(h.nqaload.getcol("x"))
                yHand = array(h.nqaload.getcol("y"))
                xElbow = array(h.nqaload.getcol("ex"))
                yElbow = array(h.nqaload.getcol("ey"))
                shAng = array(h.nqaload.getcol("ang0"))
                elAng = array(h.nqaload.getcol("ang1"))
                shExtMusLength = array(h.nqaload.getcol("ML0"))
                shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                elExtMusLength = array(h.nqaload.getcol("ML2"))
                elFlexMusLength = array(h.nqaload.getcol("ML3"))
                
                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_train-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)

                h('spikesload = new NQS(%s)' % outfilestem)
                h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                tspike=array(h.spikesload.getcol("t"))
                cellid=array(h.spikesload.getcol("id"))
                celltype=array(h.spikesload.getcol("type"))
                muscleid=array(h.spikesload.getcol("mid"))
                
                #save to matlab
                #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                scipy.io.savemat(('%s_target%d_i%d_w%d_spk.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    
                            
                scipy.io.savemat(('%s_target%d_i%d_w%d_arm.mat'%(simdatadir, iparam2, iiseed, iwseed)), \
                mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    

def save2MatlabMistPerturb(simdatadir, t1Arg, t2Arg):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456, 398115, 534031, 796321, 895199]
    iseedvals = [1235, 2837, 3955, 4506, 6789]

    param1_range = arange(0,192,1) # mistCell
    param2_range= arange(200, 600, 200) # mistStart
    param3_range = arange(200,201) # mistDuration
    param4_range = arange(250,501,250) # mistRate
    param5_range = arange(1) # only four targets
    wseed_range = arange(1)
    inseed_range = arange(1)

    param1_range = arange(0,192)#192,1) # mistCell
    param2_range= arange(0,1)# 600, 200) # mistStart
    param3_range = arange(1,2) # mistDuration
    param4_range = arange(1,2)#(501,250) # mistRate
    param5_range = arange(1) # only four targets
    wseed_range = arange(1)
    inseed_range = arange(1)

    #with open('%s/params'% (simdatadir)) as f:
    #    param1_range, param2_range, param3_range, param4_range, param5_range, wseed_range, inseed_range = pickle.load(f)

    param1_range = arange(20,192,1) # mistCell
        
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    # Loop over param1 vals
    for param1 in param1_range:
        # Loop over param2 vals
        for param2 in param2_range:
            # Loop over param3 vals
            for param3 in param3_range:
                # Loop over param4 vals
                for param4 in param4_range:
                    # Loop over wiring seed...
                    for wseed in wseed_range:
                        # Loop over input seed...
                        for iseed in inseed_range:
                            # set filename
                            #outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                            outfilestem = '%s/cell-%d_start-%d_dur-%d_rate-%d_test-nqa.nqs' % (simdatadir, param1, param2, param3, param4)
                            
                            # get fields from nqs
                            if os.path.isfile(outfilestem):
                                outfilestem = '"%s/cell-%d_start-%d_dur-%d_rate-%d_test-nqa.nqs"' % (simdatadir, param1, param2, param3, param4)
                                h('nqaload = new NQS(%s)'%outfilestem)
                                h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                t=array(h.nqaload.getcol("t"))
                                xHand = array(h.nqaload.getcol("x"))
                                yHand = array(h.nqaload.getcol("y"))
                                xElbow = array(h.nqaload.getcol("ex"))
                                yElbow = array(h.nqaload.getcol("ey"))
                                shAng = array(h.nqaload.getcol("ang0"))
                                elAng = array(h.nqaload.getcol("ang1"))
                                shExtMusLength = array(h.nqaload.getcol("ML0"))
                                shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                                elExtMusLength = array(h.nqaload.getcol("ML2"))
                                elFlexMusLength = array(h.nqaload.getcol("ML3"))

                                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                                outfilestem = '"%s/cell-%d_start-%d_dur-%d_rate-%d_test-spk.nqs"' % (simdatadir, param1, param2, param3, param4)

                                h('spikesload = new NQS(%s)' % outfilestem)
                                h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                tspike=array(h.spikesload.getcol("t"))
                                cellid=array(h.spikesload.getcol("id"))
                                celltype=array(h.spikesload.getcol("type"))
                                muscleid=array(h.spikesload.getcol("mid"))

                                #save to matlab
                                #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                                #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                                scipy.io.savemat(('%s/cell-%d_start-%d_dur-%d_rate-%d_spk.mat'%(simdatadir, param1, param2, param3, param4)), \
                                mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    

                                scipy.io.savemat(('%s/cell-%d_start-%d_dur-%d_rate-%d_arm.mat'%(simdatadir, param1, param2, param3, param4)), \
                                mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    

def save2MatlabMistProbing(simdatadir, t1Arg, t2Arg):
    target_range = [1,3] #arange(4) # only four targets
    param1_range = [0,1] # kill cells or synapses
    param2_range = [0,5,10] # kill 5% or 10%
    param3_range = arange(0,99)#arange(0,191,1) # mistCell
    param4_range = [0]#[400] # mistStart
    param5_range = [0]#[200] # mistDuration
    param6_range = [0]#[250] # mistRate

    param1_range = [0,1] #kill cells or synapses
    param2_range = [5,10] # kill 5% or 10%
    param3_range = arange(0,191)#arange(0,191,1) # mistCell
    param4_range = [100]#[400] # mistStart
    param5_range = [200]#[200] # mistDuration
    param6_range = [250]#[250] # mistRate

    h('objref nqaload')
    h('objref spikesload')

    # Loop over target vals
    for target in target_range:
        # Loop over param1 vals
        for param1 in param1_range:
            # Loop over param2 vals
            for param2 in param2_range:
                # Loop over param3 vals
                for param3 in param3_range:
                    # Loop over param4 vals
                    for param4 in param4_range:
                        # Loop over param5 vals
                        for param5 in param5_range:
                            # Loop over param6 vals
                            for param6 in param6_range:

                                # set filename
                                outfilestem = '%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_test-nqa.nqs' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)
                                
                                # get fields from nqs
                                if os.path.isfile(outfilestem):
                                    outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_test-nqa.nqs"' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)
                                    h('nqaload = new NQS(%s)'%outfilestem)
                                    h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                    t=array(h.nqaload.getcol("t"))
                                    xHand = array(h.nqaload.getcol("x"))
                                    yHand = array(h.nqaload.getcol("y"))
                                    xElbow = array(h.nqaload.getcol("ex"))
                                    yElbow = array(h.nqaload.getcol("ey"))
                                    shAng = array(h.nqaload.getcol("ang0"))
                                    elAng = array(h.nqaload.getcol("ang1"))
                                    shExtMusLength = array(h.nqaload.getcol("ML0"))
                                    shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                                    elExtMusLength = array(h.nqaload.getcol("ML2"))
                                    elFlexMusLength = array(h.nqaload.getcol("ML3"))

                                    #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                                    outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_test-spk.nqs"' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)

                                    h('spikesload = new NQS(%s)' % outfilestem)
                                    h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                    tspike=array(h.spikesload.getcol("t"))
                                    cellid=array(h.spikesload.getcol("id"))
                                    celltype=array(h.spikesload.getcol("type"))
                                    muscleid=array(h.spikesload.getcol("mid"))

                                    #save to matlab
                                    #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                                    #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                                    scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_spk.mat' % (simdatadir, target, param1, param2, param3, param4, param5, param6)), \
                                    mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    

                                    scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_arm.mat' % (simdatadir, target, param1, param2, param3, param4, param5, param6)), \
                                    mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    

def save2MatlabMistProbingMulti(simdatadir, t1Arg, t2Arg):
    target_range = [1,3] #arange(4) # only four targets
    param1_range = [0,1] # kill cells or synapses
    param2_range = [0,5,10] # kill 5% or 10%
    param3_range = arange(0,99)#arange(0,191,1) # mistCell
    param4_range = [0]#[400] # mistStart
    param5_range = [0]#[200] # mistDuration
    param6_range = [0]#[250] # mistRate

    param1_range = [0,1] #kill cells or synapses
    param2_range = [5,10] # kill 5% or 10%
    param3_range = arange(0,191)#arange(0,191,1) # mistCell
    param4_range = [100]#[400] # mistStart
    param5_range = [200]#[200] # mistDuration
    param6_range = [250]#[250] # mistRate

    h('objref nqaload')
    h('objref spikesload')

    # Loop over target vals
    for target in target_range:
        # Loop over param1 vals
        for param1 in param1_range:
            # Loop over param2 vals
            for param2 in param2_range:
                # Loop over param3 vals
                for param3 in param3_range:
                    # Loop over param4 vals
                    for param4 in param4_range:
                        # Loop over param5 vals
                        for param5 in param5_range:
                            # Loop over param6 vals
                            for param6 in param6_range:

                                # set filename
                                outfilestem = '%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_multi_test-nqa.nqs' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)
                                
                                # get fields from nqs
                                if os.path.isfile(outfilestem):
                                    outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_multi_test-nqa.nqs"' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)
                                    h('nqaload = new NQS(%s)'%outfilestem)
                                    h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                    t=array(h.nqaload.getcol("t"))
                                    xHand = array(h.nqaload.getcol("x"))
                                    yHand = array(h.nqaload.getcol("y"))
                                    xElbow = array(h.nqaload.getcol("ex"))
                                    yElbow = array(h.nqaload.getcol("ey"))
                                    shAng = array(h.nqaload.getcol("ang0"))
                                    elAng = array(h.nqaload.getcol("ang1"))
                                    shExtMusLength = array(h.nqaload.getcol("ML0"))
                                    shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                                    elExtMusLength = array(h.nqaload.getcol("ML2"))
                                    elFlexMusLength = array(h.nqaload.getcol("ML3"))

                                    #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                                    outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_multi_test-spk.nqs"' % \
                                (simdatadir, target, param1, param2, param3, param4, param5, param6)

                                    h('spikesload = new NQS(%s)' % outfilestem)
                                    h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                                    tspike=array(h.spikesload.getcol("t"))
                                    cellid=array(h.spikesload.getcol("id"))
                                    celltype=array(h.spikesload.getcol("type"))
                                    muscleid=array(h.spikesload.getcol("mid"))

                                    #save to matlab
                                    #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                                    #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                                    scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_multi_spk.mat' % (simdatadir, target, param1, param2, param3, param4, param5, param6)), \
                                    mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    

                                    scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_cell-%d_start-%d_dur-%d_rate-%d_multi_arm.mat' % (simdatadir, target, param1, param2, param3, param4, param5, param6)), \
                                    mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    
                 
def save2MatlabMistRepair(simdatadir, t1Arg, t2Arg):
    target_range = [1,3] #arange(4) # only four targets
    param1_range = [0,1] # kill cells or synapses
    param2_range = [0,5,10] # kill 5% or 10%
    param3_range = ['perturb', 'repair']

    h('objref nqaload')
    h('objref spikesload')

    # Loop over target vals
    for target in target_range:
        # Loop over param1 vals
        for param1 in param1_range:
            # Loop over param2 vals
            for param2 in param2_range:
                # Loop over param3 vals
                for param3 in param3_range:

                    # set filename
                    outfilestem = '%s/target-%d_ptype-%d_pperc-%d_%s_test-nqa.nqs' % \
                    (simdatadir, target, param1, param2, param3)
                    
                    # get fields from nqs
                    if os.path.isfile(outfilestem):
                        outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_%s_test-nqa.nqs"' % \
                    (simdatadir, target, param1, param2, param3)
                        h('nqaload = new NQS(%s)'%outfilestem)
                        h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        t=array(h.nqaload.getcol("t"))
                        xHand = array(h.nqaload.getcol("x"))
                        yHand = array(h.nqaload.getcol("y"))
                        xElbow = array(h.nqaload.getcol("ex"))
                        yElbow = array(h.nqaload.getcol("ey"))
                        shAng = array(h.nqaload.getcol("ang0"))
                        elAng = array(h.nqaload.getcol("ang1"))
                        shExtMusLength = array(h.nqaload.getcol("ML0"))
                        shFlexMusLength =  array(h.nqaload.getcol("ML1"))
                        elExtMusLength = array(h.nqaload.getcol("ML2"))
                        elFlexMusLength = array(h.nqaload.getcol("ML3"))

                        #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                        outfilestem = '"%s/target-%d_ptype-%d_pperc-%d_%s_test-spk.nqs"' % \
                    (simdatadir, target, param1, param2, param3)

                        h('spikesload = new NQS(%s)' % outfilestem)
                        h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
                        tspike=array(h.spikesload.getcol("t"))
                        cellid=array(h.spikesload.getcol("id"))
                        celltype=array(h.spikesload.getcol("type"))
                        muscleid=array(h.spikesload.getcol("mid"))

                        #save to matlab
                        #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
                        #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

                        scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_%s_spk.mat' % (simdatadir, target, param1, param2, param3)), \
                        mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    

                        scipy.io.savemat(('%s/target-%d_ptype-%d_pperc-%d_%s_arm.mat' % (simdatadir, target, param1, param2, param3)), \
                        mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    


def save2MatlabMistSingle(simdatadir, t1Arg, t2Arg):
    ###########################
    # read and save spike and trajectory data (t,x,y from nqa) 
    ########################
    
    
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    h('objref spikesload')
    
    # set filename
    #outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
    outfilestem = '%s_test-nqa.nqs' % (simdatadir)
                            
    # get fields from nqs
    if os.path.isfile(outfilestem):
        outfilestem = '"%s_test-nqa.nqs"' % (simdatadir)
        h('nqaload = new NQS(%s)'%outfilestem)
        h('nqaload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
        t=array(h.nqaload.getcol("t"))
        xHand = array(h.nqaload.getcol("x"))
        yHand = array(h.nqaload.getcol("y"))
        xElbow = array(h.nqaload.getcol("ex"))
        yElbow = array(h.nqaload.getcol("ey"))
        shAng = array(h.nqaload.getcol("ang0"))
        elAng = array(h.nqaload.getcol("ang1"))
        shExtMusLength = array(h.nqaload.getcol("ML0"))
        shFlexMusLength =  array(h.nqaload.getcol("ML1"))
        elExtMusLength = array(h.nqaload.getcol("ML2"))
        elFlexMusLength = array(h.nqaload.getcol("ML3"))

        #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_test-spk.nqs"' % (simdatadir, param1, param2, iseed, wseed)
        outfilestem = '"%s_test-spk.nqs"' % (simdatadir)

        h('spikesload = new NQS(%s)' % outfilestem)
        h('spikesload.select("t","[]", %d, %d)'%(t1Arg,t2Arg))
        tspike=array(h.spikesload.getcol("t"))
        cellid=array(h.spikesload.getcol("id"))
        celltype=array(h.spikesload.getcol("type"))
        muscleid=array(h.spikesload.getcol("mid"))

        #save to matlab
        #scipy.io.savemat(('bmm_%s_p1-%.2f_p2-%d_i-%d_w-%d_test_spk.mat'%simdatadir, param1, param2, iseed, wseed), \
        #mdict={'t': t, 'cellId': 'cellId, 'celltype': celltype, 'muscleid': muscleid})    

        scipy.io.savemat(('%s_spk.mat'%(simdatadir)), \
    mdict={'tspike': tspike, 'cellid': cellid, 'celltype': celltype, 'muscleid': muscleid})    

        scipy.io.savemat(('%s_arm.mat'%(simdatadir)), \
    mdict={'t': t, 'xhand': xHand, 'yHand': yHand, 'xElbow': xElbow,  'yElbow': yElbow, 'shAng': shAng, 'elAng': elAng, 'shExtMusLength': shExtMusLength, 'shFlexMusLength': shFlexMusLength,  'elExtMusLength': elExtMusLength, 'elFlexMusLength': elFlexMusLength })    


# plot center out reaching trajectory for single target
def plotNCM14CenterOutTraj(simdatadir, trajStartArg, trajStopArg):

    #######################
    ##### obtain target locations
    #########################
    targetsAng = calculateTargetsCenterOut()
    # convert to cartesian position
    varmLen=[0.4634 - 0.173, 0.7169 - 0.4634]
    varmCenter = angles2pos(0.62, 1.53, varmLen[0], varmLen[1])
    elbowCenter = [varmLen[0] * cos(0.62), varmLen[0] * sin(0.62)]
    targetsPos = zeros((8,2))
    for i in range(len(targetsAng)):
        targetsPos[i]=angles2pos(targetsAng[i][0], targetsAng[i][1], varmLen[0], varmLen[1])
    print targetsPos
    print varmCenter
    
    
    ###########################
    # read and plot trajectory data (t,x,y from nqa) 
    ########################
    wseedvals =[120456]#, 398115]#, 534031, 796321, 895199]
    iseedvals = arange(1235,  1235+(17*100), 17)#, 2837, 3955, 4506, 6789, 1236, 2838, 3956, 4507, 6790, 1237, 2839, 3957, 4508, 6791, 1238, 2840, 3958, 4509, 6792]
#[1235]#, 2837]#, 3955, 4506, 6789]
    
    # parameter values
    param1_range = arange(1,2,1)#  arange(0.8,1.24,0.04)#arange(200,375,25)#arange(20,180,20)##arange(10,110,10)# [0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15, 0.175, 0.2]##arange(0.8,1.24,0.04)# arange(1,9)#arange(50,550,50)# arange(10,110,10)## #[0.01, 0.025, 0.05, 0.075, 0.1, 0.125, 0.15]#[10,20,30,40,50,60,70,80,90,100]#[50, 100,250,500,750,1000]#arange(20,180,20)##
    param2_range = [0,1,2,3] # [0,1,2,3,4,5,6,7]

    # Visualization options
    figsize =[800,800]# [1300,500] # Figure size in pixels
    targetColor = 'green';#array([(1,0.4,0) , (0,0.2,0.8)]) # Define excitatory and inhibitory colors -- orange and turquoise
    targetLine = 2*2 #??
    targetSize = 30 *2
    targetMarker = 'x'
    colorlist = ['black','blue', 'red', 'green', 'brown', 'cyan','darkgrey', 'magenta', 'orange','yellow']
    trajLine = 1
    trajSize = 1
    trajMarker = '.'    
    trajStart = trajStartArg#2 # sample number to start from (each sample = 10 ms)
    trajStop = trajStopArg#400 # sample number to finish on (max 400=4sec)
        
    # create subplots for each of the param values tested
    ion()
    fig = figure(figsize=(figsize[0]/100,figsize[1]/100),dpi=100)
    fig.subplots_adjust(left=0.02) # Less space on left
    fig.subplots_adjust(right=0.93) # Less space on right
    fig.subplots_adjust(bottom=0.08) # Less space on bottom
    subplotsx = 1#5
    subplotsy = 1#2
    maxSubplots = subplotsx*subplotsy
    param1Subplots =  [None] * min(len(param1_range), maxSubplots)
    border = 0.1
    for i in range(min(len(param1_range), maxSubplots)):
        param1Subplots[i] = subplot(subplotsy, subplotsx, i+1)
        #setp(param1Subplots[i].get_xticklabels(), visible=False) # hide x and y ticks
        #setp(param1Subplots[i].get_yticklabels(), visible=False)
        param1Subplots[i].set_xlim([targetsPos[1][0]-border*1.5, targetsPos[0][0]+border]) # set x-y lims
        param1Subplots[i].set_ylim([targetsPos[3][1]-border*2.5, targetsPos[2][1]+border])
    
    # plot data for each of the 4 starting conditions in same color; different color for each target
    # Loop over param1 values
    h('objref nqaload')
    iparam1=0
    param1Subplots[iparam1].set_title('Trajectories for multiple input seeds')
    # Loop over param2 values
    iparam2=-1
    for param2 in param2_range:
        iparam2 = iparam2 + 1
        iwseed = -1

        if param2==0:
            simdatadir="data/14feb07_iseeds7"
            trajStop=200
        elif param2==1:
            simdatadir="data/14mar10_iseeds8"
            trajStop=200
        else:
            simdatadir="data/14mar24_iseeds6"
            trajStop=trajStopArg
        # plot target location
        param1Subplots[iparam1].scatter(targetsPos[iparam2][0], targetsPos[iparam2][1],  c=colorlist[iparam2], marker=targetMarker, linewidth=targetLine, s=targetSize)
        # Loop over wiring seed...
        for wseed in wseedvals:
            iwseed = iwseed + 1
            iiseed = -1
            # Loop over input seed...
            for iseed in iseedvals:
                iiseed = iiseed + 1
                # set filename
                #outfilestem = '"%s/p1-%.2f_p2-%d_i-%d_w-%d_train-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                #outfilestem = '"%s/p1-%d_p2-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param1, param2, iseed, wseed)
                outfilestem = '"%s/target-%d_i-%d_w-%d_test-nqa.nqs"' % (simdatadir, param2, iseed, wseed)
                filename = '%s/target-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, param2, iseed, wseed)
                if os.path.isfile(filename):
                    print outfilestem
                    # get errors from nqs
                    h('nqaload = new NQS(%s)'%outfilestem)
                    t=array(h.nqaload.getcol("t"))
                    x=array(h.nqaload.getcol("x"))
                    y=array(h.nqaload.getcol("y"))
                    param1Subplots[iparam1].scatter(x[trajStart:trajStop], y[trajStart:trajStop],  edgecolor=colorlist[iparam2], marker=trajMarker, linewidth=trajLine, s=trajSize) 
                
    # draw arm
    plot([0,elbowCenter[0]], [0, elbowCenter[1]],'k-',lw=4) #elbow
    plot([elbowCenter[0], varmCenter[0]], [elbowCenter[1], varmCenter[1]],'k-',lw=4)
    tight_layout()
    fig.savefig('gif/%s_traj_%d.png' % (simdatadir[5:], trajStop))

# save EMG data to matlab
def save2MatlabEMGs(simfolder):
    for file in os.listdir(simfolder): 
        if file.endswith("muscles.p"):
            with open(simfolder+'/'+file) as f:
                [musExcSeq, musActSeq, musForcesSeq]=pickle.load(f)
            #print [musExcSeq, musActSeq, musForcesSeq]
            print "Saving %s .mat" % (file)
            scipy.io.savemat(('%s/%s.mat'%(simfolder,file[0:-2])), \
            mdict={'musExcSeq': musExcSeq, 'musActSeq': musActSeq, 'musForcesSeq': musForcesSeq})    

# generate pdf figure of fitness evolution
def figureFitEvol(outfilestem):    
    with open('data/%s/fitevol'% (outfilestem), 'r') as f:
        [stat_gens, stat_avgfits, stat_worstfits, stat_bestfits] = pickle.load(f) 
        
        fontsiz=14
        genMax = 100;
        figure()        
        plot(stat_gens[0:genMax], stat_avgfits[0:genMax], 'b-')
        plot(stat_gens[0:genMax], stat_worstfits[0:genMax], 'r--')
        plot(stat_gens[0:genMax], stat_bestfits[0:genMax], 'g--')
        xlabel('Generation number',fontsize=fontsiz)
        ylabel('Fitness',fontsize=fontsiz)
        savefig("fitevol.pdf",format='pdf') 

# find best fitness solutions and generate data for analysis
def bestFitness(folder,islands=10, top=50, runsims=1, iseeds=1, wseeds=5, calculateWeights=0):  
    #folder='14dec18_evol'
    #islands=10
    #top=50
    timeout = 500
    wseedvals = [120456, 398115, 534031, 796321, 895199]
    wseedvals = wseedvals[0:wseeds]
    iseedvals = arange(1235,  1235+(17*100), 17)
    iseedvals = iseedvals[0:iseeds]
        
    print "loading data..."
    # load data from fitness evaluations
    dataFrom = 'fitness'  # 'individuals'
    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 \
    = analyse_funcs.loadData(folder,islands,dataFrom)
    
    print "calculating best solutions..."
    # combine data from all islands
    all_fits = all_gens = all_cands = all_params = all_islands = []
    for isl in fits_sort_isl: 
        all_fits = all_fits + isl  # get fits from all islands
        all_islands = all_islands + [fits_sort_isl.index(isl)] * len(isl) # create vector with island
    for isl in gens_sort_isl: 
        all_gens = all_gens + isl # get gens from all islands
    for isl in cands_sort_isl: 
        all_cands = all_cands + isl # get gens from all islands
    for isl in params_sort_isl: 
        all_params = all_params + isl # get params from all islands
    
    # unique params and sort fits 
    all_params, unique_indices = analyse_funcs.uniqueList(all_params) 
    all_fits, all_gens, all_cands, all_islands = zip(*[(all_fits[i],all_gens[i],all_cands[i],all_islands[i]) for i in unique_indices])
    sort_indices = np.argsort(all_fits) # sort fitness
    all_fits, all_gens, all_cands, all_params, all_islands  = zip(*[(all_fits[i],all_gens[i],all_cands[i],all_params[i],all_islands[i]) for i in sort_indices])

    if runsims: 
        print "running sims ..."
        # for each of the top solutions run sim to get data
        if iseeds==1:
            for i in range(top):
                # call runpbatchpbs to run sims saving muscle, weight and LFP data  
                command = 'python runbatchpbs_iseeds_16params.py %s_best data/%s_island_%s/gen_%s_cand_%s %d 0 -1' % (folder,folder,all_islands[i], all_gens[i], all_cands[i], iseeds)
                print command   
                os.system(command)  
        elif iseeds>1:
            i = top
            # call runpbatchpbs to run sims saving muscle, weight and LFP data  
            command = 'python runbatchpbs_iseeds_16params.py %s_best data/%s_island_%s/gen_%s_cand_%s %d 0 -1' % (folder,folder,all_islands[i], all_gens[i], all_cands[i], iseeds)
            print command   
            os.system(command) 
        
        print 'Sleeping for 10 mins to wait for results...'
        time.sleep(100)
        
    # for each of the top solutions run sim to get data
    if top > 1 and iseeds == 1:
        total_jobs = top
        converted = [None]*total_jobs
        num_iters = 0
        jobs_completed=0
        while jobs_completed < total_jobs:
            print str(num_iters)+" iterations, "+str(jobs_completed)+" / "+str(total_jobs)+" jobs completed"
            unfinished = [i for i, x in enumerate(converted) if x is None]
            print "unfinished:" + str(unfinished)
            for i in unfinished:
                # check if file exists
                subfolder = 'island_%s_gen_%s_cand_%s' % (all_islands[i], all_gens[i], all_cands[i])
                simdatadir = 'data/%s/%s' % (folder+'_best', subfolder)
                for iseed in iseedvals:         
                    for wseed in wseedvals:
                        try:
                            for itarget in range(4):
                                filename = '%s/target-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, itarget, iseed, wseed)    
                                with open(filename) as f:
                                    print 'Reading file: ' + f.name
                            save2Matlab_iseeds(simdatadir, 20, 1000) 
                            save2MatlabEMGs(simdatadir)
                            if calculateWeights:                        
                                popMuscles(simdatadir)   
                            converted[i] = 1
                            jobs_completed+=1
                        except:
                            #print 'Not found: '+filename
                            pass
            num_iters += 1
            if num_iters >= timeout: 
                print "max iterations reached without completing all jobs"
                for j in unfinished:
                    converted[j] = 2
                    jobs_completed += 1
            time.sleep(30)
    # for each of the iseeds -- NEEDS FIXING, SOMETIMES NOT ALL .MAT FILES ARE GENERATED
    elif top >= 1 and iseeds > 1:
        total_jobs = iseeds
        converted = [None]*total_jobs
        num_iters = 0
        jobs_completed=0
        while jobs_completed < total_jobs:
            print str(num_iters)+" iterations, "+str(jobs_completed)+" / "+str(total_jobs)+" jobs completed"
            unfinished = [i for i, x in enumerate(converted) if x is None]
            print "unfinished:" + str(unfinished)
            isol = top # top solution only
            # check if file exists
            subfolder = 'island_%s_gen_%s_cand_%s' % (all_islands[isol], all_gens[isol], all_cands[isol])
            simdatadir = 'data/%s/%s' % (folder+'_best', subfolder)
            for iseed,iseedval in enumerate(iseedvals):         
                for wseed in wseedvals:
                    try:
                        for itarget in range(4):
                            filename = '%s/target-%d_i-%d_w-%d_test-nqa.nqs' % (simdatadir, itarget, iseedval, wseed)    
                            with open(filename) as f:
                                print 'Reading file: ' + f.name  
                        converted[iseed] = 1
                        jobs_completed+=1
                    except:
                        #print 'Not found: '+filename
                        pass
            num_iters += 1
            if num_iters >= timeout: 
                print "max iterations reached without completing all jobs"
                for j in unfinished:
                    converted[j] = 2
                    jobs_completed += 1
            time.sleep(30)
        # after all iseeds have finished - convert all together
        save2Matlab_iseeds(simdatadir, 20, 1000) 
        save2MatlabEMGs(simdatadir)
        if calculateWeights:                        
            popMuscles(simdatadir) 
    else:
        print "Have to choose either top>1 and iseeds=1 or top=1 and iseeds>1"


# Script code (run always)
#


# Load default parameters and initialize the network.
# COMMENT OUT WHEN RUNNING SIM.PY !!!
#h.xopen("/usr/site/nrniv/simctrl/hoc/setup.hoc")
#h.xopen("/usr/site/nrniv/simctrl/hoc/nrnoc.hoc")

#h.load_file("syncode.hoc")
#h.load_file("decnqs.hoc")
#h.load_file("analysis.hoc")

#bestFitness('14dec18_evol')