From 628138cb9426ec4c1dcf98d3ef8a5108d0ba783b Mon Sep 17 00:00:00 2001
From: Dilawar Singh <dilawars@ncbs.res.in>
Date: Wed, 15 Feb 2017 12:09:49 +0530
Subject: [PATCH] Squashed 'moose-core/' changes from f3d0477..434c478
434c478 qMerge branch 'master' of https://github.com/BhallaLab/moose-core
425d8c8 writeKKit: checked if matplotlib is install before importing, moose.py:removed all unwanted check's on genesis_support and sbml_support
0791eb5 Fixes to SeqSynHandler
51a3a30 Moved Pool::increment, decrement and nIn to PoolBase for uniformity with Solvers. Dummy functions except in Pool and BufPool
1c9a9df group info is added to plot
03dbc2c While connecting table for plotting, a check is made to see if its pool
fdc4924 Cleaner way to replace the string
4f92820 libsbml spell corrected
22a276d unwanted print statement cleanup
2c48b1c reading group information from the SBML if exist into moose
e68a509 spell correction
8deaed8 switched from new to old to old to new while providing the files for mergering
ef9edd1 clenup while loading model and solvers are deleted
c2383e2 python3 compatibility with print statement
4ac79b1 Function cleanup while mergering and path cleanup
6fd5ed6 chemUtil and merge folders also needed to copied while using "make install" while using cmake
caaafc9 Fixed test script for python3. Removed VERSION file. Pull request BhallaLab/moose#206 .
a9a8344 Merge commit '7c71305865364216fd1a00b4a2498da44f0cbb0b'
6d907ee Script testing DifShells
0dd5cd6 Merge commit 'a5122cfb6dd82f1a831f64d1fb1a29a7d1392334'
4fc5760 Fixed docstring on example.
git-subtree-dir: moose-core
git-subtree-split: 434c47885637c72d820f92914660a679b34c0301
---
VERSION | 2 +-
kinetics/Pool.cpp | 58 +------
kinetics/Pool.h | 6 +-
kinetics/PoolBase.cpp | 45 ++++++
kinetics/PoolBase.h | 6 +
python/moose/SBML/readSBML.py | 56 +++++--
python/moose/SBML/writeSBML.py | 13 +-
python/moose/genesis/writeKkit.py | 184 +++++++++++++----------
python/moose/merge/merge.py | 241 +++++++++++++++---------------
python/moose/moose.py | 47 +-----
python/setup.py | 2 +
synapse/RollingMatrix.cpp | 14 ++
synapse/SeqSynHandler.cpp | 125 +++++++++-------
synapse/SeqSynHandler.h | 29 ++--
synapse/testSynapse.cpp | 43 ++++++
tests/python/chan_proto.py | 140 +++++++++++++++++
tests/python/param_chan.py | 80 ++++++++++
tests/python/params.py | 26 ++++
tests/python/test_all.sh | 3 +-
tests/python/test_difshells.py | 190 +++++++++++++++++++++++
20 files changed, 922 insertions(+), 388 deletions(-)
create mode 100644 tests/python/chan_proto.py
create mode 100644 tests/python/param_chan.py
create mode 100644 tests/python/params.py
create mode 100644 tests/python/test_difshells.py
diff --git a/VERSION b/VERSION
index ffbc5f42..5b9fc17b 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-3.1.1-79-gc1dec94
\ No newline at end of file
+3.1.1-100-g51a3a30
\ No newline at end of file
diff --git a/kinetics/Pool.cpp b/kinetics/Pool.cpp
index 8b8b6561..7ebc22f9 100644
--- a/kinetics/Pool.cpp
+++ b/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/kinetics/Pool.h b/kinetics/Pool.h
index dec4dfa4..9d0cca01 100644
--- a/kinetics/Pool.h
+++ b/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/kinetics/PoolBase.cpp b/kinetics/PoolBase.cpp
index 5dd4ea24..59c66ba5 100644
--- a/kinetics/PoolBase.cpp
+++ b/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/kinetics/PoolBase.h b/kinetics/PoolBase.h
index ace7ae86..8173c445 100644
--- a/kinetics/PoolBase.h
+++ b/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/python/moose/SBML/readSBML.py b/python/moose/SBML/readSBML.py
index 801d3681..6a627054 100644
--- a/python/moose/SBML/readSBML.py
+++ b/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/python/moose/SBML/writeSBML.py b/python/moose/SBML/writeSBML.py
index bcfaf2eb..90400306 100644
--- a/python/moose/SBML/writeSBML.py
+++ b/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/python/moose/genesis/writeKkit.py b/python/moose/genesis/writeKkit.py
index ba716031..2bd33e21 100644
--- a/python/moose/genesis/writeKkit.py
+++ b/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/python/moose/merge/merge.py b/python/moose/merge/merge.py
index f76ca494..02ecc805 100644
--- a/python/moose/merge/merge.py
+++ b/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/python/moose/moose.py b/python/moose/moose.py
index 904d6c32..1309c68f 100644
--- a/python/moose/moose.py
+++ b/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/python/setup.py b/python/setup.py
index 47b5f3aa..a6c3852c 100644
--- a/python/setup.py
+++ b/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/synapse/RollingMatrix.cpp b/synapse/RollingMatrix.cpp
index 5fd84a9b..f4f5193f 100644
--- a/synapse/RollingMatrix.cpp
+++ b/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/synapse/SeqSynHandler.cpp b/synapse/SeqSynHandler.cpp
index 833b3402..e6ca94f8 100644
--- a/synapse/SeqSynHandler.cpp
+++ b/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/synapse/SeqSynHandler.h b/synapse/SeqSynHandler.h
index 55fd0118..73578b04 100644
--- a/synapse/SeqSynHandler.h
+++ b/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/synapse/testSynapse.cpp b/synapse/testSynapse.cpp
index 2edd4165..db96a2aa 100644
--- a/synapse/testSynapse.cpp
+++ b/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/tests/python/chan_proto.py b/tests/python/chan_proto.py
new file mode 100644
index 00000000..ae37024f
--- /dev/null
+++ b/tests/python/chan_proto.py
@@ -0,0 +1,140 @@
+"""\
+Create general Channel Proto, pass in name, x and y power, and params
+
+Also, create the library of channels
+Might need a few other chan_proto types, such as
+ inf-tau channels
+ Ca dep channels
+chan_proto quires alpha and beta params for both activation and inactivation
+If no inactivation, just send in empty Yparam array.
+"""
+
+from __future__ import print_function, division
+import moose
+import numpy as np
+from util import NamedList
+
+SSTauChannelParams = NamedList('SSTauChannelParams', '''
+ Arate
+ A_B
+ A_C
+ Avhalf
+ Avslope
+ taumin
+ tauVdep
+ tauPow
+ tauVhalf
+ tauVslope''')
+
+AlphaBetaChannelParams = NamedList('AlphaBetaChannelParams', '''
+ A_rate
+ A_B
+ A_C
+ Avhalf
+ A_vslope
+ B_rate
+ B_B
+ B_C
+ Bvhalf
+ B_vslope''')
+
+ZChannelParams = NamedList('ZChannelParams', 'Kd power tau')
+
+
+ChannelSettings = NamedList('ChannelSettings', 'Xpow Ypow Zpow Erev name')
+
+def interpolate_values_in_table(tabA,V_0,l=40):
+ import param_chan
+ '''This function interpolates values in the table
+ around tabA[V_0]. '''
+ V = np.linspace(param_chan.VMIN, param_chan.VMAX, len(tabA))
+ idx = abs(V-V_0).argmin()
+ A_min = tabA[idx-l]
+ V_min = V[idx-l]
+ A_max = tabA[idx+l]
+ V_max = V[idx+l]
+ a = (A_max-A_min)/(V_max-V_min)
+ b = A_max - a*V_max
+ tabA[idx-l:idx+l] = V[idx-l:idx+l]*a+b
+ return tabA
+
+def fix_singularities(Params,Gate):
+ import param_chan
+
+ if Params.A_C < 0:
+
+ V_0 = Params.A_vslope*np.log(-Params.A_C)-Params.Avhalf
+
+ if V_0 > param_chan.VMIN and V_0 < param_chan.VMAX:
+ #change values in tableA and tableB, because tableB contains sum of alpha and beta
+ tabA = interpolate_values_in_table(Gate.tableA,V_0)
+ tabB = interpolate_values_in_table(Gate.tableB,V_0)
+ Gate.tableA = tabA
+ Gate.tableB = tabB
+
+ if Params.B_C < 0:
+
+ V_0 = Params.B_vslope*np.log(-Params.B_C)-Params.Bvhalf
+
+ if V_0 > param_chan.VMIN and V_0 < param_chan.VMAX:
+ #change values in tableB
+ tabB = interpolate_values_in_table(Gate.tableB,V_0)
+ Gate.tableB = tabB
+
+ return Gate
+
+#may need a CaV channel if X gate uses alpha,beta and Ygate uses inf tau
+#Or, have Y form an option - if in tau, do something like NaF
+def chan_proto(chanpath, params):
+ import param_chan
+
+ chan = moose.HHChannel(chanpath)
+ chan.Xpower = params.channel.Xpow
+ if params.channel.Xpow > 0:
+ xGate = moose.HHGate(chan.path + '/gateX')
+ xGate.setupAlpha(params.X + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ xGate = fix_singularities(params.X, xGate)
+
+ chan.Ypower = params.channel.Ypow
+ if params.channel.Ypow > 0:
+ yGate = moose.HHGate(chan.path + '/gateY')
+ yGate.setupAlpha(params.Y + [param_chan.VDIVS, param_chan.VMIN, param_chan.VMAX])
+ yGate = fix_singularities(params.Y, yGate)
+ if params.channel.Zpow > 0:
+ chan.Zpower = params.channel.Zpow
+ zgate = moose.HHGate(chan.path + '/gateZ')
+ ca_array = np.linspace(param_chan.CAMIN, param_chan.CAMAX, param_chan.CADIVS)
+ zgate.min = param_chan.CAMIN
+ zgate.max = param_chan.CAMAX
+ caterm = (ca_array/params.Z.Kd) ** params.Z.power
+ inf_z = caterm / (1 + caterm)
+ tau_z = params.Z.tau * np.ones(len(ca_array))
+ zgate.tableA = inf_z / tau_z
+ zgate.tableB = 1 / tau_z
+ chan.useConcentration = True
+ chan.Ek = params.channel.Erev
+ return chan
+
+
+TypicalOneDalpha = NamedList('TypicalOneDalpha',
+ '''channel X Y Z=[] calciumPermeable=False calciumPermeable2=False''')
+
+_FUNCTIONS = {
+ TypicalOneDalpha: chan_proto,
+
+}
+
+#*params... passes the set of values not as a list but as individuals
+def make_channel(chanpath, params):
+ func = _FUNCTIONS[params.__class__]
+ return func(chanpath, params)
+
+def chanlib():
+ import param_chan
+ if not moose.exists('/library'):
+ moose.Neutral('/library')
+ #Adding all the channels to the library. *list removes list elements from the list,
+ #so they can be used as function arguments
+ chan = [make_channel('/library/'+key, value) for key, value in param_chan.ChanDict.items()]
+
+ return chan
diff --git a/tests/python/param_chan.py b/tests/python/param_chan.py
new file mode 100644
index 00000000..0db3f697
--- /dev/null
+++ b/tests/python/param_chan.py
@@ -0,0 +1,80 @@
+# -*- coding: utf-8 -*-
+
+from util import NamedDict
+from chan_proto import (
+ SSTauChannelParams,
+ AlphaBetaChannelParams,
+ ZChannelParams,
+ ChannelSettings,
+ TypicalOneDalpha,)
+
+#chanDictSP.py
+#contains all gating parameters and reversal potentials
+# Gate equations have the form:
+#
+# y(x) = (rate + B * x) / (C + exp((x + vhalf) / vslope))
+#
+# OR
+# y(x) = tau_min + tau_vdep / (1 + exp((x + vhalf) / vslope))
+#
+# where x is membrane voltage and y is the rate constant
+#KDr params used by Sriram, RE paper1, Krp params used by RE paper 2
+#Parameters for Ca channels may need to be shifted - see Dorman model
+krev=-87e-3
+narev=50e-3
+carev=140e-3 #assumes CaExt=2 mM and CaIn=50e-3
+ZpowCDI=0
+
+VMIN = -100e-3
+VMAX = 50e-3
+VDIVS = 3000 #0.5 mV steps
+CAMIN=0.01e-3 #10 nM
+CAMAX=40e-3 #40 uM, might want to go up to 100 uM with spines
+CADIVS=3000 #10 nM steps
+
+#mtau: Ogata fig 5, no qfactor accounted in mtau, 1.2 will improve spike shape
+#activation minf fits Ogata 1990 figure 3C (which is cubed root)
+#inactivation hinf fits Ogata 1990 figure 6B
+#htau fits the main -50 through -10 slope of Ogata figure 9 (log tau), but a qfact of 2 is already taken into account.
+
+
+CaTparam = ChannelSettings(Xpow=3, Ypow=1, Zpow=1, Erev=carev, name='CaT')
+qfactCaT = 1
+
+#Original inactivation ws too slow compared to activation, made closder the alpha1G
+CaT_X_params = AlphaBetaChannelParams(A_rate = 5342.4*qfactCaT,
+ A_B = 2100.*qfactCaT,
+ A_C = 11.9,
+ Avhalf = 1e-3,
+ A_vslope = -12e-3,
+ B_rate = 289.7*qfactCaT,
+ B_B = 0.*qfactCaT,
+ B_C = 1.,
+ Bvhalf = 0.0969,
+ B_vslope = 0.0141)
+#CaT_Y_params = CaT_X_params
+CaT_Y_params = AlphaBetaChannelParams(A_rate = 0*qfactCaT,
+ A_B = -74.*qfactCaT,
+ A_C = 1.,
+ Avhalf = 9e-2,
+ A_vslope = 5e-3,
+ B_rate = 15*qfactCaT,
+ B_B = -1.5*qfactCaT,
+ B_C = 1.5,
+ Bvhalf = 0.05,
+ B_vslope = -15e-3)
+
+#These CDI params can be used with every channel, make ZpowCDI=2
+#If ZpowCDI=0 the CDI will not be used, power=-4 is to transform
+#(Ca/Kd)^pow/(1+(Ca/Kd)^pow) to 1/(1+(ca/Kd)^-pow)
+CDI_Z_params = ZChannelParams(Kd = 0.12e-3,
+ power = -4,
+ tau = 142e-3)
+
+#Dictionary of "standard" channels, to create channels using a loop
+#NaF doesn't fit since it uses different prototype form
+#will need separate dictionary for BK
+
+ChanDict = NamedDict(
+ 'ChannelParams',
+ CaT = TypicalOneDalpha(CaTparam, CaT_X_params, CaT_Y_params, CDI_Z_params, calciumPermeable=True))
diff --git a/tests/python/params.py b/tests/python/params.py
new file mode 100644
index 00000000..40696573
--- /dev/null
+++ b/tests/python/params.py
@@ -0,0 +1,26 @@
+t_stop = 10
+dend_diameter = 2.2627398e-6
+dend_length = 1.131369936e-6
+Cm = 4.021231698e-12
+Rm = 1865100032
+Em = -0.07100000232
+Vm_0 = -0.0705
+dt = 50e-6
+spines_no = 0
+difshell_no = 2
+difshell_name = "Ca_shell"
+Ca_basal = 50e-6
+Ca_initial = Ca_basal*200
+dca = 200.0e-12
+difbuff_no = 1
+difbuff_name = "Buff"
+btotal = 80.0e-3
+kf = 0.028e6
+kb = 19.6
+d = 66e-12
+inject = 0.1e-9
+gbar = 1
+#MMpump
+km = 0.3e-3
+kcat = 85e-22
+pumps = 1
diff --git a/tests/python/test_all.sh b/tests/python/test_all.sh
index 9ec99b4e..338a20cb 100755
--- a/tests/python/test_all.sh
+++ b/tests/python/test_all.sh
@@ -6,7 +6,8 @@ test_scripts="./test_function.py
test_mumbl.py
test_pymoose.py
test_synchan.py
- test_vec.py"
+ test_vec.py
+ test_difshells.py"
for testFile in $test_scripts; do
echo "Executing $testFile"
diff --git a/tests/python/test_difshells.py b/tests/python/test_difshells.py
new file mode 100644
index 00000000..669f46a8
--- /dev/null
+++ b/tests/python/test_difshells.py
@@ -0,0 +1,190 @@
+import moose
+import numpy as np
+import matplotlib.pyplot as plt
+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)
+
+ #dif.C = Ca_initial
+ dif.Ceq = Ca_basal
+ dif.D = dca
+ dif.valence = 2
+ dif.leak = 0
+ dif.shapeMode = 0
+ dif.length = comp.length
+ 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)
+ buf = moose.DifBuffer(new_name)
+ buf.bTot = btotal
+ buf.shapeMode = 0
+ buf.length = comp.length
+ 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.
+ difbuffer = []
+ for i in range(difshell_no):
+ 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 difbuff_no > 0:
+ difbuffer.append([])
+ for j in range(difbuff_no):
+ 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)
+ 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")
+ return chan
+
+
+if __name__ == '__main__':
+
+ for tick in range(0, 7):
+ moose.setClock(tick, dt)
+ moose.setClock(8, 0.005) # set output clock
+ model = moose.Neutral('/model')
+ dend = moose.Compartment('/model/dend')
+ pulse = moose.PulseGen('/model/pulse')
+ data = moose.Neutral('/data')
+ vmtab = moose.Table('/data/dend_Vm')
+ gktab = moose.Table('/data/CaT_Gk')
+ iktab = moose.Table('/data/CaT_Ik')
+
+ dend.Cm = Cm
+ dend.Rm = Rm
+ dend.Em = Em
+ 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)
+
+ moose.connect(gktab, 'requestOut', chan, 'getGk')
+ moose.connect(iktab, 'requestOut', chan, 'getIk')
+ diftab = []
+ buftab = []
+ difs, difb = add_difshells_and_buffers(dend)
+ if pumps:
+ pump = moose.MMPump('/model/dend/pump')
+ pump.Vmax = kcat
+ pump.Kd = km
+ 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))
+ diftab.append(res_dif)
+ 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))
+ buftab[i].append(res_buf)
+ moose.connect(buftab[i][j], 'requestOut', buf, 'getBBound')
+
+ moose.reinit()
+ if not gbar:
+ for i, dif in enumerate(difs):
+ if i == 0:
+ dif.C = Ca_initial
+ else:
+ dif.C = 0
+ 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 )
+ 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
+
+ for i in range(difshell_no):
+ 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='')
+
+ # plt.show()
--
GitLab