#
# MSN model
#
# Robert Lindroos <robert.lindroos at ki.se>
# 
# Original version by 
# Alexander Kozlov <akozlov at kth.se>
# Kai Du <kai.du at ki.se>
#


from __future__ import print_function, division
import mpi4py as MPI
import sys
import json
import numpy as np
from neuron import h
from math import exp
import pickle
import os
import glob
import matplotlib.pyplot as plt



mod        = "./mod/"
defparams  = "./params-rob.json"
morphology = "./morphology/"



#h.nrn_load_dll('/pdc/vol/neuron/7.4-py27/x86_64/.libs/libnrnmech.so')
h.load_file('stdlib.hoc')
h.load_file('import3d.hoc')



# global result dict
RES = {}


# Distributions:
'''
T-type Ca: (1.0/(1+np.exp((v-70)/-4.5)))
naf (den): (0.1 + 0.9/(1 + np.exp((v-60.0)/10.0)))

'''

#                        1,      0.8, 2.5,60,-25
#calculate_distribution(d3, dist, a4, a5, a6, a7, g8)
def calculate_distribution(d3, dist, a4, a5, a6, a7, g8):
    '''
    Used for setting the maximal conductance of a segment.
    Scales the maximal conductance based on somatic distance and distribution type.
    
    Parameters:
    d3   = distribution type:
         0 linear, 
         1 sigmoidal, 
         2 exponential
         3 step function
    dist = somatic distance of segment
    a4-7 = distribution parameters 
    g8   = maximal conductance of channel
    
    '''
    
    if   d3 == 0: 
        value = a4 + a5*dist
    elif d3 == 1: 
        value = a4 + a5/(1 + exp((dist-a6)/a7) )
    elif d3 == 2: 
        value = a4 + a5*exp((dist-a6)/a7)
    elif d3 == 3:
        if (dist > a6) and (dist < a7):
            value = a4
        else:
            value = a5
            
    if value < 0:
        value = 0
        
    value = value*g8
    return value 


  
    
def alpha(tstart, gmax, tau):
    # calc and returns a "magnitude" using an alpha function -> used for [DA] in cascade
    
    v = 1 - (h.t - tstart) / tau
    e = exp(v)
    mag = gmax * (h.t - tstart) / tau * e
    
    return mag
    
    
    
def save_vector(t, v, outfile):
    
    with open(outfile, "w") as out:
        for time, y in zip(t, v):
            out.write("%g %g\n" % (time, y))
            
            
            
def calc_rand_Modulation(mod_list, range_list=False, distribution='uniform'):
    
    '''
    uses numpy to draws random modulation factors in range [0,2] from a given distribution:
        uniform          (same probability for all values)
        centered         (highest probability for values at center of interval)
        inverse centered (higest probability for values close to limits)
    for each channel in mod_list.
    The factors can also be mapped to an arbitrary interval. 
    This is done if a range_list is given.
    
    mod_list     = list of channels to be modulated
    range_list   = list of [min, max] values to be used in modulation
    distribution = distribution to draw factors from
    '''
    
    mod_factors = []
    
    A=0
    B=2
    
    for i,channel in enumerate(mod_list):
        
        if distribution=='centered':
            factor = 1.0 + ( np.random.uniform() - np.random.uniform() )
        elif distribution=='inv_centered':
            factor = 1.0 + ( np.random.uniform() - np.random.uniform() )
            if factor <= 1:
                factor = factor + 1.0
            else:
                factor = factor - 1.0
        elif distribution=='uniform':
            factor = 2.0 * np.random.uniform()
        else:
            print('Error in MF distribution--line ~121')
            eegsjsd
        
        
        if range_list:
            
            a       = range_list[i][0]
            b       = range_list[i][1]
            
            factor = (b-a) / (B-A) * (factor-A) + a
       
        mod_factors.append(factor)
        
    return mod_factors 
    
    


def save_obj(obj, name ):
    with open('save/'+ name + '.pkl', 'wb') as f:
        pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)

def load_obj(name ):
    with open('save/' + name + '.pkl', 'rb') as f:
        return pickle.load(f) 
        
        
        
def getSpikedata_x_y(x,y):
    
    ''' 
    There's probably a Neuron function for this--use instead?
    
    getSpikedata_x_y(x,y) -> return count
    
    Extracts and returns the number of spikes from spike trace data.
    
    # arguments
    x = time vector
    y = vm vector 
    
    # returns
    count = number of spikes in trace (int)
    
    # extraction algorithm
    -threshold y and store list containing index for all points larger than 0 V 
    -sorts out and counts the index that are the first one(s) crossing the threshold, i.e. 
        the first index of each spike. This is done by looping over all index and check if 
        the index is equal to the previous index + 1. If not it is the first index of a 
        spike.
        
        If no point is above threshold in the trace the function returns 0.
        
    '''
    
    count = 0
    spikes = []
    
    # pick out index for all points above zero potential for potential trace
    spikeData = [i for i,v in enumerate(y) if v > 0]

    # if no point above 0
    if len(spikeData) == 0:
        
        return spikes
	
    else:
        # pick first point of each individual transient (spike)...
        for j in range(0, len(spikeData)-1):
            if j==0:
                
                count += 1
                spikes.append(x[spikeData[j]])
		
            # ...by checking above stated criteria
            elif not spikeData[j] == spikeData[j-1]+1:
                count += 1
                spikes.append(x[spikeData[j]])
            
    return spikes 
    

        
        

# ======================= the MSN class ==================================================

class MSN:
    def __init__(self, params=defparams, factors=None):
        Import = h.Import3d_SWC_read()
        Import.input(morphology + 'latest_WT-P270-20-14ak.swc')
        imprt = h.Import3d_GUI(Import, 0)
        imprt.instantiate(None)
        h.define_shape()
        # h.cao0_ca_ion = 2  # default in nrn
        h.celsius = 35
        self._create_sectionlists()
        self._set_nsegs()
        self.v_init = -80
        for sec in self.allseclist:
            sec.Ra = 150
            sec.cm = 1.0
            sec.insert('pas')
            #sec.g_pas = 1e-5 # set using json file
            sec.e_pas = -70 # -73
        for sec in self.somalist:
            sec.insert('naf')
            sec.insert('kaf')
            sec.insert('kas')
            sec.insert('kdr')
            sec.insert('kir')
            sec.ena = 50
            sec.ek = -85 # -90
            sec.insert('cal12')
            sec.insert('cal13')
            sec.insert('car')
            sec.insert('cadyn')
            sec.insert('caldyn')
            sec.insert('sk')
            sec.insert('bk')
            sec.insert('can')
	    #sec.kb_cadyn = 200.
        for sec in self.axonlist:
            sec.insert('naf')
            #sec.insert('kaf')
            sec.insert('kas')
            #sec.insert('kdr')
            #sec.insert('kir')
            sec.ena = 50
            sec.ek = -85 # -90
        for sec in self.dendlist:
            sec.insert('naf')
            sec.insert('kaf')
            sec.insert('kas')
            sec.insert('kdr')
            sec.insert('kir')
            sec.ena = 50
            sec.ek = -85 # -90
            sec.insert('cal12')
            sec.insert('cal13')
            sec.insert('car')
            sec.insert('cadyn')
            sec.insert('caldyn')
            sec.insert('sk')
            sec.insert('bk')
            sec.insert('cat32')
            sec.insert('cat33')

        with open(params) as file:
            par = json.load(file)

        self.distribute_channels("soma", "g_pas", 0, 1, 0, 0, 0, float(par['g_pas_all']['Value']))
        self.distribute_channels("axon", "g_pas", 0, 1, 0, 0, 0, float(par['g_pas_all']['Value']))
        self.distribute_channels("dend", "g_pas", 0, 1, 0, 0, 0, float(par['g_pas_all']['Value']))

        self.distribute_channels("soma", "gbar_naf", 0, 1, 0, 0, 0, float(par['gbar_naf_somatic']['Value']),factors=factors)
        self.distribute_channels("soma", "gbar_kaf", 0, 1, 0, 0, 0, float(par['gbar_kaf_somatic']['Value']))
        self.distribute_channels("soma", "gbar_kas", 0, 1, 0, 0, 0, float(par['gbar_kas_somatic']['Value']))
        self.distribute_channels("soma", "gbar_kdr", 0, 1, 0, 0, 0, float(par['gbar_kdr_somatic']['Value']))
        self.distribute_channels("soma", "gbar_kir", 0, 1, 0, 0, 0, float(par['gbar_kir_somatic']['Value']))
        self.distribute_channels("soma", "gbar_sk", 0, 1, 0, 0, 0, float(par['gbar_sk_somatic']['Value']))
        self.distribute_channels("soma", "gbar_bk", 0, 1, 0, 0, 0, float(par['gbar_bk_somatic']['Value']))
        
        #self.distribute_channels("axon", "gbar_naf", 0, 1, 0, 0, 0, float(par['gbar_naf_somatic']['Value']),factors=factors)
        self.distribute_channels("axon", "gbar_naf", 3, 1, 1.1, 30, 500, float(par['gbar_naf_axonal']['Value']),factors=factors)
        #self.distribute_channels("axon", "gbar_naf", 3, 1, 1.1, 20, 500, float(par['gbar_naf_axonal']['Value']))
        #self.distribute_channels("dend", "gbar_naf", 1, 1,  1.2, 30, -5, float(par['gbar_naf_axonal']['Value']))
        self.distribute_channels("axon", "gbar_kas", 0, 1, 0, 0, 0, float(par['gbar_kas_axonal']['Value']))
        
        self.distribute_channels("dend", "gbar_naf", 1, 0.1, 0.9, 60.0, 10.0, float(par['gbar_naf_basal']['Value']),factors=factors)
        self.distribute_channels("dend", "gbar_kaf", 1, 1,  0.5, 120, -30, float(par['gbar_kaf_basal']['Value']))
        #self.distribute_channels("dend", "gbar_naf", 0, 1, -0.0072, 0, 0, float(par['gbar_naf_basal']['Value']))
        #self.distribute_channels("dend", "gbar_kaf", 0, 1,  0.0167, 0, 0, float(par['gbar_kaf_basal']['Value']))
        self.distribute_channels("dend", "gbar_kas", 2, 1, 9.0, 0.0, -5.0, float(par['gbar_kas_basal']['Value']))
        self.distribute_channels("dend", "gbar_kdr", 0, 1, 0, 0, 0, float(par['gbar_kdr_basal']['Value']))
        self.distribute_channels("dend", "gbar_kir", 0, 1, 0, 0, 0, float(par['gbar_kir_basal']['Value']))
        self.distribute_channels("dend", "gbar_sk", 0, 1, 0, 0, 0, float(par['gbar_sk_basal']['Value']))
        self.distribute_channels("dend", "gbar_bk", 0, 1, 0, 0, 0, float(par['gbar_bk_basal']['Value']))

        self.distribute_channels("soma", "pbar_cal12", 0, 1, 0, 0, 0, 1e-5)
        self.distribute_channels("soma", "pbar_cal13", 0, 1, 0, 0, 0, 1e-6)
        self.distribute_channels("soma", "pbar_car", 0, 1, 0, 0, 0, 1e-4)
        self.distribute_channels("soma", "pbar_can", 0, 1, 0, 0, 0, 3e-5)
        #self.distribute_channels("soma", "kb_cadyn", 0, 1, 0, 0, 0, 200.0)
        self.distribute_channels("dend", "pbar_cal12", 0, 1, 0, 0, 0, 1e-5)
        self.distribute_channels("dend", "pbar_cal13", 0, 1, 0, 0, 0, 1e-6)
        self.distribute_channels("dend", "pbar_car", 0, 1, 0, 0, 0, 1e-4)
        self.distribute_channels("dend", "pbar_cat32", 1, 0, 1.0, 70.0, -4.5, 1e-7)
        self.distribute_channels("dend", "pbar_cat33", 1, 0, 1.0, 70.0, -4.5, 1e-8)

    def _create_sectionlists(self):
        self.allsecnames = []
        self.allseclist = h.SectionList()
        for sec in h.allsec():
            self.allsecnames.append(sec.name())
            self.allseclist.append(sec=sec)
        self.nsomasec = 0
        self.somalist = h.SectionList()
        for sec in h.allsec():
            if sec.name().find('soma') >= 0:
                self.somalist.append(sec=sec)
                if self.nsomasec == 0:
                    self.soma = sec
                self.nsomasec += 1
        self.axonlist = h.SectionList()
        for sec in h.allsec():
            if sec.name().find('axon') >= 0:
                self.axonlist.append(sec=sec)
        self.dendlist = h.SectionList()
        for sec in h.allsec():
            if sec.name().find('dend') >= 0:
                self.dendlist.append(sec=sec)

    def _set_nsegs(self):
        for sec in self.allseclist:
            sec.nseg = 2*int(sec.L/40.0)+1
        for sec in self.axonlist:
            sec.nseg = 2  # two segments in axon initial segment

    def _max_dist(self, axon_excluding=True):
        h.distance(sec=self.soma)
        dmax = 0
        for sec in self.allseclist:
	        if axon_excluding and sec.name().find('axon') == 0: continue
                dmax = max(dmax, h.distance(1, sec=sec))
        return dmax

    def distribute_channels(self, as1, as2, d3, a4, a5, a6, a7, g8, factors=None):
        h.distance(sec=self.soma)
        dmax = self._max_dist()
        for sec in self.allseclist:
            if sec.name().find(as1) >= 0:
                for seg in sec:
                    dist = h.distance(seg.x, sec=sec)
                    val = calculate_distribution(d3, dist, a4, a5, a6, a7, g8)
                    cmd = 'seg.%s = %g' % (as2, val)
                    exec(cmd)
                    
                    #names = ['mVhalf_naf', 'hVhalf_naf', 'mSlope_naf', 'hSlope_naf']
                    names  = ['taum_naf', 'tauh_naf', 'taun_naf']
                    
                    if factors:
                        for i,factor in enumerate(factors):
                            
                            print('updating factors! ', factors, factor, names[i])
                            
                            # add minus here if running with slope and half activation; cmd = 'seg.%s = %g' % (names[i], -factor)
                            cmd = 'seg.%s = %g' % (names[i], factor)
                            exec(cmd)
                            
 
                  
   
#=========================================================================================



def main(par="./params-msn.json", \
                            sim='vm',       \
                            amp=0.265,      \
                            run=None,       \
                            modulation=1,   \
                            simDur=7000,    \
                            stimDur=900,    \
                            factors=None,   \
                            section=None,   \
                            randMod=None,   \
                            chan2mod=['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can'] ): 
    
    
    
    
    # initiate cell
    cell = MSN(params=par, factors=factors)
    
    # set cascade ---- move to MSN def?
    casc = h.D1_reduced_cascade2_0(0.5, sec=cell.soma) # other cascades also possible...
    
    pointer = casc._ref_Target1p       #totalActivePKA    (if full cascade used)
    
    # set edge of soma as reference for distance 
    h.distance(1, sec=h.soma[0])
    
    # set current injection
    stim = h.IClamp(0.5, sec=cell.soma)
    stim.amp = amp  
    stim.delay = 100
    stim.dur = stimDur            # 2ms 2nA to elicit single AP, following Day et al 2008 Ca dyn    
    
    # record vectors
    tm = h.Vector()
    tm.record(h._ref_t)
    vm = h.Vector()
    vm.record(cell.soma(0.5)._ref_v)
    
    pka = h.Vector()
    pka.record(pointer)
    
    # peak n dipp parameters
    da_peak   = 1500   # concentration [nM]
    da_tstart = 500    # stimulation time [ms]
    da_tau    = 500    # time constant [ms]
    
    
    tstop = simDur               # [ms]
    
    
    # all channels to modulate
    mod_list = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can' ]
    
    
    not2mod = [] #['kaf']
    
    
    # find channels that should not be modulated
    for chan in mod_list:
        
        if chan not in chan2mod:
            
            not2mod.append(chan)
    
    
    # calc modulation factors--------------------------------------------------------------
    base_mod    = casc.init_Target1p
    
    # for random modulation: modValues = np.arange(0.1, 2.0, 0.1) 
    if randMod == 1:
    
        if amp == 0.32:
        
            mod_fact = calc_rand_Modulation(mod_list, range_list=[[0.60,0.80],    \
                                                                  [0.65,0.85],  \
                                                                  [0.75,0.85],  \
                                                                  [0.85,1.25],  \
                                                                  [1.0,2.0],    \
                                                                  [1.0,2.0],    \
                                                                  [0.0,1.0]],
                                                                  distribution='uniform'  )
            
        else:
            
            mod_fact = RES[run]['factors']
    
    else:
        mod_fact = [ 0.8, 0.8, 0.8, 1.25, 2.0, 2.0, 0.5  ]
           
    
    factors     = []
    for i,mech in enumerate(mod_list):
        
        # normalization of factors to substrate range. Only for dynamical mod?
        factor  = (mod_fact[i] - 1) / (2317.1 - base_mod) #2317.1
        
        factors.append(factor)
        
        #print(mech, mod_fact[i], factor) # --------------------------------------------------------
            
    
    # set pointers 
    for sec in h.allsec():
        
        for seg in sec:
            
            # naf and kas is in all sections
            h.setpointer(pointer, 'pka', seg.kas )
            h.setpointer(pointer, 'pka', seg.naf )
            
            if sec.name().find('axon') < 0:    
                
                # these channels are not in the axon section
                
                h.setpointer(pointer, 'pka', seg.kaf )
                h.setpointer(pointer, 'pka', seg.cal12 )
                h.setpointer(pointer, 'pka', seg.cal13 )
                h.setpointer(pointer, 'pka', seg.kir )
                #h.setpointer(pointerc, 'pka', seg.car )
                
                if sec.name().find('soma') >= 0:
                    
                    # can is only distributed to the soma section
                    h.setpointer(pointer, 'pka', seg.can )
            
    
    # dynamical modulation ---------------------------------------------------------------
    if sim == 'modulation':
        
        print('inne ', sim)
        
        for sec in h.allsec():
            
            for seg in sec:
                
                for mech in seg:
                    
                    # if this first statement is active the axon will not be modulated
                    '''if sec.name().find('axon') >= 0     \
                            and mech.name() in mod_list:
                            
                        mech.factor = 0.0
                        print(sec.name(), seg.x, mech.name() )'''
                        
                    if mech.name() in not2mod:
                        
                        mech.factor = 0.0
                        print(mech.name(), 'and channel:', not2mod, mech.factor, sec.name())
                        
                    elif mech.name() in mod_list:
                    
                        mech.base       = base_mod
                        index           = mod_list.index( mech.name() )
                        mech.factor     = factors[index]
                    
                    
                        
                    
    
    # static modulation
    elif sim == 'directMod':
        
        #print('inne ', sim)
    
        for sec in h.allsec():
            
            for seg in sec:
                
                for mech in seg:
                
                    if mech.name() in mod_list: 
                        
                        if sec.name().find('axon') < 0: # if 0: no axon modulated; if 10: all sections modulated
                        
                            factor = mod_fact[mod_list.index(mech.name() )]
                            
                            if mech.name() in not2mod:
                                mech.factor = 0.0
                            elif mech.name()[0] == 'c':
                                pbar = mech.pbar
                                mech.pbar = pbar * factor
                                #print(''.join(['setting pbar ', mech.name(), ' ', str(factor) ]) )
                            else:
                                gbar = mech.gbar
                                mech.gbar = gbar * factor
                                #if seg.x < 0.2:
                                #print(''.join(['setting gbar ', mech.name(), ' ', str(factor), ' ', str(gbar), ' ', str(factor*gbar) ]) )
                        
                        else:
                            
                            print(sec.name(), seg.x, sec.name().find('axon'))
    
    
    
    # set ampa and nmda epsp's--what compartments to use?
    elif sim == 'plateau':
        
        print('inne ', sim)
        
        dend_name = 'dend[' + str(int(section)) + ']'
        
        for sec in h.allsec():
            
            if sec.name() == dend_name:
                
                x = 0.5
                
                ampa = h.ampa(x, sec=sec)
                ampa.onset = 100
                ampa.gmax  = 5e-3
                #h.setpointer(pointer, 'pka', seg.kaf )
                
                nmda = h.nmda(x, sec=sec)
                nmda.onset = 100
                nmda.gmax  = 10e-2
                #h.setpointer(pointer, 'pka', seg.kaf )
                
                vmL = h.Vector()
                vmL.record(sec(x)._ref_v)
                
                d2soma = int(h.distance(x, sec=sec))
                
                #print(sec.name(), h.distance(seg.x, sec=sec))
    
            
    
    # record Ca dynamics!
    elif sim == 'ca':
        
        print('inne ', sim)
        
        for i,sec in enumerate(h.allsec()):
            
            if sec.name().find('axon') < 0: # don't record in axon
            
                for j,seg in enumerate(sec):
                    
                    sName = sec.name().split('[')[0]
                    
                    # N, P/Q, R Ca pool
                    cmd = 'ca_%s%s_%s = h.Vector()' % (sName, str(i), str(j))
                    exec(cmd)
                    
                    cmd = 'ca_%s%s_%s.record(seg._ref_cai)' % (sName, str(i), str(j))
                    exec(cmd)   
                    
                    # the L-type Ca
                    cmd = 'cal_%s%s_%s = h.Vector()' % (sName, str(i), str(j))
                    exec(cmd)
                    
                    cmd = 'cal_%s%s_%s.record(seg._ref_cali)' % (sName, str(i), str(j))
                    exec(cmd)   
                    
                    # uncomment here if testing kaf blocking effect on bAP
                    #gbar = seg.kaf.gbar
                    #seg.kaf.gbar = 0.8 * gbar
    
    
              
    # solver------------------------------------------------------------------------------            
    cvode = h.CVode()
    
    h.finitialize(cell.v_init)
    
    # run simulation
    while h.t < tstop:
    
        if modulation == 1:
        
            if h.t > da_tstart: 
                
                # set DA and ACh values (using alpha function)
                casc.DA = alpha(da_tstart, da_peak, da_tau) 
                #casc.ACh = ach_base - alpha(ach_tstart, ach_base, ach_tau)
                
        h.fadvance()
        
    
    
    # save output
    if sim in ['vm', 'directMod', 'modulation']:
        
        '''s = ''
        for chan in not2mod:
            s = chan + s'''
        
        save_vector(tm, vm, ''.join(['./vm_', sim, '_', str(int(amp*1e3)), '.out']) )
        '''
        spikes      = getSpikedata_x_y(tm,vm) 
        amp         = int(amp*1e3) 
        
        
        if amp == 320:
            
            RES[run]                = {}
            RES[run]['factors']     = mod_fact
            
            if run == 0:
            
                RES['channels'] = mod_list
        
        RES[run][amp]   = {'spikes': spikes}'''
        
        
    elif sim == 'plateau':
    
        save_vector(tm, vm, ''.join(['./vm_', sim, str(d2soma), '_dend', str(int(section)), '.out']) )
        save_vector(tm, vmL, ''.join(['./vmL_', sim, str(d2soma), '_dend', str(int(section)), '.out']) )
    
    
    elif sim == 'ca':
        
        # vm
        save_vector(tm, vm, ''.join(['./vm_', sim, '_', str(int(amp*1e3)), '.out']) )
        
        # Ca
        for i,sec in enumerate(h.allsec()):
        
            if sec.name().find('axon') < 0:
            
                for j,seg in enumerate(sec):
                    
                    sName   = sec.name().split('[')[0]
                    
                    vName   = 'ca_%s%s_%s' %  (sName, str(i), str(j)    )
                    v2Name   = 'cal_%s%s_%s' %  (sName, str(i), str(j)    )
                    fName   = 'Ca/Org/ca_%s_%s.out'  %  (str(int(np.round(h.distance(seg.x)))), vName )
                    
                    cmd     = 'save_vector(tm, np.add(%s, %s), %s)' % (vName, v2Name, 'fName' ) 
                    
                    exec(cmd)
                



# if run from terminal...   ===============================================================
if __name__ == "__main__":
    
    sys.argv = ['a', 'ca']
    print('starting sim')
    
    arguments = None
    
    # define modulation vectors--used to investigate m and h gate differences, i.e. shift gates
    '''
    kaf_m_vhalf = [8, 9, 10, 11, 12]
    kaf_h_vhalf = [73.6, 74.6, 75.6, 76.6, 77.6]
    kaf_m_slope = [26, 26.5, 27, 27.5, 28]
    kaf_h_slope = [-8, -8.5, -9, -9.5, -10]
    
    # 
    naf_m_vhalf = [23, 24, 25, 26, 27]
    naf_h_vhalf = [60, 61, 62, 63, 64]
    naf_m_slope = [8.2, 8.7, 9.2, 9.7, 10.2]
    naf_h_slope = [-5, -5.5, -6, -6.5, -7]
    '''
    
    # tau sim 1
    #taum = [0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15]
    #tauh = [0.14, 0.19, 0.24, 0.29, 0.34, 0.39, 0.44, 0.49, 0.54]
    
    # tau sim 2 vhalf
    #taum = [52, 54, 56, 58, 60, 62, 64]
    #tauh = [26, 28, 30, 32, 34, 36, 38]
    
    # tau sim 3 slope
    #taum = [4, 6, 8, 10, 12]
    #taun = [-23, -28, -33, -38, -43]
    #tauh = [2.5, 3.5, 4.5, 6.5, 7.5]
    
    factors = [1] #list(itertools.product(taum, tauh, taun))
    
    if sys.argv[1] == 'ca':     # single run
    
        current = 2000
        main( par="./params-rob.json",          \
                    amp=current*1e-3,           \
                    modulation=0,               \
                    simDur=200,                 \
                    stimDur=2,                  \
                    sim='ca'                    )
                    # dur increased from 2 to 25
                    # amp decreased from 2000 to 1000
                                                    
                                                    
                                                    
                                                    
    elif sys.argv[1] in ['vm', 'directMod']:  #for current in currents:
    
        mod_list = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can']
        
        randMod = 0
        
        currents = np.arange(310,335,10)
        
        RES = {}
        
        #for n in range(400):
        
        files = glob.glob('vm*')
        
        plot = 1
        if plot:
            for f in files:
                
                plt.plot(*np.loadtxt(f, unpack=True))
                
            plt.show()
        else:
        
            for current in currents:
                main( par="./params-rob.json",      \
                        amp=current*1e-3,           \
                        run=1,                      \
                        modulation=0,               \
                        simDur=1000,                \
                        stimDur=900,                \
                        sim=sys.argv[1],            \
                        randMod=randMod,            \
                        chan2mod=mod_list           )
                
                
                if randMod == 1 and n % 5 == 0:
                    
                    print('inside save loop', n)
                
                    if n == 0:
                        
                        mod_fact = RES[0]['factors']
                    
                        ID = ''
                        
                        for i,mech in enumerate(mod_list):
                            
                            ID = ID + mech + str( int(mod_fact[i]*100) )
                
                    save_obj(RES, ''.join(['./StatUnikaf-noAxon-', ID]) )
                
                
            
                        
                        
                        
                        
                        
                        
    elif sys.argv[1] == 'plateau':  #for sec in sections (dendritic):
        
        sections = np.arange(1,57) #[4, 8, 18, 40, 57]
        for sec in sections:
            main( par="./params-rob.json",      \
                    amp=0,                      \
                    modulation=0,               \
                    simDur=500,                 \
                    sim=sys.argv[1],            \
                    stimDur=800,                \
                    section=sec                 )
                        
                        
                        
                        
    elif sys.argv[1] == 'modulation':     # multiple currents/modulations?
        
        mod_list = ['naf', 'kas', 'kaf', 'kir', 'cal12', 'cal13', 'can']
        
        currents = np.arange(300,355,10)
        
        randMod = 1
        
        for n in range(80):
            
            RES = {}
        
            for i,current in enumerate(currents):
                 
                main( par="./params-rob.json",      \
                        amp=current*1e-3,           \
                        run=i,                      \
                        modulation=1,               \
                        simDur=3000,                \
                        stimDur=3000,               \
                        sim=sys.argv[1],            \
                        randMod=randMod,            \
                        chan2mod=mod_list           )
                                
            '''
                 chan2mod=[item for item in mod_list if not item == chan ] \

            '''
        
            if randMod == 1:
            
                mod_fact = RES['factors']
                
                ID = ''
                
                for i,mech in enumerate(mod_list):
                    
                    ID = ID + mech + str( int(mod_fact[i]*100) )
            
                save_obj(RES, ''.join(['RMD-', ID]) )