diff --git a/moose-core/VERSION b/moose-core/VERSION
index 498dab1f0ef5ca79fb07efae9f6df24223b7108e..90dd2fedc8c7a43f0b6f4c6181ab599c8c86e23c 100644
--- a/moose-core/VERSION
+++ b/moose-core/VERSION
@@ -1 +1 @@
-3.1.1-126-g0d5fb79
+3.1.1-100-g51a3a30
diff --git a/moose-core/kinetics/Pool.cpp b/moose-core/kinetics/Pool.cpp
index 8b8b656148f78c6a0390f22471e05eb29fea80a4..7ebc22f93905867695599f10bb0efeb7598c2694 100644
--- a/moose-core/kinetics/Pool.cpp
+++ b/moose-core/kinetics/Pool.cpp
@@ -32,23 +32,8 @@ const Cinfo* Pool::initCinfo()
);
//////////////////////////////////////////////////////////////
- // MsgDest Definitions: All but increment and decrement inherited
+ // MsgDest Definitions: All inherited
//////////////////////////////////////////////////////////////
- static DestFinfo increment( "increment",
- "Increments mol numbers by specified amount. Can be +ve or -ve",
- new OpFunc1< Pool, double >( &Pool::increment )
- );
-
- static DestFinfo decrement( "decrement",
- "Decrements mol numbers by specified amount. Can be +ve or -ve",
- new OpFunc1< Pool, double >( &Pool::decrement )
- );
-
- static DestFinfo nIn( "nIn",
- "Set the number of molecules by specified amount",
- new OpFunc1< Pool, double >( &Pool::nIn )
- );
-
//////////////////////////////////////////////////////////////
// SrcFinfo Definitions: All inherited.
//////////////////////////////////////////////////////////////
@@ -57,9 +42,6 @@ const Cinfo* Pool::initCinfo()
//////////////////////////////////////////////////////////////
static Finfo* poolFinfos[] = {
&isBuffered, // ElementValueFinfo
- &increment, // DestFinfo
- &decrement, // DestFinfo
- &nIn, // DestFinfo
};
static Dinfo< Pool > dinfo;
@@ -125,14 +107,6 @@ void Pool::vProcess( const Eref& e, ProcPtr p )
{
// double A = e.sumBuf( aSlot );
// double B = e.sumBuf( bSlot );
- /*
- if ( n_ < 0 )
- cout << "nugh" << e.objId().path() << endl;
- if ( B_ < 0 )
- cout << "bugh" << e.objId().path() << endl;
- if ( p->dt < 0 )
- cout << "tugh" << e.objId().path() << endl;
- */
if ( n_ > EPSILON && B_ > EPSILON ) {
double C = exp( -B_ * p->dt / n_ );
@@ -162,7 +136,7 @@ void Pool::vReac( double A, double B )
B_ += B;
}
-void Pool::increment( double val )
+void Pool::vIncrement( double val )
{
if ( val > 0 )
A_ += val;
@@ -170,7 +144,7 @@ void Pool::increment( double val )
B_ -= val;
}
-void Pool::decrement( double val )
+void Pool::vDecrement( double val )
{
if ( val < 0 )
A_ -= val;
@@ -178,37 +152,13 @@ void Pool::decrement( double val )
B_ += val;
}
-void Pool::nIn( double val)
+void Pool::vnIn( double val)
{
n_ = val;
B_ = 0;
A_ = 0;
}
-/*
-void Pool::vRemesh( const Eref& e,
- double oldvol,
- unsigned int numTotalEntries, unsigned int startEntry,
- const vector< unsigned int >& localIndices,
- const vector< double >& vols )
-{
- if ( e.dataIndex() != 0 )
- return;
- Neutral* n = reinterpret_cast< Neutral* >( e.data() );
- assert( vols.size() > 0 );
- double concInit = nInit_ / ( NA * oldvol );
- if ( vols.size() != e.element()->dataHandler()->localEntries() )
- n->setLastDimension( e, q, vols.size() );
- // Note that at this point the Pool pointer may be invalid!
- // But we need to update the concs anyway.
- assert( e.element()->dataHandler()->localEntries() == vols.size() );
- Pool* pooldata = reinterpret_cast< Pool* >( e.data() );
- for ( unsigned int i = 0; i < vols.size(); ++i ) {
- pooldata[i].nInit_ = pooldata[i].n_ = concInit * vols[i] * NA;
- }
-}
-*/
-
void Pool::vHandleMolWt( const Eref& e, double v )
{
; // Here I should update DiffConst too.
diff --git a/moose-core/kinetics/Pool.h b/moose-core/kinetics/Pool.h
index dec4dfa4f136057f79545254448ff89fa2d53227..9d0cca01c0c2e07ffaf0e2dc239a55410cc9d3d7 100644
--- a/moose-core/kinetics/Pool.h
+++ b/moose-core/kinetics/Pool.h
@@ -69,13 +69,13 @@ class Pool: public PoolBase
void vProcess( const Eref& e, ProcPtr p );
void vReinit( const Eref& e, ProcPtr p );
void vReac( double A, double B );
+ void vIncrement( double val );
+ void vDecrement( double val );
+ void vnIn( double val );
//////////////////////////////////////////////////////////////////
// Novel Dest funcs not present in Pool base class.
//////////////////////////////////////////////////////////////////
- void increment( double val );
- void decrement( double val );
- void nIn( double val );
//////////////////////////////////////////////////////////////////
static const Cinfo* initCinfo();
diff --git a/moose-core/kinetics/PoolBase.cpp b/moose-core/kinetics/PoolBase.cpp
index 5dd4ea248fedb3ad8930f7c84035ae745b42851e..59c66ba55c726acdb7c5da89ecf5807f87ba8139 100644
--- a/moose-core/kinetics/PoolBase.cpp
+++ b/moose-core/kinetics/PoolBase.cpp
@@ -104,6 +104,24 @@ const Cinfo* PoolBase::initCinfo()
"Should only be used in SharedMsg with species.",
new EpFunc1< PoolBase, double >( &PoolBase::handleMolWt )
);
+ //////////////////////////////////////////////////////////////
+ // MsgDest Definitions: These three are used for non-reaction
+ // calculations involving algebraically defined rate terms.
+ //////////////////////////////////////////////////////////////
+ static DestFinfo increment( "increment",
+ "Increments mol numbers by specified amount. Can be +ve or -ve",
+ new OpFunc1< PoolBase, double >( &PoolBase::increment )
+ );
+
+ static DestFinfo decrement( "decrement",
+ "Decrements mol numbers by specified amount. Can be +ve or -ve",
+ new OpFunc1< PoolBase, double >( &PoolBase::decrement )
+ );
+
+ static DestFinfo nIn( "nIn",
+ "Assigns the number of molecules in Pool to specified value",
+ new OpFunc1< PoolBase, double >( &PoolBase::nIn )
+ );
//////////////////////////////////////////////////////////////
// SrcFinfo Definitions
@@ -155,6 +173,9 @@ const Cinfo* PoolBase::initCinfo()
&concInit, // Value
&volume, // Readonly Value
&speciesId, // Value
+ &increment, // DestFinfo
+ &decrement, // DestFinfo
+ &nIn, // DestFinfo
&reac, // SharedFinfo
&proc, // SharedFinfo
&species, // SharedFinfo
@@ -208,6 +229,21 @@ void PoolBase::reinit( const Eref& e, ProcPtr p )
vReinit( e, p );
}
+void PoolBase::increment( double val )
+{
+ vIncrement(val);
+}
+
+void PoolBase::decrement( double val )
+{
+ vDecrement( val );
+}
+
+void PoolBase::nIn( double val)
+{
+ vnIn(val);
+}
+
void PoolBase::reac( double A, double B )
{
vReac( A, B );
@@ -234,6 +270,15 @@ void PoolBase::vReac( double A, double B )
void PoolBase::vHandleMolWt( const Eref& e, double v )
{;}
+void PoolBase::vIncrement( double val )
+{;}
+
+void PoolBase::vDecrement( double val )
+{;}
+
+void PoolBase::vnIn( double val)
+{;}
+
//////////////////////////////////////////////////////////////
// Field Definitions
//////////////////////////////////////////////////////////////
diff --git a/moose-core/kinetics/PoolBase.h b/moose-core/kinetics/PoolBase.h
index ace7ae86924b2b0680a837501d6aea5ff2942591..8173c4450c7fe5a989872bf870a31b831d92ee09 100644
--- a/moose-core/kinetics/PoolBase.h
+++ b/moose-core/kinetics/PoolBase.h
@@ -126,6 +126,9 @@ class PoolBase
void reinit( const Eref& e, ProcPtr p );
void reac( double A, double B );
void handleMolWt( const Eref& e, double v );
+ void increment( double val );
+ void decrement( double val );
+ void nIn( double val );
//////////////////////////////////////////////////////////////////
// Virtual Dest funcs. Most of these have a generic do-nothing
@@ -136,6 +139,9 @@ class PoolBase
virtual void vReinit( const Eref& e, ProcPtr p );
virtual void vReac( double A, double B );
virtual void vHandleMolWt( const Eref& e, double v);
+ virtual void vIncrement( double val );
+ virtual void vDecrement( double val );
+ virtual void vnIn( double val );
//////////////////////////////////////////////////////////////////
static const Cinfo* initCinfo();
diff --git a/moose-core/python/moose/SBML/readSBML.py b/moose-core/python/moose/SBML/readSBML.py
index 801d3681ec23536b60dce134f2faa67480c70571..6a627054181e6d06e26ed7a4e2db2317cd962b25 100644
--- a/moose-core/python/moose/SBML/readSBML.py
+++ b/moose-core/python/moose/SBML/readSBML.py
@@ -13,7 +13,7 @@
** copyright (C) 2003-2017 Upinder S. Bhalla. and NCBS
Created : Thu May 12 10:19:00 2016(+0530)
Version
-Last-Updated: Wed Jan 11 2017
+Last-Updated: Tue Feb 10 2017
By:
**********************************************************************/
@@ -24,7 +24,7 @@ import os.path
import collections
import moose
from validation import validateModel
-
+import re
'''
TODO in
-Compartment
@@ -62,7 +62,8 @@ def mooseReadSBML(filepath, loadpath, solver="ee"):
if not foundLibSBML_:
print('No python-libsbml found.'
'\nThis module can be installed by following command in terminal:'
- '\n\t easy_install python-libsbl'
+ '\n\t easy_install python-libsbml'
+ '\n\t apt-get install python-libsbml'
)
return None
@@ -116,10 +117,10 @@ def mooseReadSBML(filepath, loadpath, solver="ee"):
return moose.element('/')
else:
baseId = moose.Neutral(loadpath)
+ basePath = baseId
# All the model will be created under model as
# a thumbrule
- basePath = moose.Neutral(
- baseId.path + '/model')
+ basePath = moose.Neutral(baseId.path + '/model')
# Map Compartment's SBML id as key and value is
# list of[ Moose ID and SpatialDimensions ]
global comptSbmlidMooseIdMap
@@ -139,6 +140,7 @@ def mooseReadSBML(filepath, loadpath, solver="ee"):
errorFlag,warning = createSpecies(
basePath, model, comptSbmlidMooseIdMap, specInfoMap, modelAnnotaInfo)
#print(errorFlag,warning)
+
if errorFlag:
errorFlag = createRules(
model, specInfoMap, globparameterIdValue)
@@ -147,6 +149,7 @@ def mooseReadSBML(filepath, loadpath, solver="ee"):
model, specInfoMap, modelAnnotaInfo, globparameterIdValue)
getModelAnnotation(
model, baseId, basePath)
+
if not errorFlag:
print(msg)
# Any time in the middle if SBML does not read then I
@@ -170,7 +173,6 @@ def setupEnzymaticReaction(enz, groupName, enzName,
cplx = (modelAnnotaInfo[groupName]["complex"])
cplx = str(idBeginWith(cplx))
complx = moose.element(specInfoMap[cplx]["Mpath"].path)
-
enzyme_ = moose.Enz(enzParent.path + '/' + enzName)
moose.move(complx, enzyme_)
moose.connect(enzyme_, "cplx", complx, "reac")
@@ -286,7 +288,7 @@ def getModelAnnotation(obj, baseId, basepath):
for plots in plotlist:
plots = plots.replace(" ", "")
plotorg = plots
- if moose.exists(basepath.path + plotorg):
+ if( moose.exists(basepath.path + plotorg) and isinstance(moose.element(basepath.path+plotorg),moose.PoolBase)) :
plotSId = moose.element(
basepath.path + plotorg)
# plotorg = convertSpecialChar(plotorg)
@@ -301,8 +303,7 @@ def getModelAnnotation(obj, baseId, basepath):
if not fullPath in tablelistname:
tab = moose.Table2(fullPath)
tablelistname.append(fullPath)
- moose.connect(
- tab, "requestOut", plotSId, "getConc")
+ moose.connect(tab, "requestOut", plotSId, "getConc")
def getObjAnnotation(obj, modelAnnotationInfo):
@@ -335,6 +336,8 @@ def getObjAnnotation(obj, modelAnnotationInfo):
annotateMap[nodeName] = nodeValue
if nodeName == "textColor":
annotateMap[nodeName] = nodeValue
+ if nodeName == "Group":
+ annotateMap[nodeName] = nodeValue
return annotateMap
@@ -459,7 +462,14 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue):
for ritem in range(0, model.getNumReactions()):
reactionCreated = False
groupName = ""
+
reac = model.getReaction(ritem)
+ group = ""
+ reacAnnoInfo = {}
+ reacAnnoInfo = getObjAnnotation(reac, modelAnnotaInfo)
+ if "Group" in reacAnnoInfo:
+ group = reacAnnoInfo["Group"]
+
if (reac.isSetId()):
rId = reac.getId()
if (reac.isSetName()):
@@ -529,6 +539,11 @@ def createReaction(model, specInfoMap, modelAnnotaInfo, globparameterIdValue):
sp = react.getSpecies()
sp = str(idBeginWith(sp))
speCompt = specInfoMap[sp]["comptId"].path
+ if group:
+ if moose.exists(speCompt+'/'+group):
+ speCompt = speCompt+'/'+group
+ else:
+ speCompt = (moose.Neutral(speCompt+'/'+group)).path
reaction_ = moose.Reac(speCompt + '/' + rName)
reactionCreated = True
reactSBMLIdMooseId[rName] = {
@@ -758,10 +773,12 @@ def createRules(model, specInfoMap, globparameterIdValue):
exp = rule.getFormula()
for mem in ruleMemlist:
if (mem in specInfoMap):
- exp1 = exp.replace(mem, str(speFunXterm[mem]))
+ #exp1 = exp.replace(mem, str(speFunXterm[mem]))
+ exp1 = re.sub(r'\b%s\b'% (mem), speFunXterm[mem], exp)
exp = exp1
elif(mem in globparameterIdValue):
- exp1 = exp.replace(mem, str(globparameterIdValue[mem]))
+ #exp1 = exp.replace(mem, str(globparameterIdValue[mem]))
+ exp1 = re.sub(r'\b%s\b'% (mem), globparameterIdValue[mem], exp)
exp = exp1
else:
print("Math expression need to be checked")
@@ -810,6 +827,12 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap,
else:
for sindex in range(0, model.getNumSpecies()):
spe = model.getSpecies(sindex)
+ group = ""
+ specAnnoInfo = {}
+ specAnnoInfo = getObjAnnotation(spe, modelAnnotaInfo)
+ if "Group" in specAnnoInfo:
+ group = specAnnoInfo["Group"]
+
sName = None
sId = spe.getId()
if spe.isSetName():
@@ -828,7 +851,11 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap,
hasonlySubUnit = spe.getHasOnlySubstanceUnits()
# "false": is {unit of amount}/{unit of size} (i.e., concentration or density).
# "true": then the value is interpreted as having a unit of amount only.
-
+ if group:
+ if moose.exists(comptEl+'/'+group):
+ comptEl = comptEl+'/'+group
+ else:
+ comptEl = (moose.Neutral(comptEl+'/'+group)).path
if (boundaryCondition):
poolId = moose.BufPool(comptEl + '/' + sName)
else:
@@ -836,13 +863,13 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap,
if (spe.isSetNotes):
pullnotes(spe, poolId)
- specAnnoInfo = {}
- specAnnoInfo = getObjAnnotation(spe, modelAnnotaInfo)
+
if specAnnoInfo:
if not moose.exists(poolId.path + '/info'):
poolInfo = moose.Annotator(poolId.path + '/info')
else:
poolInfo = moose.element(poolId.path + '/info')
+
for k, v in list(specAnnoInfo.items()):
if k == 'xCord':
poolInfo.x = float(v)
@@ -904,6 +931,7 @@ def createSpecies(basePath, model, comptSbmlidMooseIdMap,
if (rule_variable == sId):
found = True
break
+
if not (found):
print(
"Invalid SBML: Either initialConcentration or initialAmount must be set or it should be found in assignmentRule but non happening for ",
diff --git a/moose-core/python/moose/SBML/writeSBML.py b/moose-core/python/moose/SBML/writeSBML.py
index bcfaf2eb2824b52337fd8b058f7b2dd5ba4d5fc8..90400306877ebe14cb5fccc4e7c42211d648ae20 100644
--- a/moose-core/python/moose/SBML/writeSBML.py
+++ b/moose-core/python/moose/SBML/writeSBML.py
@@ -47,9 +47,10 @@ def mooseWriteSBML(modelpath, filename, sceneitems={}):
if not foundLibSBML_:
print('No python-libsbml found.'
'\nThis module can be installed by following command in terminal:'
- '\n\t easy_install python-libsbml'
+ '\n\t easy_install python-libsbml or'
+ '\n\t apt-get install python-libsbml'
)
- return -1, msg
+ return -2, "Could not save the model in to SBML file. \nThis module can be installed by following command in terminal: \n\t easy_install python-libsbml or \n\t apt-get install python-libsbml",''
sbmlDoc = SBMLDocument(3, 1)
filepath, filenameExt = os.path.split(filename)
@@ -876,12 +877,12 @@ def writeSimulationAnnotation(modelpath):
graphSpefound = True
if graphSpefound:
if not plots:
- #plots = ori[ori.find(q.name)-1:len(ori)]
- plots = '/' + q.name + '/' + name
+ plots = ori[ori.find(q.name)-1:len(ori)]
+ #plots = '/' + q.name + '/' + name
else:
- #plots = plots + "; "+ori[ori.find(q.name)-1:len(ori)]
- plots = plots + "; /" + q.name + '/' + name
+ 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 ba71603148d5e5096a41e14499027a7a11b442c2..2bd33e212590cad4216ca54e733b45bde7ad60b3 100644
--- a/moose-core/python/moose/genesis/writeKkit.py
+++ b/moose-core/python/moose/genesis/writeKkit.py
@@ -8,15 +8,21 @@ __version__ = "1.0.0"
__maintainer__ = "Harsha Rani"
__email__ = "hrani@ncbs.res.in"
__status__ = "Development"
-
+__updated__ = "Feb 13 2017"
import sys
import random
import re
-import matplotlib
import moose
from moose.chemUtil.chemConnectUtil import *
from moose.chemUtil.graphUtils import *
+foundmatplotlib_ = False
+try:
+ import matplotlib
+ foundmatplotlib_ = True
+except Exception as e:
+ pass
+
GENESIS_COLOR_SEQUENCE = ((248, 0, 255), (240, 0, 255), (232, 0, 255), (224, 0, 255), (216, 0, 255), (208, 0, 255),
(200, 0, 255), (192, 0, 255), (184, 0, 255), (176, 0, 255), (168, 0, 255), (160, 0, 255), (152, 0, 255), (144, 0, 255),
(136, 0, 255), (128, 0, 255), (120, 0, 255), (112, 0, 255), (104, 0, 255), (96, 0, 255), (88, 0, 255), (80, 0, 255),
@@ -38,78 +44,91 @@ GENESIS_COLOR_SEQUENCE = ((248, 0, 255), (240, 0, 255), (232, 0, 255), (224, 0,
# --StimulusTable
def mooseWriteKkit( modelpath, filename, sceneitems={}):
- if filename.rfind('.') != -1:
- filename = filename[:filename.rfind('.')]
- else:
- filename = filename[:len(filename)]
- filename = filename+'.g'
- global NA
- NA = 6.0221415e23
- global cmin,cmax,xmin,xmax,ymin,ymax
- cmin, xmin, ymin = 0, 0, 0
- cmax, xmax, ymax = 1, 1, 1
-
- compt = moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]')
- maxVol = estimateDefaultVol(compt)
- positionInfoExist = True
- if compt:
- if bool(sceneitems):
- cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems)
- elif not bool(sceneitems):
- srcdesConnection = {}
- setupItem(modelpath,srcdesConnection)
- meshEntry,xmin,xmax,ymin,ymax,positionInfoExist,sceneitems = setupMeshObj(modelpath)
- if not positionInfoExist:
- #cmin,cmax,sceneitems = autoCoordinates(meshEntry,srcdesConnection)
- sceneitems = autoCoordinates(meshEntry,srcdesConnection)
-
- if not positionInfoExist:
- # if position are not from kkit, then zoom factor is applied while
- # writing to genesis. Like if position is from pyqtSceneItem or auto-coordinates
- cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems)
- for k,v in list(sceneitems.items()):
- anno = moose.element(k.path+'/info')
- x1 = calPrime(v['x'])
- y1 = calPrime(v['y'])
- sceneitems[k]['x'] = x1
- sceneitems[k]['y'] = y1
-
- f = open(filename, 'w')
- writeHeader (f,maxVol)
-
- gtId_vol = writeCompartment(modelpath,compt,f)
- writePool(modelpath,f,gtId_vol,sceneitems)
- reacList = writeReac(modelpath,f,sceneitems)
- enzList = writeEnz(modelpath,f,sceneitems)
- 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]')
- first, second = " ", " "
- if tgraphs:
- first,second = writeplot(tgraphs,f)
- if first:
- f.write(first)
- f.write("simundump xgraph /moregraphs/conc3 0 0 100 0 1 0\n"
- "simundump xgraph /moregraphs/conc4 0 0 100 0 1 0\n")
- if second:
- f.write(second)
- f.write("simundump xcoredraw /edit/draw 0 -6 4 -2 6\n"
- "simundump xtree /edit/draw/tree 0 \\\n"
- " /kinetics/#[],/kinetics/#[]/#[],/kinetics/#[]/#[]/#[][TYPE!=proto],/kinetics/#[]/#[]/#[][TYPE!=linkinfo]/##[] \"edit_elm.D <v>; drag_from_edit.w <d> <S> <x> <y> <z>\" auto 0.6\n"
- "simundump xtext /file/notes 0 1\n")
- storeReacMsg(reacList,f)
- storeEnzMsg(enzList,f)
- if tgraphs:
- storePlotMsgs(tgraphs,f)
- writeFooter1(f)
- writeNotes(modelpath,f)
- writeFooter2(f)
- print('Written to file '+filename)
- return True
- else:
- print(("Warning: writeKkit:: No model found on " , modelpath))
+ global foundmatplotlib_
+ if not foundmatplotlib_:
+ print('No maplotlib found.'
+ '\nThis module can be installed by following command in terminal:'
+ '\n\t sudo apt install python-maplotlib', "")
return False
+ else:
+
+ ignoreColor= ["mistyrose","antiquewhite","aliceblue","azure","bisque","black","blanchedalmond","blue","cornsilk","darkolivegreen","darkslategray","dimgray","floralwhite","gainsboro","ghostwhite","honeydew","ivory","lavender","lavenderblush","lemonchiffon","lightcyan","lightgoldenrodyellow","lightgray","lightyellow","linen","mediumblue","mintcream","navy","oldlace","papayawhip","saddlebrown","seashell","snow","wheat","white","whitesmoke","aquamarine","lightsalmon","moccasin","limegreen","snow","sienna","beige","dimgrey","lightsage"]
+ matplotcolor = {}
+ for name,hexno in matplotlib.colors.cnames.items():
+ matplotcolor[name]=hexno
+
+ if filename.rfind('.') != -1:
+ filename = filename[:filename.rfind('.')]
+ else:
+ filename = filename[:len(filename)]
+ filename = filename+'.g'
+ global NA
+ NA = 6.0221415e23
+ global cmin,cmax,xmin,xmax,ymin,ymax
+ cmin, xmin, ymin = 0, 0, 0
+ cmax, xmax, ymax = 1, 1, 1
+
+ compt = moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]')
+ maxVol = estimateDefaultVol(compt)
+ positionInfoExist = True
+ if compt:
+ if bool(sceneitems):
+ cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems)
+ elif not bool(sceneitems):
+ srcdesConnection = {}
+ setupItem(modelpath,srcdesConnection)
+ meshEntry,xmin,xmax,ymin,ymax,positionInfoExist,sceneitems = setupMeshObj(modelpath)
+ if not positionInfoExist:
+ #cmin,cmax,sceneitems = autoCoordinates(meshEntry,srcdesConnection)
+ sceneitems = autoCoordinates(meshEntry,srcdesConnection)
+
+ if not positionInfoExist:
+ # if position are not from kkit, then zoom factor is applied while
+ # writing to genesis. Like if position is from pyqtSceneItem or auto-coordinates
+ cmin,cmax,xmin1,xmax1,ymin1,ymax1 = findMinMax(sceneitems)
+ for k,v in list(sceneitems.items()):
+ anno = moose.element(k.path+'/info')
+ x1 = calPrime(v['x'])
+ y1 = calPrime(v['y'])
+ sceneitems[k]['x'] = x1
+ sceneitems[k]['y'] = y1
+
+ f = open(filename, 'w')
+ writeHeader (f,maxVol)
+
+ gtId_vol = writeCompartment(modelpath,compt,f)
+ writePool(modelpath,f,gtId_vol,sceneitems)
+ reacList = writeReac(modelpath,f,sceneitems)
+ enzList = writeEnz(modelpath,f,sceneitems)
+ 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]')
+ first, second = " ", " "
+ if tgraphs:
+ first,second = writeplot(tgraphs,f)
+ if first:
+ f.write(first)
+ f.write("simundump xgraph /moregraphs/conc3 0 0 100 0 1 0\n"
+ "simundump xgraph /moregraphs/conc4 0 0 100 0 1 0\n")
+ if second:
+ f.write(second)
+ f.write("simundump xcoredraw /edit/draw 0 -6 4 -2 6\n"
+ "simundump xtree /edit/draw/tree 0 \\\n"
+ " /kinetics/#[],/kinetics/#[]/#[],/kinetics/#[]/#[]/#[][TYPE!=proto],/kinetics/#[]/#[]/#[][TYPE!=linkinfo]/##[] \"edit_elm.D <v>; drag_from_edit.w <d> <S> <x> <y> <z>\" auto 0.6\n"
+ "simundump xtext /file/notes 0 1\n")
+ storeReacMsg(reacList,f)
+ storeEnzMsg(enzList,f)
+ if tgraphs:
+ storePlotMsgs(tgraphs,f)
+ writeFooter1(f)
+ writeNotes(modelpath,f)
+ writeFooter2(f)
+ print('Written to file '+filename)
+ return True
+ else:
+ print(("Warning: writeKkit:: No model found on " , modelpath))
+ return False
def findMinMax(sceneitems):
cmin = 0.0
@@ -484,11 +503,6 @@ def getColorCheck(color,GENESIS_COLOR_SEQUENCE):
else:
raise Exception("Invalid Color Value!")
-ignoreColor= ["mistyrose","antiquewhite","aliceblue","azure","bisque","black","blanchedalmond","blue","cornsilk","darkolivegreen","darkslategray","dimgray","floralwhite","gainsboro","ghostwhite","honeydew","ivory","lavender","lavenderblush","lemonchiffon","lightcyan","lightgoldenrodyellow","lightgray","lightyellow","linen","mediumblue","mintcream","navy","oldlace","papayawhip","saddlebrown","seashell","snow","wheat","white","whitesmoke","aquamarine","lightsalmon","moccasin","limegreen","snow","sienna","beige","dimgrey","lightsage"]
-matplotcolor = {}
-for name,hexno in matplotlib.colors.cnames.items():
- matplotcolor[name]=hexno
-
def getRandColor():
k = random.choice(list(matplotcolor.keys()))
if k in ignoreColor:
@@ -617,13 +631,19 @@ def writeFooter2(f):
if __name__ == "__main__":
import sys
-
+ import os
filename = sys.argv[1]
- modelpath = filename[0:filename.find('.')]
+ filepath, filenameWithext = os.path.split(filename)
+ if filenameWithext.find('.') != -1:
+ modelpath = filenameWithext[:filenameWithext.find('.')]
+ else:
+ modelpath = filenameWithext
+
moose.loadModel(filename,'/'+modelpath,"gsl")
output = modelpath+"_.g"
- written = write('/'+modelpath,output)
+ written = mooseWriteKkit('/'+modelpath,output)
+
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/merge/merge.py b/moose-core/python/moose/merge/merge.py
index f76ca494b3daffccc1fa6318dd9783564c3a09d1..02ecc805775f0fc54f4a1788b8cf8a88585f432c 100644
--- a/moose-core/python/moose/merge/merge.py
+++ b/moose-core/python/moose/merge/merge.py
@@ -33,40 +33,55 @@ from moose.chemUtil.graphUtils import *
def mergeChemModel(A,B):
""" Merges two model or the path """
-
- modelA,loadedA = loadModels(A)
- modelB,loadedB = loadModels(B)
- if not loadedA or not loadedB:
- if not loadedA:
- modelB = moose.Shell('/')
- if not loadedB:
- modelA = moose.Shell('/')
+ loadedA = False
+ loadedB = False
+
+ if os.path.isfile(A):
+ modelA,loadedA = loadModels(A)
+ elif moose.exists(A):
+ modelA = A
+ loadedA = True
+ else:
+ print ("%s path or file doesnot exists. Mergering will exist" % (A))
+ exit(0)
+
+ if os.path.isfile(B):
+ modelB,loadedB = loadModels(B)
+ elif moose.exists(B):
+ modelB = B
+ loadedB = True
else:
- directory, bfname = os.path.split(B)
- global grpNotcopiedyet,poolListina
+ print ("%s path or file doesnot exists. Mergering will exist " % (B))
+ exit(0)
+
+ if loadedA and loadedB:
+ ## yet deleteSolver is called to make sure all the moose object are off from solver
+ deleteSolver(modelA)
+ deleteSolver(modelB)
+
+ global poolListina
poolListina = {}
grpNotcopiedyet = []
dictComptA = dict( [ (i.name,i) for i in moose.wildcardFind(modelA+'/##[ISA=ChemCompt]') ] )
dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] )
-
poolNotcopiedyet = []
- for key in list(dictComptB.keys()):
- if key not in dictComptA:
+ for key in list(dictComptA.keys()):
+ if key not in dictComptB:
# if compartmentname from modelB does not exist in modelA, then copy
- copy = moose.copy(dictComptB[key],moose.element(modelA))
+ copy = moose.copy(dictComptA[key],moose.element(modelB))
else:
#if compartmentname from modelB exist in modelA,
#volume is not same, then change volume of ModelB same as ModelA
- if abs(dictComptA[key].volume - dictComptB[key].volume):
+ if abs(dictComptB[key].volume - dictComptA[key].volume):
#hack for now
- while (abs(dictComptA[key].volume - dictComptB[key].volume) != 0.0):
- dictComptB[key].volume = float(dictComptA[key].volume)
- dictComptA = dict( [ (i.name,i) for i in moose.wildcardFind(modelA+'/##[ISA=ChemCompt]') ] )
+ while (abs(dictComptB[key].volume - dictComptA[key].volume) != 0.0):
+ dictComptA[key].volume = float(dictComptB[key].volume)
+ dictComptB = dict( [ (i.name,i) for i in moose.wildcardFind(modelB+'/##[ISA=ChemCompt]') ] )
#Mergering pool
- poolMerge(dictComptA[key],dictComptB[key],poolNotcopiedyet)
-
+ poolMerge(dictComptB[key],dictComptA[key],poolNotcopiedyet)
+
if grpNotcopiedyet:
# objA = moose.element(comptA).parent.name
# if not moose.exists(objA+'/'+comptB.name+'/'+bpath.name):
@@ -74,29 +89,30 @@ def mergeChemModel(A,B):
# moose.copy(bpath,moose.element(objA+'/'+comptB.name))
pass
- comptAdict = comptList(modelA)
- poolListina = {}
- poolListina = updatePoolList(comptAdict)
- funcNotallowed = []
+ comptBdict = comptList(modelB)
+ poolListinb = {}
+ poolListinb = updatePoolList(comptBdict)
R_Duplicated, R_Notcopiedyet,R_Daggling = [], [], []
E_Duplicated, E_Notcopiedyet,E_Daggling = [], [], []
- for key in list(dictComptB.keys()):
- funcNotallowed = []
- funcNotallowed = functionMerge(dictComptA,dictComptB,key)
+ for key in list(dictComptA.keys()):
+ funcExist, funcNotallowed = [], []
+ funcExist,funcNotallowed = functionMerge(dictComptB,dictComptA,key)
- poolListina = updatePoolList(dictComptA)
- R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptA,dictComptB,key,poolListina)
-
- poolListina = updatePoolList(dictComptA)
- E_Duplicated,E_Notcopiedyet,E_Daggling = enzymeMerge(dictComptA,dictComptB,key,poolListina)
-
- print("\n Model is merged to %s" %modelA)
+ poolListinb = updatePoolList(dictComptB)
+ R_Duplicated,R_Notcopiedyet,R_Daggling = reacMerge(dictComptB,dictComptA,key,poolListinb)
+
+ poolListinb = updatePoolList(dictComptB)
+ E_Duplicated,E_Notcopiedyet,E_Daggling = enzymeMerge(dictComptB,dictComptA,key,poolListinb)
- if funcNotallowed:
+ print("\n Model is merged to %s" %modelB)
+ if funcExist:
print( "\nPool already connected to a function, this function is not to connect to same pool, since no two function are allowed to connect to same pool:")
+ for fl in list(funcExist):
+ print("\t [Pool]: %s [Function]: %s \n" %(str(fl.parent.name), str(fl.path)))
+ if funcNotallowed:
+ print( "\nThis function is not to copied, since pool connected to function input are from different compartment:")
for fl in list(funcNotallowed):
print("\t [Pool]: %s [Function]: %s \n" %(str(fl.parent.name), str(fl.path)))
-
if R_Duplicated or E_Duplicated:
print ("Reaction / Enzyme are Duplicate"
"\n 1. The once whoes substrate / product names are different for a give reaction name "
@@ -134,62 +150,50 @@ def mergeChemModel(A,B):
print ("Enzyme:")
for ed in list(E_Daggling):
print ("%s " %str(ed.name))
-
+
def functionMerge(comptA,comptB,key):
- funcNotallowed = []
+ funcNotallowed, funcExist = [], []
comptApath = moose.element(comptA[key]).path
comptBpath = moose.element(comptB[key]).path
- funcListina = moose.wildcardFind(comptApath+'/##[ISA=PoolBase]')
- funcListinb = moose.wildcardFind(comptBpath+'/##[ISA=Function]')
objA = moose.element(comptApath).parent.name
objB = moose.element(comptBpath).parent.name
- #For function mergering pool name is taken as reference
- funcNotcopiedyet = []
+
+ #This will give us all the function which exist in modelB
+ funcListinb = moose.wildcardFind(comptBpath+'/##[ISA=Function]')
for fb in funcListinb:
-
- if fb.parent.className in ['ZombiePool','Pool','ZombieBufPool','BufPool']:
- objA = moose.element(comptApath).parent.name
- fbpath = fb.path
- #funcpath = fbpath[fbpath.find(findCompartment(fb).name)-1:len(fbpath)]
- funcparentB = fb.parent.path
- funcpath = fbpath.replace(objB,objA)
- funcparentA = funcparentB.replace(objB,objA)
- tt = moose.element(funcparentA).neighbors['setN']
- if tt:
- funcNotallowed.append(fb)
- else:
- if len(moose.element(fb.path+'/x').neighbors["input"]):
- #inputB = moose.element(fb.path+'/x').neighbors["input"]
- inputB = subprdList(moose.element(fb.path+'/x'),"input")
- inputB_expr = fb.expr
- if moose.exists(funcpath):
- #inputA = moose.element(objA+funcpath+'/x').neighbors["input"]
- inputA = subprdList(moose.element(funcpath+'/x'),"input")
- inputA_expr = moose.element(funcpath).expr
- hassameExpr = False
- if inputA_expr == inputB_expr:
- hassameExpr = True
- hassameLen,hassameS,hassameVols = same_len_name_vol(inputA,inputB)
- if not all((hassameLen,hassameS,hassameVols,hassameExpr)):
- fb.name = fb.name+'_duplicatedF'
- createFunction(fb,inputB,objB,objA)
+ #This will give us all the pools that its connected to, for this function
+ fvalueOut = moose.element(fb).neighbors['valueOut']
+ for poolinB in fvalueOut:
+ poolinBpath = poolinB.path
+ poolinA = poolinBpath.replace(objB,objA)
+ connectionexist = []
+ if moose.exists(poolinA):
+ #This is give us if pool which is to be connected already exist any connection
+ connectionexist = moose.element(poolinA).neighbors['setN']+moose.element(poolinA).neighbors['setConc']+ moose.element(poolinA).neighbors['increment']
+ if len(connectionexist) == 0:
+ #This pool in model A doesnot exist with any function
+ inputs = moose.element(fb.path+'/x').neighbors['input']
+ volumes = []
+ for ins in inputs:
+ volumes.append((findCompartment(moose.element(ins))).volume)
+ if len(set(volumes)) == 1:
+ # If all the input connected belongs to one compartment then copy
+ createFunction(fb,poolinA,objB,objA)
else:
- #function doesnot exist then copy
- if len(inputB):
- volinput = []
- for inb in inputB:
- volinput.append(findCompartment(moose.element(inb)).volume)
- if len(set(volinput)) == 1:
- # If all the input connected belongs to one compartment then copy
- createFunction(fb,inputB,objB,objA)
- else:
- # moose doesn't allow function input to come from different compartment
- funcNotallowed.append(fb)
- return funcNotallowed
-
-def createFunction(fb,inputB,objB,objA):
+ # moose doesn't allow function's input to come from different compartment
+ funcNotallowed.append(fb)
+ else:
+ #Pool in model 'A' already exist function "
+ funcExist.append(fb)
+ else:
+ print(" Path in model A doesn't exists %s" %(poolinA))
+
+ return funcExist,funcNotallowed
+
+def createFunction(fb,setpool,objB,objA):
fapath1 = fb.path.replace(objB,objA)
fapath = fapath1.replace('[0]','')
+
if not moose.exists(fapath):
# if fb.parent.className in ['CubeMesh','CyclMesh']:
# des = moose.Function('/'+objA+'/'+fb.parent.name+'/'+fb.name)
@@ -198,7 +202,13 @@ def createFunction(fb,inputB,objB,objA):
# if fb.parent.name == akey.name:
# des = moose.Function(akey.path+'/'+fb.name)
des = moose.Function(fapath)
- moose.connect(des, 'valueOut', moose.element(fapath).parent,'setN' )
+ else:
+ des = moose.element(fapath)
+ inputB = moose.element(fb.path+'/x').neighbors["input"]
+ moose.connect(des, 'valueOut', moose.element(setpool),'setN' )
+ inputA = []
+ inputA = moose.element(fapath+'/x').neighbors["input"]
+ if not inputA:
for src in inputB:
pool = ((src.path).replace(objB,objA)).replace('[0]','')
numVariables = des.numVars
@@ -209,46 +219,45 @@ def createFunction(fb,inputB,objB,objA):
des.expr = expr
moose.connect( pool, 'nOut', des.x[numVariables], 'input' )
- #if fb.expr != des.expr:
- # print "Function ",des, " which is duplicated from modelB, expression is different, this is tricky in moose to know what those constants are connected to "
- # print "ModelB ", fb, fb.expr, "\nModelA ",des, des.expr
-
def comptList(modelpath):
comptdict = {}
for ca in moose.wildcardFind(modelpath+'/##[ISA=ChemCompt]'):
comptdict[ca.name] = ca
return comptdict
-def loadModels(filename):
+def loadModels(filepath):
""" load models into moose if file, if moosepath itself it passes back the path and
delete solver if exist """
modelpath = '/'
loaded = False
- if os.path.isfile(filename) :
- modelpath = filename[filename.rfind('/'): filename.rfind('.')]
- ext = os.path.splitext(filename)[1]
- filename = filename.strip()
- modeltype = mtypes.getType(filename)
- subtype = mtypes.getSubtype(filename, modeltype)
+ if os.path.isfile(filepath) :
+ fpath, filename = os.path.split(filepath)
+ # print " head and tail ",head, " ",tail
+ # modelpath = filename[filename.rfind('/'): filename.rfind('.')]
+ # print "modelpath ",modelpath
+ # ext = os.path.splitext(filename)[1]
+ # filename = filename.strip()
+ modelpath = '/'+filename[:filename.rfind('.')]
+ modeltype = mtypes.getType(filepath)
+ subtype = mtypes.getSubtype(filepath, modeltype)
+
if subtype == 'kkit' or modeltype == "cspace":
- moose.loadModel(filename,modelpath)
+ moose.loadModel(filepath,modelpath)
loaded = True
elif subtype == 'sbml':
- #moose.ReadSBML()
+ #moose.mooseReadSBML(filename,modelpath)
+ #loaded = True
pass
else:
print("This file is not supported for mergering")
modelpath = moose.Shell('/')
- elif moose.exists(filename):
- modelpath = filename
+
+ elif moose.exists(filepath):
+ modelpath = filepath
loaded = True
- ## default is 'ee' solver while loading the model using moose.loadModel,
- ## yet deleteSolver is called just to be assured
- if loaded:
- deleteSolver(modelpath)
return modelpath,loaded
@@ -300,15 +309,15 @@ def copy_deleteUnlyingPoolObj(pool,path):
# which will automatically copie's the pool
copied = False
- if pool.parent.className not in ["Enz","ZombieEnz"]:
+ if pool.parent.className not in ["Enz","ZombieEnz","MMenz","ZombieMMenz"]:
poolcopied = moose.copy(pool,path)
copied = True
# deleting function and enzyme which gets copied if exist under pool
# This is done to ensure daggling function / enzyme not copied.
funclist = []
for types in ['setConc','setN','increment']:
-
funclist.extend(moose.element(poolcopied).neighbors[types])
+
for fl in funclist:
moose.delete(fl)
enzlist = moose.element(poolcopied).neighbors['reac']
@@ -364,20 +373,15 @@ def enzymeMerge(comptA,comptB,key,poolListina):
allclean = True
else:
# didn't find sub or prd for this Enzyme
- #print ("didn't find sub or prd for this reaction" )
RE_Notcopiedyet.append(eb)
else:
# -- it is dagging reaction
RE_Daggling.append(eb)
- #print ("This reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied")
- #war_msg = war_msg+"\nThis reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied"
else:
#Same Enzyme name
# -- Same substrate and product including same volume then don't copy
# -- different substrate/product or if sub/prd's volume is different then DUPLICATE the Enzyme
allclean = False
- #ea = moose.element('/'+obj+'/'+enzcompartment.name+'/'+enzparent.name+'/'+eb.name)
- #ea = moose.element(pA.path+'/'+eb.name)
ea = moose.element(eb.path.replace(objB,objA))
eAsubname = subprdList(ea,"sub")
eBsubname = subprdList(eb,"sub")
@@ -407,7 +411,6 @@ def enzymeMerge(comptA,comptB,key,poolListina):
moose.connect(moose.element(enz).parent,"nOut",moose.element(enz),"enzDest")
#moose.connect(moose.element(enz),"enz",moose.element(enz).parent,"reac")
- #moose.connect( cplxItem, 'reac', enz, 'cplx' )
connectObj(enz,eBsubname,"sub",comptA,war_msg)
connectObj(enz,eBprdname,"prd",comptA,war_msg)
RE_Duplicated.append(enz)
@@ -422,12 +425,9 @@ def enzymeMerge(comptA,comptB,key,poolListina):
# -- it may be connected Enzyme cplx
if eBsubname and eBprdname:
RE_Notcopiedyet.append(eb)
- #print ("This Enzyme \""+eb.path+"\" has no substrate/product must be connect to cplx")
- #war_msg = war_msg+ "\nThis Enzyme \""+rb.path+"\" has no substrate/product must be connect to cplx"
else:
RE_Daggling.append(eb)
- #print ("This enzyme \""+eb.path+"\" has no substrate/product daggling reaction are not copied")
- #war_msg = war_msg+"\nThis reaction \""+eb.path+"\" has no substrate/product daggling reaction are not copied"
+
return RE_Duplicated,RE_Notcopiedyet,RE_Daggling
def reacMerge(comptA,comptB,key,poolListina):
@@ -511,17 +511,12 @@ def reacMerge(comptA,comptB,key,poolListina):
# -- it may be connected Enzyme cplx
if rBsubname and rBprdname:
RE_Notcopiedyet.append(rb)
- #print ("This reaction \""+rb.path+"\" has no substrate/product must be connect to cplx")
- #war_msg = war_msg+ "\nThis reaction \""+rb.path+"\" has no substrate/product must be connect to cplx"
else:
RE_Daggling.append(rb)
- #print ("This reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied")
- #war_msg = war_msg+"\nThis reaction \""+rb.path+"\" has no substrate/product daggling reaction are not copied"
return RE_Duplicated,RE_Notcopiedyet,RE_Daggling
def subprdList(reac,subprd):
- #print "Reac ",reac
rtype = moose.element(reac).neighbors[subprd]
rname = []
for rs in rtype:
@@ -572,8 +567,6 @@ def connectObj(reac,spList,spType,comptA,war_msg):
allclean = True
else:
#It should not come here unless the sub/prd is connected to enzyme cplx pool
- #print ("This pool \""+rsp.name+"\" doesnot exists in this "+comptName+" compartment to connect to this reaction \""+reac.name+"\"")
- #war_msg = war_msg+ "This pool \""+rsp.name+"\" doesnot exists in this "+comptName+" compartment to connect to this reaction \""+reac.name+"\""
allclean = False
return allclean
@@ -604,6 +597,6 @@ def mooseIsInstance(element, classNames):
if __name__ == "__main__":
- modelA = '/home/harsha/genesis_files/gfile/acc92.g'
- modelB = '/home/harsha/genesis_files/gfile/acc50.g'
+ modelA = 'acc50.g'
+ modelB = 'acc92.g'
mergered = mergeChemModel(modelA,modelB)
diff --git a/moose-core/python/moose/moose.py b/moose-core/python/moose/moose.py
index 904d6c32e19f45d0cefe06e07d6f3854ca7699d0..1309c68fc838d4b225e75cad82c3347232c932fe 100644
--- a/moose-core/python/moose/moose.py
+++ b/moose-core/python/moose/moose.py
@@ -17,30 +17,9 @@ from collections import defaultdict
from . import _moose
from ._moose import *
import __main__ as main
-
-sbmlSupport_, genesisSupport_ = True, True
-try:
- import SBML.readSBML
- import SBML.writeSBML
-except Exception as e:
- print( 'MOOSE could not load SBML support due to %s' % e )
- sbmlSupport_ = False
-
-try:
- import genesis.writeKkit
-except Exception as e:
- print('MOOSE could not load GENESIS support')
- print('\tError was %s' % e)
- genesisSupport_ = False
-
-chemUtilSupport_ = True
-try:
- import chemUtil.add_Delete_ChemicalSolver
-except Exception as e:
- chemUtilSupport_ = False
- print( 'Failed to import utility module chemUtil' )
- print( '\tError was %s' % e )
-
+import genesis.writeKkit
+import SBML.readSBML
+import SBML.writeSBML
sequence_types = ['vector<double>',
'vector<int>',
'vector<long>',
@@ -77,17 +56,10 @@ def mooseReadSBML(filepath, loadpath, solver='ee'):
solver -- Solver to use (default 'ee' ) \n
"""
- global sbmlSupport_
- if not sbmlSupport_:
- print('SBML support was not loaded')
- return None
- if not os.path.isfile(filepath):
- raise UserWarning('File %s not found' % filepath)
-
return SBML.readSBML.mooseReadSBML( filepath, loadpath, solver )
-def mooseWriteSBML(modelpath, filenpath, sceneitems={}):
+def mooseWriteSBML(modelpath, filepath, sceneitems={}):
"""Writes loaded model under modelpath to a file in SBML format.
keyword arguments:\n
@@ -103,12 +75,6 @@ def mooseWriteSBML(modelpath, filenpath, sceneitems={}):
--- else, auto-coordinates is used for layout position and passed
"""
-
- global sbmlSupport_
- if not sbmlSupport_:
- print('SBML support was not loaded')
- return None
-
return SBML.writeSBML.mooseWriteSBML(modelpath, filepath, sceneitems)
@@ -120,11 +86,6 @@ def mooseWriteKkit(modelpath, filepath,sceneitems={}):
modelpath -- model path in moose \n
filepath -- Path of output file.
"""
- global genesisSupport_
- if not genesisSupport_:
- print('GENESIS(kkit) support was not loaded')
- return None
-
return genesis.writeKkit.mooseWriteKkit(modelpath, filepath,sceneitems)
diff --git a/moose-core/python/setup.py b/moose-core/python/setup.py
index 47b5f3aa0d06da2ac630645b00ce7e9cc005e5c5..a6c3852c787d09e3d59d8b03ba82c359759e8f11 100644
--- a/moose-core/python/setup.py
+++ b/moose-core/python/setup.py
@@ -62,6 +62,8 @@ setup(
, 'moose.SBML'
, 'moose.neuroml'
, 'moose.genesis'
+ , 'moose.chemUtil'
+ , 'moose.merge'
],
package_dir = {
'moose' : 'moose'
diff --git a/moose-core/synapse/RollingMatrix.cpp b/moose-core/synapse/RollingMatrix.cpp
index 5fd84a9bfe83c4ab42a239731bc2d96a827759de..f4f5193feab6dc2746ec65c94bdaf8e324b55b15 100644
--- a/moose-core/synapse/RollingMatrix.cpp
+++ b/moose-core/synapse/RollingMatrix.cpp
@@ -11,6 +11,7 @@
using namespace std;
#include "RollingMatrix.h"
#include <cassert>
+// #include <iostream>
RollingMatrix::RollingMatrix()
: nrows_(0), ncolumns_(0), currentStartRow_(0)
@@ -67,9 +68,21 @@ void RollingMatrix::sumIntoRow( const vector< double >& input, unsigned int row
double RollingMatrix::dotProduct( const vector< double >& input,
unsigned int row, unsigned int startColumn ) const
{
+ /// startColumn is the middle of the kernel.
unsigned int index = (row + currentStartRow_) % nrows_;
const SparseVector& sv = rows_[index];
+ unsigned int i2 = input.size()/2;
+ unsigned int istart = (startColumn <= i2) ? 0 : startColumn - i2;
+ unsigned int colstart = (startColumn <= i2) ? 0 : startColumn - i2;
+ unsigned int iend = (sv.size()-startColumn > i2 ) ? input.size() :
+ i2 - startColumn + sv.size();
+ // if ( iend >= istart ) cout << startColumn << i2 << istart << iend << colstart << "\n";
+ double ret = 0;
+ for (unsigned int i = istart, j = 0; i < iend; ++i, ++j )
+ ret += sv[j + colstart] * input[i];
+
+ /*
double ret = 0;
if ( input.size() + startColumn <= sv.size() ) {
for (unsigned int i = 0; i < input.size(); ++i )
@@ -79,6 +92,7 @@ double RollingMatrix::dotProduct( const vector< double >& input,
for (unsigned int i = 0; i < end; ++i )
ret += sv[i + startColumn] * input[i];
}
+ */
return ret;
}
diff --git a/moose-core/synapse/SeqSynHandler.cpp b/moose-core/synapse/SeqSynHandler.cpp
index 833b340221fb080dcb957d0889e9468c376c82a0..e6ca94f8e8bccadd54293ee9267d136bce0cfcb9 100644
--- a/moose-core/synapse/SeqSynHandler.cpp
+++ b/moose-core/synapse/SeqSynHandler.cpp
@@ -42,17 +42,22 @@ const Cinfo* SeqSynHandler::initCinfo()
"kernel with the history[time][synapse#] matrix."
"\nThe local response can affect the synapse in three ways: "
"1. It can sum the entire response vector, scale by the "
- "*responseScale* term, and send to the synapse as a steady "
+ "*sequenceScale* term, and send to the synapse as a steady "
"activation. Consider this a cell-wide immediate response to "
"a sequence that it likes.\n"
"2. It do an instantaneous scaling of the weight of each "
"individual synapse by the corresponding entry in the response "
- "vector. It uses the *weightScale* term to do this. Consider "
+ "vector. It uses the *plasticityScale* term to do this. "
+ "Consider "
"this a short-term plasticity effect on specific synapses. \n"
"3. It can do long-term plasticity of each individual synapse "
"using the matched local entries in the response vector and "
"individual synapse history as inputs to the learning rule. "
"This is not yet implemented.\n"
+ "In addition to all these, the SeqSynHandler can act just like "
+ "a regular synapse, where it responds to individual synaptic "
+ "input according to the weight of the synapse. The size of "
+ "this component of the output is scaled by *baseScale*\n"
};
static FieldElementFinfo< SynHandlerBase, Synapse > synFinfo(
@@ -88,22 +93,28 @@ const Cinfo* SeqSynHandler::initCinfo()
&SeqSynHandler::setHistoryTime,
&SeqSynHandler::getHistoryTime
);
- static ValueFinfo< SeqSynHandler, double > responseScale(
- "responseScale",
+ static ValueFinfo< SeqSynHandler, double > baseScale(
+ "baseScale",
+ "Basal scaling factor for regular synaptic activation.",
+ &SeqSynHandler::setBaseScale,
+ &SeqSynHandler::getBaseScale
+ );
+ static ValueFinfo< SeqSynHandler, double > sequenceScale(
+ "sequenceScale",
"Scaling factor for sustained activation of synapse by seq",
- &SeqSynHandler::setResponseScale,
- &SeqSynHandler::getResponseScale
+ &SeqSynHandler::setSequenceScale,
+ &SeqSynHandler::getSequenceScale
);
static ReadOnlyValueFinfo< SeqSynHandler, double > seqActivation(
"seqActivation",
"Reports summed activation of synaptic channel by sequence",
&SeqSynHandler::getSeqActivation
);
- static ValueFinfo< SeqSynHandler, double > weightScale(
- "weightScale",
- "Scaling factor for weight of each synapse by response vector",
- &SeqSynHandler::setWeightScale,
- &SeqSynHandler::getWeightScale
+ static ValueFinfo< SeqSynHandler, double > plasticityScale(
+ "plasticityScale",
+ "Scaling factor for doing plasticity by scaling each synapse by response vector",
+ &SeqSynHandler::setPlasticityScale,
+ &SeqSynHandler::getPlasticityScale
);
static ReadOnlyValueFinfo< SeqSynHandler, vector< double > >
weightScaleVec(
@@ -128,9 +139,10 @@ const Cinfo* SeqSynHandler::initCinfo()
&kernelWidth, // Field
&seqDt, // Field
&historyTime, // Field
- &responseScale, // Field
+ &sequenceScale, // Field
+ &baseScale, // Field
&seqActivation, // Field
- &weightScale, // Field
+ &plasticityScale, // Field
&weightScaleVec, // Field
&kernel, // ReadOnlyField
&history // ReadOnlyField
@@ -160,12 +172,12 @@ SeqSynHandler::SeqSynHandler()
kernelWidth_( 5 ),
historyTime_( 2.0 ),
seqDt_ ( 1.0 ),
- responseScale_( 1.0 ),
- weightScale_( 0.0 ),
+ sequenceScale_( 1.0 ),
+ baseScale_( 0.0 ),
+ plasticityScale_( 0.0 ),
seqActivation_( 0.0 )
{
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
- history_.resize( numHistory, 0 );
+ history_.resize( numHistory(), 0 );
}
SeqSynHandler::~SeqSynHandler()
@@ -193,8 +205,7 @@ void SeqSynHandler::vSetNumSynapses( const unsigned int v )
for ( unsigned int i = prevSize; i < v; ++i )
synapses_[i].setHandler( this );
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
- history_.resize( numHistory, v );
+ history_.resize( numHistory(), v );
latestSpikes_.resize( v, 0.0 );
weightScaleVec_.resize( v, 0.0 );
updateKernel();
@@ -229,9 +240,9 @@ void SeqSynHandler::updateKernel()
p.DefineConst(_T("e"), (mu::value_type)M_E);
p.SetExpr( kernelEquation_ );
kernel_.clear();
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
- kernel_.resize( numHistory );
- for ( int i = 0; i < numHistory; ++i ) {
+ int nh = numHistory();
+ kernel_.resize( nh );
+ for ( int i = 0; i < nh; ++i ) {
kernel_[i].resize( kernelWidth_ );
t = i * seqDt_;
for ( unsigned int j = 0; j < kernelWidth_; ++j ) {
@@ -268,8 +279,7 @@ void SeqSynHandler::setSeqDt( double v )
{
seqDt_ = v;
updateKernel();
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
- history_.resize( numHistory, vGetNumSynapses() );
+ history_.resize( numHistory(), vGetNumSynapses() );
}
double SeqSynHandler::getSeqDt() const
@@ -280,8 +290,7 @@ double SeqSynHandler::getSeqDt() const
void SeqSynHandler::setHistoryTime( double v )
{
historyTime_ = v;
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
- history_.resize( numHistory, vGetNumSynapses() );
+ history_.resize( numHistory(), vGetNumSynapses() );
updateKernel();
}
@@ -290,14 +299,24 @@ double SeqSynHandler::getHistoryTime() const
return historyTime_;
}
-void SeqSynHandler::setResponseScale( double v )
+void SeqSynHandler::setBaseScale( double v )
+{
+ baseScale_ = v;
+}
+
+double SeqSynHandler::getBaseScale() const
+{
+ return baseScale_;
+}
+
+void SeqSynHandler::setSequenceScale( double v )
{
- responseScale_ = v;
+ sequenceScale_ = v;
}
-double SeqSynHandler::getResponseScale() const
+double SeqSynHandler::getSequenceScale() const
{
- return responseScale_;
+ return sequenceScale_;
}
double SeqSynHandler::getSeqActivation() const
@@ -305,26 +324,26 @@ double SeqSynHandler::getSeqActivation() const
return seqActivation_;
}
-double SeqSynHandler::getWeightScale() const
+double SeqSynHandler::getPlasticityScale() const
{
- return weightScale_;
+ return plasticityScale_;
}
-vector< double >SeqSynHandler::getWeightScaleVec() const
+void SeqSynHandler::setPlasticityScale( double v )
{
- return weightScaleVec_;
+ plasticityScale_ = v;
}
-void SeqSynHandler::setWeightScale( double v )
+vector< double >SeqSynHandler::getWeightScaleVec() const
{
- weightScale_ = v;
+ return weightScaleVec_;
}
vector< double > SeqSynHandler::getKernel() const
{
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
+ int nh = numHistory();
vector< double > ret;
- for ( int i = 0; i < numHistory; ++i ) {
+ for ( int i = 0; i < nh; ++i ) {
ret.insert( ret.end(), kernel_[i].begin(), kernel_[i].end() );
}
return ret;
@@ -332,11 +351,11 @@ vector< double > SeqSynHandler::getKernel() const
vector< double > SeqSynHandler::getHistory() const
{
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
+ int nh = numHistory();
int numX = vGetNumSynapses();
- vector< double > ret( numX * numHistory, 0.0 );
+ vector< double > ret( numX * nh, 0.0 );
vector< double >::iterator k = ret.begin();
- for ( int i = 0; i < numHistory; ++i ) {
+ for ( int i = 0; i < nh; ++i ) {
for ( int j = 0; j < numX; ++j )
*k++ = history_.get( i, j );
}
@@ -375,10 +394,10 @@ void SeqSynHandler::dropSynapse( unsigned int msgLookup )
void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
{
// Here we look at the correlations and do something with them.
- int numHistory = static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
+ int nh = numHistory();
// Check if we need to do correlations at all.
- if ( numHistory > 0 && kernel_.size() > 0 ) {
+ if ( nh > 0 && kernel_.size() > 0 ) {
// Check if timestep rolls over a seqDt boundary
if ( static_cast< int >( p->currTime / seqDt_ ) >
static_cast< int >( (p->currTime - p->dt) / seqDt_ ) ) {
@@ -388,22 +407,22 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
// Build up the sum of correlations over time
vector< double > correlVec( vGetNumSynapses(), 0.0 );
- for ( int i = 0; i < numHistory; ++i )
+ for ( int i = 0; i < nh; ++i )
history_.correl( correlVec, kernel_[i], i );
- if ( responseScale_ > 0.0 ) { // Sum all responses, send to chan
+ if ( sequenceScale_ > 0.0 ) { // Sum all responses, send to chan
seqActivation_ = 0.0;
for ( vector< double >::iterator y = correlVec.begin();
y != correlVec.end(); ++y )
seqActivation_ += *y;
// We'll use the seqActivation_ to send a special msg.
- seqActivation_ *= responseScale_;
+ seqActivation_ *= sequenceScale_;
}
- if ( weightScale_ > 0.0 ) { // Short term changes in individual wts
+ if ( plasticityScale_ > 0.0 ) { // Short term changes in individual wts
weightScaleVec_ = correlVec;
for ( vector< double >::iterator y=weightScaleVec_.begin();
y != weightScaleVec_.end(); ++y )
- *y *= weightScale_;
+ *y *= plasticityScale_;
}
}
}
@@ -412,15 +431,15 @@ void SeqSynHandler::vProcess( const Eref& e, ProcPtr p )
// We can't leave it to the base class vProcess, because we need
// to scale the weights individually in some cases.
double activation = seqActivation_; // Start with seq activation
- if ( weightScale_ > 0.0 ) {
+ if ( plasticityScale_ > 0.0 ) {
while( !events_.empty() && events_.top().time <= p->currTime ) {
- activation += events_.top().weight *
+ activation += events_.top().weight * baseScale_ *
weightScaleVec_[ events_.top().synIndex ] / p->dt;
events_.pop();
}
} else {
while( !events_.empty() && events_.top().time <= p->currTime ) {
- activation += events_.top().weight / p->dt;
+ activation += baseScale_ * events_.top().weight / p->dt;
events_.pop();
}
}
@@ -435,3 +454,7 @@ void SeqSynHandler::vReinit( const Eref& e, ProcPtr p )
events_.pop();
}
+int SeqSynHandler::numHistory() const
+{
+ return static_cast< int >( 1.0 + floor( historyTime_ * (1.0 - 1e-6 ) / seqDt_ ) );
+}
diff --git a/moose-core/synapse/SeqSynHandler.h b/moose-core/synapse/SeqSynHandler.h
index 55fd01181bcf868b610c0192456fc04f5c05b71d..73578b0420bbec7c877726bf7aa0831725b5b455 100644
--- a/moose-core/synapse/SeqSynHandler.h
+++ b/moose-core/synapse/SeqSynHandler.h
@@ -55,15 +55,20 @@ class SeqSynHandler: public SynHandlerBase
double getSeqDt() const;
void setHistoryTime( double v );
double getHistoryTime() const;
- void setResponseScale( double v );
- double getResponseScale() const;
+ void setBaseScale( double v );
+ double getBaseScale() const;
+ void setSequenceScale( double v );
+ double getSequenceScale() const;
double getSeqActivation() const; // summed activation of syn chan
- void setWeightScale( double v );
- double getWeightScale() const;
+ void setPlasticityScale( double v );
+ double getPlasticityScale() const;
vector< double > getWeightScaleVec() const;
vector< double > getKernel() const;
vector< double > getHistory() const;
+ ////////////////////////////////////////////////////////////////
+ // Utility func
+ int numHistory() const;
////////////////////////////////////////////////////////////////
static const Cinfo* initCinfo();
private:
@@ -84,13 +89,19 @@ class SeqSynHandler: public SynHandlerBase
string kernelEquation_;
unsigned int kernelWidth_; // Width in terms of number of synapses
- // Time to store history. KernelDt defines num of rows
+ /// Time to store history. KernelDt defines num of rows
double historyTime_;
double seqDt_; // Time step for successive entries in kernel
- // Scaling factor for sustained activation of synapse from response
- double responseScale_;
- // Scaling factor for weight changes in each synapse from response
- double weightScale_;
+ /// Scaling factor for baseline synaptic responses.
+ double baseScale_;
+ /// Scaling factor for sequence recognition responses.
+ double sequenceScale_;
+
+ /**
+ * Scaling factor for short-term plastic weight changes in each
+ * synapse arising from sequential input.
+ */
+ double plasticityScale_;
///////////////////////////////////////////
// Some readonly fields
diff --git a/moose-core/synapse/testSynapse.cpp b/moose-core/synapse/testSynapse.cpp
index 2edd4165c1c5866effe687f2e08e5bc943e792fe..db96a2aacd8ef5b116f35e2cefd83c193113d749 100644
--- a/moose-core/synapse/testSynapse.cpp
+++ b/moose-core/synapse/testSynapse.cpp
@@ -20,6 +20,48 @@
#include "../shell/Shell.h"
#include "../randnum/randnum.h"
+double doCorrel( RollingMatrix& rm, vector< vector< double >> & kernel )
+{
+ int nr = kernel.size();
+ vector< double > correlVec( nr, 0.0 );
+ for ( int i = 0; i < nr; ++i )
+ rm.correl( correlVec, kernel[i], i );
+ double seqActivation = 0.0;
+ for ( int i = 0; i < nr; ++i )
+ seqActivation += correlVec[i];
+ return seqActivation;
+}
+
+void testRollingMatrix2()
+{
+ int nr = 5;
+ RollingMatrix rm;
+ rm.resize( nr, nr );
+ vector< vector< double > > kernel( nr );
+ for ( int i = 0; i < nr; ++i ) {
+ kernel[i].resize( nr, 0.0 );
+ rm.zeroOutRow( i );
+ for ( int j = 0; j < nr; ++j ) {
+ kernel[i][j] = 16 - (i-j)*(i-j); // symmetric, forward
+ rm.sumIntoEntry( (i==j), i, j );
+ }
+ }
+ double ret1 = doCorrel( rm, kernel );
+
+ for ( int i = 0; i < nr; ++i ) {
+ kernel[i].clear();
+ kernel[i].resize( nr, 0.0 );
+ rm.zeroOutRow( i );
+ for ( int j = 0; j < nr; ++j ) {
+ int k = nr-i-1;
+ kernel[i][j] = 16 - (k-j)*(k-j); // symmetric, backwards
+ rm.sumIntoEntry( (k==j), i, j );
+ }
+ }
+ double ret2 = doCorrel( rm, kernel );
+ assert( doubleEq( ret1, ret2 ) );
+}
+
void testRollingMatrix()
{
int nr = 5;
@@ -150,6 +192,7 @@ void testSynapse()
{
#ifdef DO_UNIT_TESTS
testRollingMatrix();
+ testRollingMatrix2();
testSeqSynapse();
#endif // DO_UNIT_TESTS
}
diff --git a/moose-core/tests/python/test_difshells.py b/moose-core/tests/python/test_difshells.py
index 8e6da7c8ccb8850ba98f529d8c690b8f21cad281..367b953a8a55d3bcb8d6563b4fa57f7688a42450 100644
--- a/moose-core/tests/python/test_difshells.py
+++ b/moose-core/tests/python/test_difshells.py
@@ -5,17 +5,19 @@ import chan_proto
from params import *
-def linoid(x,param):
- #plt.figure()
- #plt.plot(x,(param[2]+np.exp((V+param[3])/param[4])))
- den = (param[2]+np.exp((V+param[3])/param[4]))
- nom = (param[0]+param[1]*V)
- return nom/den
-
-def add_difshell(comp,i,shell_thickness,shell_radius):
- new_name = comp.path.split('[')[0]+'/'+comp.name+'/'+difshell_name+str(i)
- dif = moose.DifShell(new_name)
+def linoid(x, param):
+ # plt.figure()
+ # plt.plot(x,(param[2]+np.exp((V+param[3])/param[4])))
+ den = (param[2] + np.exp((V + param[3]) / param[4]))
+ nom = (param[0] + param[1] * V)
+ return nom / den
+
+
+def add_difshell(comp, i, shell_thickness, shell_radius):
+ new_name = comp.path.split('[')[0] + '/' + \
+ comp.name + '/' + difshell_name + str(i)
+ dif = moose.DifShell(new_name)
#dif.C = Ca_initial
dif.Ceq = Ca_basal
@@ -24,71 +26,81 @@ def add_difshell(comp,i,shell_thickness,shell_radius):
dif.leak = 0
dif.shapeMode = 0
dif.length = comp.length
- dif.diameter = 2*shell_radius
+ dif.diameter = 2 * shell_radius
dif.thickness = shell_thickness
-
+
return dif
-def add_difbuffer_to_dishell(comp,i,j,shell_thickness,shell_radius):
- new_name = comp.path.split('[')[0]+'/'+comp.name+'/'+difshell_name+str(i)+'_'+difbuff_name+str(j)
+
+def add_difbuffer_to_dishell(comp, i, j, shell_thickness, shell_radius):
+ new_name = comp.path.split('[')[0] + '/' + comp.name + \
+ '/' + difshell_name + str(i) + '_' + difbuff_name + str(j)
buf = moose.DifBuffer(new_name)
buf.bTot = btotal
buf.shapeMode = 0
buf.length = comp.length
- buf.diameter = 2*shell_radius
+ buf.diameter = 2 * shell_radius
buf.thickness = shell_thickness
buf.kf = kf
buf.kb = kb
buf.D = d
return buf
+
def add_difshells_and_buffers(comp):
-
+
if difshell_no < 1:
return [], []
-
- difshell= []
- shell_thickness = dend.diameter/difshell_no/2.
+
+ difshell = []
+ shell_thickness = dend.diameter / difshell_no / 2.
difbuffer = []
for i in range(difshell_no):
- shell_radius = comp.diameter/2-i*shell_thickness
- dif = add_difshell(comp,i,shell_thickness,shell_radius)
+ shell_radius = comp.diameter / 2 - i * shell_thickness
+ dif = add_difshell(comp, i, shell_thickness, shell_radius)
difshell.append(dif)
-
- if i > 0 :
- moose.connect(difshell[i-1],"outerDifSourceOut",difshell[i],"fluxFromOut")
- moose.connect(difshell[i],"innerDifSourceOut",difshell[i-1],"fluxFromIn")
-
+
+ if i > 0:
+ moose.connect(
+ difshell[i - 1], "outerDifSourceOut", difshell[i], "fluxFromOut")
+ moose.connect(difshell[i], "innerDifSourceOut",
+ difshell[i - 1], "fluxFromIn")
+
if difbuff_no > 0:
difbuffer.append([])
for j in range(difbuff_no):
- buf = add_difbuffer_to_dishell(comp,i,j,shell_thickness,shell_radius)
+ buf = add_difbuffer_to_dishell(
+ comp, i, j, shell_thickness, shell_radius)
difbuffer[i].append(buf)
- moose.connect(difshell[i],"concentrationOut",buf,"concentration")
- moose.connect(buf,"reactionOut",difshell[i],"reaction")
- if i>0:
- moose.connect(difbuffer[i-1][j],"outerDifSourceOut",difbuffer[i][j],"fluxFromOut")
- moose.connect(difbuffer[i][j],"innerDifSourceOut",difbuffer[i-1][j],"fluxFromIn")
-
- return difshell,difbuffer
-
-def addOneChan(chanpath,gbar,comp):
-
- SA = np.pi * comp.length*comp.diameter
- proto = moose.element('/library/'+chanpath)
+ moose.connect(
+ difshell[i], "concentrationOut", buf, "concentration")
+ moose.connect(buf, "reactionOut", difshell[i], "reaction")
+ if i > 0:
+ moose.connect(
+ difbuffer[i - 1][j], "outerDifSourceOut", difbuffer[i][j], "fluxFromOut")
+ moose.connect(difbuffer[i][j], "innerDifSourceOut", difbuffer[
+ i - 1][j], "fluxFromIn")
+
+ return difshell, difbuffer
+
+
+def addOneChan(chanpath, gbar, comp):
+
+ SA = np.pi * comp.length * comp.diameter
+ proto = moose.element('/library/' + chanpath)
chan = moose.copy(proto, comp, chanpath)
chan.Gbar = gbar * SA
- #If we are using GHK AND it is a calcium channel, connect it to GHK
- moose.connect(comp,'VmOut',chan,'Vm')
- moose.connect(chan,"channelOut",comp,"handleChannel")
+ # If we are using GHK AND it is a calcium channel, connect it to GHK
+ moose.connect(comp, 'VmOut', chan, 'Vm')
+ moose.connect(chan, "channelOut", comp, "handleChannel")
return chan
if __name__ == '__main__':
-
+
for tick in range(0, 7):
moose.setClock(tick, dt)
- moose.setClock(8,0.005) #set output clock
+ moose.setClock(8, 0.005) # set output clock
model = moose.Neutral('/model')
dend = moose.Compartment('/model/dend')
pulse = moose.PulseGen('/model/pulse')
@@ -103,20 +115,18 @@ if __name__ == '__main__':
dend.initVm = Vm_0
dend.diameter = dend_diameter
dend.length = dend_length
-
pulse.delay[0] = 8.
pulse.width[0] = 500e-3
pulse.level[0] = inject
pulse.delay[1] = 1e9
m = moose.connect(pulse, 'output', dend, 'injectMsg')
-
+
moose.connect(vmtab, 'requestOut', dend, 'getVm')
chan_proto.chanlib()
- chan = addOneChan('CaT',gbar,dend)
+ chan = addOneChan('CaT', gbar, dend)
-
moose.connect(gktab, 'requestOut', chan, 'getGk')
moose.connect(iktab, 'requestOut', chan, 'getIk')
diftab = []
@@ -126,57 +136,54 @@ if __name__ == '__main__':
pump = moose.MMPump('/model/dend/pump')
pump.Vmax = kcat
pump.Kd = km
- moose.connect(pump,"PumpOut",difs[0],"mmPump")
+ moose.connect(pump, "PumpOut", difs[0], "mmPump")
if difs:
- moose.connect(chan,"IkOut",difs[0],"influx")
- moose.connect(difs[0],'concentrationOut',chan,'concen')
-
- for i,dif in enumerate(difs):
- res_dif = moose.Table('/data/'+difshell_name+str(i))
+ moose.connect(chan, "IkOut", difs[0], "influx")
+ moose.connect(difs[0], 'concentrationOut', chan, 'concen')
+
+ for i, dif in enumerate(difs):
+ res_dif = moose.Table('/data/' + difshell_name + str(i))
diftab.append(res_dif)
- moose.connect(diftab[i],'requestOut',dif,'getC')
+ moose.connect(diftab[i], 'requestOut', dif, 'getC')
if (difbuff_no):
buftab.append([])
- for j,buf in enumerate(difb[i]):
- res_buf = moose.Table('/data/'+difshell_name+str(i)+'_'+difbuff_name+str(j))
+ for j, buf in enumerate(difb[i]):
+ res_buf = moose.Table(
+ '/data/' + difshell_name + str(i) + '_' + difbuff_name + str(j))
buftab[i].append(res_buf)
- moose.connect(buftab[i][j],'requestOut',buf,'getBBound')
-
-
+ moose.connect(buftab[i][j], 'requestOut', buf, 'getBBound')
+
moose.reinit()
if not gbar:
- for i,dif in enumerate(difs):
+ for i, dif in enumerate(difs):
if i == 0:
dif.C = Ca_initial
else:
dif.C = 0
- for j,dbuf in enumerate(difb[i]):
+ for j, dbuf in enumerate(difb[i]):
dbuf.bFree = dbuf.bTot
-
+
moose.start(t_stop)
-
+
t = np.linspace(0, t_stop, len(vmtab.vector))
- fname = 'moose_results_difshell_no_'+str(difshell_no)+'_difbuffer_no_'+str(difbuff_no)+'_pump_'+str(pumps)+'_gbar_'+str(gbar)+'.txt'
- print fname
+ fname = 'moose_results_difshell_no_' + str(difshell_no) + '_difbuffer_no_' + str(
+ difbuff_no) + '_pump_' + str(pumps) + '_gbar_' + str(gbar) + '.txt'
+ print( fname )
header = 'time Vm Ik Gk'
- number = 4+difshell_no*(difbuff_no+1)
- res = np.zeros((len(t),number))
- res[:,0] = t
- res[:,1] = vmtab.vector
- res[:,2] = iktab.vector
- res[:,3] = gktab.vector
-
+ number = 4 + difshell_no * (difbuff_no + 1)
+ res = np.zeros((len(t), number))
+ res[:, 0] = t
+ res[:, 1] = vmtab.vector
+ res[:, 2] = iktab.vector
+ res[:, 3] = gktab.vector
+
for i in range(difshell_no):
- header += ' difshell_'+str(i)
- res[:,4+i*(difbuff_no+1)] = diftab[i].vector
+ header += ' difshell_' + str(i)
+ res[:, 4 + i * (difbuff_no + 1)] = diftab[i].vector
if (difbuff_no):
- for j,buf in enumerate(buftab[i]):
- res[:,4+i*(difbuff_no+1)+j+1] = buf.vector
- header += ' difshell_'+str(i)+'_difbuff_'+str(j)
- np.savetxt(fname,res,header=header, comments='')
-
-
-
-
+ for j, buf in enumerate(buftab[i]):
+ res[:, 4 + i * (difbuff_no + 1) + j + 1] = buf.vector
+ header += ' difshell_' + str(i) + '_difbuff_' + str(j)
+ np.savetxt(fname, res, header=header, comments='')
# plt.show()