diff --git a/modcc/module.cpp b/modcc/module.cpp
index 0e5c4de5eda7f42df867146f96dee8e86fb6bacb..df34279f5953ab90f6ef6c526563c7e984d1e2c5 100644
--- a/modcc/module.cpp
+++ b/modcc/module.cpp
@@ -375,7 +375,21 @@ bool Module::semantic() {
auto deriv = solve_expression->procedure();
if (deriv->kind()==procedureKind::kinetic) {
- kinetic_rewrite(deriv->body())->accept(solver.get());
+ auto rewrite_body = kinetic_rewrite(deriv->body());
+ bool linear_kinetic = true;
+
+ for (auto& s: rewrite_body->is_block()->statements()) {
+ if(s->is_assignment() && !state_vars.empty()) {
+ linear_test_result r = linear_test(s->is_assignment()->rhs(), state_vars);
+ linear_kinetic &= r.is_linear;
+ }
+ }
+
+ if (!linear_kinetic) {
+ solver = std::make_unique<SparseNonlinearSolverVisitor>();
+ }
+
+ rewrite_body->accept(solver.get());
}
else if (deriv->kind()==procedureKind::linear) {
solver = std::make_unique<LinearSolverVisitor>(state_vars);
diff --git a/modcc/solvers.cpp b/modcc/solvers.cpp
index 4a53c0d9eb4081c82f4ca0644f0f283840794853..19e3700ae872df4eea6bd53a2bcbd715df869150 100644
--- a/modcc/solvers.cpp
+++ b/modcc/solvers.cpp
@@ -206,7 +206,8 @@ void SparseSolverVisitor::visit(CompartmentExpression *e) {
return;
}
auto idx = it - dvars_.begin();
- scale_factor_[idx] = e->scale_factor()->clone();
+ scale_factor_[idx] = make_expression<DivBinaryExpression>(
+ loc, make_expression<NumberExpression>(loc, 1.0), e->scale_factor()->clone());
}
}
@@ -272,7 +273,7 @@ void SparseSolverVisitor::visit(AssignmentExpression *e) {
make_expression<MulBinaryExpression>(loc, r.coef[dvars_[j]]->clone(), dt_expr->clone());
if (scale_factor_[j]) {
- expr = make_expression<DivBinaryExpression>(loc, std::move(expr), scale_factor_[j]->clone());
+ expr = make_expression<MulBinaryExpression>(loc, std::move(expr), scale_factor_[j]->clone());
}
}
@@ -328,8 +329,8 @@ void SparseSolverVisitor::visit(ConserveExpression *e) {
}
}
else {
- error({"ICE: coefficient in state variable is not an identifier", loc});
- return;
+ error({"ICE: coefficient in state variable is not an identifier", loc});
+ return;
}
// Replace that row with the conserve statement
@@ -345,7 +346,7 @@ void SparseSolverVisitor::visit(ConserveExpression *e) {
if (it != terms.end()) {
auto expr = (*it)->is_stoich_term()->coeff()->clone();
if (scale_factor_[j]) {
- expr = make_expression<MulBinaryExpression>(loc, scale_factor_[j]->clone(), std::move(expr));
+ expr = make_expression<DivBinaryExpression>(loc, std::move(expr), scale_factor_[j]->clone());
}
auto local_a_term = make_unique_local_assign(scope, expr.get(), "a_");
@@ -477,6 +478,7 @@ void LinearSolverVisitor::visit(LinearExpression *e) {
rhs_.push_back(symtbl_.define(e->rhs()->is_identifier()->spelling()));
++deq_index_;
}
+
void LinearSolverVisitor::finalize() {
A_.augment(rhs_);
@@ -528,6 +530,273 @@ void LinearSolverVisitor::finalize() {
BlockRewriterBase::finalize();
}
+void SparseNonlinearSolverVisitor::visit(BlockExpression* e) {
+ // Do a first pass to initialize some local variables and extract state variables
+
+ for (auto& stmt: e->statements()) {
+ if (stmt && stmt->is_assignment() && stmt->is_assignment()->lhs()->is_derivative()) {
+ auto id = stmt->is_assignment()->lhs()->is_derivative();
+
+ // Save the state variables
+ dvars_.push_back(id->name());
+
+ // Create identifiers out of the state variables, they will be used to initialize some values
+ auto dvar_ident = make_expression<IdentifierExpression>(e->location(), id->name());
+
+ // Create two sets of local variables and assign them to the initial values of the state variables
+ // The first set doesn't change across iterations
+ auto init_dvar_term = make_unique_local_assign(e->scope(), dvar_ident.get(), "p_");
+ auto init_dvar = init_dvar_term.id->is_identifier()->spelling();
+
+ // The second set is updated in every iteration of Newton's method
+ auto temp_dvar_term = make_unique_local_assign(e->scope(), dvar_ident.get(), "t_");
+ auto temp_dvar = temp_dvar_term.id->is_identifier()->spelling();
+
+ // Save the variable names
+ dvar_init_.push_back(init_dvar);
+ dvar_temp_.push_back(temp_dvar);
+
+ // Add local declarations and assignment statements
+ statements_.push_back(std::move(init_dvar_term.local_decl));
+ statements_.push_back(std::move(init_dvar_term.assignment));
+
+ statements_.push_back(std::move(temp_dvar_term.local_decl));
+ statements_.push_back(std::move(temp_dvar_term.assignment));
+ }
+ }
+ scale_factor_.resize(dvars_.size());
+
+ BlockRewriterBase::visit(e);
+}
+
+void SparseNonlinearSolverVisitor::visit(CompartmentExpression *e) {
+ auto loc = e->location();
+
+ for (auto& s: e->is_compartment()->state_vars()) {
+ auto it = std::find(dvars_.begin(), dvars_.end(), s->is_identifier()->spelling());
+ if (it == dvars_.end()) {
+ error({"COMPARTMENT variable is not used", loc});
+ return;
+ }
+ auto idx = it - dvars_.begin();
+
+ auto scale_inv = make_expression<DivBinaryExpression>(
+ loc, make_expression<NumberExpression>(loc, 1.0), e->scale_factor()->clone());
+ auto local_s_term = make_unique_local_assign(e->scope(), scale_inv.get(), "s_");
+
+ statements_.push_back(std::move(local_s_term.local_decl));
+ statements_.push_back(std::move(local_s_term.assignment));
+ scale_factor_[idx] = std::move(local_s_term.id);
+ }
+}
+
+void SparseNonlinearSolverVisitor::visit(AssignmentExpression *e) {
+ if (A_.empty()) {
+ unsigned n = dvars_.size();
+ A_ = symge::sym_matrix(n, n);
+ }
+
+ auto loc = e->location();
+ scope_ptr scope = e->scope();
+
+ auto lhs = e->lhs();
+ auto rhs = e->rhs();
+ auto deriv = lhs->is_derivative();
+
+ if (!deriv) {
+ statements_.push_back(e->clone());
+
+ auto id = lhs->is_identifier();
+ if (id) {
+ auto expand = substitute(rhs, local_expr_);
+ if (involves_identifier(expand, dvars_)) {
+ local_expr_[id->spelling()] = std::move(expand);
+ }
+ }
+ return;
+ }
+
+ auto s = deriv->name();
+ auto expanded_rhs = substitute(rhs, local_expr_);
+
+ // Populate sparse symbolic matrix for GE.
+ if (s!=dvars_[deq_index_]) {
+ error({"ICE: inconsistent ordering of derivative assignments", loc});
+ return;
+ }
+
+ auto dt_expr = make_expression<IdentifierExpression>(loc, "dt");
+ auto one_expr = make_expression<NumberExpression>(loc, 1.0);
+
+ // Form and save F(y) = y - x(t) - dt G(y)
+ // y are stored in dvar_temp_ and are updated at every iteration of Newton's method
+ // x(t) are stored in dvar_init_ and are constant across iterations of Newton's method
+ // G(y) is the rhs of the derivative assignment expression
+ // Newton's method is used to find x(t+1) = y s.t. F(y) = 0.
+
+ // `scale_factors` multiply the lhs of a differential equation.
+ // After applying Backward Euler and Newton's method, the new F(y) is
+ // F(y) = y - x(t) - dt * diag(s^-1) * G(y)
+
+ expression_ptr F_x;
+ F_x = make_expression<MulBinaryExpression>(loc, expanded_rhs->clone(), dt_expr->clone());
+ if (scale_factor_[deq_index_]) {
+ F_x = make_expression<MulBinaryExpression>(loc, std::move(F_x), scale_factor_[deq_index_]->clone());
+ }
+ F_x = make_expression<AddBinaryExpression>(loc, make_expression<IdentifierExpression>(loc, dvar_init_[deq_index_]), std::move(F_x));
+ F_x = make_expression<SubBinaryExpression>(loc, make_expression<IdentifierExpression>(loc, dvar_temp_[deq_index_]), std::move(F_x));
+
+ for (unsigned k = 0; k < dvars_.size(); ++k) {
+ F_x = substitute(F_x, dvars_[k], make_expression<IdentifierExpression>(loc, dvar_temp_[k]));
+ }
+
+ auto local_f_term = make_unique_local_assign(scope, F_x.get(), "f_");
+ auto a_ = local_f_term.id->is_identifier()->spelling();
+
+ statements_.push_back(std::move(local_f_term.local_decl));
+ F_.push_back(std::move(local_f_term.assignment));
+
+ // Form and save the Jacobian J(x) of F(x)
+ // J(x) = I - dt * ∂G(x)/∂x
+ // If scale_factor[x] exists
+ // J(x) = I - dt * diag(s^-1) * ∂G(x)/∂x
+
+ linear_test_result r = linear_test(expanded_rhs, dvars_);
+ for (unsigned j = 0; j<dvars_.size(); ++j) {
+ expression_ptr J_x;
+
+ // For zero coefficient and diagonal element, the matrix entry is 1.
+ // For non-zero coefficient c and diagonal element, the entry is 1-c*dt.
+ // Otherwise, for non-zero coefficient c, the entry is -c*dt.
+
+ if (r.coef.count(dvars_[j])) {
+ J_x = make_expression<MulBinaryExpression>(loc,
+ r.coef[dvars_[j]]->clone(),
+ dt_expr->clone());
+
+ if (scale_factor_[deq_index_]) {
+ J_x = make_expression<MulBinaryExpression>(loc, std::move(J_x), scale_factor_[deq_index_]->clone());
+ }
+ }
+
+ if (j==deq_index_) {
+ if (J_x) {
+ J_x = make_expression<SubBinaryExpression>(loc, one_expr->clone(), std::move(J_x));
+ }
+ else {
+ J_x = one_expr->clone();
+ }
+ }
+ else if (J_x) {
+ J_x = make_expression<NegUnaryExpression>(loc, std::move(J_x));
+ }
+
+ if (J_x) {
+ for (unsigned k = 0; k < dvars_.size(); ++k) {
+ J_x = substitute(J_x, dvars_[k], make_expression<IdentifierExpression>(loc, dvar_temp_[k]));
+ }
+ }
+
+ if (!J_x) continue;
+
+ auto local_j_term = make_unique_local_assign(scope, J_x.get(), "j_");
+ auto j_ = local_j_term.id->is_identifier()->spelling();
+
+ statements_.push_back(std::move(local_j_term.local_decl));
+ J_.push_back(std::move(local_j_term.assignment));
+
+ A_[deq_index_].push_back({j, symtbl_.define(j_)});
+ }
+ ++deq_index_;
+}
+
+void SparseNonlinearSolverVisitor::finalize() {
+
+ // Create rhs of A_
+ std::vector<symge::symbol> rhs;
+ for (const auto& var: F_) {
+ auto id = var->is_assignment()->lhs()->is_identifier()->spelling();
+ rhs.push_back(symtbl_.define(id));
+ }
+
+ A_.augment(rhs);
+ symge::gj_reduce(A_, symtbl_);
+
+ // Create and assign intermediate variables.
+ // Save gaussian elimination solution statements in S_
+ std::vector<expression_ptr> S_;
+ for (unsigned i = 0; i < symtbl_.size(); ++i) {
+ symge::symbol s = symtbl_[i];
+
+ if (primitive(s)) continue;
+
+ auto expr = as_expression(definition(s));
+ auto local_t_term = make_unique_local_assign(block_scope_, expr.get(), "s_");
+ auto t_ = local_t_term.id->is_identifier()->spelling();
+ symtbl_.name(s, t_);
+
+ statements_.push_back(std::move(local_t_term.local_decl));
+ S_.push_back(std::move(local_t_term.assignment));
+ }
+
+ // Create the statements that update the temporary state variables
+ // (dvar_temp) after a Newton's iteration and save them in U_
+ std::vector<expression_ptr> U_;
+
+ Location loc;
+ for (unsigned i = 0; i < A_.nrow(); ++i) {
+ const symge::sym_row &row = A_[i];
+ unsigned rhs_col = A_.augcol();
+ int lhs_col = -1;
+ for (unsigned r = 0; r < A_.nrow(); r++) {
+ if (row[r]) {
+ lhs_col = r;
+ break;
+ }
+ }
+
+ if (lhs_col == -1) {
+ error({"Something went wrong in solving the sparse matrix", loc});
+ return;
+ }
+
+ auto var = make_expression<IdentifierExpression>(loc, dvar_temp_[lhs_col]);
+ auto expr = make_expression<AssignmentExpression>(loc, var->clone(),
+ make_expression<SubBinaryExpression>(loc, var->clone(),
+ make_expression<DivBinaryExpression>(loc, make_expression<IdentifierExpression>(loc, symge::name(A_[i][rhs_col])),
+ make_expression<IdentifierExpression>(loc, symge::name(A_[i][lhs_col])))));
+
+ U_.push_back(std::move(expr));
+ }
+
+ // Do 3 Newton iterations
+ for (unsigned n = 0; n < 3; n++) {
+ // Print out the statements that calulate F(xn), J(xn), solve J(xn)^-1*F(xn), update xn -> xn+1
+ for (auto &s: F_) {
+ statements_.push_back(s->clone());
+ }
+ for (auto &s: J_) {
+ statements_.push_back(s->clone());
+ }
+ for (auto &s: S_) {
+ statements_.push_back(s->clone());
+ }
+ for (auto &s: U_) {
+ statements_.push_back(s->clone());
+ }
+ }
+
+ // Finally update the state variables, dvars_ = dvar_temp_
+ for (unsigned i = 0; i < dvars_.size(); ++i) {
+ auto expr = make_expression<AssignmentExpression>(loc,
+ make_expression<IdentifierExpression>(loc, dvars_[i]),
+ make_expression<IdentifierExpression>(loc, dvar_temp_[i]));
+ statements_.push_back(std::move(expr));
+ }
+
+ BlockRewriterBase::finalize();
+}
+
// Implementation for `remove_unused_locals`: uses two visitors,
// `UnusedVisitor` and `RemoveVariableVisitor` below.
@@ -599,7 +868,7 @@ public:
std::set<std::string> unused_locals() {
if (!computed_) {
for (auto& id: used_ids) {
- remove_deps_from_unused(id);
+ remove_deps_from_unused(id, {});
}
computed_ = true;
}
@@ -616,11 +885,12 @@ public:
private:
bool computed_ = false;
- void remove_deps_from_unused(const std::string& id) {
+ void remove_deps_from_unused(const std::string& id, std::set<std::string> visited) {
auto range = deps.equal_range(id);
for (auto i = range.first; i != range.second; ++i) {
- if (unused_ids.count(i->second)) {
- remove_deps_from_unused(i->second);
+ if (unused_ids.count(i->second) && visited.find(i->second) == visited.end()) {
+ visited.insert(i->second);
+ remove_deps_from_unused(i->second, visited);
}
}
unused_ids.erase(id);
diff --git a/modcc/solvers.hpp b/modcc/solvers.hpp
index 01bdba33a3499a61c5d6dac8bc28c0e624bc68fd..b375edc9b88f93dfb98c50236b6c006b16783b13 100644
--- a/modcc/solvers.hpp
+++ b/modcc/solvers.hpp
@@ -121,6 +121,55 @@ public:
}
};
+class SparseNonlinearSolverVisitor : public SolverVisitorBase {
+protected:
+ // 'Current' differential equation is for variable with this
+ // index in `dvars`.
+ unsigned deq_index_ = 0;
+
+ // Expanded local assignments that need to be substituted in for derivative
+ // calculations.
+ substitute_map local_expr_;
+
+ // F(x) and the Jacobian J(x) for every state variable
+ // Needed for Newton's method
+ std::vector<expression_ptr> F_;
+ std::vector<expression_ptr> J_;
+
+ std::vector<std::string> dvar_temp_;
+ std::vector<std::string> dvar_init_;
+
+ // Symbolic matrix for backwards Euler step.
+ symge::sym_matrix A_;
+
+ // 'Symbol table' for symbolic manipulation.
+ symge::symbol_table symtbl_;
+
+ // State variable multiplier/divider
+ std::vector<expression_ptr> scale_factor_;
+
+public:
+ using SolverVisitorBase::visit;
+
+ SparseNonlinearSolverVisitor() {}
+ SparseNonlinearSolverVisitor(scope_ptr enclosing): SolverVisitorBase(enclosing) {}
+
+ virtual void visit(BlockExpression* e) override;
+ virtual void visit(AssignmentExpression *e) override;
+ virtual void visit(CompartmentExpression *e) override;
+ virtual void visit(ConserveExpression *e) override {};
+ virtual void finalize() override;
+ virtual void reset() override {
+ deq_index_ = 0;
+ local_expr_.clear();
+ F_.clear();
+ J_.clear();
+ symtbl_.clear();
+ scale_factor_.clear();
+ SolverVisitorBase::reset();
+ }
+};
+
class LinearSolverVisitor : public SolverVisitorBase {
protected:
// 'Current' differential equation is for variable with this
diff --git a/test/unit/CMakeLists.txt b/test/unit/CMakeLists.txt
index 29e8591e08876e5c1c946e17f46e89d546a53ad5..6791ccabe9b9de6c12e0a218787edeff5c31020d 100644
--- a/test/unit/CMakeLists.txt
+++ b/test/unit/CMakeLists.txt
@@ -13,6 +13,9 @@ set(test_mechanisms
test1_kin_diff
test1_kin_conserve
test1_kin_compartment
+ test2_kin_diff
+ test3_kin_diff
+ test4_kin_compartment
test1_kin_steadystate
fixed_ica_current
point_ica_current
diff --git a/test/unit/mod/test2_kin_diff.mod b/test/unit/mod/test2_kin_diff.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6b2ef3b88817ec8a44f62521b3edc9fd68ad76df
--- /dev/null
+++ b/test/unit/mod/test2_kin_diff.mod
@@ -0,0 +1,25 @@
+NEURON {
+ SUFFIX test2_kin_diff
+}
+
+STATE {
+ a b c
+}
+
+BREAKPOINT {
+ SOLVE state METHOD sparse
+}
+
+KINETIC state {
+ LOCAL f, r
+ f = 2
+ r = 1
+
+ ~ 2a + b <-> c (f, r)
+}
+
+INITIAL {
+ a = 0.2
+ b = 0.3
+ c = 0.5
+}
diff --git a/test/unit/mod/test3_kin_diff.mod b/test/unit/mod/test3_kin_diff.mod
new file mode 100644
index 0000000000000000000000000000000000000000..286def99ce9b635fadd1fdada126cb5f1eb63d1d
--- /dev/null
+++ b/test/unit/mod/test3_kin_diff.mod
@@ -0,0 +1,28 @@
+NEURON {
+ SUFFIX test3_kin_diff
+}
+
+STATE {
+ a b c
+}
+
+BREAKPOINT {
+ SOLVE state METHOD sparse
+}
+
+KINETIC state {
+ LOCAL f0, f1, r0, r1
+ f0 = 2
+ r0 = 1
+ f1 = 3
+ r1 = 0
+
+ ~ a + b <-> c (f0, r0)
+ ~ c <-> b (f1, r1)
+}
+
+INITIAL {
+ a = 0.2
+ b = 0.3
+ c = 0.5
+}
diff --git a/test/unit/mod/test4_kin_compartment.mod b/test/unit/mod/test4_kin_compartment.mod
new file mode 100644
index 0000000000000000000000000000000000000000..28bc29897eced1dcd00828639eee5623f3262150
--- /dev/null
+++ b/test/unit/mod/test4_kin_compartment.mod
@@ -0,0 +1,41 @@
+NEURON {
+ SUFFIX test4_kin_compartment
+}
+
+STATE {
+ A
+ B
+ C
+ d
+ e
+}
+
+PARAMETER {
+ x = 4
+ y = 0.4
+ z = 1.8
+ w = 21
+ s0 = 27
+ s1 = 5.9
+}
+
+BREAKPOINT {
+ SOLVE state METHOD sparse
+}
+
+INITIAL {
+ A = 4.5
+ B = 6.6
+ C = 0.28
+ d = 2
+ e = 0
+
+}
+
+KINETIC state {
+ COMPARTMENT s0 {A B C}
+ COMPARTMENT s1 {d e}
+
+ ~ A + B <-> C ( x, y )
+ ~ A + d <-> e ( z, w )
+}
\ No newline at end of file
diff --git a/test/unit/test_kinetic_linear.cpp b/test/unit/test_kinetic_linear.cpp
index e4bd84035be2b504be96dc69468edef99b389218..5a9d00965d1e6740da9c89cf3974af686df87f51 100644
--- a/test/unit/test_kinetic_linear.cpp
+++ b/test/unit/test_kinetic_linear.cpp
@@ -30,7 +30,8 @@ void run_test(std::string mech_name,
std::vector<std::string> state_variables,
std::unordered_map<std::string, fvm_value_type> assigned_variables,
std::vector<fvm_value_type> t0_values,
- std::vector<fvm_value_type> t1_values) {
+ std::vector<fvm_value_type> t1_values,
+ fvm_value_type dt) {
auto cat = make_unit_test_catalogue();
@@ -75,7 +76,7 @@ void run_test(std::string mech_name,
}
}
- shared_state->update_time_to(0.5, 0.5);
+ shared_state->update_time_to(dt, dt);
shared_state->set_dt();
test->nrn_state();
@@ -89,89 +90,135 @@ void run_test(std::string mech_name,
}
}
-TEST(mech_kinetic, kintetic_scaled) {
+TEST(mech_kinetic, kintetic_linear_scaled) {
std::vector<std::string> state_variables = {"s", "h", "d"};
std::vector<fvm_value_type> t0_values = {0.5, 0.2, 0.3};
std::vector<fvm_value_type> t1_0_values = {0.373297, 0.591621, 0.0350817};
std::vector<fvm_value_type> t1_1_values = {0.329897, 0.537371, 0.132732};
- run_test<multicore::backend>("test0_kin_compartment", state_variables, {}, t0_values, t1_0_values);
- run_test<multicore::backend>("test1_kin_compartment", state_variables, {}, t0_values, t1_1_values);
+ run_test<multicore::backend>("test0_kin_compartment", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<multicore::backend>("test1_kin_compartment", state_variables, {}, t0_values, t1_1_values, 0.5);
+
}
-TEST(mech_kinetic, kintetic_1_conserve) {
+TEST(mech_kinetic, kintetic_linear_1_conserve) {
std::vector<std::string> state_variables = {"s", "h", "d"};
std::vector<fvm_value_type> t0_values = {0.5, 0.2, 0.3};
std::vector<fvm_value_type> t1_0_values = {0.380338, 0.446414, 0.173247};
std::vector<fvm_value_type> t1_1_values = {0.218978, 0.729927, 0.0510949};
- run_test<multicore::backend>("test0_kin_diff", state_variables, {}, t0_values, t1_0_values);
- run_test<multicore::backend>("test0_kin_conserve", state_variables, {}, t0_values, t1_0_values);
- run_test<multicore::backend>("test0_kin_steadystate", state_variables, {}, t0_values, t1_1_values);
+ run_test<multicore::backend>("test0_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<multicore::backend>("test0_kin_conserve", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<multicore::backend>("test0_kin_steadystate", state_variables, {}, t0_values, t1_1_values, 0.5);
}
-TEST(mech_kinetic, kintetic_2_conserve) {
+TEST(mech_kinetic, kintetic_linear_2_conserve) {
std::vector<std::string> state_variables = {"a", "b", "x", "y"};
std::vector<fvm_value_type> t0_values = {0.2, 0.8, 0.6, 0.4};
std::vector<fvm_value_type> t1_0_values = {0.217391304, 0.782608696, 0.33333333, 0.66666666};
std::vector<fvm_value_type> t1_1_values = {0.230769, 0.769231, 0.189189, 0.810811};
- run_test<multicore::backend>("test1_kin_diff", state_variables, {}, t0_values, t1_0_values);
- run_test<multicore::backend>("test1_kin_conserve", state_variables, {}, t0_values, t1_0_values);
- run_test<multicore::backend>("test1_kin_steadystate", state_variables, {}, t0_values, t1_1_values);
+ run_test<multicore::backend>("test1_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<multicore::backend>("test1_kin_conserve", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<multicore::backend>("test1_kin_steadystate", state_variables, {}, t0_values, t1_1_values, 0.5);
+}
+
+TEST(mech_kinetic, kintetic_nonlinear) {
+ std::vector<std::string> state_variables = {"a", "b", "c"};
+ std::vector<fvm_value_type> t0_values = {0.2, 0.3, 0.5};
+ std::vector<fvm_value_type> t1_0_values = {0.222881, 0.31144, 0.48856};
+ std::vector<fvm_value_type> t1_1_values = {0.2078873133, 0.34222075, 0.45777925};
+
+ run_test<multicore::backend>("test2_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.025);
+ run_test<multicore::backend>("test3_kin_diff", state_variables, {}, t0_values, t1_1_values, 0.025);
+
+}
+
+TEST(mech_kinetic, kintetic_nonlinear_scaled) {
+ std::vector<std::string> state_variables = {"A", "B", "C", "d", "e"};
+ std::vector<fvm_value_type> t0_values = {4.5, 6.6, 0.28, 2, 0};
+ std::vector<fvm_value_type> t1_values = {4.087281958014442,
+ 6.224088678118931,
+ 0.6559113218810689,
+ 1.8315624742412617,
+ 0.16843752575873824};
+
+ run_test<multicore::backend>("test4_kin_compartment", state_variables, {}, t0_values, t1_values, 0.1);
}
-TEST(mech_linear, linear) {
+TEST(mech_linear, linear_system) {
std::vector<std::string> state_variables = {"h", "s", "d"};
std::vector<fvm_value_type> values = {0.5, 0.2, 0.3};
std::unordered_map<std::string, fvm_value_type> assigned_variables = {{"a0", 2.5}, {"a1",0.5}, {"a2",3}, {"a3",2.3}};
- run_test<multicore::backend>("test_linear_state", state_variables, assigned_variables, {}, values);
- run_test<multicore::backend>("test_linear_init", state_variables, assigned_variables, values, {});
- run_test<multicore::backend>("test_linear_init_shuffle", state_variables, assigned_variables, values, {});
+ run_test<multicore::backend>("test_linear_state", state_variables, assigned_variables, {}, values, 0.5);
+ run_test<multicore::backend>("test_linear_init", state_variables, assigned_variables, values, {}, 0.5);
+ run_test<multicore::backend>("test_linear_init_shuffle", state_variables, assigned_variables, values, {}, 0.5);
}
#ifdef ARB_GPU_ENABLED
-TEST(mech_kinetic_gpu, kintetic_scaled) {
+TEST(mech_kinetic_gpu, kintetic_linear_scaled) {
std::vector<std::string> state_variables = {"s", "h", "d"};
std::vector<fvm_value_type> t0_values = {0.5, 0.2, 0.3};
std::vector<fvm_value_type> t1_0_values = {0.373297, 0.591621, 0.0350817};
std::vector<fvm_value_type> t1_1_values = {0.329897, 0.537371, 0.132732};
- run_test<gpu::backend>("test0_kin_compartment", state_variables, {}, t0_values, t1_0_values);
- run_test<gpu::backend>("test1_kin_compartment", state_variables, {}, t0_values, t1_1_values);
+ run_test<gpu::backend>("test0_kin_compartment", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<gpu::backend>("test1_kin_compartment", state_variables, {}, t0_values, t1_1_values, 0.5);
}
-TEST(mech_kinetic_gpu, kintetic_1_conserve) {
+TEST(mech_kinetic_gpu, kintetic_linear_1_conserve) {
std::vector<std::string> state_variables = {"s", "h", "d"};
std::vector<fvm_value_type> t0_values = {0.5, 0.2, 0.3};
std::vector<fvm_value_type> t1_0_values = {0.380338, 0.446414, 0.173247};
std::vector<fvm_value_type> t1_1_values = {0.218978, 0.729927, 0.0510949};
- run_test<gpu::backend>("test0_kin_diff", state_variables, {}, t0_values, t1_0_values);
- run_test<gpu::backend>("test0_kin_conserve", state_variables, {}, t0_values, t1_0_values);
- run_test<gpu::backend>("test0_kin_steadystate", state_variables, {}, t0_values, t1_1_values);
+ run_test<gpu::backend>("test0_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<gpu::backend>("test0_kin_conserve", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<gpu::backend>("test0_kin_steadystate", state_variables, {}, t0_values, t1_1_values, 0.5);
}
-TEST(mech_kinetic_gpu, kintetic_2_conserve) {
+TEST(mech_kinetic_gpu, kintetic_linear_2_conserve) {
std::vector<std::string> state_variables = {"a", "b", "x", "y"};
std::vector<fvm_value_type> t0_values = {0.2, 0.8, 0.6, 0.4};
std::vector<fvm_value_type> t1_0_values = {0.217391304, 0.782608696, 0.33333333, 0.66666666};
std::vector<fvm_value_type> t1_1_values = {0.230769, 0.769231, 0.189189, 0.810811};
- run_test<gpu::backend>("test1_kin_diff", state_variables, {}, t0_values, t1_0_values);
- run_test<gpu::backend>("test1_kin_conserve", state_variables, {}, t0_values, t1_0_values);
- run_test<gpu ::backend>("test1_kin_steadystate", state_variables, {}, t0_values, t1_1_values);
+ run_test<gpu::backend>("test1_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<gpu::backend>("test1_kin_conserve", state_variables, {}, t0_values, t1_0_values, 0.5);
+ run_test<gpu::backend>("test1_kin_steadystate", state_variables, {}, t0_values, t1_1_values, 0.5);
+}
+
+TEST(mech_kinetic, kintetic_nonlinear) {
+ std::vector<std::string> state_variables = {"a", "b", "c"};
+ std::vector<fvm_value_type> t0_values = {0.2, 0.3, 0.5};
+ std::vector<fvm_value_type> t1_0_values = {0.222881, 0.31144, 0.48856};
+ std::vector<fvm_value_type> t1_1_values = {0.2078873133, 0.34222075, 0.45777925};
+
+ run_test<gpu::backend>("test2_kin_diff", state_variables, {}, t0_values, t1_0_values, 0.025);
+ run_test<gpu::backend>("test3_kin_diff", state_variables, {}, t0_values, t1_1_values, 0.025);
+}
+
+TEST(mech_kinetic, kintetic_nonlinear_scaled) {
+ std::vector<std::string> state_variables = {"A", "B", "C", "d", "e"};
+ std::vector<fvm_value_type> t0_values = {4.5, 6.6, 0.28, 2, 0};
+ std::vector<fvm_value_type> t1_values = {4.087281958014442,
+ 6.224088678118931,
+ 0.6559113218810689,
+ 1.8315624742412617,
+ 0.16843752575873824};
+
+ run_test<gpu::backend>("test4_kin_compartment", state_variables, {}, t0_values, t1_values, 0.1);
}
-TEST(mech_linear_gpu, linear) {
+TEST(mech_linear_gpu, linear_system) {
std::vector<std::string> state_variables = {"h", "s", "d"};
std::vector<fvm_value_type> values = {0.5, 0.2, 0.3};
std::unordered_map<std::string, fvm_value_type> assigned_variables = {{"a0", 2.5},{"a1",0.5},{"a2",3},{"a3",2.3}};
- run_test<gpu::backend>("test_linear_state", state_variables, assigned_variables, {}, values);
- run_test<gpu::backend>("test_linear_init", state_variables, assigned_variables, values, {});
- run_test<gpu::backend>("test_linear_init_shuffle", state_variables, assigned_variables, values, {});
+ run_test<gpu::backend>("test_linear_state", state_variables, assigned_variables, {}, values, 0.5);
+ run_test<gpu::backend>("test_linear_init", state_variables, assigned_variables, values, {}, 0.5);
+ run_test<gpu::backend>("test_linear_init_shuffle", state_variables, assigned_variables, values, {}, 0.5);
}
#endif
diff --git a/test/unit/unit_test_catalogue.cpp b/test/unit/unit_test_catalogue.cpp
index f86219ac8efb48e8ba71973c0eea89e2a2396d33..e68bbdfbc056b62a83d861adadffb8fe32e8066d 100644
--- a/test/unit/unit_test_catalogue.cpp
+++ b/test/unit/unit_test_catalogue.cpp
@@ -19,6 +19,9 @@
#include "mechanisms/test1_kin_compartment.hpp"
#include "mechanisms/test1_kin_diff.hpp"
#include "mechanisms/test1_kin_conserve.hpp"
+#include "mechanisms/test2_kin_diff.hpp"
+#include "mechanisms/test3_kin_diff.hpp"
+#include "mechanisms/test4_kin_compartment.hpp"
#include "mechanisms/test1_kin_steadystate.hpp"
#include "mechanisms/fixed_ica_current.hpp"
#include "mechanisms/point_ica_current.hpp"
@@ -58,8 +61,11 @@ mechanism_catalogue make_unit_test_catalogue() {
ADD_MECH(cat, test0_kin_compartment)
ADD_MECH(cat, test1_kin_diff)
ADD_MECH(cat, test1_kin_conserve)
+ ADD_MECH(cat, test2_kin_diff)
+ ADD_MECH(cat, test3_kin_diff)
ADD_MECH(cat, test1_kin_steadystate)
ADD_MECH(cat, test1_kin_compartment)
+ ADD_MECH(cat, test4_kin_compartment)
ADD_MECH(cat, fixed_ica_current)
ADD_MECH(cat, point_ica_current)
ADD_MECH(cat, linear_ca_conc)