diff --git a/moose-core/external/muparser/include/muParserBase.h b/moose-core/external/muparser/include/muParserBase.h index 004fe20e9f1221ba1554fb0b7046f84c2e9f9a3f..774330cce8fc914920b754846fb141d574e7858a 100644 --- a/moose-core/external/muparser/include/muParserBase.h +++ b/moose-core/external/muparser/include/muParserBase.h @@ -288,7 +288,7 @@ private: mutable stringbuf_type m_vStringBuf; ///< String buffer, used for storing string function arguments stringbuf_type m_vStringVarBuf; - std::auto_ptr<token_reader_type> m_pTokenReader; ///< Managed pointer to the token reader object. + std::unique_ptr<token_reader_type> m_pTokenReader; ///< Managed pointer to the token reader object. funmap_type m_FunDef; ///< Map of function names and pointers. funmap_type m_PostOprtDef; ///< Postfix operator callbacks diff --git a/moose-core/external/muparser/include/muParserToken.h b/moose-core/external/muparser/include/muParserToken.h index b3a879a7fc1c45379fcb44380781ff5093126d94..25a40163002feb831b21efa7ac280736bf6c05a8 100644 --- a/moose-core/external/muparser/include/muParserToken.h +++ b/moose-core/external/muparser/include/muParserToken.h @@ -69,7 +69,7 @@ namespace mu TString m_strTok; ///< Token string TString m_strVal; ///< Value for string variables value_type m_fVal; ///< the value - std::auto_ptr<ParserCallback> m_pCallback; + std::unique_ptr<ParserCallback> m_pCallback; public: diff --git a/moose-core/external/muparser/src/muParserTest.cpp b/moose-core/external/muparser/src/muParserTest.cpp index c03624e89de990f43b0f65332ff2bcd7f65d0330..3b1dec6c8a2e13a826c374aaf8eac0b686bf46c5 100644 --- a/moose-core/external/muparser/src/muParserTest.cpp +++ b/moose-core/external/muparser/src/muParserTest.cpp @@ -1258,7 +1258,7 @@ namespace mu try { - std::auto_ptr<Parser> p1; + std::unique_ptr<Parser> p1; Parser p2, p3; // three parser objects // they will be used for testing copy and assignment operators // p1 is a pointer since i'm going to delete it in order to test if diff --git a/moose-core/python/rdesigneur/__init__.py b/moose-core/python/rdesigneur/__init__.py index 5f50f095337268b66d6ab731acb68013df524db7..188e877d28b4f897865afd7f15a6332c191559a1 100644 --- a/moose-core/python/rdesigneur/__init__.py +++ b/moose-core/python/rdesigneur/__init__.py @@ -2,3 +2,4 @@ from __future__ import print_function, absolute_import from rdesigneur.rdesigneur import * +from rdesigneur.rdesigneur import rdesigneur diff --git a/moose-core/python/rdesigneur/rdesigneur.py b/moose-core/python/rdesigneur/rdesigneur.py index 4220f3763ecc8e5ba84062b1d82cedcbe775148d..5d1de41ef95d72dc860cb02fb503c3eae6507048 100644 --- a/moose-core/python/rdesigneur/rdesigneur.py +++ b/moose-core/python/rdesigneur/rdesigneur.py @@ -16,39 +16,27 @@ ## latter in the former, including mapping entities like calcium and ## channel conductances, between them. ########################################################################## -from __future__ import print_function +from __future__ import print_function, absolute_import + import imp import os import moose import numpy as np import pylab import math -import rmoogli -#import rdesigneurProtos -from rdesigneurProtos import * + +import rdesigneur.rmoogli +from rdesigneur.rdesigneurProtos import * + from moose.neuroml.NeuroML import NeuroML from moose.neuroml.ChannelML import ChannelML +# In python3, cElementTree is deprecated. We do not plan to support python <2.7 +# in future, so other imports have been removed. try: from lxml import etree except ImportError: - try: - # Python 2.5 - import xml.etree.cElementTree as etree - except ImportError: - try: - # Python 2.5 - import xml.etree.ElementTree as etree - except ImportError: - try: - # normal cElementTree install - import cElementTree as etree - except ImportError: - try: - # normal ElementTree install - import elementtree.ElementTree as etree - except ImportError: - print("Failed to import ElementTree from any known place") + import xml.etree.ElementTree as etree import csv @@ -85,8 +73,8 @@ class rdesigneur: diffusionLength= 2e-6, meshLambda = -1.0, #This is a backward compatibility hack temperature = 32, - chemDt= 0.001, # Much finer than MOOSE, for multiscale - diffDt= 0.001, # 10x finer than MOOSE, for multiscale + chemDt= 0.1, # Much finer than MOOSE, for multiscale + diffDt= 0.01, # 10x finer than MOOSE, for multiscale elecDt= 50e-6, # Same default as from MOOSE chemPlotDt = 1.0, # Same default as from MOOSE elecPlotDt = 0.1e-3, # Same default as from MOOSE diff --git a/moose-core/tests/python/abstrModelEqns9.py b/moose-core/tests/python/abstrModelEqns9.py new file mode 100644 index 0000000000000000000000000000000000000000..ae57be3f67667d603f4bda398264e1d6990da81d --- /dev/null +++ b/moose-core/tests/python/abstrModelEqns9.py @@ -0,0 +1,240 @@ +# -*- coding: utf-8 -*- +import re +import moose + +# Equations here are: +# Adot = 1 -6A + 5A^2 - A^3, or spread out as: +# Adot = k0a + k1a.A + k2a.A.A + k3a.A.A.A + k4a.Ca.A/(1+A+10*B) - k5a.A.B +# Bdot = k1b.A - k2b.B + + +def parseExpr( expr, params, hasCa ): + if hasCa: + expr = expr.replace( 'Ca', 'x0' ) + expr = expr.replace( 'A', 'x1' ) + expr = expr.replace( 'B', 'x2' ) + else: + expr = expr.replace( 'Ca', 'x0' ) # happens in the negFF model + expr = expr.replace( 'A', 'x0' ) # This is the usual case. + expr = expr.replace( 'B', 'x1' ) + + parts = re.split( 'k', expr ) + ret = parts[0] + for i in parts[1:]: + ret += str( params[ 'k' + i[:2] ] ) + ret += i[2:] + + if hasCa: + return 'x3*( ' + ret + ')' + else: + return 'x2*( ' + ret + ')' + +def makeChemProto( name, Aexpr, Bexpr, params ): + sw = params['stimWidth'] + diffLength = params['diffusionLength'] + dca = params['diffConstA'] * diffLength * diffLength + dcb = params['diffConstB'] * diffLength * diffLength + + # Objects + chem = moose.Neutral( '/library/' + name ) + compt = moose.CubeMesh( '/library/' + name + '/' + name ) + A = moose.Pool( compt.path + '/A' ) + B = moose.Pool( compt.path + '/B' ) + Z = moose.BufPool( compt.path + '/Z' ) + Ca = moose.BufPool( compt.path + '/Ca' ) + phase = moose.BufPool( compt.path + '/phase' ) + vel = moose.BufPool( compt.path + '/vel' ) + ampl = moose.BufPool( compt.path + '/ampl' ) + Adot = moose.Function( A.path + '/Adot' ) + Bdot = moose.Function( B.path + '/Bdot' ) + CaStim = moose.Function( Ca.path + '/CaStim' ) + A.diffConst = dca + B.diffConst = dcb + + # Equations + + Adot.expr = parseExpr( Aexpr, params, True ) + Bdot.expr = parseExpr( Bexpr, params, False ) + CaStim.expr = 'x2 * exp( -((x0 - t)^2)/(2* ' + str(sw*sw) + ') )' + + #print Adot.expr + #print Bdot.expr + #print CaStim.expr + + # Connections + Adot.x.num = 4 + moose.connect( Ca, 'nOut', Adot.x[0], 'input' ) + moose.connect( A, 'nOut', Adot.x[1], 'input' ) + moose.connect( B, 'nOut', Adot.x[2], 'input' ) + moose.connect( Z, 'nOut', Adot.x[3], 'input' ) + moose.connect( Adot, 'valueOut', A, 'increment' ) + + Bdot.x.num = 3 + if name[:5] == 'negFF': + moose.connect( Ca, 'nOut', Bdot.x[0], 'input' ) + print('Doing special msg') + else: + moose.connect( A, 'nOut', Bdot.x[0], 'input' ) + moose.connect( B, 'nOut', Bdot.x[1], 'input' ) + moose.connect( Z, 'nOut', Bdot.x[2], 'input' ) + moose.connect( Bdot, 'valueOut', B, 'increment' ) + + CaStim.x.num = 3 + moose.connect( phase, 'nOut', CaStim.x[0], 'input' ) + moose.connect( vel, 'nOut', CaStim.x[1], 'input' ) + moose.connect( ampl, 'nOut', CaStim.x[2], 'input' ) + moose.connect( CaStim, 'valueOut', Ca, 'setN' ) + + return compt + + +def makeBis( args ): + params = { + 'k0a':0.1, # Constant + 'k1a':-5.0, # Coeff for A + 'k2a':5.0, # Coeff for A^2 + 'k3a':-1.0, # Coeff for A^3 + 'k4a':10.0, # turnon of A by A and Ca + 'k5a':-5.0, # Turnoff of A by B + 'k1b':0.01, # turnon of B by A + 'k2b':-0.01, # Decay rate of B + 'diffusionLength':1.0e-6, # Diffusion characteristic length, used as voxel length too. + 'dendDiameter': 10e-6, # Diameter of section of dendrite in model + 'dendLength': 100e-6, # Length of section of dendrite in model + 'diffConstA':5.0, # Diffusion constant of A + 'diffConstB':2.0, # Diffusion constant of B + 'stimWidth' :1.0, # Stimulus width in seconds + 'stimAmplitude':1.0, # Stimulus amplitude, arb units. From FHN review + 'blankVoxelsAtEnd':10, # of voxels to leave blank at end of cylinder + 'preStimTime':10.0, # Time to run before turning on stimulus. + 'postStimTime':40.0, # Time to run after stimulus. ~3x decay time + 'settleTime':20.0, # Settling time of response, after stimulus. + # To include in analysis of total response over + # entire dendrite. + 'fnumber':1, # Number to append to fname + } + + for i in args: + params[i] = args[i] + + makeChemProto( 'bis', + 'k0a + k1a*A + k2a*A*A + k3a*A*A*A + k4a*Ca*A/(1+A+10*B) + k5a*A*B', + 'k1b*A*A + k2b*B', + params ) + return params + +def makeFHN( args ): + params = { + 'k_t':2.5, # Time-const. + 'k_a':0.7, # Coeff1 + 'k_b':0.8, # Coeff2 + 'kxa': 2.0, # Offset for V for FHN eqns. + 'kxb': 0.8, # Offset for W for FHN eqns. + 'diffusionLength':1.0e-6, # Diffusion characteristic length, used as voxel length too. + 'dendDiameter': 10e-6, # Diameter of section of dendrite in model + 'dendLength': 100e-6, # Length of section of dendrite in model + 'diffConstA':7.5, # Diffusion constant of A + 'diffConstB':5.0, # Diffusion constant of B + 'stimWidth' :1.0, # Stimulus width in seconds + 'stimAmplitude':0.4, # Stimulus amplitude, arb units. From FHN review + 'blankVoxelsAtEnd':10, # of voxels to leave blank at end of cylinder + 'preStimTime':10.0, # Time to run before turning on stimulus. + 'postStimTime':40.0, # Time to run after stimulus. ~3x decay time + 'settleTime':20.0, # Settling time of response, after stimulus. + # To include in analysis of total response over + # entire dendrite. + 'fnumber': 1, # Number to append to fname + } + + for i in args: + params[i] = args[i] + + makeChemProto( 'fhn', + '5.0*(A - kxa - (A-kxa)^3/3 - (B-kxb) + Ca)', + '(A-kxa + k_a - k_b*(B-kxb))/k_t', + params ) + # We do this to get the system to settle at the start. + B = moose.element( '/library/fhn/fhn/B' ) + B.nInit = 0.2 + return params + + +def makeNegFB( args ): + params = { + 'k1a':-0.1, # Coeff for decay of A, slow. + 'k2a':-0.2, # Coeff for turnoff of A by B, medium. + 'k3a':10.0, # turnon of A by Ca, fast. + 'k1b':0.2, # turnon of B by A, medium + 'k2b':-0.1, # Decay rate of B, slow + 'diffusionLength':1.0e-6, # Diffusion characteristic length, used as voxel length too. + 'dendDiameter': 10e-6, # Diameter of section of dendrite in model + 'dendLength': 100e-6, # Length of section of dendrite in model + 'diffConstA':0.5, # Diffusion constant of A + 'diffConstB':1.0, # Diffusion constant of B + 'stimWidth' :1.0, # Stimulus width in seconds + 'stimAmplitude':1.0, # Stimulus amplitude, arb units. From FHN review + 'blankVoxelsAtEnd':10, #of voxels to leave blank at end of cylinder + 'preStimTime':10.0, # Time to run before turning on stimulus. + 'postStimTime':40.0, # Time to run after stimulus. ~3x decay time + 'settleTime':20.0, # Settling time of response, after stimulus. + # To include in analysis of total response over + # entire dendrite. + 'fnumber':1, # Number to append to fname + } + + for i in args: + params[i] = args[i] + + makeChemProto( 'negFB', + 'k1a*A + k2a*A*B + k3a*Ca', + 'k1b*A + k2b*B', + params ) + return params + +# Was negFF2 in earlier versions of abstrModelEqns +def makeNegFF( args ): + params = { + 'k1a':-0.1, # Coeff for decay of A, slow. + 'k2a':-0.01, # Coeff for turnoff of A by B, medium. + 'k3a':10.0, # turnon of A by Ca, fast. + 'k4a':40.0, # amount of B inhibition of turnon of A by Ca. + 'k1b':2.0, # turnon of B by Ca, medium + 'k2b':-0.05, # Decay rate of B, slow + 'diffusionLength':1.0e-6, # Diffusion characteristic length, used as voxel length too. + 'dendDiameter': 10e-6, # Diameter of section of dendrite in model + 'dendLength': 100e-6, # Length of section of dendrite in model + 'diffConstA':0.02, # Diffusion constant of A + 'diffConstB':0.4, # Diffusion constant of B + 'stimWidth' :1.0, # Stimulus width in seconds + 'stimAmplitude':10.0, # Stimulus amplitude, arb units. From FHN review + 'blankVoxelsAtEnd':10, # of voxels to leave blank at end of cylinder + 'preStimTime':10.0, # Time to run before turning on stimulus. + 'postStimTime':40.0, # Time to run after stimulus. ~3x decay time + 'settleTime':20.0, # Settling time of response, after stimulus. + # To include in analysis of total response over + # entire dendrite. + 'fnumber':1, # Number to append to fname + } + + for i in args: + params[i] = args[i] + + makeChemProto( 'negFF', + 'k1a*A + k2a*A*B + k3a*Ca/(1+k4a*B*B)', + 'k1b*Ca + k2b*B', + params ) + + return params + +if __name__ == '__main__': + moose.Neutral( '/library' ) + print("Making Bistable model") + makeBis() + print("Making FHN model") + makeFHN() + print("Making Negative Feedback model") + makeNegFB() + print("Making Negative Feedforward models") + makeNegFF() + + diff --git a/moose-core/tests/python/test_muparser.py b/moose-core/tests/python/test_muparser.py new file mode 100644 index 0000000000000000000000000000000000000000..d889df4680ff8403ee369beac57212b8fab3cdd0 --- /dev/null +++ b/moose-core/tests/python/test_muparser.py @@ -0,0 +1,167 @@ +# -*- coding: utf-8 -*- +"""test_muparser.py: + +""" + +__author__ = "Dilawar Singh" +__copyright__ = "Copyright 2017-, Dilawar Singh" +__version__ = "1.0.0" +__maintainer__ = "Dilawar Singh" +__email__ = "dilawars@ncbs.res.in" +__status__ = "Development" + +import sys +import os +import numpy as np +import sys +import numpy as np +import moose +import abstrModelEqns9 as ame +import rdesigneur as rd + + +def singleCompt( name, params ): + mod = moose.copy( '/library/' + name + '/' + name, '/model' ) + A = moose.element( mod.path + '/A' ) + Z = moose.element( mod.path + '/Z' ) + Z.nInit = 1 + Ca = moose.element( mod.path + '/Ca' ) + CaStim = moose.element( Ca.path + '/CaStim' ) + runtime = params['preStimTime'] + params['postStimTime'] + steptime = 50 + + CaStim.expr += ' + x2 * (t > 100+' + str( runtime ) + ' ) * ( t < 100+' + str( runtime + steptime ) + ' )' + print(CaStim.expr) + tab = moose.Table2( '/model/' + name + '/Atab' ) + ampl = moose.element( mod.path + '/ampl' ) + phase = moose.element( mod.path + '/phase' ) + moose.connect( tab, 'requestOut', A, 'getN' ) + ampl.nInit = params['stimAmplitude'] * 1 + phase.nInit = params['preStimTime'] + + ksolve = moose.Ksolve( mod.path + '/ksolve' ) + stoich = moose.Stoich( mod.path + '/stoich' ) + stoich.compartment = mod + stoich.ksolve = ksolve + stoich.path = mod.path + '/##' + + moose.reinit() + runtime += 100 + steptime*2 + moose.start( runtime ) + t = np.arange( 0, runtime + 1e-9, tab.dt ) + return name, t, tab.vector + + +def plotPanelC(): + panelC = [] + panelCticks = [] + panelC.append( singleCompt( 'negFB', ame.makeNegFB( [] ) ) ) + panelC.append( singleCompt( 'negFF', ame.makeNegFF( [] ) ) ) + panelC.append( singleCompt( 'fhn', ame.makeFHN( [] ) ) ) + panelC.append( singleCompt( 'bis', ame.makeBis( [] ) ) ) + + panelCticks.append( np.arange( 0, 15.00001, 5 ) ) + panelCticks.append( np.arange( 0, 1.50001, 0.5 ) ) + panelCticks.append( np.arange( 0, 5.00002, 1 ) ) + panelCticks.append( np.arange( 0, 5.00002, 1 ) ) + moose.delete( '/model' ) + for i in zip( panelC, panelCticks, list(range( len( panelC ))) ): + plotPos = i[2] + 5 + doty = i[1][-1] * 0.95 + print('doty', doty ) + +def runPanelDEFG( name, dist, seqDt, numSpine, seq, stimAmpl ): + preStim = 10.0 + blanks = 20 + rdes = rd.rdesigneur( + useGssa = False, + turnOffElec = True, + chemPlotDt = 0.1, + diffusionLength = 1e-6, + cellProto = [['cell', 'soma']], + chemProto = [['dend', name]], + chemDistrib = [['dend', 'soma', 'install', '1' ]], + plotList = [['soma', '1', 'dend' + '/A', 'n', '# of A']], + ) + rdes.buildModel() + A = moose.vec( '/model/chem/dend/A' ) + Z = moose.vec( '/model/chem/dend/Z' ) + print(moose.element( '/model/chem/dend/A/Adot' ).expr) + print(moose.element( '/model/chem/dend/B/Bdot' ).expr) + print(moose.element( '/model/chem/dend/Ca/CaStim' ).expr) + phase = moose.vec( '/model/chem/dend/phase' ) + ampl = moose.vec( '/model/chem/dend/ampl' ) + vel = moose.vec( '/model/chem/dend/vel' ) + vel.nInit = 1e-6 * seqDt + ampl.nInit = stimAmpl + stride = int( dist ) / numSpine + phase.nInit = 10000 + Z.nInit = 0 + for j in range( numSpine ): + k = int( blanks + j * stride ) + Z[k].nInit = 1 + phase[k].nInit = preStim + seq[j] * seqDt + + moose.reinit() + runtime = 50 + snapshot = preStim + seqDt * (numSpine - 0.8) + print(snapshot) + moose.start( snapshot ) + avec = moose.vec( '/model/chem/dend/A' ).n + moose.start( runtime - snapshot ) + tvec = [] + for i in range( 5 ): + tab = moose.element( '/model/graphs/plot0[' + str( blanks + i * stride ) + ']' ) + dt = tab.dt + tvec.append( tab.vector ) + moose.delete( '/model' ) + return dt, tvec, avec + +def makePassiveSoma( name, length, diameter ): + elecid = moose.Neuron( '/library/' + name ) + dend = moose.Compartment( elecid.path + '/soma' ) + dend.diameter = diameter + dend.length = length + dend.x = length + return elecid + +def plotOnePanel( tLabel, dt, tplot, numSyn, plotRange, tick ): + t = np.arange( 0, len( tplot[0] ), 1.0 ) * dt + for i in range( 5 ): + print( tplot[i] ) + + +def plotPanelDEFG( seq, row ): + makePassiveSoma( 'cell', 100e-6, 10e-6 ) + start = (row -1) * 4 + tLabel = chr( ord( 'B' ) + row - 1 ) + xLabel = chr( ord( 'D' ) + row - 1 ) + xplot = [] + + + dt, tplot, avec = runPanelDEFG( 'negFB', 5.0, 2.0, 5, seq, 1.0 ) + xplot.append( avec ) + t = np.arange( 0, len( tplot[0] ), 1.0 ) * dt + + dt, tplot, avec = runPanelDEFG( 'negFF', 10.0, 1.0, 5, seq, 1.0 ) + xplot.append( avec ) + t = np.arange( 0, len( tplot[0] ), 1.0 ) * dt + dt, tplot, avec = runPanelDEFG( 'fhn', 5.0, 1.5, 5, seq, 0.4 ) + xplot.append( avec ) + t = np.arange( 0, len( tplot[0] ), 1.0 ) * dt + + dt, tplot, avec = runPanelDEFG( 'bis', 15.0, 2.0, 5, seq, 1.0 ) + xplot.append( avec ) + t = np.arange( 0, len( tplot[0] ), 1.0 ) * dt + ref = np.array( [0.04840247926323106, 0.060446860947786119 + , 0.047612129439079991, 0.081329604404641223, 0.050365470686926379] ) + res = np.array( [ np.mean( x ) for x in tplot ] ) + assert np.isclose( ref.all(), res.all() ), "Expected %s got %s" % (ref,res) + +if __name__ == '__main__': + moose.Neutral( '/library' ) + moose.Neutral( '/model' ) + plotPanelC() + plotPanelDEFG( [0,1,2,3,4], 3 ) + plotPanelDEFG( [4,1,0,3,2], 4 ) + print( 'All done' ) diff --git a/moose-core/tests/python/test_rdesigneur.py b/moose-core/tests/python/test_rdesigneur.py new file mode 100644 index 0000000000000000000000000000000000000000..2972a3c06f71ae4f1312fbdae1282e84fe136bae --- /dev/null +++ b/moose-core/tests/python/test_rdesigneur.py @@ -0,0 +1,40 @@ +"""test_rdesigneur.py: + +""" + +__author__ = "Dilawar Singh" +__copyright__ = "Copyright 2017-, Dilawar Singh" +__version__ = "1.0.0" +__maintainer__ = "Dilawar Singh" +__email__ = "dilawars@ncbs.res.in" +__status__ = "Development" + +import sys +import os +import numpy as np +import moose +import rdesigneur as rd + +def test2( ): + if moose.exists( '/model' ): + moose.delete( '/model' ) + + rdes = rd.rdesigneur( + stimList = [['soma', '1', '.', 'inject', '(t>0.1 && t<0.2) * 2e-8' ]], + plotList = [['soma', '1', '.', 'Vm', 'Soma membrane potential']] + ) + rdes.buildModel() + moose.reinit() + moose.start( 0.3 ) + +def test1( ): + rdes = rd.rdesigneur() + rdes.buildModel() + moose.showfields( rdes.soma ) + +def main( ): + test1() + test2( ) + +if __name__ == '__main__': + main()