fitting.py

import random
import csv
import math
from functools import partial
import time
import sys
import os
import glob
import shutil
import subprocess
#import matplotlib
#matplotlib.use('Agg')
#import matplotlib.pyplot as plt
import readconffile as rcf
import readexpfile as ref

global nefun
global seedinitvalue

def fitting(configuration_file, experimental_file, mod_file, all_traces, singletrace, demo, singletrace_number):  
    '''
        Accepts:
            configuration_file:
                type: text/plain
                description: text file, specifying parameters of the simulation.
            experimental_file:
                type: text/plain
                description: text file, specifying experimental traces.
            mod_file:
                type: video/mpeg
                description: mod file, specifying the model. 
            all_traces:
                type: bool
                description: If True, all traces are fitted 100 times [default=False].
            singletrace:
                type: bool
                description: If True, only one trace is fitted 100 times [default=False].
            demo:
                type: bool
                description: If True, one trace is fitted 10 times [default=False].
            singletrace_number:
                type: double
                description: number of the trace to be fitted [default=3].
        Returns:
            res: text/csv
    '''
    
    #subprocess.call(['nrnivmodl', ])
    
    import neuron
    start_time = time.time()

    singletrace_number = int(singletrace_number)
    rcf.filename=configuration_file
    ref.filename2=experimental_file

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

    seedinitvalue = seedinitvaluef

    neuron.h.use_mcell_ran4(1)
    neuron.h.mcell_ran4_init(seedinitvalue)

    synapse_rng = neuron.h.Random()
    synapse_rng.MCellRan4(12345)
    synapse_rng.uniform(0, 1)

    trace_number = 0

    randnum = neuron.h.Random()
    randnum.MCellRan4(12345)

    vecparamforfit = fixed(nrparamsfit,paraminitval)

    listofvecs = []
    for fitnr in range(100):
        vecparamforfit = fixed(nrparamsfit,paraminitval)
        vecparamforfit2 = vecparamforfit
        for i in range(nrparamsfit):
            vecparamforfit2[i] = vecparamforfit[i]+randnum.uniform(-vecparamforfit[i], vecparamforfit[i])
        listofvecs.append(vecparamforfit2)

    import fitness
    fitness.filename3=mod_file
    if demo=='True':
        twoargs = [(num, fitnr) for num in range(singletrace_number, singletrace_number+1) for fitnr in range(5)]
    if all_traces=='True':
        twoargs = [(num, fitnr) for num in range(1,nrtraces+1) for fitnr in range(100)]
    if singletrace=='True':
        twoargs = [(num, fitnr) for num in range(singletrace_number, singletrace_number+1) for fitnr in range(20)] 
    pc = neuron.h.ParallelContext()
    pc.runworker()
    info = runsim(twoargs,pc,seedinitvalue,listofvecs,nrparamsfit)
    nums = info[0]
    fitnrs = info[1]
    errorrs = info[2]
    soglias = info[3]
    minvals = info[4]
    vecparamfitteds = info[5]
    sizeofsws = info[6]
    maxofsws = info[7]
    vec5ss = info[8]
    timevecss = info[9]
    cutsinss = info[10]
    vec5realss = info[11]
    with open("test.csv", "a") as myfile2:
        for k in range(len(nums)):
            myfile2.write("%i\t" % nums[k])
            myfile2.write("%s\t" % fitnrs[k])
            myfile2.write("%s\t" % errorrs[k])
            vecps = vecparamfitteds[k]
            for numparams in range(nrparamsfit):
                myfile2.write("%s\t" % vecps[numparams])
            myfile2.write("%s\t" % soglias[k])
            myfile2.write("%s\n" % minvals[k])
    myfile2.close()
    if (singletrace or demo):
        vecerrorssingletrace=[]
        for k in range(len(nums)):
            vecerrorssingletrace.append(errorrs[k])
        if len(vecerrorssingletrace)>0:
            indexplot=[n for n,i in enumerate(vecerrorssingletrace) if i==min(vecerrorssingletrace) ][0]
            print("indexplot ", indexplot)
            print("error ", errorrs[indexplot])
            if (len(nums)==0):
                for k in range(len(timevecss)):
                    vec5pss = vec5ss[k]
                    timevecpss = timevecss[k]
                #plt.plot(timevecpss,vec5pss,'g')
            else:
                sizeofsws2 = sizeofsws[indexplot]
                maxofsws2 = maxofsws[indexplot]
                vec5s2 = vec5ss[indexplot]
                timevecs2 =timevecss[indexplot]
                cutsinss2 = cutsinss[indexplot]
                vec5realss2 = vec5realss[indexplot]
                timevecreal = []
                vec5s2real = []
                for k in range(cutsinss2,cutsinss2+len(vec5realss2)):
                    timevecreal.append(timevecs2[k])
                    vec5s2real.append(vec5s2[k])
                vec5s2final = []
                for k in range(len(vec5s2)):
                    vec5s2final.append(vec5s2[k]-max(vec5s2real))      
                if (sizeofsws2<maxofsws2):
                    for ss in paramname:
                        if 'netstim.start' in ss:
                            startipoz=paramname.index(ss)
                    starti=vecps[startipoz]
                    vecps[startipoz]=starti+timevecs2[cutsinss2]
                    model_current = fitness.run_model(vecps,time_trace=timevecs2);
                    errorverif=0
                    for k in range(cutsinss2,cutsinss2+len(vec5realss2)):
                        errorverif=errorverif+(model_current[k]-vec5s2final[k])*(model_current[k]-vec5s2final[k])
                    errorverif=errorverif/len(vec5s2real)
                    #plt.plot(timevecs2,vec5s2final,'g',timevecs2,model_current,'b')
            imagename='tracefit.png'
            #plt.savefig(imagename, dpi=300)      
    print("done")
    pc.done()
    neuron.h.quit()
    
def fixed(nrparamsfit,paraminitval):
    """Return fixed initialisation"""
    vecinit= []
    for i in range(nrparamsfit):
        vecinit.append(paraminitval[i])
    return vecinit

def optim(twoargss,seedinitvalue,listofvecs,nrparamsfit):
    import neuron
    import fitness
    (num,fitnr)=twoargss
    #print num
    #print fitnr
    #global nefun
    #print "fitnr ", fitnr
    #num=3
    #print "Start Process", os.getpid(), "with args", num, fitnr
    #print "args", num, fitnr
    neuron.h.attr_praxis(1e-4, .5, 0)
    #print seedinitvalue
    neuron.h.attr_praxis(seedinitvalue)
    vecparamforfit2=listofvecs[fitnr]
    #print "initial values ", vecparamforfit2
    vecparamforfit3=[]
    for i in range(int(nrparamsfit)):
        vecparamforfit3.append(math.log(vecparamforfit2[i]))
    vec=neuron.h.Vector(nrparamsfit)
    for i in range(int(nrparamsfit)):
        vec.x[i]=vecparamforfit3[i]
    fitness.nefun=0
    fitness.nquad=0
    neuron.h.stop_praxis(0)
    [sizeofsw,maxofsw,vec5,timevec,cutsin]=fitness.finaltrace(trace_number=num);
    flagsw=0
    soglia=(0.1*min(vec5))**2
    #print "soglia", soglia
    minval=min(vec5)
    [timevecprov,vecallprov]=ref.readexpfile(num=num)
    vecall = []
    for i in range(len(vecallprov)):
        vecall.append(vecallprov[i])
    timevecall = []
    for i in range(len(vecall)):
        timevecall.append(timevecprov[i]-timevecprov[0])
    if (sizeofsw<=maxofsw):
        flagsw=1
        error = neuron.h.fit_praxis(partial(fitness.migliore_eval, trace_number=num, timevec=timevec, vec5=vec5),vec)
        vecparamfitted=[]
        for i in range(int(nrparamsfit)):
            vecparamfitted.append(math.exp(vec.x[i]))
        return (num, fitnr, error, soglia, minval, vecparamfitted, fitness.nefun, sizeofsw, maxofsw, vecall, timevecall, cutsin, vec5)
    else:
        return (num, fitnr, 1e6, soglia, minval, vecparamforfit3, 0, sizeofsw, maxofsw, vecall, timevecall, cutsin, vec5)

def runsim(twoargs,pc,seedinitvalue,listofvecs,nrparamsfit):
    for items in twoargs:
        pc.submit(optim,items,seedinitvalue,listofvecs,nrparamsfit)
    num2s=[]
    fitnr2s=[]
    error2s=[]
    minval2s=[]
    vecparamfitted2s=[]
    soglia2s=[]
    sizeofswrs=[]
    maxofswrs=[]
    vec5rs=[]
    timevecrs=[]
    cutsinrs=[]
    vec5reals=[]
    while (pc.working()):
        num2, fitnr2, error2, soglia2, minval2, vecparamfitted2, nefun2, sizeofswr, maxofswr, vec5r, timevecr, cutsinr, vec5real = pc.pyret()
        print(num2,fitnr2,error2,soglia2,minval2,nefun2)
        if (error2<soglia2):
            num2s.append(num2)
            fitnr2s.append(fitnr2)
            error2s.append(error2)
            soglia2s.append(soglia2)
            minval2s.append(minval2)
        vecparamfitted2s.append(vecparamfitted2)
        sizeofswrs.append(sizeofswr)
        maxofswrs.append(maxofswr)
        vec5rs.append(vec5r)
        timevecrs.append(timevecr)
        cutsinrs.append(cutsinr)
        vec5reals.append(vec5real)
    return (num2s, fitnr2s, error2s, soglia2s, minval2s, vecparamfitted2s, sizeofswrs, maxofswrs, vec5rs, timevecrs, cutsinrs, vec5reals)

if __name__ == "__main__":
    fitting(sys.argv[1],sys.argv[2],sys.argv[3],sys.argv[4],sys.argv[5],sys.argv[6],sys.argv[7])