fitness.py

import neuron
import math
import copy
import readconffile as rcf
import readexpfile as ref

filename3=''

neuron.h.load_file("nrngui.hoc")
[inputfilename,modfilename,parametersfilename,flagdata,flagcut,nrtraces,Vrestf,esynf,nrparamsfit,paramnr,paramname,paraminitval,paramsconstraints,nrdepnotfit,depnotfit,nrdepfit,depfit,seedinitvaluef]=rcf.readconffile()


def cuttrace( trace_number, sizeofsw ):
    #print "resizing trace: ", trace_number
    #trace_number=0
    #sizeofsw=10
    perccut=10
    flagcut=1
    [timevecprov,vec5]=ref.readexpfile(num=trace_number)
    timevec = []
    for i in range(len(vec5)):
        timevec.append(timevecprov[i])
                        
    # CUT TRACE
    vec5forsliding = copy.copy(vec5)
    
    derivativessliding = []
    idxderivativessliding = []
    for i in range(0,len(vec5)-sizeofsw,sizeofsw):
        vector1ms = []
        for j in range(i,i+sizeofsw+1):
            vector1ms.append(vec5forsliding[j])
        if (vector1ms.index(max(vector1ms))<vector1ms.index(min(vector1ms))):
            derivativessliding.append(max(vector1ms)-min(vector1ms))
            idxderivativessliding.append(i)
    a=0
    maxes = []
    while ((a<=len(derivativessliding)-1) and (len(maxes)<2)):
        if (derivativessliding[a]==max(derivativessliding)):
            if ((max(derivativessliding)/(max(vec5)-min(vec5)))*100>=perccut):
                maxes.append(idxderivativessliding[a])
            derivativessliding[a]=0
        a=a+1
    #print "length maxes: ", len(maxes)
    
    if (len(maxes)==1):
        maxes.append(len(vec5)-1)

    vec5forsliding=[]
    for i in range(maxes[0],maxes[1]+1):
            vec5forsliding.append(vec5[i])

    timevecaftersliding = []
    for i in range(maxes[0],maxes[1]+1):
        timevecaftersliding.append(timevec[i])

    timevec = []
    for i in range(0,len(timevecaftersliding)):
        timevec.append(timevecaftersliding[i]-timevecaftersliding[0])

    vec3 = []
    for i in range(0,len(vec5forsliding)):
        vec3.append(vec5forsliding[i])
    #print "trace ", trace_number , " cutted"
    
    #REMOVE HOLDING CURRENT
    vec5 = []
    if (vec3[0]>=max(vec3)):
        for i in range(0,len(vec3)):
            vec5.append(vec3[i]-vec3[0])
    else:
        for i in range(0,len(vec3)):
            vec5.append(vec3[i]-max(vec3))

    for i in range(0,len(vec5)):
        if (vec5[i]==min(vec5)):
            imin=i
            i=len(vec5)-1

    vecsin = []
    for i in range(0,imin):
        vecsin.append(vec5[i])
    ides=len(vec5)
    
    a=len(vecsin)+1
    while (a<=len(vec5)-1):
        if (vec5[a]>max(vecsin)):
            ides=a
            a=len(vec5)-1
        a=a+1

    vecfin = []
    for i in range(0,ides):
        vecfin.append(vec5[i])

    vecfinfin = []
    for i in range(0,ides):
        vecfinfin.append(vecfin[i]-max(vecfin))

    vec5 = []
    for i in range(0,ides):
        vec5.append(vecfinfin[i])

    timevectemp = []
    for i in range(0,ides):
        timevectemp.append(timevec[i])
    timevec = []
    for i in range(0,ides):
        timevec.append(timevectemp[i])
    #print "removed holding current "
    
    if (flagcut==1):
        timemin=vec5.index(min(vec5)) 
        flag=0
        i=timemin
        while (flag<1 and i<len(vec5)):
            if (vec5[i]>0.2*min(vec5)):
                flag=1
            i=i+1
        i=i-1
        if (i<len(vec5) and timevec[i]>4):
            #print("resizing trace %d to %g ms, %d\n", trace_number, timevec[i], i)
            vec5prov=[]
            for j in range(i):
                vec5prov.append(vec5[j])
            vec5=[]
            for pp in range(len(vec5prov)):
                vec5.append(vec5prov[pp])
            timevecprov=[]
            for j in range(i):
                timevecprov.append(timevec[j])
            timevec=[]
            for pp in range(len(timevecprov)):
                timevec.append(timevecprov[pp])
    vecfinfin = []
    for i in range(0,len(vec5)):
        vecfinfin.append(vec5[i])
    vec5 = []
    for i in range(0,len(vecfinfin)):
        vec5.append(vecfinfin[i])
    timevecfinfin = []
    for i in range(0,len(timevec)):
        timevecfinfin.append(timevec[i])
    timevec = []
    for i in range(0,len(timevecfinfin)):
        timevec.append(timevecfinfin[i]) 
    return [vec5, timevec, maxes[0]];
    

def finaltrace(trace_number=0):
    [timevecprov,vecprov]=ref.readexpfile(num=trace_number)
    exp_current = []
    exp_times = []
    i=0
    while vecprov[i]<1000:
        exp_current.append(vecprov[i])
        exp_times.append(timevecprov[i])
        i=i+1
        if i==len(vecprov): break
        
    vec5 = exp_current
    vec52 = iter(vec5)
    valm = min(vec52)    
    maxofsw = vec5.index(valm)
    #print "maxim sizeofsw ", maxofsw
    
    sizeofsw = int(maxofsw/2)
    #print "sizeofsw ", sizeofsw
    
    perccut=10
    [vec5,timevec,cutsin]=cuttrace(trace_number,sizeofsw);
    
    vec52=iter(vec5)
    valm=min(vec52)
    indvalm=vec5.index(valm)
    while ((indvalm<=20) and (sizeofsw<=maxofsw)):
        sizeofsw=sizeofsw+1
        [vec5,timevec,cutsin]=cuttrace(trace_number,sizeofsw);
        vec52=iter(vec5)
        valm=min(vec52)
        indvalm=vec5.index(valm)
    
    while ((((len(timevec)-indvalm)<=40)) and (sizeofsw<=maxofsw)):
        sizeofsw=sizeofsw+1
        [vec5,timevec,cutsin]=cuttrace(trace_number,sizeofsw);
        vec52=iter(vec5)
        valm=min(vec52)
        indvalm=vec5.index(valm)
    #print "final sizeofsw of trace ", trace_number, "sizeofsw ", sizeofsw
    return [sizeofsw,maxofsw,vec5,timevec, cutsin]

nefun=0
def migliore_eval( vec, timevec, vec5, trace_number=0):
    global nefun
    """Evaluation function of synaptic optimiser"""
    exp_current=vec5
    vecparams=[]
    for i in range(len(vec)):
        vecparams.append(math.exp(vec.x[i]))
    nefun +=1
    #print "nefun ", nefun
    if (nefun>=2000):
        neuron.h.stop_praxis()
    model_current = run_model(
    vecparams,
    time_trace=timevec)
        
    if model_current is None:
        model_error=1e6
    else:
        model_error=0
        for i in range(len(model_current)):
            model_error=model_error+(model_current[i]-exp_current[i])*(model_current[i]-exp_current[i])
        model_error=model_error/len(model_current)
    return model_error
   

def run_model(parameters, time_trace=None):
    """Run the model with the specified set of parameters"""
    import cellprop
    tstop = 100
    e_syn = esynf
    Vrest = Vrestf
    netstim = neuron.h.NetStims(0.5, sec=cellprop.soma)
    netstim.freqhz = 18.0
    netstim.q = 0.0
    netstim.prob = 2.0
    netstim.noise = 1.0
    netstim.number = 1.0
    
    vclamp = neuron.h.VClamp(0.5, sec=cellprop.soma)
    vclamp.dur[0] = tstop
    vclamp.amp[0] = Vrest
    
    with open(filename3) as ff:
        searchlines=ff.readlines()
    for kk, line in enumerate(searchlines):
        if "POINT_PROCESS" in line:
            break
    l2=line.split()
    synapse = neuron.h.__getattribute__(l2[1])(.5)
    synapse.verboseLevel = 0
    synapse.Use = 1.0
    synapse.u0 = 1.0
    synapse.e_GABAA = e_syn
    synapse.setRNG(cellprop.synapse_rng)

    netcon = neuron.h.NetCon(netstim, synapse )
    netcon.delay = 0.0
    netcon.threshold = 0.0

    vclamp_i = neuron.h.Vector()
    timevec = neuron.h.Vector()
    timevec.from_python(time_trace)
  
    #print "nrparamsfit ", nrparamsfit
    neuron.h('''nrparamsfit=0''')
    neuron.h.nrparamsfit=nrparamsfit
    neuron.h('''objref paramnamenrn[nrparamsfit]''')
    for i in range(nrparamsfit):
        neuron.h.paramnamenrn[i]=neuron.h.String()
        neuron.h.paramnamenrn[i].s=paramname[i]   
    neuron.h('''objref parametersnrn''')
    neuron.h('''parametersnrn =new Vector()''')
    neuron.h.parametersnrn.from_python(parameters)
    for i in range(nrparamsfit):
        #print i
        #cmd=paramname[i]+"="+str(parameters[i])
        #print cmd
        #exec cmd
        #print neuron.h.paramnamenrn[i].s
        #print neuron.h.parametersnrn.x[i]
        neuron.h('strdef cmdstr')
        neuron.h('a=0')
        neuron.h.a=i
        neuron.h.execute('sprint(cmdstr,"%s = %g", paramnamenrn[a].s, parametersnrn.x[a])')
        #print neuron.h.cmdstr
        #neuron.h(neuron.h.cmdstr)
        exec(neuron.h.cmdstr)
    #print "geph", synapse.geph
    for i in range(nrdepnotfit):
        cmd=depnotfit[i]
        exec(cmd)
    #print "nhalf", synapse.nhalf

    neuron.h.tstop = tstop
    
    exc=0
    for i in range(nrdepfit):
        cmd=depfit[i]
        if eval(cmd):
            exc=exc or 1
    
    paramnr=0
    for row in paramsconstraints:
        low=row[0]
        high=row[1]
        if (parameters[paramnr]<low or parameters[paramnr]>high):
            exc=exc or 1      
        paramnr=paramnr+1
    
    vclamp_i.record(synapse._ref_i, timevec)
    if (exc==1):
        return None
    else:
        try:
            neuron.h.run()
        except RuntimeError:
            return None
    del netcon
    vclamp_i.mul(1000.0)
    #print cellprop.soma.g_pas
    return vclamp_i.to_python()