diff --git a/biophysics/Neuron.cpp b/biophysics/Neuron.cpp index 4d5993f167..fc72bf01b1 100644 --- a/biophysics/Neuron.cpp +++ b/biophysics/Neuron.cpp @@ -19,13 +19,13 @@ #include "Neuron.h" #include "../basecode/global.h" -#include "../external/muparser/include/muParser.h" +#include "../builtins/MooseParser.h" -class nuParser: public mu::Parser +class nuParser: public MooseParser { public: nuParser( const string& expr ): - mu::Parser(), + MooseParser(), p(0.0), // geometrical path distance wound through dendrite g(0.0), // geometrical path distance direct from soma. L(0.0), // electrical distance arg @@ -109,40 +109,40 @@ const Cinfo* Neuron::initCinfo() // ValueFinfos ///////////////////////////////////////////////////////////////////// static ValueFinfo< Neuron, double > RM( "RM", - "Membrane resistivity, in ohm.m^2. Default value is 1.0.", - &Neuron::setRM, - &Neuron::getRM - ); + "Membrane resistivity, in ohm.m^2. Default value is 1.0.", + &Neuron::setRM, + &Neuron::getRM + ); static ValueFinfo< Neuron, double > RA( "RA", - "Axial resistivity of cytoplasm, in ohm.m. Default value is 1.0.", - &Neuron::setRA, - &Neuron::getRA - ); + "Axial resistivity of cytoplasm, in ohm.m. Default value is 1.0.", + &Neuron::setRA, + &Neuron::getRA + ); static ValueFinfo< Neuron, double > CM( "CM", - "Membrane Capacitance, in F/m^2. Default value is 0.01", - &Neuron::setCM, - &Neuron::getCM - ); + "Membrane Capacitance, in F/m^2. Default value is 0.01", + &Neuron::setCM, + &Neuron::getCM + ); static ValueFinfo< Neuron, double > Em( "Em", - "Resting membrane potential of compartments, in Volts. " - "Default value is -0.065.", - &Neuron::setEm, - &Neuron::getEm - ); + "Resting membrane potential of compartments, in Volts. " + "Default value is -0.065.", + &Neuron::setEm, + &Neuron::getEm + ); static ValueFinfo< Neuron, double > theta( "theta", "Angle to rotate cell geometry, around long axis of neuron. " "Think Longitude. Units are radians. " "Default value is zero, which means no rotation. ", &Neuron::setTheta, &Neuron::getTheta - ); + ); static ValueFinfo< Neuron, double > phi( "phi", "Angle to rotate cell geometry, around elevation of neuron. " "Think Latitude. Units are radians. " "Default value is zero, which means no rotation. ", &Neuron::setPhi, &Neuron::getPhi - ); + ); static ValueFinfo< Neuron, string > sourceFile( "sourceFile", "Name of source file from which to load a model. " @@ -151,317 +151,317 @@ const Cinfo* Neuron::initCinfo() "definitions should have been loaded into /library. ", &Neuron::setSourceFile, &Neuron::getSourceFile - ); + ); static ValueFinfo< Neuron, double > compartmentLengthInLambdas( - "compartmentLengthInLambdas", - "Units: meters (SI). \n" - "Electrotonic length to use for the largest compartment in the " - "model. Used to define subdivision of branches into compartments. " - "For example, if we set *compartmentLengthInLambdas* to 0.1, " - "and *lambda* (electrotonic length) is 250 microns, then it " - "sets the compartment length to 25 microns. Thus a dendritic " - "branch of 500 microns is subdivided into 20 commpartments. " - "If the branch is shorter than *compartmentLengthInLambdas*, " - "then it is not subdivided. " - "If *compartmentLengthInLambdas* is set to 0 then the original " - "compartmental structure of the model is preserved. " - " Note that this routine does NOT merge branches, even if " - "*compartmentLengthInLambdas* is bigger than the branch. " - "While all this subdivision is being done, the Neuron class " - "preserves as detailed a geometry as it can, so it can rebuild " - "the more detailed version if needed. " - "Default value of *compartmentLengthInLambdas* is 0. ", - &Neuron::setCompartmentLengthInLambdas, - &Neuron::getCompartmentLengthInLambdas - ); + "compartmentLengthInLambdas", + "Units: meters (SI). \n" + "Electrotonic length to use for the largest compartment in the " + "model. Used to define subdivision of branches into compartments. " + "For example, if we set *compartmentLengthInLambdas* to 0.1, " + "and *lambda* (electrotonic length) is 250 microns, then it " + "sets the compartment length to 25 microns. Thus a dendritic " + "branch of 500 microns is subdivided into 20 commpartments. " + "If the branch is shorter than *compartmentLengthInLambdas*, " + "then it is not subdivided. " + "If *compartmentLengthInLambdas* is set to 0 then the original " + "compartmental structure of the model is preserved. " + " Note that this routine does NOT merge branches, even if " + "*compartmentLengthInLambdas* is bigger than the branch. " + "While all this subdivision is being done, the Neuron class " + "preserves as detailed a geometry as it can, so it can rebuild " + "the more detailed version if needed. " + "Default value of *compartmentLengthInLambdas* is 0. ", + &Neuron::setCompartmentLengthInLambdas, + &Neuron::getCompartmentLengthInLambdas + ); static ElementValueFinfo< Neuron, vector< string > > - channelDistribution( - "channelDistribution", - "Specification for distribution of channels, CaConcens and " - "any other model components that are defined as prototypes and " - "have to be placed on the electrical compartments.\n" - "Arguments: proto path field expr [field expr]...\n" - " Each entry is terminated with an empty string. " - "The prototype is any object created in */library*, " - "If a channel matching the prototype name already exists, then " - "all subsequent operations are applied to the extant channel and " - "a new one is not created. " - "The paired arguments are as follows: \n" - "The *field* argument specifies the name of the parameter " - "that is to be assigned by the expression.\n" - "The *expression* argument is a mathematical expression in " - "the muparser framework, which permits most operations including " - "trig and transcendental ones. Of course it also handles simple " - "numerical values like 1.0, 1e-10 and so on. " - "Available arguments for muParser are:\n" - " p, g, L, len, dia, maxP, maxG, maxL \n" - " p: path distance from soma, measured along dendrite, in metres.\n" - " g: geometrical distance from soma, in metres.\n" - " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" - " len: length of compartment, in metres.\n" - " dia: for diameter of compartment, in metres.\n" - " maxP: Maximum value of *p* for this neuron. \n" - " maxG: Maximum value of *g* for this neuron. \n" - " maxL: Maximum value of *L* for this neuron.\n" - "The expression for the first field must evaluate to > 0 " - "for the channel to be installed. For example, for " - "channels, if Field == Gbar, and func( r, L, len, dia) < 0, \n" - "then the channel is not installed. This feature is typically " - "used with the sign() or Heaviside H() function to limit range: " - "for example: H(1 - L) will only put channels closer than " - "one length constant from the soma, and zero elsewhere. \n" - "Available fields are: \n" - "Channels: Gbar (install), Ek \n" - "CaConcen: shellDia (install), shellFrac (install), tau, min\n" - "Unless otherwise noted, all fields are scaled appropriately by " - "the dimensions of their compartment. Thus the channel " - "maximal conductance Gbar is automatically scaled by the area " - "of the compartment, and the user does not need to insert this " - "scaling into the calculations.\n" - "All parameters are expressed in SI units. Conductance, for " - "example, is Siemens/sq metre. " - "\n\n" - "Some example function forms might be for a channel Gbar: \n" - " p < 10e-6 ? 400 : 0.0 \n" - " equivalently, \n" - " H(10e-6 - p) * 400 \n" - " equivalently, \n" - " ( sign(10e-6 - p) + 1) * 200 \n" - "Each of these forms instruct the function to " - "set channel Gbar to 400 S/m^2 only within 10 microns path " - "distance of soma\n" - "\n" - " L < 1.0 ? 100 * exp( -L ) : 0.0 \n" - " ->Set channel Gbar to an exponentially falling function of " - "electrotonic distance from soma, provided L is under " - "1.0 lambdas. \n", + channelDistribution( + "channelDistribution", + "Specification for distribution of channels, CaConcens and " + "any other model components that are defined as prototypes and " + "have to be placed on the electrical compartments.\n" + "Arguments: proto path field expr [field expr]...\n" + " Each entry is terminated with an empty string. " + "The prototype is any object created in */library*, " + "If a channel matching the prototype name already exists, then " + "all subsequent operations are applied to the extant channel and " + "a new one is not created. " + "The paired arguments are as follows: \n" + "The *field* argument specifies the name of the parameter " + "that is to be assigned by the expression.\n" + "The *expression* argument is a mathematical expression in " + "the muparser framework, which permits most operations including " + "trig and transcendental ones. Of course it also handles simple " + "numerical values like 1.0, 1e-10 and so on. " + "Available arguments for muParser are:\n" + " p, g, L, len, dia, maxP, maxG, maxL \n" + " p: path distance from soma, measured along dendrite, in metres.\n" + " g: geometrical distance from soma, in metres.\n" + " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" + " len: length of compartment, in metres.\n" + " dia: for diameter of compartment, in metres.\n" + " maxP: Maximum value of *p* for this neuron. \n" + " maxG: Maximum value of *g* for this neuron. \n" + " maxL: Maximum value of *L* for this neuron.\n" + "The expression for the first field must evaluate to > 0 " + "for the channel to be installed. For example, for " + "channels, if Field == Gbar, and func( r, L, len, dia) < 0, \n" + "then the channel is not installed. This feature is typically " + "used with the sign() or Heaviside H() function to limit range: " + "for example: H(1 - L) will only put channels closer than " + "one length constant from the soma, and zero elsewhere. \n" + "Available fields are: \n" + "Channels: Gbar (install), Ek \n" + "CaConcen: shellDia (install), shellFrac (install), tau, min\n" + "Unless otherwise noted, all fields are scaled appropriately by " + "the dimensions of their compartment. Thus the channel " + "maximal conductance Gbar is automatically scaled by the area " + "of the compartment, and the user does not need to insert this " + "scaling into the calculations.\n" + "All parameters are expressed in SI units. Conductance, for " + "example, is Siemens/sq metre. " + "\n\n" + "Some example function forms might be for a channel Gbar: \n" + " p < 10e-6 ? 400 : 0.0 \n" + " equivalently, \n" + " H(10e-6 - p) * 400 \n" + " equivalently, \n" + " ( sign(10e-6 - p) + 1) * 200 \n" + "Each of these forms instruct the function to " + "set channel Gbar to 400 S/m^2 only within 10 microns path " + "distance of soma\n" + "\n" + " L < 1.0 ? 100 * exp( -L ) : 0.0 \n" + " ->Set channel Gbar to an exponentially falling function of " + "electrotonic distance from soma, provided L is under " + "1.0 lambdas. \n", &Neuron::setChannelDistribution, &Neuron::getChannelDistribution - ); + ); static ElementValueFinfo< Neuron, vector< string > > - passiveDistribution( - "passiveDistribution", - "Specification for distribution of passive properties of cell.\n" - "Arguments: . path field expr [field expr]...\n" - "Note that the arguments list starts with a period. " - " Each entry is terminated with an empty string. " - "The paired arguments are as follows: \n" - "The *field* argument specifies the name of the parameter " - "that is to be assigned by the expression.\n" - "The *expression* argument is a mathematical expression in " - "the muparser framework, which permits most operations including " - "trig and transcendental ones. Of course it also handles simple " - "numerical values like 1.0, 1e-10 and so on. " - "Available arguments for muParser are:\n" - " p, g, L, len, dia, maxP, maxG, maxL \n" - " p: path distance from soma, measured along dendrite, in metres.\n" - " g: geometrical distance from soma, in metres.\n" - " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" - " len: length of compartment, in metres.\n" - " dia: for diameter of compartment, in metres.\n" - " maxP: Maximum value of *p* for this neuron. \n" - " maxG: Maximum value of *g* for this neuron. \n" - " maxL: Maximum value of *L* for this neuron.\n" - "Available fields are: \n" - "RM, RA, CM, Rm, Ra, Cm, Em, initVm \n" - "The first three fields are scaled appropriately by " - "the dimensions of their compartment. Thus the membrane " - "resistivity RM (ohms.m^2) is automatically scaled by the area " - "of the compartment, and the user does not need to insert this " - "scaling into the calculations to compute Rm." - "Using the Rm field lets the user directly assign the " - "membrane resistance (in ohms), presumably using len and dia.\n" - "Similarly, RA (ohms.m) and CM (Farads/m^2) are specific units " - "and the actual values for each compartment are assigned by " - "scaling by length and diameter. Ra (ohms) and Cm (Farads) " - "require explicit evaluation of the expression. " - "All parameters are expressed in SI units. Conductance, for " - "example, is Siemens/sq metre.\n" - "Note that time these calculations do NOT currently include spines\n", + passiveDistribution( + "passiveDistribution", + "Specification for distribution of passive properties of cell.\n" + "Arguments: . path field expr [field expr]...\n" + "Note that the arguments list starts with a period. " + " Each entry is terminated with an empty string. " + "The paired arguments are as follows: \n" + "The *field* argument specifies the name of the parameter " + "that is to be assigned by the expression.\n" + "The *expression* argument is a mathematical expression in " + "the muparser framework, which permits most operations including " + "trig and transcendental ones. Of course it also handles simple " + "numerical values like 1.0, 1e-10 and so on. " + "Available arguments for muParser are:\n" + " p, g, L, len, dia, maxP, maxG, maxL \n" + " p: path distance from soma, measured along dendrite, in metres.\n" + " g: geometrical distance from soma, in metres.\n" + " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" + " len: length of compartment, in metres.\n" + " dia: for diameter of compartment, in metres.\n" + " maxP: Maximum value of *p* for this neuron. \n" + " maxG: Maximum value of *g* for this neuron. \n" + " maxL: Maximum value of *L* for this neuron.\n" + "Available fields are: \n" + "RM, RA, CM, Rm, Ra, Cm, Em, initVm \n" + "The first three fields are scaled appropriately by " + "the dimensions of their compartment. Thus the membrane " + "resistivity RM (ohms.m^2) is automatically scaled by the area " + "of the compartment, and the user does not need to insert this " + "scaling into the calculations to compute Rm." + "Using the Rm field lets the user directly assign the " + "membrane resistance (in ohms), presumably using len and dia.\n" + "Similarly, RA (ohms.m) and CM (Farads/m^2) are specific units " + "and the actual values for each compartment are assigned by " + "scaling by length and diameter. Ra (ohms) and Cm (Farads) " + "require explicit evaluation of the expression. " + "All parameters are expressed in SI units. Conductance, for " + "example, is Siemens/sq metre.\n" + "Note that time these calculations do NOT currently include spines\n", &Neuron::setPassiveDistribution, &Neuron::getPassiveDistribution - ); + ); static ElementValueFinfo< Neuron, vector< string > >spineDistribution( - "spineDistribution", - "Specification for distribution of spines on dendrite. \n" - "Arguments: proto path spacing expr [field expr]...\n" - " Each entry is terminated with an empty string. " - "The *prototype* is any spine object created in */library*, \n" - "The *path* is the wildcard path of compartments on which to " - "place the spine.\n" - "The *spacing* is the spacing of spines, in metres. \n" - "The *expression* argument is a mathematical expression in " - "the muparser framework, which permits most operations including " - "trig and transcendental ones. Of course it also handles simple " - "numerical values like 1.0, 1e-10 and so on. " - "The paired arguments are as follows: \n" - "The *field* argument specifies the name of the parameter " - "that is to be assigned by the expression.\n" - "The *expression* argument is a mathematical expression as above. " - "Available arguments for muParser are:\n" - " p, g, L, len, dia, maxP, maxG, maxL \n" - " p: path distance from soma, measured along dendrite, in metres.\n" - " g: geometrical distance from soma, in metres.\n" - " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" - " len: length of compartment, in metres.\n" - " dia: for diameter of compartment, in metres.\n" - " maxP: Maximum value of *p* for this neuron. \n" - " maxG: Maximum value of *g* for this neuron. \n" - " maxL: Maximum value of *L* for this neuron.\n" - "The expression for the *spacing* field must evaluate to > 0 for " - "the spine to be installed. For example, if the expresssion is\n" - " H(1 - L) \n" - "then the system will only put spines closer than " - "one length constant from the soma, and zero elsewhere. \n" - "Available spine parameters are: \n" - "spacing, minSpacing, size, sizeDistrib " - "angle, angleDistrib \n" - "minSpacing sets the granularity of sampling (typically about 0.1*" - "spacing) for the usual case where spines are spaced randomly. " - "If minSpacing < 0 then the spines are spaced equally at " - "'spacing', unless the dendritic segment length is smaller than " - "'spacing'. In that case it falls back to the regular random " - "placement method.", + "spineDistribution", + "Specification for distribution of spines on dendrite. \n" + "Arguments: proto path spacing expr [field expr]...\n" + " Each entry is terminated with an empty string. " + "The *prototype* is any spine object created in */library*, \n" + "The *path* is the wildcard path of compartments on which to " + "place the spine.\n" + "The *spacing* is the spacing of spines, in metres. \n" + "The *expression* argument is a mathematical expression in " + "the muparser framework, which permits most operations including " + "trig and transcendental ones. Of course it also handles simple " + "numerical values like 1.0, 1e-10 and so on. " + "The paired arguments are as follows: \n" + "The *field* argument specifies the name of the parameter " + "that is to be assigned by the expression.\n" + "The *expression* argument is a mathematical expression as above. " + "Available arguments for muParser are:\n" + " p, g, L, len, dia, maxP, maxG, maxL \n" + " p: path distance from soma, measured along dendrite, in metres.\n" + " g: geometrical distance from soma, in metres.\n" + " L: electrotonic distance (# of lambdas) from soma, along dend. No units.\n" + " len: length of compartment, in metres.\n" + " dia: for diameter of compartment, in metres.\n" + " maxP: Maximum value of *p* for this neuron. \n" + " maxG: Maximum value of *g* for this neuron. \n" + " maxL: Maximum value of *L* for this neuron.\n" + "The expression for the *spacing* field must evaluate to > 0 for " + "the spine to be installed. For example, if the expresssion is\n" + " H(1 - L) \n" + "then the system will only put spines closer than " + "one length constant from the soma, and zero elsewhere. \n" + "Available spine parameters are: \n" + "spacing, minSpacing, size, sizeDistrib " + "angle, angleDistrib \n" + "minSpacing sets the granularity of sampling (typically about 0.1*" + "spacing) for the usual case where spines are spaced randomly. " + "If minSpacing < 0 then the spines are spaced equally at " + "'spacing', unless the dendritic segment length is smaller than " + "'spacing'. In that case it falls back to the regular random " + "placement method.", &Neuron::setSpineDistribution, &Neuron::getSpineDistribution - ); + ); static ReadOnlyValueFinfo< Neuron, unsigned int > numCompartments( - "numCompartments", - "Number of electrical compartments in model. ", - &Neuron::getNumCompartments - ); + "numCompartments", + "Number of electrical compartments in model. ", + &Neuron::getNumCompartments + ); static ReadOnlyValueFinfo< Neuron, unsigned int > numSpines( - "numSpines", - "Number of dendritic spines in model. ", - &Neuron::getNumSpines - ); + "numSpines", + "Number of dendritic spines in model. ", + &Neuron::getNumSpines + ); static ReadOnlyValueFinfo< Neuron, unsigned int > numBranches( - "numBranches", - "Number of branches in dendrites. ", - &Neuron::getNumBranches - ); + "numBranches", + "Number of branches in dendrites. ", + &Neuron::getNumBranches + ); static ReadOnlyValueFinfo< Neuron, vector< double > > pathDistFromSoma( - "pathDistanceFromSoma", - "geometrical path distance of each segment from soma, measured by " - "threading along the dendrite.", - &Neuron::getPathDistFromSoma - ); + "pathDistanceFromSoma", + "geometrical path distance of each segment from soma, measured by " + "threading along the dendrite.", + &Neuron::getPathDistFromSoma + ); static ReadOnlyValueFinfo< Neuron, vector< double > > geomDistFromSoma( - "geometricalDistanceFromSoma", - "geometrical distance of each segment from soma.", - &Neuron::getGeomDistFromSoma - ); + "geometricalDistanceFromSoma", + "geometrical distance of each segment from soma.", + &Neuron::getGeomDistFromSoma + ); static ReadOnlyValueFinfo< Neuron, vector< double > > elecDistFromSoma( - "electrotonicDistanceFromSoma", - "geometrical distance of each segment from soma, as measured along " - "the dendrite.", - &Neuron::getElecDistFromSoma - ); + "electrotonicDistanceFromSoma", + "geometrical distance of each segment from soma, as measured along " + "the dendrite.", + &Neuron::getElecDistFromSoma + ); static ReadOnlyValueFinfo< Neuron, vector< ObjId > > compartments( - "compartments", - "Vector of ObjIds of electrical compartments. Order matches order " - "of segments, and also matches the order of the electrotonic and " - "geometricalDistanceFromSoma vectors. ", - &Neuron::getCompartments - ); + "compartments", + "Vector of ObjIds of electrical compartments. Order matches order " + "of segments, and also matches the order of the electrotonic and " + "geometricalDistanceFromSoma vectors. ", + &Neuron::getCompartments + ); static ReadOnlyLookupElementValueFinfo< Neuron, string, vector< ObjId > > - compartmentsFromExpression( - "compartmentsFromExpression", - "Vector of ObjIds of electrical compartments that match the " - "'path expression' pair in the argument string.", - &Neuron::getExprElist - ); + compartmentsFromExpression( + "compartmentsFromExpression", + "Vector of ObjIds of electrical compartments that match the " + "'path expression' pair in the argument string.", + &Neuron::getExprElist + ); static ReadOnlyLookupElementValueFinfo< Neuron, string, vector< double > > - valuesFromExpression( - "valuesFromExpression", - "Vector of values computed for each electrical compartment that " - "matches the 'path expression' pair in the argument string." - "This has 13 times the number of entries as # of compartments." - "For each compartment the entries are: \n" - "val, p, g, L, len, dia, maxP, maxG, maxL, x, y, z, 0", - &Neuron::getExprVal - ); + valuesFromExpression( + "valuesFromExpression", + "Vector of values computed for each electrical compartment that " + "matches the 'path expression' pair in the argument string." + "This has 13 times the number of entries as # of compartments." + "For each compartment the entries are: \n" + "val, p, g, L, len, dia, maxP, maxG, maxL, x, y, z, 0", + &Neuron::getExprVal + ); static ReadOnlyLookupElementValueFinfo< Neuron, string, vector< ObjId > > - spinesFromExpression( - "spinesFromExpression", - //"Vector of ObjIds of spines/heads sitting on the electrical " - //"compartments that match the 'path expression' pair in the " - //"argument string.", - "Vector of ObjIds of compartments comprising spines/heads " - "that match the 'path expression' pair in the " - "argument string.", - &Neuron::getSpinesFromExpression - ); + spinesFromExpression( + "spinesFromExpression", + //"Vector of ObjIds of spines/heads sitting on the electrical " + //"compartments that match the 'path expression' pair in the " + //"argument string.", + "Vector of ObjIds of compartments comprising spines/heads " + "that match the 'path expression' pair in the " + "argument string.", + &Neuron::getSpinesFromExpression + ); static ReadOnlyLookupElementValueFinfo< Neuron, ObjId,vector< ObjId > > - spinesOnCompartment( - "spinesOnCompartment", - "Vector of ObjIds of spines shafts/heads sitting on the specified " - "electrical compartment. If each spine has a shaft and a head," - "and there are 10 spines on the compartment, there will be 20 " - "entries in the returned vector, ordered " - "shaft0, head0, shaft1, head1, ... ", - &Neuron::getSpinesOnCompartment - ); + spinesOnCompartment( + "spinesOnCompartment", + "Vector of ObjIds of spines shafts/heads sitting on the specified " + "electrical compartment. If each spine has a shaft and a head," + "and there are 10 spines on the compartment, there will be 20 " + "entries in the returned vector, ordered " + "shaft0, head0, shaft1, head1, ... ", + &Neuron::getSpinesOnCompartment + ); static ReadOnlyLookupElementValueFinfo< Neuron, ObjId, ObjId > - parentCompartmentOfSpine( - "parentCompartmentOfSpine", - "Returns parent compartment of specified spine compartment." - "Both the spine head or its shaft will return the same parent.", - &Neuron::getParentCompartmentOfSpine - ); + parentCompartmentOfSpine( + "parentCompartmentOfSpine", + "Returns parent compartment of specified spine compartment." + "Both the spine head or its shaft will return the same parent.", + &Neuron::getParentCompartmentOfSpine + ); static ReadOnlyLookupElementValueFinfo< Neuron, vector< ObjId >, vector< ObjId > > - spineIdsFromCompartmentIds( - "spineIdsFromCompartmentIds", - "Vector of ObjIds of spine entries (FieldElements on this Neuron, " - "used for scaling) that map to the the specified " - "electrical compartments. If a bad compartment Id is given, the" - "corresponding spine entry is the root Id.", - &Neuron::getSpineIdsFromCompartmentIds - ); + spineIdsFromCompartmentIds( + "spineIdsFromCompartmentIds", + "Vector of ObjIds of spine entries (FieldElements on this Neuron, " + "used for scaling) that map to the the specified " + "electrical compartments. If a bad compartment Id is given, the" + "corresponding spine entry is the root Id.", + &Neuron::getSpineIdsFromCompartmentIds + ); ///////////////////////////////////////////////////////////////////// // DestFinfos ///////////////////////////////////////////////////////////////////// static DestFinfo buildSegmentTree( "buildSegmentTree", - "Build the reference segment tree structure using the child " - "compartments of the current Neuron. Fills in all the coords and " - "length constant information into the segments, for later use " - "when we build reduced compartment trees and channel " - "distributions. Should only be called once, since subsequent use " - "on a reduced model will lose the original full cell geometry. ", - new EpFunc0< Neuron >( &Neuron::buildSegmentTree ) - ); + "Build the reference segment tree structure using the child " + "compartments of the current Neuron. Fills in all the coords and " + "length constant information into the segments, for later use " + "when we build reduced compartment trees and channel " + "distributions. Should only be called once, since subsequent use " + "on a reduced model will lose the original full cell geometry. ", + new EpFunc0< Neuron >( &Neuron::buildSegmentTree ) + ); static DestFinfo setSpineAndPsdMesh( "setSpineAndPsdMesh", - "Assigns the spine and psd mesh to the Neuron. This is used " - "to build up a mapping from Spine entries on the Neuron to " - "chem spines and PSDs, so that volume change operations from " - "the Spine can propagate to the chem systems.", - new OpFunc2< Neuron, Id, Id >( &Neuron::setSpineAndPsdMesh ) - ); + "Assigns the spine and psd mesh to the Neuron. This is used " + "to build up a mapping from Spine entries on the Neuron to " + "chem spines and PSDs, so that volume change operations from " + "the Spine can propagate to the chem systems.", + new OpFunc2< Neuron, Id, Id >( &Neuron::setSpineAndPsdMesh ) + ); static DestFinfo setSpineAndPsdDsolve( "setSpineAndPsdDsolve", - "Assigns the Dsolves used by spine and PSD to the Neuron. " - "This is used " - "to handle the rescaling of diffusion rates when spines are " - "resized. ", - new OpFunc2< Neuron, Id, Id >( &Neuron::setSpineAndPsdDsolve ) - ); + "Assigns the Dsolves used by spine and PSD to the Neuron. " + "This is used " + "to handle the rescaling of diffusion rates when spines are " + "resized. ", + new OpFunc2< Neuron, Id, Id >( &Neuron::setSpineAndPsdDsolve ) + ); /* static DestFinfo rotateInSpace( "rotateInSpace", diff --git a/builtins/CMakeLists.txt b/builtins/CMakeLists.txt index 8b96a98ede..bcba8eae04 100644 --- a/builtins/CMakeLists.txt +++ b/builtins/CMakeLists.txt @@ -51,7 +51,8 @@ set(SRCS Arith.cpp Group.cpp Mstring.cpp - Func.cpp + # Func.cpp + MooseParser.cpp Function.cpp Variable.cpp InputVariable.cpp diff --git a/builtins/Function.cpp b/builtins/Function.cpp index d034a2c96e..1948e441ac 100644 --- a/builtins/Function.cpp +++ b/builtins/Function.cpp @@ -59,29 +59,29 @@ static const double TriggerThreshold = 0.0; static SrcFinfo1 *valueOut() { static SrcFinfo1 valueOut("valueOut", - "Evaluated value of the function for the current variable values."); + "Evaluated value of the function for the current variable values."); return &valueOut; } static SrcFinfo1< double > *derivativeOut() { static SrcFinfo1< double > derivativeOut("derivativeOut", - "Value of derivative of the function for the current variable values"); + "Value of derivative of the function for the current variable values"); return &derivativeOut; } static SrcFinfo1< double > *rateOut() { static SrcFinfo1< double > rateOut("rateOut", - "Value of time-derivative of the function for the current variable values"); + "Value of time-derivative of the function for the current variable values"); return &rateOut; } static SrcFinfo1< vector < double > *> *requestOut() { static SrcFinfo1< vector < double > * > requestOut( - "requestOut", - "Sends request for input variable from a field on target object"); + "requestOut", + "Sends request for input variable from a field on target object"); return &requestOut; } @@ -168,8 +168,9 @@ const Cinfo * Function::initCinfo() "max var. max of all arguments\n" "sum var. sum of all arguments\n" "avg var. mean value of all arguments\n" + "rnd 0 rand(), random float between 0 and 1, honors global moose.seed.\n" "rand 1 rand(seed), random float between 0 and 1, \n" - " if seed = -1, then a 'random' seed is created.\n" + " if seed = -1, then a 'random' seed is used.\n" "rand2 3 rand(a, b, seed), random float between a and b, \n" " if seed = -1, a 'random' seed is created using either\n" " by random_device or by reading system clock\n" @@ -235,89 +236,89 @@ const Cinfo * Function::initCinfo() // Shared messages /////////////////////////////////////////////////////////////////// static DestFinfo process( "process", - "Handles process call, updates internal time stamp.", - new ProcOpFunc< Function >( &Function::process ) ); + "Handles process call, updates internal time stamp.", + new ProcOpFunc< Function >( &Function::process ) ); static DestFinfo reinit( "reinit", - "Handles reinit call.", - new ProcOpFunc< Function >( &Function::reinit ) ); + "Handles reinit call.", + new ProcOpFunc< Function >( &Function::reinit ) ); static Finfo* processShared[] = { &process, &reinit }; static SharedFinfo proc( "proc", - "This is a shared message to receive Process messages " - "from the scheduler objects." - "The first entry in the shared msg is a MsgDest " - "for the Process operation. It has a single argument, " - "ProcInfo, which holds lots of information about current " - "time, thread, dt and so on. The second entry is a MsgDest " - "for the Reinit operation. It also uses ProcInfo. ", - processShared, sizeof( processShared ) / sizeof( Finfo* ) - ); - /* - static DestFinfo trigger( "trigger", - "Handles trigger input. Argument is timestamp of event. This is " - "compatible with spike events as well as chemical ones. ", - new OpFunc1< Function, double >( &Function::trigger ) ); - */ + "This is a shared message to receive Process messages " + "from the scheduler objects." + "The first entry in the shared msg is a MsgDest " + "for the Process operation. It has a single argument, " + "ProcInfo, which holds lots of information about current " + "time, thread, dt and so on. The second entry is a MsgDest " + "for the Reinit operation. It also uses ProcInfo. ", + processShared, sizeof( processShared ) / sizeof( Finfo* ) + ); + /* + static DestFinfo trigger( "trigger", + "Handles trigger input. Argument is timestamp of event. This is " + "compatible with spike events as well as chemical ones. ", + new OpFunc1< Function, double >( &Function::trigger ) ); + */ static Finfo *functionFinfos[] = - { - &value, - &rate, - &derivative, - &mode, - &useTrigger, - &doEvalAtReinit, - &expr, - &numVars, - &inputs, - &constants, - &independent, - &proc, - requestOut(), - valueOut(), - rateOut(), - derivativeOut(), - }; + { + &value, + &rate, + &derivative, + &mode, + &useTrigger, + &doEvalAtReinit, + &expr, + &numVars, + &inputs, + &constants, + &independent, + &proc, + requestOut(), + valueOut(), + rateOut(), + derivativeOut(), + }; static string doc[] = - { - "Name", "Function", - "Author", "Subhasis Ray", - "Description", - "General purpose function calculator using real numbers.\n" - "It can parse mathematical expression defining a function and evaluate" - " it and/or its derivative for specified variable values." - "You can assign expressions of the form::\n" - "\n" - "f(c0, c1, ..., cM, x0, x1, ..., xN, y0,..., yP ) \n" - "\n" - " where `ci`'s are constants and `xi`'s and `yi`'s are variables." - - "The constants must be defined before setting the expression and" - " variables are connected via messages. The constants can have any" - " name, but the variable names must be of the form x{i} or y{i}" - " where i is increasing integer starting from 0.\n" - " The variables can be input from other moose objects." - " Such variables must be named `x{i}` in the expression and the source" - " field is connected to Function.x[i]'s `input` destination field.\n" - " In case the input variable is not available as a source field, but is" - " a value field, then the value can be requested by connecting the" - " `requestOut` message to the `get{Field}` destination on the target" - " object. Such variables must be specified in the expression as y{i}" - " and connecting the messages should happen in the same order as the" - " y indices.\n" - " This class handles only real numbers (C-double). Predefined constants" - " are: pi=3.141592..., e=2.718281..." - }; + { + "Name", "Function", + "Author", "Subhasis Ray", + "Description", + "General purpose function calculator using real numbers.\n" + "It can parse mathematical expression defining a function and evaluate" + " it and/or its derivative for specified variable values." + "You can assign expressions of the form::\n" + "\n" + "f(c0, c1, ..., cM, x0, x1, ..., xN, y0,..., yP ) \n" + "\n" + " where `ci`'s are constants and `xi`'s and `yi`'s are variables." + + "The constants must be defined before setting the expression and" + " variables are connected via messages. The constants can have any" + " name, but the variable names must be of the form x{i} or y{i}" + " where i is increasing integer starting from 0.\n" + " The variables can be input from other moose objects." + " Such variables must be named `x{i}` in the expression and the source" + " field is connected to Function.x[i]'s `input` destination field.\n" + " In case the input variable is not available as a source field, but is" + " a value field, then the value can be requested by connecting the" + " `requestOut` message to the `get{Field}` destination on the target" + " object. Such variables must be specified in the expression as y{i}" + " and connecting the messages should happen in the same order as the" + " y indices.\n" + " This class handles only real numbers (C-double). Predefined constants" + " are: pi=3.141592..., e=2.718281..." + }; static Dinfo< Function > dinfo; static Cinfo functionCinfo("Function", - Neutral::initCinfo(), - functionFinfos, - sizeof(functionFinfos) / sizeof(Finfo*), - &dinfo, - doc, - sizeof(doc)/sizeof(string)); + Neutral::initCinfo(), + functionFinfos, + sizeof(functionFinfos) / sizeof(Finfo*), + &dinfo, + doc, + sizeof(doc)/sizeof(string)); return &functionCinfo; } diff --git a/builtins/Function.h b/builtins/Function.h index 1976103d8a..082e78c20a 100644 --- a/builtins/Function.h +++ b/builtins/Function.h @@ -12,43 +12,12 @@ // URL: // Keywords: // Compatibility: -// -// - -// Commentary: -// -// A new version of Func with FieldElements to collect data. -// -// - -// Change log: -// -// -// -// -// This program is free software; you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 3, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program; see the file COPYING. If not, write to -// the Free Software Foundation, Inc., 51 Franklin Street, Fifth -// Floor, Boston, MA 02110-1301, USA. -// -// -// Code: #ifndef _MOOSE_FUNCTION_H_ #define _MOOSE_FUNCTION_H_ -#include "../external/muparser/include/muParser.h" +#include "MooseParser.h" /** Simple function parser and evaluator for MOOSE. This can take a mathematical @@ -158,7 +127,7 @@ class Function map< string, double *> _constbuf; // for constants string _independent; // index of independent variable - mu::Parser _parser; + MooseParser _parser; void _clearBuffer(); void _showError(mu::Parser::exception_type &e) const; diff --git a/builtins/MooseParser.cpp b/builtins/MooseParser.cpp new file mode 100644 index 0000000000..b8cdcb74fa --- /dev/null +++ b/builtins/MooseParser.cpp @@ -0,0 +1,18 @@ +/*** + * Description: Class MooseParser. + * + * Created: 2019-05-30 + + * Author: Dilawar Singh + * Organization: NCBS Bangalore + * License: MIT License + */ + +#include "MooseParser.h" + +MooseParser::MooseParser() +{} + +MooseParser::~MooseParser() +{} + diff --git a/builtins/MooseParser.h b/builtins/MooseParser.h new file mode 100644 index 0000000000..9d67b9709e --- /dev/null +++ b/builtins/MooseParser.h @@ -0,0 +1,26 @@ +/*** + * Description: MooseParser class. + * + * Created: 2019-05-30 + + * Author: Dilawar Singh + * Organization: NCBS Bangalore + * License: MIT License + */ + +#ifndef MOOSEPARSER_H +#define MOOSEPARSER_H + +#include "../external/muparser/include/muParser.h" + +class MooseParser : public mu::Parser +{ + public: + MooseParser(); + ~MooseParser(); + + private: + /* data */ +}; + +#endif /* end of include guard: MOOSEPARSER_H */ diff --git a/external/muparser/include/muParser.h b/external/muparser/include/muParser.h index 04b38018bf..c1baf79846 100644 --- a/external/muparser/include/muParser.h +++ b/external/muparser/include/muParser.h @@ -107,7 +107,8 @@ namespace mu static value_type Rand2(value_type, value_type, value_type); // Random number between 0 and 1, non-deterministic seed. - static value_type Rand( value_type ); + static value_type Rand( value_type seed ); + static value_type Rnd( ); // Prefix operators // !!! Unary Minus is a MUST if you want to use negative signs !!! diff --git a/external/muparser/src/muParser.cpp b/external/muparser/src/muParser.cpp index 4ef2416058..2f9321ee14 100644 --- a/external/muparser/src/muParser.cpp +++ b/external/muparser/src/muParser.cpp @@ -25,6 +25,7 @@ */ #include "../include/muParser.h" #include "../include/muParserTemplateMagic.h" +#include "../../../basecode/global.h" //--- Standard includes ------------------------------------------------------------------------ #include @@ -112,10 +113,24 @@ namespace mu value_type Parser::Fmod(value_type v1, value_type v2) { return fmod(v1, v2); } value_type Parser::Quot(value_type v1, value_type v2) { return (int)(v1 / v2); } + value_type Parser::Rand( value_type seed = -1 ) + { + static bool isSeedSet_ = false; + + if( ! isSeedSet_ ) + { + mu::rng.setSeed( seed ); + isSeedSet_ = true; + } + return rng.uniform( ); /* Between 0 and 1 */ + } + // If no seed is given, - value_type Parser::Rand( value_type seed ) + value_type Parser::Rnd( ) { static bool isSeedSet_ = false; + // check if global seed is set + size_t seed = moose::getGlobalSeed(); if( ! isSeedSet_ ) { mu::rng.setSeed( seed ); @@ -330,6 +345,7 @@ namespace mu DefineFun(_T("rint"), Rint); DefineFun(_T("abs"), Abs); DefineFun(_T("fmod"), Fmod); + DefineFun(_T("rnd"), Rnd); DefineFun(_T("rand"), Rand); DefineFun(_T("rand2"), Rand2); // Functions with variable number of arguments diff --git a/ksolve/FuncTerm.h b/ksolve/FuncTerm.h index 2b83686bbd..71bf2de725 100644 --- a/ksolve/FuncTerm.h +++ b/ksolve/FuncTerm.h @@ -10,7 +10,7 @@ #ifndef _FUNC_TERM_H #define _FUNC_TERM_H -#include "../external/muparser/include/muParser.h" +#include "../builtins/MooseParser.h" class FuncTerm { @@ -44,7 +44,7 @@ class FuncTerm double* args_; // Look up reactants in the S vec. vector< unsigned int > reactantIndex_; - mu::Parser parser_; + MooseParser parser_; string expr_; /** * Scale factor to account for pool volume if we are assigning conc diff --git a/tests/alpha/ex7.6_func_controls_reac_rate.py b/tests/alpha/ex7.6_func_controls_reac_rate.py new file mode 100644 index 0000000000..b8c8fb4c64 --- /dev/null +++ b/tests/alpha/ex7.6_func_controls_reac_rate.py @@ -0,0 +1,67 @@ +######################################################################## +# This example illustrates how a function can be used to control a reaction +# rate. This kind of calculation is appropriate when we need to link +# different kinds of physical processses with chemical reactions, for +# example, membrane curvature with molecule accumulation. The use of +# functions to modify reaction rates should be avoided in purely chemical +# systems since they obscure the underlying chemistry, and do not map +# cleanly to stochastic calculations. +# +# In this example we simply have a molecule C that controls the forward +# rate of a reaction that converts A to B. C is a function of location +# on the cylinder, and is fixed. In more elaborate computations we could +# have a function of multiple molecules, some of which could be changing and +# others could be buffered. +# +# Copyright (C) Upinder S. Bhalla NCBS 2018 +# Released under the terms of the GNU Public License V3. +######################################################################## + +import numpy as np +import moose +import pylab +import rdesigneur as rd + + +def makeFuncRate(): + model = moose.Neutral( '/library' ) + model = moose.Neutral( '/library/chem' ) + compt = moose.CubeMesh( '/library/chem/compt' ) + compt.volume = 1e-15 + A = moose.Pool( '/library/chem/compt/A' ) + B = moose.Pool( '/library/chem/compt/B' ) + C = moose.Pool( '/library/chem/compt/C' ) + reac = moose.Reac( '/library/chem/compt/reac' ) + func = moose.Function( '/library/chem/compt/reac/func' ) + func.x.num = 1 + func.expr = "(x0/1e8)^2" + moose.connect( C, 'nOut', func.x[0], 'input' ) + moose.connect( func, 'valueOut', reac, 'setNumKf' ) + moose.connect( reac, 'sub', A, 'reac' ) + moose.connect( reac, 'prd', B, 'reac' ) + + A.concInit = 1 + B.concInit = 0 + C.concInit = 0 + reac.Kb = 1 + + +makeFuncRate() + +rdes = rd.rdesigneur( + turnOffElec = True, + #This subdivides the 50-micron cylinder into 2 micron voxels + diffusionLength = 2e-6, + cellProto = [['somaProto', 'soma', 5e-6, 50e-6]], + chemProto = [['chem', 'chem']], + chemDistrib = [['chem', 'soma', 'install', '1' ]], + plotList = [['soma', '1', 'dend/A', 'conc', 'A conc', 'wave'], + ['soma', '1', 'dend/C', 'conc', 'C conc', 'wave']], +) +rdes.buildModel() + +C = moose.element( '/model/chem/dend/C' ) +C.vec.concInit = [ 1+np.sin(x/5.0) for x in range( len(C.vec) ) ] +moose.reinit() +moose.start(10) +rdes.display() diff --git a/tests/alpha/function.py b/tests/alpha/function.py new file mode 100644 index 0000000000..1dc982bf35 --- /dev/null +++ b/tests/alpha/function.py @@ -0,0 +1,160 @@ +# function.py --- + +import numpy as np +import sys +import moose + +simtime = 1.0 +plot_ = False +if plot_: + import matplotlib.pyplot as plt + +def example(): + """Function objects can be used to evaluate expressions with arbitrary + number of variables and constants. We can assign expression of the + form:: + + f(c0, c1, ..., cM, x0, x1, ..., xN, y0,..., yP ) + + where `c_i`'s are constants and `x_i`'s and `y_i`'s are variables. + + The constants must be defined before setting the expression and + variables are connected via messages. The constants can have any + name, but the variable names must be of the form x{i} or y{i} + where i is increasing integer starting from 0. + + The `x_i`'s are field elements and you have to set their number + first (function.x.num = N). Then you can connect any source field + sending out double to the 'input' destination field of the + `x[i]`. + + The `y_i`'s are useful when the required variable is a value field + and is not available as a source field. In that case you connect + the `requestOut` source field of the function element to the + `get{Field}` destination field on the target element. The `y_i`'s + are automatically added on connecting. Thus, if you call:: + + moose.connect(function, 'requestOut', a, 'getSomeField') + moose.connect(function, 'requestOut', b, 'getSomeField') + + then ``a.someField`` will be assigned to ``y0`` and + ``b.someField`` will be assigned to ``y1``. + + In this example we evaluate the expression: ``z = c0 * exp(c1 * + x0) * cos(y0)`` + + with x0 ranging from -1 to +1 and y0 ranging from -pi to + +pi. These values are stored in two stimulus tables called xtab + and ytab respectively, so that at each timestep the next values of + x0 and y0 are assigned to the function. + + Along with the value of the expression itself we also compute its + derivative with respect to y0 and its derivative with respect to + time (rate). The former uses a five-point stencil for the + numerical differentiation and has a glitch at y=0. The latter uses + backward difference divided by dt. + + Unlike Func class, the number of variables and constants are + unlimited in Function and you can set all the variables via + messages. + + """ + demo = moose.Neutral('/model') + function = moose.Function('/model/function') + function.c['c0'] = 1.0 + function.c['c1'] = 2.0 + function.x.num = 1 + function.expr = 'c0 * exp(c1*x0) * cos(y0) + sin(t)' + # mode 0 - evaluate function value, derivative and rate + # mode 1 - just evaluate function value, + # mode 2 - evaluate derivative, + # mode 3 - evaluate rate + function.mode = 0 + function.independent = 'y0' + nsteps = 10000 + xarr = np.linspace(0.0, 1.0, nsteps) + # Stimulus tables allow you to store sequences of numbers which + # are delivered via the 'output' message at each time step. This + # is a placeholder and in real scenario you will be using any + # sourceFinfo that sends out a double value. + input_x = moose.StimulusTable('/xtab') + input_x.vector = xarr + input_x.startTime = 0.0 + input_x.stepPosition = xarr[0] + input_x.stopTime = simtime + moose.connect(input_x, 'output', function.x[0], 'input') + + yarr = np.linspace(-np.pi, np.pi, nsteps) + input_y = moose.StimulusTable('/ytab') + input_y.vector = yarr + input_y.startTime = 0.0 + input_y.stepPosition = yarr[0] + input_y.stopTime = simtime + moose.connect(function, 'requestOut', input_y, 'getOutputValue') + + # data recording + result = moose.Table('/ztab') + moose.connect(result, 'requestOut', function, 'getValue') + derivative = moose.Table('/zprime') + moose.connect(derivative, 'requestOut', function, 'getDerivative') + rate = moose.Table('/dz_by_dt') + moose.connect(rate, 'requestOut', function, 'getRate') + x_rec = moose.Table('/xrec') + moose.connect(x_rec, 'requestOut', input_x, 'getOutputValue') + y_rec = moose.Table('/yrec') + moose.connect(y_rec, 'requestOut', input_y, 'getOutputValue') + + dt = simtime/nsteps + for ii in range(32): + moose.setClock(ii, dt) + moose.reinit() + moose.start(simtime) + + # Uncomment the following lines and the import matplotlib.pyplot as plt on top + # of this file to display the plot. + if plot_: + plt.plot(x_rec.vector, result.vector, 'r-', label='z = {}'.format(function.expr)) + plt.subplot(4,1,1) + z = function.c['c0'] * np.exp(function.c['c1'] * xarr) * np.cos(yarr) + np.sin(np.arange(len(xarr)) * dt) + zz = result.vector[1:] + err = z[1:] - zz[1:] + assert (np.abs(err) <= 0.05).all(), err + assert (np.mean(err) <= 0.001), np.mean(err) + + if plot_: + plt.plot(xarr, z, 'b--', label='numpy computed') + plt.xlabel('x') + plt.ylabel('z') + plt.legend() + + if plot_: + plt.subplot(4,1,2) + plt.plot(y_rec.vector, derivative.vector, 'r-', label='dz/dy0') + # derivatives computed by putting x values in the analytical formula + dzdy = function.c['c0'] * np.exp(function.c['c1'] * xarr) * (- np.sin(yarr)) + err = np.abs(dzdy - derivative.vector[1:]) + assert (err < 1e-2).all(), err + assert (np.mean(err) < 1e-3), np.mean(err) + + if plot_: + plt.plot(yarr, dzdy, 'b--', label='numpy computed') + plt.xlabel('y') + plt.ylabel('dz/dy') + plt.legend() + + if plot_: + plt.subplot(4,1,3) + # *** BEWARE *** The first two entries are spurious. Entry 0 is + # *** from reinit sending out the defaults. Entry 2 is because + # *** there is no lastValue for computing real forward difference. + plt.plot(np.arange(2, len(rate.vector), 1) * dt, rate.vector[2:], 'r-', label='dz/dt') + dzdt = np.diff(z)/dt + plt.plot(np.arange(0, len(dzdt), 1.0) * dt, dzdt, 'b--', label='numpy computed') + plt.xlabel('t') + plt.ylabel('dz/dt') + plt.legend() + plt.tight_layout() + plt.show() + +if __name__ == '__main__': + example() diff --git a/tests/python/test_Xreacs2.py b/tests/python/test_Xreacs2.py index 2fb36cfb16..c02a2cc3f3 100644 --- a/tests/python/test_Xreacs2.py +++ b/tests/python/test_Xreacs2.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- import os -import sys import moose +print( "[INFO ] Using moose from %s" % moose.__file__ ) import numpy as np import moose.fixXreacs as fixXreacs diff --git a/tests/python/test_expr_parser.py b/tests/python/test_expr_parser.py new file mode 100644 index 0000000000..fabbc28e40 --- /dev/null +++ b/tests/python/test_expr_parser.py @@ -0,0 +1,346 @@ +# -*- coding: utf-8 -*- +"""test_muparser.py: + +odified from https://elifesciences.org/articles/25827 +Fig4.py +""" +import sys +import re +import numpy as np +import moose +import rdesigneur as rd +print( "[INFO ] Using moose from %s" % moose.__file__ ) + +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 = None ): + 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 = None ): + 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 = None ): + 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 = None ): + 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 + + +def singleCompt( name, params ): + print('=============') + print('[INFO] Making compartment %s' % name) + 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' ) + + print( '\n\n[INFO] CaStim %s' % CaStim.path ) + runtime = params['preStimTime'] + params['postStimTime'] + steptime = 50 + CaStim.expr += '+x2*(t>100+'+str(runtime)+')*(t<100+'+str(runtime+steptime)+ ')' + print("[INFO] CaStim.expr = %s" % 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 + '/##' + + print( 'REINIT AND START' ) + moose.reinit() + runtime += 100 + steptime*2 + moose.start( runtime ) + t = np.arange( 0, runtime + 1e-9, tab.dt ) + return name, t, tab.vector + + +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 run(): + panelC = [] + panelCticks = [] + panelC.append( singleCompt( 'negFB', makeNegFB( [] ) ) ) + panelC.append( singleCompt( 'negFF', makeNegFF( [] ) ) ) + panelC.append( singleCompt( 'fhn', makeFHN( [] ) ) ) + panelC.append( singleCompt( 'bis', 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 ) + print( '[INFO] Run is over') + +if __name__ == '__main__': + moose.Neutral( '/library' ) + moose.Neutral( '/model' ) + run() + print( 'All done' ) diff --git a/tests/python/test_function.py b/tests/python/test_function.py index 9066775d57..79e9453424 100644 --- a/tests/python/test_function.py +++ b/tests/python/test_function.py @@ -1,70 +1,37 @@ # -*- coding: utf-8 -*- # test_function.py --- -# # Filename: test_function.py # Description: # Author: subha -# Maintainer: +# Maintainer: Dilawar Singh # Created: Sat Mar 28 19:34:20 2015 (-0400) # Version: -# Last-Updated: -# By: -# Update #: 0 -# URL: -# Keywords: -# Compatibility: -# -# -# Commentary: -# -# -# -# - -# Change log: -# -# -# -# -# This program is free software; you can redistribute it and/or -# modify it under the terms of the GNU General Public License as -# published by the Free Software Foundation; either version 3, or -# (at your option) any later version. -# -# This program is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -# General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program; see the file COPYING. If not, write to -# the Free Software Foundation, Inc., 51 Franklin Street, Fifth -# Floor, Boston, MA 02110-1301, USA. -# -# - -# Code: -"""Check variable ordering - bug #161""" from __future__ import print_function - import numpy as np import moose print( "[INFO ] Using moose %s form %s" % (moose.version(), moose.__file__) ) +def create_func( funcname, expr ): + f = moose.Function( funcname ) + f.expr = expr + t = moose.Table( funcname + 'tab' ) + moose.connect( t, 'requestOut', f, 'getValue' ) + moose.setClock( f.tick, 0.1) + moose.setClock( t.tick, 0.1) + return f, t + def test_var_order(): """The y values are one step behind the x values because of scheduling sequences""" + moose.delete( '/' ) nsteps = 5 simtime = nsteps dt = 1.0 # fn0 = moose.Function('/fn0') - # fn0.x.num = 2 - # fn0.expr = 'x0 + x1' - # fn0.mode = 1 fn1 = moose.Function('/fn1') fn1.x.num = 2 - fn1.expr = 'y1 + y0 + x1 + x0' + fn1.expr = 'y1+y0+x1+x0' fn1.mode = 1 inputs = np.arange(0, nsteps+1, 1.0) x0 = moose.StimulusTable('/x0') @@ -91,6 +58,8 @@ def test_var_order(): y1.startTime = 0.0 y1.stopTime = simtime y1.stepPosition = 0.0 + # print( fn1, type(fn1) ) + # print( moose.showmsg(fn1.x) ) moose.connect(x0, 'output', fn1.x[0], 'input') moose.connect(x1, 'output', fn1.x[1], 'input') moose.connect(fn1, 'requestOut', y0, 'getOutputValue') @@ -101,28 +70,54 @@ def test_var_order(): moose.setClock(ii, dt) moose.reinit() moose.start(simtime) - for ii in range(len(z1.vector)): - print(ii, z1.vector[ii]) expected = [0, 1.1, 2.211, 3.322, 4.433, 5.544] assert np.allclose(z1.vector, expected), "Excepted %s, got %s" % (expected, z1.vector ) + print( 'Passed order vars' ) def test_t( ): - moose.delete( '/' ) - f = moose.Function( '/funct' ) - f.expr = 't' - t = moose.Table( '/table1' ) - moose.setClock( f.tick, 0.1) - moose.setClock( t.tick, 0.1) - moose.connect( t, 'requestOut', f, 'getValue' ) + f, t = create_func( 'funct', 't/2.0') + moose.reinit() + moose.start(1) + y = t.vector + d = np.diff( y[1:] ) + assert np.mean(d) == d[0] + print( 'Passed t/2' ) + +def test_trig( ): + f, t = create_func('func2', '(sin(t)^2+cos(t)^2)-1') moose.reinit() moose.start( 1 ) - y = t.vector[1:] - print( 'Raw', y ) - print( 'Diff', np.diff(y) ) - assert np.isclose(np.diff(y), 0.1).all(), np.diff(y) + y = t.vector + print(y) + assert np.isclose(np.mean(y), 0.0), np.mean(y) + assert np.isclose(np.std(y), 0.0), np.std(y) + print( 'Passed sin^2 x + cos^x=1' ) + +def test_rand( ): + moose.seed( 10 ) + f, t = create_func( 'random', 'rnd()') + moose.reinit() + moose.start(1) + expected = [0.49458993, 0.44301495, 0.58332174, 0.70920801, 0.26360285, + 0.68381843, 0.33607158, 0.19812181, 0.87761494, 0.54088093, + 0.41366738] + assert np.isclose(t.vector, expected ).all(), t.vector + print( 'Passed test random' ) +def test_fmod( ): + f, t = create_func( 'fmod', 'fmod(t, 2)' ) + moose.reinit() + moose.start( 20 ) + y = t.vector + assert (np.fmod(y, 2) == y).all() + assert(np.isclose(np.max(y), 1.9)), np.max(y) + assert(np.isclose(np.min(y), 0.0)), np.min(y) + print('Passed fmod(t,2)') if __name__ == '__main__': - test_t( ) test_var_order() + test_t() + test_trig() + test_rand() + test_fmod() diff --git a/tests/python/test_function_chemsys.py b/tests/python/test_function_chemsys.py new file mode 100644 index 0000000000..fabbc28e40 --- /dev/null +++ b/tests/python/test_function_chemsys.py @@ -0,0 +1,346 @@ +# -*- coding: utf-8 -*- +"""test_muparser.py: + +odified from https://elifesciences.org/articles/25827 +Fig4.py +""" +import sys +import re +import numpy as np +import moose +import rdesigneur as rd +print( "[INFO ] Using moose from %s" % moose.__file__ ) + +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 = None ): + 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 = None ): + 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 = None ): + 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 = None ): + 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 + + +def singleCompt( name, params ): + print('=============') + print('[INFO] Making compartment %s' % name) + 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' ) + + print( '\n\n[INFO] CaStim %s' % CaStim.path ) + runtime = params['preStimTime'] + params['postStimTime'] + steptime = 50 + CaStim.expr += '+x2*(t>100+'+str(runtime)+')*(t<100+'+str(runtime+steptime)+ ')' + print("[INFO] CaStim.expr = %s" % 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 + '/##' + + print( 'REINIT AND START' ) + moose.reinit() + runtime += 100 + steptime*2 + moose.start( runtime ) + t = np.arange( 0, runtime + 1e-9, tab.dt ) + return name, t, tab.vector + + +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 run(): + panelC = [] + panelCticks = [] + panelC.append( singleCompt( 'negFB', makeNegFB( [] ) ) ) + panelC.append( singleCompt( 'negFF', makeNegFF( [] ) ) ) + panelC.append( singleCompt( 'fhn', makeFHN( [] ) ) ) + panelC.append( singleCompt( 'bis', 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 ) + print( '[INFO] Run is over') + +if __name__ == '__main__': + moose.Neutral( '/library' ) + moose.Neutral( '/model' ) + run() + print( 'All done' ) diff --git a/tests/python/test_function_derivative.py b/tests/python/test_function_derivative.py new file mode 100644 index 0000000000..c2e1dc81ed --- /dev/null +++ b/tests/python/test_function_derivative.py @@ -0,0 +1,160 @@ +# function.py --- + +import numpy as np +import sys +import moose + +simtime = 1.0 +plot_ = False +if plot_: + import matplotlib.pyplot as plt + +def example(): + """Function objects can be used to evaluate expressions with arbitrary + number of variables and constants. We can assign expression of the + form:: + + f(c0, c1, ..., cM, x0, x1, ..., xN, y0,..., yP ) + + where `c_i`'s are constants and `x_i`'s and `y_i`'s are variables. + + The constants must be defined before setting the expression and + variables are connected via messages. The constants can have any + name, but the variable names must be of the form x{i} or y{i} + where i is increasing integer starting from 0. + + The `x_i`'s are field elements and you have to set their number + first (function.x.num = N). Then you can connect any source field + sending out double to the 'input' destination field of the + `x[i]`. + + The `y_i`'s are useful when the required variable is a value field + and is not available as a source field. In that case you connect + the `requestOut` source field of the function element to the + `get{Field}` destination field on the target element. The `y_i`'s + are automatically added on connecting. Thus, if you call:: + + moose.connect(function, 'requestOut', a, 'getSomeField') + moose.connect(function, 'requestOut', b, 'getSomeField') + + then ``a.someField`` will be assigned to ``y0`` and + ``b.someField`` will be assigned to ``y1``. + + In this example we evaluate the expression: ``z = c0 * exp(c1 * + x0) * cos(y0)`` + + with x0 ranging from -1 to +1 and y0 ranging from -pi to + +pi. These values are stored in two stimulus tables called xtab + and ytab respectively, so that at each timestep the next values of + x0 and y0 are assigned to the function. + + Along with the value of the expression itself we also compute its + derivative with respect to y0 and its derivative with respect to + time (rate). The former uses a five-point stencil for the + numerical differentiation and has a glitch at y=0. The latter uses + backward difference divided by dt. + + Unlike Func class, the number of variables and constants are + unlimited in Function and you can set all the variables via + messages. + + """ + demo = moose.Neutral('/model') + function = moose.Function('/model/function') + function.c['c0'] = 1.0 + function.c['c1'] = 2.0 + function.x.num = 1 + function.expr = 'c0 * exp(c1*x0) * cos(y0) + sin(t)' + # mode 0 - evaluate function value, derivative and rate + # mode 1 - just evaluate function value, + # mode 2 - evaluate derivative, + # mode 3 - evaluate rate + function.mode = 0 + function.independent = 'y0' + nsteps = 1000 + xarr = np.linspace(0.0, 1.0, nsteps) + # Stimulus tables allow you to store sequences of numbers which + # are delivered via the 'output' message at each time step. This + # is a placeholder and in real scenario you will be using any + # sourceFinfo that sends out a double value. + input_x = moose.StimulusTable('/xtab') + input_x.vector = xarr + input_x.startTime = 0.0 + input_x.stepPosition = xarr[0] + input_x.stopTime = simtime + moose.connect(input_x, 'output', function.x[0], 'input') + + yarr = np.linspace(-np.pi, np.pi, nsteps) + input_y = moose.StimulusTable('/ytab') + input_y.vector = yarr + input_y.startTime = 0.0 + input_y.stepPosition = yarr[0] + input_y.stopTime = simtime + moose.connect(function, 'requestOut', input_y, 'getOutputValue') + + # data recording + result = moose.Table('/ztab') + moose.connect(result, 'requestOut', function, 'getValue') + derivative = moose.Table('/zprime') + moose.connect(derivative, 'requestOut', function, 'getDerivative') + rate = moose.Table('/dz_by_dt') + moose.connect(rate, 'requestOut', function, 'getRate') + x_rec = moose.Table('/xrec') + moose.connect(x_rec, 'requestOut', input_x, 'getOutputValue') + y_rec = moose.Table('/yrec') + moose.connect(y_rec, 'requestOut', input_y, 'getOutputValue') + + dt = simtime/nsteps + for ii in range(32): + moose.setClock(ii, dt) + moose.reinit() + moose.start(simtime) + + # Uncomment the following lines and the import matplotlib.pyplot as plt on top + # of this file to display the plot. + if plot_: + plt.plot(x_rec.vector, result.vector, 'r-', label='z = {}'.format(function.expr)) + plt.subplot(4,1,1) + z = function.c['c0'] * np.exp(function.c['c1'] * xarr) * np.cos(yarr) + np.sin(np.arange(len(xarr)) * dt) + zz = result.vector[1:] + err = z[1:] - zz[1:] + assert (np.abs(err) <= 0.05).all(), err[err > 0.05] + assert (np.mean(err) <= 0.001), np.mean(err) + + if plot_: + plt.plot(xarr, z, 'b--', label='numpy computed') + plt.xlabel('x') + plt.ylabel('z') + plt.legend() + + if plot_: + plt.subplot(4,1,2) + plt.plot(y_rec.vector, derivative.vector, 'r-', label='dz/dy0') + # derivatives computed by putting x values in the analytical formula + dzdy = function.c['c0'] * np.exp(function.c['c1'] * xarr) * (- np.sin(yarr)) + err = np.abs(dzdy - derivative.vector[1:]) + assert (err < 0.05).all(), (err[err>1e-2]) + assert (np.mean(err) < 1e-2), np.mean(err) + + if plot_: + plt.plot(yarr, dzdy, 'b--', label='numpy computed') + plt.xlabel('y') + plt.ylabel('dz/dy') + plt.legend() + + if plot_: + plt.subplot(4,1,3) + # *** BEWARE *** The first two entries are spurious. Entry 0 is + # *** from reinit sending out the defaults. Entry 2 is because + # *** there is no lastValue for computing real forward difference. + plt.plot(np.arange(2, len(rate.vector), 1) * dt, rate.vector[2:], 'r-', label='dz/dt') + dzdt = np.diff(z)/dt + plt.plot(np.arange(0, len(dzdt), 1.0) * dt, dzdt, 'b--', label='numpy computed') + plt.xlabel('t') + plt.ylabel('dz/dt') + plt.legend() + plt.tight_layout() + plt.show() + +if __name__ == '__main__': + example() diff --git a/tests/python/test_gsolve_parallel.py b/tests/python/test_gsolve_parallel.py index 13c4ad8a03..056762b47b 100644 --- a/tests/python/test_gsolve_parallel.py +++ b/tests/python/test_gsolve_parallel.py @@ -2,8 +2,14 @@ import moose print( '[INFO] Using moose from %s' % moose.__file__ ) import time + +# Does not guarantee thread determinism in multithreaded Gsolve/Ksolve. moose.seed( 10 ) +def printCompt(compt): + print( 'x0=%s, x1=%s, diffLength=%s, numDiffCompt=%d' % (compt.x0, compt.x1, + compt.diffLength, compt.numDiffCompts)) + def main( nT ): """ This example implements a reaction-diffusion like system which is @@ -14,25 +20,22 @@ def main( nT ): uses a reaction plus a function object to control its rates. """ print( 'Using %d threads' % nT ) - dt = 0.1 - # define the geometry compt = moose.CylMesh( '/cylinder' ) compt.r0 = compt.r1 = 100e-9 - compt.x1 = 200e-9 + compt.x1 = 200e-09 compt.diffLength = 0.2e-9 - assert( compt.numDiffCompts == compt.x1/compt.diffLength ) + assert( compt.numDiffCompts == compt.x1/compt.diffLength) #define the molecule. Its geometry is defined by its parent volume, cylinder c = moose.Pool( '/cylinder/pool' ) c.diffConst = 1e-13 # define diffusion constant - # Here we set up a function calculation func = moose.Function( '/cylinder/pool/func' ) - func.expr = "(-x0 * (30e-9 - x0) * (100e-9 - x0))*0.0001" - func.x.num = 1 #specify number of input variables. + func.expr = "(-x0*(30e-9-x0)*(100e-9-x0))*0.0001" + # func.x.num = 1 #specify number of input variables. #Connect the molecules to the func moose.connect( c, 'nOut', func.x[0], 'input' ) @@ -45,6 +48,7 @@ def main( nT ): ksolve.numThreads = nT except Exception as e: print( 'OLD MOOSE. Does not support multithreading' ) + dsolve = moose.Dsolve( '/cylinder/dsolve' ) stoich = moose.Stoich( '/cylinder/stoich' ) stoich.compartment = compt @@ -62,16 +66,32 @@ def main( nT ): runtime = updateDt * 10 t1 = time.time() res = [] + clk = moose.element( '/clock' ) for t in range( 0, runtime-1, updateDt ): - moose.start( updateDt ) y = c.vec.n - res.append( (np.mean(y), np.std(y)) ) + s = np.sum(y) + v = (np.mean(y), np.max(y), np.min(y), s) + print(v) + res.append(v) + moose.start( updateDt ) + currTime = clk.currentTime - expected = [(11.0, 0.0), (6.0, 0.0), (5.0, 0.0), (5.0, 0.0), (3.0, 0.0), - (3.0, 0.0), (3.0, 0.0), (3.0, 0.0), (3.0, 0.0), (3.0, 0.0)] + # One molecule here and there because thread launching has undeterministic + # characteristics. Even after setting moose.seed; we may not see same + # numbers on all platfroms. + expected = [ + (1000.0, 1000.0, 1000.0, 1000000.0) + , (9.908, 10.0, 8.0, 9908.0) + , (6.869, 7.0, 6.0, 6869.0) + , (5.354, 6.0, 5.0, 5354.0) + , (4.562, 5.0, 4.0, 4562.0) + , (3.483, 4.0, 3.0, 3483.0) + , (3.043, 4.0, 3.0, 3043.0) + , (2.261, 3.0, 2.0, 2261.0) + , (1.967, 2.0, 1.0, 1967.0) + , (1.997, 2.0, 1.0, 1997.0) ] print("Time = ", time.time() - t1) - print( res ) - assert np.isclose(res, expected, 0.1, 0.1).all(), (res, expected) + assert np.isclose(res, expected, atol=1, rtol=1).all(), "Got %s, expected %s" % (res, expected) if __name__ == '__main__': import sys diff --git a/tests/python/test_moose_attribs.py b/tests/python/test_moose_attribs.py index b10a2b30a4..3c8746b773 100644 --- a/tests/python/test_moose_attribs.py +++ b/tests/python/test_moose_attribs.py @@ -18,7 +18,7 @@ 'Compartment', 'CompartmentBase', 'ConcChan', 'CplxEnzBase', 'CubeMesh', 'CylMesh', 'DestField', 'DiagonalMsg', 'DifBuffer', 'DifBufferBase', 'DifShell', 'DifShellBase', 'DiffAmp', 'Dsolve', 'ElementField', - 'EndoMesh', 'Enz', 'EnzBase', 'ExIF', 'Finfo', 'Func', 'Function', + 'EndoMesh', 'Enz', 'EnzBase', 'ExIF', 'Finfo', 'Function', 'GapJunction', 'GraupnerBrunel2012CaPlasticitySynHandler', 'Group', 'Gsolve', 'HHChannel', 'HHChannel2D', 'HHChannelBase', 'HHGate', 'HHGate2D', 'HSolve', 'INFINITE', 'IntFire', 'IntFireBase', 'Interpol',