From 66ef7141d6cd59877caf8984e938523c1b96d7c6 Mon Sep 17 00:00:00 2001 From: Alan Garny Date: Mon, 4 Nov 2024 18:35:55 +1300 Subject: [PATCH] Generator: generated some wrong indices when untracking some variables. --- src/generator.cpp | 56 ++++++++++++++++--- src/generator_p.h | 8 ++- ...r.tracking.untracked.algebraic.variables.c | 48 ++++++++-------- ....tracking.untracked.algebraic.variables.py | 48 ++++++++-------- ...del.dae.for.tracking.untracked.variables.c | 48 ++++++++-------- ...el.dae.for.tracking.untracked.variables.py | 48 ++++++++-------- ...acked.algebraic.variables.with.externals.c | 6 +- ...cked.algebraic.variables.with.externals.py | 6 +- ...racked.computed.constants.with.externals.c | 6 +- ...acked.computed.constants.with.externals.py | 6 +- ...model.untracked.variables.with.externals.c | 10 ++-- ...odel.untracked.variables.with.externals.py | 10 ++-- 12 files changed, 173 insertions(+), 127 deletions(-) diff --git a/src/generator.cpp b/src/generator.cpp index e348d4b1e..434324b1f 100644 --- a/src/generator.cpp +++ b/src/generator.cpp @@ -44,6 +44,48 @@ void Generator::GeneratorImpl::reset() mCode = {}; } +std::string Generator::GeneratorImpl::doVariableIndexString(const AnalyserModelPtr &model, + const AnalyserVariablePtr &variable, + const std::vector &variables) +{ + // Determine the actual index of the variable in the list of variables by accounting for the fact that some + // variables may be untracked. + + size_t i = MAX_SIZE_T; + size_t res = MAX_SIZE_T; + + for (;;) { + auto var = variables[++i]; + + if (doIsTrackedVariable(model, var)) { + ++res; + } + + if (variable == var) { + break; + } + } + + return convertToString(res); +} + +std::string Generator::GeneratorImpl::variableIndexString(const AnalyserModelPtr &model, + const AnalyserVariablePtr &variable) +{ + switch (variable->type()) { + case AnalyserVariable::Type::CONSTANT: + return doVariableIndexString(model, variable, model->constants()); + case AnalyserVariable::Type::COMPUTED_CONSTANT: + return doVariableIndexString(model, variable, model->computedConstants()); + case AnalyserVariable::Type::ALGEBRAIC: + return doVariableIndexString(model, variable, model->algebraic()); + default: + break; + } + + return convertToString(variable->index()); +} + bool Generator::GeneratorImpl::doIsTrackedEquation(const AnalyserEquationPtr &equation, bool tracked) { switch (equation->type()) { @@ -1177,7 +1219,7 @@ void Generator::GeneratorImpl::addNlaSystemsCode(const AnalyserModelPtr &model) mProfile->algebraicArrayString(); methodBody += mProfile->indentString() - + arrayString + mProfile->openArrayString() + convertToString(variable->index()) + mProfile->closeArrayString() + + arrayString + mProfile->openArrayString() + variableIndexString(model, variable) + mProfile->closeArrayString() + mProfile->equalityString() + mProfile->uArrayString() + mProfile->openArrayString() + convertToString(++i) + mProfile->closeArrayString() + mProfile->commandSeparatorString() + "\n"; @@ -1255,7 +1297,7 @@ void Generator::GeneratorImpl::addNlaSystemsCode(const AnalyserModelPtr &model) methodBody += mProfile->indentString() + mProfile->uArrayString() + mProfile->openArrayString() + convertToString(++i) + mProfile->closeArrayString() + mProfile->equalityString() - + arrayString + mProfile->openArrayString() + convertToString(variable->index()) + mProfile->closeArrayString() + + arrayString + mProfile->openArrayString() + variableIndexString(model, variable) + mProfile->closeArrayString() + mProfile->commandSeparatorString() + "\n"; } @@ -1281,7 +1323,7 @@ void Generator::GeneratorImpl::addNlaSystemsCode(const AnalyserModelPtr &model) mProfile->algebraicArrayString(); methodBody += mProfile->indentString() - + arrayString + mProfile->openArrayString() + convertToString(variable->index()) + mProfile->closeArrayString() + + arrayString + mProfile->openArrayString() + variableIndexString(model, variable) + mProfile->closeArrayString() + mProfile->equalityString() + mProfile->uArrayString() + mProfile->openArrayString() + convertToString(++i) + mProfile->closeArrayString() + mProfile->commandSeparatorString() + "\n"; @@ -1322,7 +1364,7 @@ std::string generateDoubleCode(const std::string &value) } std::string Generator::GeneratorImpl::generateDoubleOrConstantVariableNameCode(const AnalyserModelPtr &model, - const VariablePtr &variable) const + const VariablePtr &variable) { if (isCellMLReal(variable->initialValue())) { return generateDoubleCode(variable->initialValue()); @@ -1331,7 +1373,7 @@ std::string Generator::GeneratorImpl::generateDoubleOrConstantVariableNameCode(c auto initialValueVariable = owningComponent(variable)->variable(variable->initialValue()); auto analyserInitialValueVariable = model->variable(initialValueVariable); - return mProfile->constantsArrayString() + mProfile->openArrayString() + convertToString(analyserInitialValueVariable->index()) + mProfile->closeArrayString(); + return mProfile->constantsArrayString() + mProfile->openArrayString() + variableIndexString(model, analyserInitialValueVariable) + mProfile->closeArrayString(); } std::string Generator::GeneratorImpl::generateVariableNameCode(const AnalyserModelPtr &model, @@ -1373,7 +1415,7 @@ std::string Generator::GeneratorImpl::generateVariableNameCode(const AnalyserMod arrayName = mProfile->externalArrayString(); } - return arrayName + mProfile->openArrayString() + convertToString(analyserVariable->index()) + mProfile->closeArrayString(); + return arrayName + mProfile->openArrayString() + variableIndexString(model, analyserVariable) + mProfile->closeArrayString(); } std::string Generator::GeneratorImpl::generateOperatorCode(const AnalyserModelPtr &model, const std::string &op, @@ -2252,7 +2294,7 @@ std::string Generator::GeneratorImpl::generateEquationCode(const AnalyserModelPt + generateVariableNameCode(model, variable->variable()) + mProfile->equalityString() + replace(mProfile->externalVariableMethodCallString(modelHasOdes(model)), - "[INDEX]", convertToString(variable->index())) + "[INDEX]", variableIndexString(model, variable)) + mProfile->commandSeparatorString() + "\n"; } diff --git a/src/generator_p.h b/src/generator_p.h index a88582733..3489d4786 100644 --- a/src/generator_p.h +++ b/src/generator_p.h @@ -43,12 +43,16 @@ struct Generator::GeneratorImpl: public Logger::LoggerImpl void reset(); + std::string doVariableIndexString(const AnalyserModelPtr &model, const AnalyserVariablePtr &variable, + const std::vector &variables); + std::string variableIndexString(const AnalyserModelPtr &model, const AnalyserVariablePtr &variable); + bool doIsTrackedEquation(const AnalyserEquationPtr &equation, bool tracked); bool isTrackedEquation(const AnalyserEquationPtr &equation); bool isUntrackedEquation(const AnalyserEquationPtr &equation); - bool doIsTrackedVariable(const AnalyserModelPtr &model, const AnalyserVariablePtr &variable, bool tracked); + bool doIsTrackedVariable(const AnalyserModelPtr &model, const AnalyserVariablePtr &variable, bool tracked = true); bool doIsTrackedVariable(const AnalyserVariablePtr &variable, bool tracked); bool isTrackedVariable(const AnalyserVariablePtr &variable); @@ -161,7 +165,7 @@ struct Generator::GeneratorImpl: public Logger::LoggerImpl std::string generateMethodBodyCode(const std::string &methodBody) const; std::string generateDoubleOrConstantVariableNameCode(const AnalyserModelPtr &model, - const VariablePtr &variable) const; + const VariablePtr &variable); std::string generateVariableNameCode(const AnalyserModelPtr &model, const VariablePtr &variable, bool state = true); diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.c b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.c index 841dc42ca..a6f115c0a 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.c +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.c @@ -274,9 +274,9 @@ void objectiveFunction6(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[7] = u[0]; + algebraic[6] = u[0]; - f[0] = algebraic[7]-4.0*exp(states[0]/18.0)-0.0; + f[0] = algebraic[6]-4.0*exp(states[0]/18.0)-0.0; } void findRoot6(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -284,11 +284,11 @@ void findRoot6(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[7]; + u[0] = algebraic[6]; nlaSolve(objectiveFunction6, u, 1, &rfi); - algebraic[7] = u[0]; + algebraic[6] = u[0]; } void objectiveFunction7(double *u, double *f, void *data) @@ -304,7 +304,7 @@ void objectiveFunction7(double *u, double *f, void *data) double sodium_channel_m_gate_alpha_m = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); - f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[7]*states[2])-0.0; + f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[6]*states[2])-0.0; } void findRoot7(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -328,9 +328,9 @@ void objectiveFunction8(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[8] = u[0]; + algebraic[7] = u[0]; - f[0] = algebraic[8]-0.07*exp(states[0]/20.0)-0.0; + f[0] = algebraic[7]-0.07*exp(states[0]/20.0)-0.0; } void findRoot8(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -338,11 +338,11 @@ void findRoot8(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[8]; + u[0] = algebraic[7]; nlaSolve(objectiveFunction8, u, 1, &rfi); - algebraic[8] = u[0]; + algebraic[7] = u[0]; } void objectiveFunction9(double *u, double *f, void *data) @@ -354,9 +354,9 @@ void objectiveFunction9(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[9] = u[0]; + algebraic[8] = u[0]; - f[0] = algebraic[9]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0; + f[0] = algebraic[8]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0; } void findRoot9(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -364,11 +364,11 @@ void findRoot9(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[9]; + u[0] = algebraic[8]; nlaSolve(objectiveFunction9, u, 1, &rfi); - algebraic[9] = u[0]; + algebraic[8] = u[0]; } void objectiveFunction10(double *u, double *f, void *data) @@ -382,7 +382,7 @@ void objectiveFunction10(double *u, double *f, void *data) rates[1] = u[0]; - f[0] = rates[1]-(algebraic[8]*(1.0-states[1])-algebraic[9]*states[1])-0.0; + f[0] = rates[1]-(algebraic[7]*(1.0-states[1])-algebraic[8]*states[1])-0.0; } void findRoot10(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -432,9 +432,9 @@ void objectiveFunction12(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[10] = u[0]; + algebraic[9] = u[0]; - f[0] = algebraic[10]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0; + f[0] = algebraic[9]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0; } void findRoot12(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -442,11 +442,11 @@ void findRoot12(double voi, double *states, double *rates, double *constants, do RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[10]; + u[0] = algebraic[9]; nlaSolve(objectiveFunction12, u, 1, &rfi); - algebraic[10] = u[0]; + algebraic[9] = u[0]; } void objectiveFunction13(double *u, double *f, void *data) @@ -458,9 +458,9 @@ void objectiveFunction13(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[11] = u[0]; + algebraic[10] = u[0]; - f[0] = algebraic[11]-0.125*exp(states[0]/80.0)-0.0; + f[0] = algebraic[10]-0.125*exp(states[0]/80.0)-0.0; } void findRoot13(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -468,11 +468,11 @@ void findRoot13(double voi, double *states, double *rates, double *constants, do RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[11]; + u[0] = algebraic[10]; nlaSolve(objectiveFunction13, u, 1, &rfi); - algebraic[11] = u[0]; + algebraic[10] = u[0]; } void objectiveFunction14(double *u, double *f, void *data) @@ -486,7 +486,7 @@ void objectiveFunction14(double *u, double *f, void *data) rates[3] = u[0]; - f[0] = rates[3]-(algebraic[10]*(1.0-states[3])-algebraic[11]*states[3])-0.0; + f[0] = rates[3]-(algebraic[9]*(1.0-states[3])-algebraic[10]*states[3])-0.0; } void findRoot14(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -522,11 +522,11 @@ void initialiseVariables(double *states, double *rates, double *constants, doubl algebraic[3] = 0.0; algebraic[4] = 0.0; algebraic[5] = 0.0; + algebraic[6] = 0.0; algebraic[7] = 0.0; algebraic[8] = 0.0; algebraic[9] = 0.0; algebraic[10] = 0.0; - algebraic[11] = 0.0; } void computeComputedConstants(double *constants, double *computedConstants) diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.py b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.py index a6d6678a0..bf19a353d 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.py +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.algebraic.variables.py @@ -225,19 +225,19 @@ def objective_function_6(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[7] = u[0] + algebraic[6] = u[0] - f[0] = algebraic[7]-4.0*exp(states[0]/18.0)-0.0 + f[0] = algebraic[6]-4.0*exp(states[0]/18.0)-0.0 def find_root_6(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[7] + u[0] = algebraic[6] u = nla_solve(objective_function_6, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[7] = u[0] + algebraic[6] = u[0] def objective_function_7(u, f, data): @@ -252,7 +252,7 @@ def objective_function_7(u, f, data): sodium_channel_m_gate_alpha_m = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) - f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[7]*states[2])-0.0 + f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[6]*states[2])-0.0 def find_root_7(voi, states, rates, constants, computed_constants, algebraic): @@ -273,19 +273,19 @@ def objective_function_8(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[8] = u[0] + algebraic[7] = u[0] - f[0] = algebraic[8]-0.07*exp(states[0]/20.0)-0.0 + f[0] = algebraic[7]-0.07*exp(states[0]/20.0)-0.0 def find_root_8(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[8] + u[0] = algebraic[7] u = nla_solve(objective_function_8, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[8] = u[0] + algebraic[7] = u[0] def objective_function_9(u, f, data): @@ -296,19 +296,19 @@ def objective_function_9(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[9] = u[0] + algebraic[8] = u[0] - f[0] = algebraic[9]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0 + f[0] = algebraic[8]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0 def find_root_9(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[9] + u[0] = algebraic[8] u = nla_solve(objective_function_9, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[9] = u[0] + algebraic[8] = u[0] def objective_function_10(u, f, data): @@ -321,7 +321,7 @@ def objective_function_10(u, f, data): rates[1] = u[0] - f[0] = rates[1]-(algebraic[8]*(1.0-states[1])-algebraic[9]*states[1])-0.0 + f[0] = rates[1]-(algebraic[7]*(1.0-states[1])-algebraic[8]*states[1])-0.0 def find_root_10(voi, states, rates, constants, computed_constants, algebraic): @@ -365,19 +365,19 @@ def objective_function_12(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[10] = u[0] + algebraic[9] = u[0] - f[0] = algebraic[10]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0 + f[0] = algebraic[9]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0 def find_root_12(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[10] + u[0] = algebraic[9] u = nla_solve(objective_function_12, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[10] = u[0] + algebraic[9] = u[0] def objective_function_13(u, f, data): @@ -388,19 +388,19 @@ def objective_function_13(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[11] = u[0] + algebraic[10] = u[0] - f[0] = algebraic[11]-0.125*exp(states[0]/80.0)-0.0 + f[0] = algebraic[10]-0.125*exp(states[0]/80.0)-0.0 def find_root_13(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[11] + u[0] = algebraic[10] u = nla_solve(objective_function_13, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[11] = u[0] + algebraic[10] = u[0] def objective_function_14(u, f, data): @@ -413,7 +413,7 @@ def objective_function_14(u, f, data): rates[3] = u[0] - f[0] = rates[3]-(algebraic[10]*(1.0-states[3])-algebraic[11]*states[3])-0.0 + f[0] = rates[3]-(algebraic[9]*(1.0-states[3])-algebraic[10]*states[3])-0.0 def find_root_14(voi, states, rates, constants, computed_constants, algebraic): @@ -446,11 +446,11 @@ def initialise_variables(states, rates, constants, computed_constants, algebraic algebraic[3] = 0.0 algebraic[4] = 0.0 algebraic[5] = 0.0 + algebraic[6] = 0.0 algebraic[7] = 0.0 algebraic[8] = 0.0 algebraic[9] = 0.0 algebraic[10] = 0.0 - algebraic[11] = 0.0 def compute_computed_constants(constants, computed_constants): diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.c b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.c index f7e6db7bf..3d3a0767b 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.c +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.c @@ -278,9 +278,9 @@ void objectiveFunction6(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[7] = u[0]; + algebraic[6] = u[0]; - f[0] = algebraic[7]-4.0*exp(states[0]/18.0)-0.0; + f[0] = algebraic[6]-4.0*exp(states[0]/18.0)-0.0; } void findRoot6(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -288,11 +288,11 @@ void findRoot6(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[7]; + u[0] = algebraic[6]; nlaSolve(objectiveFunction6, u, 1, &rfi); - algebraic[7] = u[0]; + algebraic[6] = u[0]; } void objectiveFunction7(double *u, double *f, void *data) @@ -308,7 +308,7 @@ void objectiveFunction7(double *u, double *f, void *data) double sodium_channel_m_gate_alpha_m = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); - f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[7]*states[2])-0.0; + f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[6]*states[2])-0.0; } void findRoot7(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -332,9 +332,9 @@ void objectiveFunction8(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[8] = u[0]; + algebraic[7] = u[0]; - f[0] = algebraic[8]-0.07*exp(states[0]/20.0)-0.0; + f[0] = algebraic[7]-0.07*exp(states[0]/20.0)-0.0; } void findRoot8(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -342,11 +342,11 @@ void findRoot8(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[8]; + u[0] = algebraic[7]; nlaSolve(objectiveFunction8, u, 1, &rfi); - algebraic[8] = u[0]; + algebraic[7] = u[0]; } void objectiveFunction9(double *u, double *f, void *data) @@ -358,9 +358,9 @@ void objectiveFunction9(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[9] = u[0]; + algebraic[8] = u[0]; - f[0] = algebraic[9]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0; + f[0] = algebraic[8]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0; } void findRoot9(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -368,11 +368,11 @@ void findRoot9(double voi, double *states, double *rates, double *constants, dou RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[9]; + u[0] = algebraic[8]; nlaSolve(objectiveFunction9, u, 1, &rfi); - algebraic[9] = u[0]; + algebraic[8] = u[0]; } void objectiveFunction10(double *u, double *f, void *data) @@ -386,7 +386,7 @@ void objectiveFunction10(double *u, double *f, void *data) rates[1] = u[0]; - f[0] = rates[1]-(algebraic[8]*(1.0-states[1])-algebraic[9]*states[1])-0.0; + f[0] = rates[1]-(algebraic[7]*(1.0-states[1])-algebraic[8]*states[1])-0.0; } void findRoot10(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -440,9 +440,9 @@ void objectiveFunction12(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[10] = u[0]; + algebraic[9] = u[0]; - f[0] = algebraic[10]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0; + f[0] = algebraic[9]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0; } void findRoot12(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -450,11 +450,11 @@ void findRoot12(double voi, double *states, double *rates, double *constants, do RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[10]; + u[0] = algebraic[9]; nlaSolve(objectiveFunction12, u, 1, &rfi); - algebraic[10] = u[0]; + algebraic[9] = u[0]; } void objectiveFunction13(double *u, double *f, void *data) @@ -466,9 +466,9 @@ void objectiveFunction13(double *u, double *f, void *data) double *computedConstants = ((RootFindingInfo *) data)->computedConstants; double *algebraic = ((RootFindingInfo *) data)->algebraic; - algebraic[11] = u[0]; + algebraic[10] = u[0]; - f[0] = algebraic[11]-0.125*exp(states[0]/80.0)-0.0; + f[0] = algebraic[10]-0.125*exp(states[0]/80.0)-0.0; } void findRoot13(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -476,11 +476,11 @@ void findRoot13(double voi, double *states, double *rates, double *constants, do RootFindingInfo rfi = { voi, states, rates, constants, computedConstants, algebraic }; double u[1]; - u[0] = algebraic[11]; + u[0] = algebraic[10]; nlaSolve(objectiveFunction13, u, 1, &rfi); - algebraic[11] = u[0]; + algebraic[10] = u[0]; } void objectiveFunction14(double *u, double *f, void *data) @@ -494,7 +494,7 @@ void objectiveFunction14(double *u, double *f, void *data) rates[3] = u[0]; - f[0] = rates[3]-(algebraic[10]*(1.0-states[3])-algebraic[11]*states[3])-0.0; + f[0] = rates[3]-(algebraic[9]*(1.0-states[3])-algebraic[10]*states[3])-0.0; } void findRoot14(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic) @@ -525,11 +525,11 @@ void initialiseVariables(double *states, double *rates, double *constants, doubl algebraic[3] = 0.0; algebraic[4] = 0.0; algebraic[5] = 0.0; + algebraic[6] = 0.0; algebraic[7] = 0.0; algebraic[8] = 0.0; algebraic[9] = 0.0; algebraic[10] = 0.0; - algebraic[11] = 0.0; } void computeComputedConstants(double *constants, double *computedConstants) diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.py b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.py index 4bf475898..ae9f7dffb 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.py +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.dae.for.tracking.untracked.variables.py @@ -229,19 +229,19 @@ def objective_function_6(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[7] = u[0] + algebraic[6] = u[0] - f[0] = algebraic[7]-4.0*exp(states[0]/18.0)-0.0 + f[0] = algebraic[6]-4.0*exp(states[0]/18.0)-0.0 def find_root_6(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[7] + u[0] = algebraic[6] u = nla_solve(objective_function_6, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[7] = u[0] + algebraic[6] = u[0] def objective_function_7(u, f, data): @@ -256,7 +256,7 @@ def objective_function_7(u, f, data): sodium_channel_m_gate_alpha_m = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) - f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[7]*states[2])-0.0 + f[0] = rates[2]-(sodium_channel_m_gate_alpha_m*(1.0-states[2])-algebraic[6]*states[2])-0.0 def find_root_7(voi, states, rates, constants, computed_constants, algebraic): @@ -277,19 +277,19 @@ def objective_function_8(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[8] = u[0] + algebraic[7] = u[0] - f[0] = algebraic[8]-0.07*exp(states[0]/20.0)-0.0 + f[0] = algebraic[7]-0.07*exp(states[0]/20.0)-0.0 def find_root_8(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[8] + u[0] = algebraic[7] u = nla_solve(objective_function_8, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[8] = u[0] + algebraic[7] = u[0] def objective_function_9(u, f, data): @@ -300,19 +300,19 @@ def objective_function_9(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[9] = u[0] + algebraic[8] = u[0] - f[0] = algebraic[9]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0 + f[0] = algebraic[8]-1.0/(exp((states[0]+30.0)/10.0)+1.0)-0.0 def find_root_9(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[9] + u[0] = algebraic[8] u = nla_solve(objective_function_9, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[9] = u[0] + algebraic[8] = u[0] def objective_function_10(u, f, data): @@ -325,7 +325,7 @@ def objective_function_10(u, f, data): rates[1] = u[0] - f[0] = rates[1]-(algebraic[8]*(1.0-states[1])-algebraic[9]*states[1])-0.0 + f[0] = rates[1]-(algebraic[7]*(1.0-states[1])-algebraic[8]*states[1])-0.0 def find_root_10(voi, states, rates, constants, computed_constants, algebraic): @@ -373,19 +373,19 @@ def objective_function_12(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[10] = u[0] + algebraic[9] = u[0] - f[0] = algebraic[10]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0 + f[0] = algebraic[9]-0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0)-0.0 def find_root_12(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[10] + u[0] = algebraic[9] u = nla_solve(objective_function_12, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[10] = u[0] + algebraic[9] = u[0] def objective_function_13(u, f, data): @@ -396,19 +396,19 @@ def objective_function_13(u, f, data): computed_constants = data[4] algebraic = data[5] - algebraic[11] = u[0] + algebraic[10] = u[0] - f[0] = algebraic[11]-0.125*exp(states[0]/80.0)-0.0 + f[0] = algebraic[10]-0.125*exp(states[0]/80.0)-0.0 def find_root_13(voi, states, rates, constants, computed_constants, algebraic): u = [nan]*1 - u[0] = algebraic[11] + u[0] = algebraic[10] u = nla_solve(objective_function_13, u, 1, [voi, states, rates, constants, computed_constants, algebraic]) - algebraic[11] = u[0] + algebraic[10] = u[0] def objective_function_14(u, f, data): @@ -421,7 +421,7 @@ def objective_function_14(u, f, data): rates[3] = u[0] - f[0] = rates[3]-(algebraic[10]*(1.0-states[3])-algebraic[11]*states[3])-0.0 + f[0] = rates[3]-(algebraic[9]*(1.0-states[3])-algebraic[10]*states[3])-0.0 def find_root_14(voi, states, rates, constants, computed_constants, algebraic): @@ -449,11 +449,11 @@ def initialise_variables(states, rates, constants, computed_constants, algebraic algebraic[3] = 0.0 algebraic[4] = 0.0 algebraic[5] = 0.0 + algebraic[6] = 0.0 algebraic[7] = 0.0 algebraic[8] = 0.0 algebraic[9] = 0.0 algebraic[10] = 0.0 - algebraic[11] = 0.0 def compute_computed_constants(constants, computed_constants): diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.c b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.c index f4e5bac15..06063c5c2 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.c +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.c @@ -130,11 +130,11 @@ void computeRates(double voi, double *states, double *rates, double *constants, double membrane_i_Stim = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0; double leakage_current_i_L = constants[2]*(states[0]-computedConstants[0]); double potassium_channel_i_K = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[2]); - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); externals[0] = externalVariable(voi, states, rates, constants, computedConstants, algebraic, externals, 0); rates[0] = -(-membrane_i_Stim+externals[0]+potassium_channel_i_K+leakage_current_i_L)/constants[0]; double sodium_channel_m_gate_beta_m = 4.0*exp(states[0]/18.0); - rates[2] = algebraic[3]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2]; + rates[2] = algebraic[0]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2]; double sodium_channel_h_gate_alpha_h = 0.07*exp(states[0]/20.0); double sodium_channel_h_gate_beta_h = 1.0/(exp((states[0]+30.0)/10.0)+1.0); rates[1] = sodium_channel_h_gate_alpha_h*(1.0-states[1])-sodium_channel_h_gate_beta_h*states[1]; @@ -145,6 +145,6 @@ void computeRates(double voi, double *states, double *rates, double *constants, void computeVariables(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic, double *externals, ExternalVariable externalVariable) { - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); externals[0] = externalVariable(voi, states, rates, constants, computedConstants, algebraic, externals, 0); } diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.py b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.py index a3ced765d..c9dd42140 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.py +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.algebraic.variables.with.externals.py @@ -99,11 +99,11 @@ def compute_rates(voi, states, rates, constants, computed_constants, algebraic, membrane_i_Stim = -20.0 if and_func(geq_func(voi, 10.0), leq_func(voi, 10.5)) else 0.0 leakage_current_i_L = constants[2]*(states[0]-computed_constants[0]) potassium_channel_i_K = constants[4]*pow(states[3], 4.0)*(states[0]-computed_constants[2]) - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) externals[0] = external_variable(voi, states, rates, constants, computed_constants, algebraic, externals, 0) rates[0] = -(-membrane_i_Stim+externals[0]+potassium_channel_i_K+leakage_current_i_L)/constants[0] sodium_channel_m_gate_beta_m = 4.0*exp(states[0]/18.0) - rates[2] = algebraic[3]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2] + rates[2] = algebraic[0]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2] sodium_channel_h_gate_alpha_h = 0.07*exp(states[0]/20.0) sodium_channel_h_gate_beta_h = 1.0/(exp((states[0]+30.0)/10.0)+1.0) rates[1] = sodium_channel_h_gate_alpha_h*(1.0-states[1])-sodium_channel_h_gate_beta_h*states[1] @@ -113,5 +113,5 @@ def compute_rates(voi, states, rates, constants, computed_constants, algebraic, def compute_variables(voi, states, rates, constants, computed_constants, algebraic, externals, external_variable): - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) externals[0] = external_variable(voi, states, rates, constants, computed_constants, algebraic, externals, 0) diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.c b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.c index 126baf57e..d064cb298 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.c +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.c @@ -126,7 +126,7 @@ void initialiseVariables(double *states, double *rates, double *constants, doubl void computeComputedConstants(double *constants, double *computedConstants) { - computedConstants[2] = constants[1]+12.0; + computedConstants[0] = constants[1]+12.0; } void computeRates(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic, double *externals, ExternalVariable externalVariable) @@ -134,7 +134,7 @@ void computeRates(double voi, double *states, double *rates, double *constants, algebraic[0] = ((voi >= 10.0) && (voi <= 10.5))?-20.0:0.0; double leakage_current_E_L = constants[1]-10.613; algebraic[1] = constants[2]*(states[0]-leakage_current_E_L); - algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[2]); + algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[0]); algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); externals[0] = externalVariable(voi, states, rates, constants, computedConstants, algebraic, externals, 0); rates[0] = -(-algebraic[0]+externals[0]+algebraic[2]+algebraic[1])/constants[0]; @@ -157,7 +157,7 @@ void computeVariables(double voi, double *states, double *rates, double *constan algebraic[4] = 4.0*exp(states[0]/18.0); algebraic[5] = 0.07*exp(states[0]/20.0); algebraic[6] = 1.0/(exp((states[0]+30.0)/10.0)+1.0); - algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[2]); + algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computedConstants[0]); algebraic[7] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0); algebraic[8] = 0.125*exp(states[0]/80.0); } diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.py b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.py index d98d9a538..a459f81d8 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.py +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.computed.constants.with.externals.py @@ -96,14 +96,14 @@ def initialise_variables(states, rates, constants, computed_constants, algebraic def compute_computed_constants(constants, computed_constants): - computed_constants[2] = constants[1]+12.0 + computed_constants[0] = constants[1]+12.0 def compute_rates(voi, states, rates, constants, computed_constants, algebraic, externals, external_variable): algebraic[0] = -20.0 if and_func(geq_func(voi, 10.0), leq_func(voi, 10.5)) else 0.0 leakage_current_E_L = constants[1]-10.613 algebraic[1] = constants[2]*(states[0]-leakage_current_E_L) - algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computed_constants[2]) + algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computed_constants[0]) algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) externals[0] = external_variable(voi, states, rates, constants, computed_constants, algebraic, externals, 0) rates[0] = -(-algebraic[0]+externals[0]+algebraic[2]+algebraic[1])/constants[0] @@ -125,6 +125,6 @@ def compute_variables(voi, states, rates, constants, computed_constants, algebra algebraic[4] = 4.0*exp(states[0]/18.0) algebraic[5] = 0.07*exp(states[0]/20.0) algebraic[6] = 1.0/(exp((states[0]+30.0)/10.0)+1.0) - algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computed_constants[2]) + algebraic[2] = constants[4]*pow(states[3], 4.0)*(states[0]-computed_constants[0]) algebraic[7] = 0.01*(states[0]+10.0)/(exp((states[0]+10.0)/10.0)-1.0) algebraic[8] = 0.125*exp(states[0]/80.0) diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.c b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.c index 5eb5595c9..661b551a8 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.c +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.c @@ -111,7 +111,7 @@ void initialiseVariables(double *states, double *rates, double *constants, doubl void computeComputedConstants(double *constants, double *computedConstants) { double membrane_E_R = 0.0; - computedConstants[2] = membrane_E_R+12.0; + computedConstants[0] = membrane_E_R+12.0; } void computeRates(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic, double *externals, ExternalVariable externalVariable) @@ -122,12 +122,12 @@ void computeRates(double voi, double *states, double *rates, double *constants, double leakage_current_E_L = membrane_E_R-10.613; double leakage_current_i_L = leakage_current_g_L*(states[0]-leakage_current_E_L); double potassium_channel_g_K = 36.0; - double potassium_channel_i_K = potassium_channel_g_K*pow(states[3], 4.0)*(states[0]-computedConstants[2]); - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); + double potassium_channel_i_K = potassium_channel_g_K*pow(states[3], 4.0)*(states[0]-computedConstants[0]); + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); externals[0] = externalVariable(voi, states, rates, constants, computedConstants, algebraic, externals, 0); rates[0] = -(-membrane_i_Stim+externals[0]+potassium_channel_i_K+leakage_current_i_L)/constants[0]; double sodium_channel_m_gate_beta_m = 4.0*exp(states[0]/18.0); - rates[2] = algebraic[3]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2]; + rates[2] = algebraic[0]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2]; double sodium_channel_h_gate_alpha_h = 0.07*exp(states[0]/20.0); double sodium_channel_h_gate_beta_h = 1.0/(exp((states[0]+30.0)/10.0)+1.0); rates[1] = sodium_channel_h_gate_alpha_h*(1.0-states[1])-sodium_channel_h_gate_beta_h*states[1]; @@ -138,6 +138,6 @@ void computeRates(double voi, double *states, double *rates, double *constants, void computeVariables(double voi, double *states, double *rates, double *constants, double *computedConstants, double *algebraic, double *externals, ExternalVariable externalVariable) { - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0); externals[0] = externalVariable(voi, states, rates, constants, computedConstants, algebraic, externals, 0); } diff --git a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.py b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.py index a8bc614fe..4bf0b335c 100644 --- a/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.py +++ b/tests/resources/generator/hodgkin_huxley_squid_axon_model_1952/model.untracked.variables.with.externals.py @@ -81,7 +81,7 @@ def initialise_variables(states, rates, constants, computed_constants, algebraic def compute_computed_constants(constants, computed_constants): membrane_E_R = 0.0 - computed_constants[2] = membrane_E_R+12.0 + computed_constants[0] = membrane_E_R+12.0 def compute_rates(voi, states, rates, constants, computed_constants, algebraic, externals, external_variable): @@ -91,12 +91,12 @@ def compute_rates(voi, states, rates, constants, computed_constants, algebraic, leakage_current_E_L = membrane_E_R-10.613 leakage_current_i_L = leakage_current_g_L*(states[0]-leakage_current_E_L) potassium_channel_g_K = 36.0 - potassium_channel_i_K = potassium_channel_g_K*pow(states[3], 4.0)*(states[0]-computed_constants[2]) - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) + potassium_channel_i_K = potassium_channel_g_K*pow(states[3], 4.0)*(states[0]-computed_constants[0]) + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) externals[0] = external_variable(voi, states, rates, constants, computed_constants, algebraic, externals, 0) rates[0] = -(-membrane_i_Stim+externals[0]+potassium_channel_i_K+leakage_current_i_L)/constants[0] sodium_channel_m_gate_beta_m = 4.0*exp(states[0]/18.0) - rates[2] = algebraic[3]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2] + rates[2] = algebraic[0]*(1.0-states[2])-sodium_channel_m_gate_beta_m*states[2] sodium_channel_h_gate_alpha_h = 0.07*exp(states[0]/20.0) sodium_channel_h_gate_beta_h = 1.0/(exp((states[0]+30.0)/10.0)+1.0) rates[1] = sodium_channel_h_gate_alpha_h*(1.0-states[1])-sodium_channel_h_gate_beta_h*states[1] @@ -106,5 +106,5 @@ def compute_rates(voi, states, rates, constants, computed_constants, algebraic, def compute_variables(voi, states, rates, constants, computed_constants, algebraic, externals, external_variable): - algebraic[3] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) + algebraic[0] = 0.1*(states[0]+25.0)/(exp((states[0]+25.0)/10.0)-1.0) externals[0] = external_variable(voi, states, rates, constants, computed_constants, algebraic, externals, 0)