diff --git a/moose-core/diffusion/Dsolve.cpp b/moose-core/diffusion/Dsolve.cpp index b90fbc267e6d419b04b912ca5070af8919436df0..9f1fc45af0e0e566eec1971ca6060223552f19ef 100644 --- a/moose-core/diffusion/Dsolve.cpp +++ b/moose-core/diffusion/Dsolve.cpp @@ -422,7 +422,9 @@ void Dsolve::calcJnChan( const DiffJunction& jn, Dsolve* other, double dt ) double lastN = myN; double otherN = otherDv.getN( j->second ); double chanN = chanDv.getN( j->first ); - double perm = myChan.permeability * chanN / NA; + // Stick in a conversion factor for the myN and otherN into + // concentrations. Note that SI is millimolar. + double perm = myChan.permeability * chanN * 1000.0 / NA; myN = integ( myN, perm * myN/j->firstVol, perm * otherN/j->secondVol, dt ); otherN += lastN - myN; // Mass consv @@ -459,7 +461,9 @@ void Dsolve::calcOtherJnChan( const DiffJunction& jn, Dsolve* other, double dt ) double lastN = myN; double otherN = otherDv.getN( j->second ); double chanN = chanDv.getN( j->second ); - double perm = otherChan.permeability * chanN / NA; + // Stick in a conversion factor for the myN and otherN into + // concentrations. Note that SI is millimolar. + double perm = otherChan.permeability * chanN * 1000.0 / NA; myN = integ( myN, perm * myN/j->firstVol, perm * otherN/j->secondVol, dt ); otherN += lastN - myN; // Mass consv @@ -1130,7 +1134,7 @@ double Dsolve::getNinit( const Eref& e ) const { return pools_[ pid ].getNinit( vox ); } - cout << "Warning: Dsolve::setNinit: Eref " << e << " out of range " << + cout << "Warning: Dsolve::getNinit: Eref " << e << " out of range " << pools_.size() << ", " << numVoxels_ << "\n"; return 0.0; } diff --git a/moose-core/kinetics/ReadKkit.cpp b/moose-core/kinetics/ReadKkit.cpp index 0b3824643f39d5650183ed4342d22a09942533f9..01254b62292d6a3cc1f2bdf92638e4c04dcadd3c 100644 --- a/moose-core/kinetics/ReadKkit.cpp +++ b/moose-core/kinetics/ReadKkit.cpp @@ -1195,6 +1195,7 @@ Id ReadKkit::buildChan( const vector< string >& args ) assert( chan != Id() ); string chanPath = clean.substr( 10 ); chanIds_[ chanPath ] = chan; + Id info = buildInfo( chan, chanMap_, args ); return chan; } diff --git a/moose-core/mesh/EndoMesh.cpp b/moose-core/mesh/EndoMesh.cpp index 7902e3bded47f9e28f2473eeadbfed264d647fc8..e1e21dfdec33ec10d8ba5c9061e42b34f09c1fd6 100644 --- a/moose-core/mesh/EndoMesh.cpp +++ b/moose-core/mesh/EndoMesh.cpp @@ -23,6 +23,9 @@ #include "NeuroMesh.h" #include "EndoMesh.h" #include "../utility/numutil.h" + +static CubeMesh defaultParent; + const Cinfo* EndoMesh::initCinfo() { ////////////////////////////////////////////////////////////// @@ -137,7 +140,7 @@ static const Cinfo* endoMeshCinfo = EndoMesh::initCinfo(); ////////////////////////////////////////////////////////////////// EndoMesh::EndoMesh() : - parent_( 0 ), + parent_( &defaultParent ), rPower_( 1.0 / 3.0 ), rScale_( 0.5 ), aPower_( 0.5 ), @@ -263,7 +266,7 @@ unsigned int EndoMesh::getMeshDimensions( unsigned int fid ) const /// Virtual function to return # of spatial dimensions of mesh unsigned int EndoMesh::innerGetDimensions() const { - return 1; + return 3; } /// Virtual function to return volume of mesh Entry. double EndoMesh::getMeshEntryVolume( unsigned int fid ) const diff --git a/moose-core/python/moose/SBML/readSBML.py b/moose-core/python/moose/SBML/readSBML.py index d9b0e90f0b531fb446123c5409c3b46e6bfd8de7..f8e74f7e13bdd9194b60dadc4f1f17a336fdd153 100644 --- a/moose-core/python/moose/SBML/readSBML.py +++ b/moose-core/python/moose/SBML/readSBML.py @@ -13,10 +13,11 @@ ** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS Created : Thu May 13 10:19:00 2016(+0530) Version -Last-Updated: Fri May 21 11:21:00 2018(+0530) +Last-Updated: Fri Oct 26 11:21:00 2018(+0530) By:HarshaRani **********************************************************************/ 2018 +Oct 26: - validator can be switchedoff by passing validate="off" while readSBML files May 18: - cleanedup and connected cplx pool to correct parent enzyme Jan 6: - only if valid model exists, then printing the no of compartment,pool,reaction etc - at reaction level a check made to see if path exist while creating a new reaction @@ -77,7 +78,7 @@ try: except ImportError: pass -def mooseReadSBML(filepath, loadpath, solver="ee"): +def mooseReadSBML(filepath, loadpath, solver="ee",validate="on"): """Load SBML model """ global foundLibSBML_ @@ -100,7 +101,10 @@ def mooseReadSBML(filepath, loadpath, solver="ee"): filep = open(filepath, "r") document = libsbml.readSBML(filepath) tobecontinue = False - tobecontinue = validateModel(document) + if validate == "on": + tobecontinue = validateModel(document) + else: + tobecontinue = True if tobecontinue: level = document.getLevel() version = document.getVersion() diff --git a/moose-core/python/moose/SBML/writeSBML.py b/moose-core/python/moose/SBML/writeSBML.py index ec9d357f24443cb458957c06ebb66f452c5b2efb..a3d898d44c2b0c02d2f8e4f55404332f6ac9c86a 100644 --- a/moose-core/python/moose/SBML/writeSBML.py +++ b/moose-core/python/moose/SBML/writeSBML.py @@ -18,8 +18,11 @@ Last-Updated: Mon 30 Apr 15:10:00 2018(+0530) **********************************************************************/ /**************************** 2018 +Oct 16: CylMesh's comparment volume is written, but zeroth volex details are populated +Oct 13: CylMesh are written to SBML with annotation field and only zeroth element/voxel (incase of cylMesh) of moose object is written +Oct 1 : corrected the spell of CyclMesh-->CylMesh, negating the yaxis for kkit is removed Apr 30: indentation corrected while writting annotation for enzymecomplex -Jan 6: for Product_formation_, k3 units depends on noofSub, prd was passed which is fixed +Jan 6 : for Product_formation_, k3 units depends on noofSub, prd was passed which is fixed 2017 Dec 15: If model path exists is checked Enz cplx is written only if enz,cplx,sub, prd exist, a clean check is made @@ -73,10 +76,12 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): if not moose.exists(modelpath): return False, "Path doesn't exist" elif moose.exists(modelpath): - mObj = moose.wildcardFind(moose.element(modelpath).path+'/##[ISA=PoolBase]'+','+ - moose.element(modelpath).path+'/##[ISA=ReacBase]'+','+ - moose.element(modelpath).path+'/##[ISA=EnzBase]'+','+ - moose.element(modelpath).path+'/##[ISA=StimulusTable]') + checkCompt = moose.wildcardFind(modelpath+'/##[0][ISA=ChemCompt]') + + mObj = moose.wildcardFind(moose.element(modelpath).path+'/##[0][ISA=PoolBase]'+','+ + moose.element(modelpath).path+'/##[0][ISA=ReacBase]'+','+ + moose.element(modelpath).path+'/##[0][ISA=EnzBase]'+','+ + moose.element(modelpath).path+'/##[0][ISA=StimulusTable]') for p in mObj: if not isinstance(moose.element(p.parent),moose.CplxEnzBase): if moose.exists(p.path+'/info'): @@ -102,7 +107,7 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): cremodel_.setLengthUnits("length") neutralNotes = "" - specieslist = moose.wildcardFind(modelpath + '/##[ISA=PoolBase]') + specieslist = moose.wildcardFind(modelpath + '/##[0][ISA=PoolBase]') if specieslist: neutralPath = getGroupinfo(specieslist[0]) if moose.exists(neutralPath.path + '/info'): @@ -117,18 +122,24 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): srcdesConnection = {} writeUnits(cremodel_) + modelAnno = writeSimulationAnnotation(modelpath) if modelAnno: cremodel_.setAnnotation(modelAnno) groupInfo = {} + compartexist, groupInfo = writeCompt(modelpath, cremodel_) + if compartexist == True: species = writeSpecies( modelpath,cremodel_,sbmlDoc,sceneitems,groupInfo) if species: writeFunc(modelpath, cremodel_) reacGroup = {} + writeReac(modelpath, cremodel_, sceneitems,groupInfo) + writeEnz(modelpath, cremodel_, sceneitems,groupInfo) + if groupInfo: for key,value in groupInfo.items(): mplugin = cremodel_.getPlugin("groups") @@ -152,6 +163,7 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): for values in value: member = group.createMember() member.setIdRef(values) + consistencyMessages = "" SBMLok = validateModel(sbmlDoc) if (SBMLok): @@ -165,9 +177,9 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}): return -1, consistencyMessages else: return False, "Atleast one compartment should exist to write SBML" - + def writeEnz(modelpath, cremodel_, sceneitems,groupInfo): - for enz in moose.wildcardFind(modelpath + '/##[ISA=EnzBase]'): + for enz in moose.wildcardFind(modelpath + '/##[0][ISA=EnzBase]'): enzannoexist = False enzGpnCorCol = " " cleanEnzname = convertSpecialChar(enz.name) @@ -630,7 +642,7 @@ def listofname(reacSub, mobjEnz): def writeReac(modelpath, cremodel_, sceneitems,reacGroup): - for reac in moose.wildcardFind(modelpath + '/##[ISA=ReacBase]'): + for reac in moose.wildcardFind(modelpath + '/##[0][ISA=ReacBase]'): reacSub = reac.neighbors["sub"] reacPrd = reac.neighbors["prd"] if (len(reacSub) != 0 and len(reacPrd) != 0): @@ -763,7 +775,7 @@ def writeReac(modelpath, cremodel_, sceneitems,reacGroup): def writeFunc(modelpath, cremodel_): - funcs = moose.wildcardFind(modelpath + '/##[ISA=Function]') + funcs = moose.wildcardFind(modelpath + '/##[0][ISA=Function]') # if func: foundFunc = False for func in funcs: @@ -824,7 +836,7 @@ def getGroupinfo(element): # if /modelpath/Compartment/Group/Group1/Pool, then I check and get Group1 # And /modelpath is also a NeutralObject,I stop till I find Compartment - while not mooseIsInstance(element, ["Neutral", "CubeMesh", "CyclMesh"]): + while not mooseIsInstance(element, ["Neutral", "CubeMesh", "CylMesh"]): element = element.parent return element @@ -836,7 +848,7 @@ def idBeginWith(name): return changedName def findGroup_compt(melement): - while not (mooseIsInstance(melement, ["Neutral","CubeMesh", "CyclMesh"])): + while not (mooseIsInstance(melement, ["Neutral","CubeMesh", "CylMesh"])): melement = melement.parent return melement @@ -865,7 +877,7 @@ def convertSpecialChar(str1): def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems,speGroup): # getting all the species - for spe in moose.wildcardFind(modelpath + '/##[ISA=PoolBase]'): + for spe in moose.wildcardFind(modelpath + '/##[0][ISA=PoolBase]'): sName = convertSpecialChar(spe.name) comptVec = findCompartment(spe) speciannoexist = False @@ -977,14 +989,25 @@ def writeSpecies(modelpath, cremodel_, sbmlDoc, sceneitems,speGroup): def writeCompt(modelpath, cremodel_): # getting all the compartments - compts = moose.wildcardFind(modelpath + '/##[ISA=ChemCompt]') + compts = moose.wildcardFind(modelpath + '/##[0][ISA=ChemCompt]') groupInfo = {} for compt in compts: + comptAnno = "" comptName = convertSpecialChar(compt.name) # converting m3 to litre - size = compt.volume * pow(10, 3) + if isinstance(compt,moose.CylMesh): + size = (compt.volume/compt.numDiffCompts)*pow(10,3) + comptAnno = "<moose:CompartmentAnnotation><moose:Mesh>" + \ + str(compt.className) + "</moose:Mesh>\n" + \ + "<moose:numDiffCompts>" + \ + str(compt.numDiffCompts)+ "</moose:numDiffCompts>\n" + \ + "</moose:CompartmentAnnotation>" + else: + size = compt.volume * pow(10, 3) ndim = compt.numDimensions c1 = cremodel_.createCompartment() + if comptAnno: + c1.setAnnotation(comptAnno) c1.setId(str(idBeginWith(comptName + "_" + str(compt.getId().value) + @@ -993,11 +1016,12 @@ def writeCompt(modelpath, cremodel_): "_"))) c1.setName(comptName) c1.setConstant(True) - c1.setSize(size) + c1.setSize(compt.volume*pow(10,3)) + #c1.setSize(size) c1.setSpatialDimensions(ndim) c1.setUnits('volume') #For each compartment get groups information along - for grp in moose.wildcardFind(compt.path+'/##[TYPE=Neutral]'): + for grp in moose.wildcardFind(compt.path+'/##[0][TYPE=Neutral]'): grp_cmpt = findGroup_compt(grp.parent) try: value = groupInfo[moose.element(grp)] @@ -1011,6 +1035,7 @@ def writeCompt(modelpath, cremodel_): else: return False,groupInfo # write Simulation runtime,simdt,plotdt + def writeSimulationAnnotation(modelpath): modelAnno = "" plots = "" @@ -1027,7 +1052,8 @@ def writeSimulationAnnotation(modelpath): modelAnno = modelAnno + "<moose:plotdt> " + \ str(mooseclock.dts[18]) + " </moose:plotdt>\n" plots = "" - graphs = moose.wildcardFind(modelpath + "/##[TYPE=Table2]") + graphs = moose.wildcardFind(modelpath + "/##[0][TYPE=Table2]") + for gphs in range(0, len(graphs)): gpath = graphs[gphs].neighbors['requestOut'] if len(gpath) != 0: @@ -1035,7 +1061,7 @@ def writeSimulationAnnotation(modelpath): ori = q.path name = convertSpecialChar(q.name) graphSpefound = False - while not(isinstance(moose.element(q), moose.CubeMesh)): + while not(isinstance(moose.element(q), moose.CubeMesh) or isinstance(moose.element(q),moose.CylMesh)): q = q.parent graphSpefound = True if graphSpefound: @@ -1046,6 +1072,7 @@ def writeSimulationAnnotation(modelpath): else: plots = plots + "; "+ori[ori.find(q.name)-1:len(ori)] #plots = plots + "; /" + q.name + '/' + name + if plots != " ": modelAnno = modelAnno + "<moose:plots> " + plots + "</moose:plots>\n" modelAnno = modelAnno + "</moose:ModelAnnotation>" diff --git a/moose-core/python/moose/genesis/writeKkit.py b/moose-core/python/moose/genesis/writeKkit.py index 8fc6356902f2a4e07f5ecb586f34327d08a6b06c..4dc03c076b34ef859a9df055c440283c945aedc6 100644 --- a/moose-core/python/moose/genesis/writeKkit.py +++ b/moose-core/python/moose/genesis/writeKkit.py @@ -9,8 +9,11 @@ __version__ = "1.0.0" __maintainer__ = "Harsha Rani" __email__ = "hrani@ncbs.res.in" __status__ = "Development" -__updated__ = "Aug 8 2017" +__updated__ = "Oct 23 2018" +# 2018 +# Oct 16: Channels are written back to genesis +# only zeroth element is taken for written back to genesis, this is true for CubeMesh and CylMesh # 2017 # Aug 8 : All the moose object which doesn't have compartment are not written. Error message are appended @@ -69,7 +72,7 @@ def mooseWriteKkit( modelpath, filename, sceneitems={}): cmin, xmin, ymin = 0, 0, 0 cmax, xmax, ymax = 1, 1, 1 - compt = moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]') + compt = moose.wildcardFind(modelpath+'/##[0][ISA=ChemCompt]') maxVol = estimateDefaultVol(compt) positionInfoExist = True if compt: @@ -108,16 +111,19 @@ def mooseWriteKkit( modelpath, filename, sceneitems={}): enzList,error = writeEnz(modelpath,f,sceneitems) errors = errors+error - error = writeSumtotal(modelpath,f) + chanList, error = writeConcChan(modelpath,f,sceneitems) errors = errors+error error = writeStimulus(modelpath,f) errors = errors+error + error = writeSumtotal(modelpath,f) + errors = errors+error + #writeSumtotal(modelpath,f) f.write("simundump xgraph /graphs/conc1 0 0 99 0.001 0.999 0\n" "simundump xgraph /graphs/conc2 0 0 100 0 1 0\n") - tgraphs = moose.wildcardFind(modelpath+'/##[ISA=Table2]') + tgraphs = moose.wildcardFind(modelpath+'/##[0][ISA=Table2]') first, second = " ", " " if tgraphs: first,second = writeplot(tgraphs,f) @@ -133,6 +139,7 @@ def mooseWriteKkit( modelpath, filename, sceneitems={}): "simundump xtext /file/notes 0 1\n") storeReacMsg(reacList,f) storeEnzMsg(enzList,f) + storeChanMsg(chanList,f) if tgraphs: storePlotMsgs(tgraphs,f) writeFooter1(f) @@ -171,7 +178,7 @@ def calPrime(x): def storeCplxEnzMsgs( enz, f ): for sub in enz.neighbors["subOut"]: - s = "addmsg /kinetics/" + trimPath( sub ) + " /kinetics/" + trimPath(enz) + " SUBSTRATE n \n"; + s = "addmsg /kinetics/" + trimPath( moose.element(sub) ) + " /kinetics/" + trimPath(enz) + " SUBSTRATE n \n"; s = s+ "addmsg /kinetics/" + trimPath( enz ) + " /kinetics/" + trimPath( sub ) + " REAC sA B \n"; f.write(s) for prd in enz.neighbors["prd"]: @@ -206,9 +213,72 @@ def storeEnzMsg( enzList, f): else: storeCplxEnzMsgs( enz, f ) +def storeChanMsg(chanList,f): + for channel in chanList: + for chanOL in channel.neighbors['out']: + eo = "addmsg /kinetics/" + trimPath(moose.element(channel)) + " /kinetics/" + trimPath(moose.element(chanOL))+ " REAC B A \n"; + eo = eo +"addmsg /kinetics/" + trimPath(moose.element(chanOL)) + " /kinetics/"+trimPath(moose.element(channel))+" PRODUCT n vol \n"; + f.write(eo) + + for chanIL in channel.neighbors['in']: + ei = "addmsg /kinetics/" + trimPath(moose.element(channel)) + " /kinetics/" + trimPath(moose.element(chanIL))+ " REAC A B \n"; + ei = ei +"addmsg /kinetics/" + trimPath(moose.element(chanIL)) + " /kinetics/"+trimPath(moose.element(channel))+" SUBSTRATE n vol \n"; + f.write(ei) + + for chanSNC in channel.neighbors['setNumChan']: + cff = "addmsg /kinetics/"+trimPath(moose.element(chanSNC))+ " /kinetics/"+trimPath(moose.element(channel))+ " NUMCHAN n \n" + f.write(cff) + +def writeConcChan(modelpath,f,sceneitems): + error = "" + concChanList = moose.wildcardFind(modelpath+'/##[0][ISA=ConcChan]') + for cChan in concChanList: + if findCompartment(cChan) == moose.element('/'): + error = error + " \n "+cChan.path+ " doesn't have compartment ignored to write to genesis" + else: + x = random.randrange(0,10) + y = random.randrange(0,10) + textcolor = "" + color = "" + if len(moose.element(cChan).neighbors['setNumChan']) == 1: + chanParent = moose.element(moose.element(cChan).neighbors['setNumChan'][0]) + + if not (isinstance(chanParent,moose.PoolBase)): + print(" raise exception Channel doesn't have pool as parent %s",moose.element(cChan).path) + return False,"raise exception Channel doesn't have pool as parent" + else: + vol = chanParent.volume * NA * 1e-3; + + cinfo = cChan.path+'/info' + if moose.exists(cinfo): + x = moose.Annotator(cinfo).getField('x') + y = moose.Annotator(cinfo).getField('y') + #x = sceneitems[cChan]['x'] + #y = sceneitems[cChan]['y'] + color = moose.Annotator(cinfo).getField('color') + color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) + + textcolor = moose.Annotator(cinfo).getField('textColor') + textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) + else: + error = error + "\n x and y co-ordinates are not specified for `" + cChan.name+ "` zero will be assigned \n " + if color == "" or color == " ": + color = getRandColor() + if textcolor == "" or textcolor == " ": + textcolor = getRandColor() + f.write("simundump kchan /kinetics/" + trimPath(cChan)+ " " + str(int(1)) + " " + str(cChan.permeability)+ " " + + str(int(0)) + " " + + str(int(0)) + " " + + str(int(0)) + " " + + str(int(0)) + " " + + str("") + " " + + str(textcolor) + " " + str(color) + " \"\"" + + " " + str(int(x)) + " " + str(int(y)) + " "+str(int(0))+"\n") + + return concChanList,error def writeEnz( modelpath,f,sceneitems): error = "" - enzList = moose.wildcardFind(modelpath+'/##[ISA=EnzBase]') + enzList = moose.wildcardFind(modelpath+'/##[0][ISA=EnzBase]') for enz in enzList: if findCompartment(enz) == moose.element('/'): error = error + " \n "+enz.path+ " doesn't have compartment ignored to write to genesis" @@ -313,7 +383,7 @@ def storeReacMsg(reacList,f): def writeReac(modelpath,f,sceneitems): error = "" - reacList = moose.wildcardFind(modelpath+'/##[ISA=ReacBase]') + reacList = moose.wildcardFind(modelpath+'/##[0][ISA=ReacBase]') for reac in reacList: if findCompartment(reac) == moose.element('/'): error = error + " \n "+reac.path+ " doesn't have compartment ignored to write to genesis" @@ -385,7 +455,7 @@ def trimPath(mobj): def writeSumtotal( modelpath,f): error = "" - funclist = moose.wildcardFind(modelpath+'/##[ISA=Function]') + funclist = moose.wildcardFind(modelpath+'/##[0][ISA=Function]') s = "" for func in funclist: fInfound = True @@ -420,7 +490,7 @@ def writeSumtotal( modelpath,f): def writeStimulus(modelpath,f): error = "" - if len(moose.wildcardFind(modelpath+'/##[ISA=StimulusTable]')): + if len(moose.wildcardFind(modelpath+'/##[0][ISA=StimulusTable]')): error = error +'\n StimulusTable is not written into genesis. This is in Todo List' return error @@ -501,7 +571,8 @@ def writePool(modelpath,f,volIndex,sceneitems): error = "" color = "" textcolor = "" - for p in moose.wildcardFind(modelpath+'/##[ISA=PoolBase]'): + + for p in moose.wildcardFind(modelpath+'/##[0][ISA=PoolBase]'): if findCompartment(p) == moose.element('/'): error = error + " \n "+p.path+ " doesn't have compartment ignored to write to genesis" else: @@ -541,9 +612,8 @@ def writePool(modelpath,f,volIndex,sceneitems): color = getColorCheck(color,GENESIS_COLOR_SEQUENCE) textcolor = moose.Annotator(pinfo).getField('textColor') textcolor = getColorCheck(textcolor,GENESIS_COLOR_SEQUENCE) - - - geometryName = volIndex[p.volume] + poolsCmpt = findCompartment(p) + geometryName = volIndex[float(poolsCmpt.volume)] volume = p.volume * NA * 1e-3 if color == "" or color == " ": color = getRandColor() @@ -616,6 +686,10 @@ def writeCompartment(modelpath,compts,f): l = len(compts) geometry = "" for compt in compts: + # if isinstance(compt,moose.CylMesh): + # print " 1 " + # size = (compt.volume/compt.numDiffCompts) + # else: size = compt.volume ndim = compt.numDimensions vecIndex = l-i-1 @@ -624,17 +698,18 @@ def writeCompartment(modelpath,compts,f): y = ymax+1 if vecIndex > 0: geometry = geometry+"simundump geometry /kinetics" + "/geometry[" + str(vecIndex) +"] 0 " + str(size) + " " + str(ndim) + " sphere " +" \"\" white black "+ str(int(x)) + " " +str(int(y)) +" 0\n"; - volIndex[size] = "/geometry["+str(vecIndex)+"]" + volIndex[float(size)] = "/geometry["+str(vecIndex)+"]" else: + geometry = geometry+"simundump geometry /kinetics" + "/geometry 0 " + str(size) + " " + str(ndim) + " sphere " +" \"\" white black " + str(int(x)) + " "+str(int(y))+ " 0\n"; - volIndex[size] = "/geometry" - f.write(geometry) + volIndex[float(size)] = "/geometry" + f.write(geometry) writeGroup(modelpath,f) return volIndex def writeGroup(modelpath,f): ignore = ["graphs","moregraphs","geometry","groups","conc1","conc2","conc3","conc4","model","data","graph_0","graph_1","graph_2","graph_3","graph_4","graph_5"] - for g in moose.wildcardFind(modelpath+'/##[TYPE=Neutral]'): + for g in moose.wildcardFind(modelpath+'/##[0][TYPE=Neutral]'): if not g.name in ignore: if trimPath(g) != None: x = xmin+1 @@ -705,12 +780,12 @@ def writeNotes(modelpath,f): notes = "" #items = moose.wildcardFind(modelpath+"/##[ISA=ChemCompt],/##[ISA=ReacBase],/##[ISA=PoolBase],/##[ISA=EnzBase],/##[ISA=Function],/##[ISA=StimulusTable]") items = [] - items = moose.wildcardFind(modelpath+"/##[ISA=ChemCompt]") +\ - moose.wildcardFind(modelpath+"/##[ISA=PoolBase]") +\ - moose.wildcardFind(modelpath+"/##[ISA=ReacBase]") +\ - moose.wildcardFind(modelpath+"/##[ISA=EnzBase]") +\ - moose.wildcardFind(modelpath+"/##[ISA=Function]") +\ - moose.wildcardFind(modelpath+"/##[ISA=StimulusTable]") + items = moose.wildcardFind(modelpath+"/##[0][ISA=ChemCompt]") +\ + moose.wildcardFind(modelpath+"/##[0][ISA=PoolBase]") +\ + moose.wildcardFind(modelpath+"/##[0][ISA=ReacBase]") +\ + moose.wildcardFind(modelpath+"/##[0][ISA=EnzBase]") +\ + moose.wildcardFind(modelpath+"/##[0][ISA=Function]") +\ + moose.wildcardFind(modelpath+"/##[0][ISA=StimulusTable]") for item in items: if moose.exists(item.path+'/info'): info = item.path+'/info' @@ -741,4 +816,4 @@ if __name__ == "__main__": if written: print(" file written to ",output) else: - print(" could be written to kkit format") + print(" could be written to kkit format") \ No newline at end of file diff --git a/moose-core/python/moose/moose.py b/moose-core/python/moose/moose.py index 069d43b9990f019dbcce29d3adc493d1b75d8ff6..a55bc7be1f1c6b86d36b838d68a7be7b9e245d17 100644 --- a/moose-core/python/moose/moose.py +++ b/moose-core/python/moose/moose.py @@ -122,7 +122,7 @@ known_types = ['void', 'melement'] + sequence_types # SBML related functions. -def mooseReadSBML(filepath, loadpath, solver='ee'): +def mooseReadSBML(filepath, loadpath, solver='ee',validate="on"): """Load SBML model. keyword arguments: \n @@ -134,7 +134,7 @@ def mooseReadSBML(filepath, loadpath, solver='ee'): """ global sbmlImport_ if sbmlImport_: - return _readSBML.mooseReadSBML( filepath, loadpath, solver ) + return _readSBML.mooseReadSBML( filepath, loadpath, solver,validate ) else: print( sbmlError_ ) return False @@ -222,7 +222,7 @@ def mergeChemModel(src, des): return False # NML2 reader and writer function. -def mooseReadNML2( modelpath ): +def mooseReadNML2( modelpath, verbose = False ): """Read NeuroML model (version 2) and return reader object. """ global nml2Import_ @@ -230,7 +230,7 @@ def mooseReadNML2( modelpath ): mu.warn( nml2ImportError_ ) raise RuntimeError( "Could not load NML2 support." ) - reader = _neuroml2.NML2Reader( ) + reader = _neuroml2.NML2Reader( verbose = verbose ) reader.read( modelpath ) return reader diff --git a/moose-core/python/moose/neuroml2/reader.py b/moose-core/python/moose/neuroml2/reader.py index 6a9dec0fafb5d0574cfa7e196ae7b5a35d7c90fa..2f2e2d85c8c35c5b227b5e436fd4c39c26521409 100644 --- a/moose-core/python/moose/neuroml2/reader.py +++ b/moose-core/python/moose/neuroml2/reader.py @@ -33,13 +33,6 @@ import moose.utils as mu from .units import SI -def _unique( ls ): - res = [ ] - for l in ls: - if l not in res: - res.append( l ) - return res - def _unique( ls ): res = [ ] for l in ls: @@ -134,6 +127,7 @@ class NML2Reader(object): self.pop_to_cell_type = {} self.seg_id_to_comp_name = {} self.paths_to_chan_elements = {} + self.network = None def read(self, filename, symmetric=True): filename = os.path.realpath( filename ) @@ -146,9 +140,9 @@ class NML2Reader(object): if len(self.doc.networks)>=1: self.network = self.doc.networks[0] - moose.celsius = self._getTemperature() + self.importConcentrationModels(self.doc) self.importIonChannels(self.doc) self.importInputs(self.doc) @@ -163,10 +157,13 @@ class NML2Reader(object): mu.info("Read all from %s"%filename) def _getTemperature(self): - if self.network.type=="networkWithTemperature": - return SI(self.network.temperature) - else: - return 0 # Why not, if there's no temp dependence in nml..? + if self.network is not None: + if self.network.type=="networkWithTemperature": + return SI(self.network.temperature) + else: + # Why not, if there's no temp dependence in nml..? + return 0 + return SI('25') def getCellInPopulation(self, pop_id, index): return self.cells_in_populations[pop_id][index] @@ -407,11 +404,8 @@ class NML2Reader(object): mu.info("Using %s to evaluate rate"%ct.name) rate = [] for v in tab: - vals = pynml.evaluate_component(ct - , req_variables = - {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()} - ) - # mu.info vals + req_vars = {'v':'%sV'%v,'vShift':vShift,'temperature':self._getTemperature()} + vals = pynml.evaluate_component(ct, req_variables = req_vars) if 'x' in vals: rate.append(vals['x']) if 't' in vals: @@ -420,6 +414,9 @@ class NML2Reader(object): rate.append(vals['r']) return np.array(rate) + print( "[WARN ] Could not determine rate: %s %s %s" %(ratefn.type,vmin,vmax)) + return np.array([]) + def importChannelsToCell(self, nmlcell, moosecell, membrane_properties): sg_to_segments = self._cell_to_sg[nmlcell] for chdens in membrane_properties.channel_densities + membrane_properties.channel_density_v_shifts: @@ -487,10 +484,11 @@ class NML2Reader(object): def createHHChannel(self, chan, vmin=-150e-3, vmax=100e-3, vdivs=5000): mchan = moose.HHChannel('%s/%s' % (self.lib.path, chan.id)) mgates = [moose.element(x) for x in [mchan.gateX, mchan.gateY, mchan.gateZ]] - assert(len(chan.gate_hh_rates) <= 3) # We handle only up to 3 gates in HHCHannel + assert len(chan.gate_hh_rates)<=3, "We handle only up to 3 gates in HHCHannel" if self.verbose: mu.info('== Creating channel: %s (%s) -> %s (%s)'%(chan.id, chan.gate_hh_rates, mchan, mgates)) + all_gates = chan.gates + chan.gate_hh_rates for ngate, mgate in zip(all_gates,mgates): if mgate.name.endswith('X'): @@ -512,9 +510,7 @@ class NML2Reader(object): # refering to tau_inf and m_inf?? fwd = ngate.forward_rate rev = ngate.reverse_rate - self.paths_to_chan_elements['%s/%s'%(chan.id,ngate.id)] = '%s/%s'%(chan.id,mgate.name) - q10_scale = 1 if ngate.q10_settings: if ngate.q10_settings.type == 'q10Fixed': @@ -536,6 +532,7 @@ class NML2Reader(object): if (fwd is not None) and (rev is not None): alpha = self.calculateRateFn(fwd, vmin, vmax, vdivs) beta = self.calculateRateFn(rev, vmin, vmax, vdivs) + mgate.tableA = q10_scale * (alpha) mgate.tableB = q10_scale * (alpha + beta) @@ -554,8 +551,9 @@ class NML2Reader(object): tau = 1 / (alpha + beta) if inf is not None: inf = self.calculateRateFn(inf, vmin, vmax, vdivs) - mgate.tableA = q10_scale * (inf / tau) - mgate.tableB = q10_scale * (1 / tau) + if len(inf) > 0: + mgate.tableA = q10_scale * (inf / tau) + mgate.tableB = q10_scale * (1 / tau) if self.verbose: mu.info('%s: Created %s for %s'%(self.filename,mchan.path,chan.id)) @@ -589,7 +587,7 @@ class NML2Reader(object): def importIonChannels(self, doc, vmin=-150e-3, vmax=100e-3, vdivs=5000): if self.verbose: - mu.info('%s : Importing the ion channels' % self.filename ) + mu.info('%s: Importing the ion channels' % self.filename ) for chan in doc.ion_channel+doc.ion_channel_hhs: if chan.type == 'ionChannelHH': @@ -603,7 +601,7 @@ class NML2Reader(object): self.nml_to_moose[chan] = mchan self.proto_chans[chan.id] = mchan if self.verbose: - mu.info( self.filename + ' : Created ion channel %s for %s %s'%( + mu.info( self.filename + ': Created ion channel %s for %s %s'%( mchan.path, chan.type, chan.id)) def importConcentrationModels(self, doc): diff --git a/moose-core/python/rdesigneur/rdesigneur.py b/moose-core/python/rdesigneur/rdesigneur.py index 492cb1802c31b6b614351578bceb0bcf58a50944..6ad93d02b89981d95b1bc215df7d6e8ecffce079 100644 --- a/moose-core/python/rdesigneur/rdesigneur.py +++ b/moose-core/python/rdesigneur/rdesigneur.py @@ -179,7 +179,7 @@ class rdesigneur: self.elecid.numCompartments, "compartments and", self.elecid.numSpines, "spines on", len( self.cellPortionElist ), "compartments.") - if hasattr( self , 'chemid' ): + if hasattr( self , 'chemid') and len( self.chemDistrib ) > 0: dmstoich = moose.element( self.dendCompt.path + '/stoich' ) print("\tChem part of model has the following compartments: ") for j in moose.wildcardFind( '/model/chem/##[ISA=ChemCompt]'): @@ -948,7 +948,7 @@ rdesigneur.rmoogli.updateMoogliViewer() } stims = moose.Neutral( self.modelPath + '/stims' ) k = 0 - # Stimlist = [path, geomExpr, relPath, field, expr_string] + # rstim class has {elecpath, geom_expr, relpath, field, expr} for i in self.stimList: pair = i.elecpath + " " + i.geom_expr dendCompts = self.elecid.compartmentsFromExpression[ pair ] @@ -1063,7 +1063,6 @@ rdesigneur.rmoogli.updateMoogliViewer() self._buildNeuroMesh() - self._configureSolvers() for i in self.adaptorList: # print(i) @@ -1326,7 +1325,7 @@ rdesigneur.rmoogli.updateMoogliViewer() ################################################################# def _configureSolvers( self ) : - if not hasattr( self, 'chemid' ): + if not hasattr( self, 'chemid' ) or len( self.chemDistrib ) == 0: return if not hasattr( self, 'dendCompt' ): raise BuildError( "configureSolvers: no chem meshes defined." ) diff --git a/moose-core/python/rdesigneur/rdesigneurProtos.py b/moose-core/python/rdesigneur/rdesigneurProtos.py index c85687c5fbe0611b97339e942c274b9f808a8cae..cade79201b388014a2719c513526a0c0eda63417 100644 --- a/moose-core/python/rdesigneur/rdesigneurProtos.py +++ b/moose-core/python/rdesigneur/rdesigneurProtos.py @@ -50,6 +50,7 @@ import math from moose import utils EREST_ACT = -0.060 +EREST_ACT_HH = -0.070 # A different value is used for the HH squid params ECA = 0.080 EK = -0.075 SOMA_A = 3.32e-9 @@ -71,7 +72,7 @@ def make_HH_Na(name = 'HH_Na', parent='/library', vmin=-110e-3, vmax=50e-3, vdiv na.Ek = 50e-3 na.Xpower = 3 na.Ypower = 1 - v = np.linspace(vmin, vmax, vdivs+1) - EREST_ACT + v = np.linspace(vmin, vmax, vdivs+1) - EREST_ACT_HH m_alpha = per_ms * (25 - v * 1e3) / (10 * (np.exp((25 - v * 1e3) / 10) - 1)) m_beta = per_ms * 4 * np.exp(- v * 1e3/ 18) m_gate = moose.element('%s/gateX' % (na.path)) @@ -100,7 +101,7 @@ def make_HH_K(name = 'HH_K', parent='/library', vmin=-120e-3, vmax=40e-3, vdivs= k = moose.HHChannel('%s/%s' % (parent, name)) k.Ek = -77e-3 k.Xpower = 4 - v = np.linspace(vmin, vmax, vdivs+1) - EREST_ACT + v = np.linspace(vmin, vmax, vdivs+1) - EREST_ACT_HH n_alpha = per_ms * (10 - v * 1e3)/(100 * (np.exp((10 - v * 1e3)/10) - 1)) n_beta = per_ms * 0.125 * np.exp(- v * 1e3 / 80) n_gate = moose.element('%s/gateX' % (k.path)) diff --git a/moose-core/tests/python/test_Xchan1.py b/moose-core/tests/python/test_Xchan1.py index af555663817579d728465bb256055821fa1567a0..d1eb35114cf02a35775b708e43f30dc985db127e 100644 --- a/moose-core/tests/python/test_Xchan1.py +++ b/moose-core/tests/python/test_Xchan1.py @@ -63,7 +63,7 @@ def makeModel(): s.concInit = 0.001 chanPool.nInit = 1000.0 # Flux (mM/s) = permeability * N * (conc_out - conc_in ) - chan.permeability = 0.1 * chanPool.volume * 6.022e23 / chanPool.nInit + chan.permeability = 0.1 * chanPool.volume * 6.022e23 *0.001 / chanPool.nInit ##################################################################### fixXreacs.fixXreacs( '/model' )