nmodl
nmodl copied to clipboard
BlueBrain
Improve Eigen solver node in the ast to allow local variables and conditional expressions
We have mod files with derivative block and derivimplicit method :
DERIVATIVE state {
LOCAL minf_VDCC, hinf_VDCC
minf_VDCC = 1/(1+exp(((vhm_VDCC-ljp_VDCC)-v)/km_VDCC))
hinf_VDCC = 1/(1+exp(((vhh_VDCC-ljp_VDCC)-v)/kh_VDCC))
A_AMPA' = -A_AMPA/tau_r_AMPA
B_AMPA' = -B_AMPA/tau_d_AMPA
A_NMDA' = -A_NMDA/tau_r_NMDA
B_NMDA' = -B_NMDA/tau_d_NMDA
m_VDCC' = (minf_VDCC-m_VDCC)/mtau_VDCC
h_VDCC' = (hinf_VDCC-h_VDCC)/htau_VDCC
cai_CR' = -(1e-09)*(ica_NMDA+ica_VDCC)*gamma_ca_CR/((1e-15)*volume_CR*2*FARADAY)-(cai_CR-min_ca_CR)/tau_ca_CR
effcai_GB' = -effcai_GB/tau_effca_GB+(cai_CR-min_ca_CR)
Rho_GB' = (-Rho_GB*(1-Rho_GB)*(rho_star_GB-Rho_GB)+potentiate_GB*gamma_p_GB*(1-Rho_GB)-depress_GB*gamma_d_GB*Rho_GB)/((1000)*tau_GB)
Use_GB' = (Use_d_GB+Rho_GB*(Use_p_GB-Use_d_GB)-Use_GB)/((1000)*tau_Use_GB)
min_ca_CR = Use_GB
}
The current code generation will produce :
void nrn_state_GluSynapse(NrnThread* nt, Memb_list* ml, int type) {
int nodecount = ml->nodecount;
int pnodecount = ml->_nodecount_padded;
const int* __restrict__ node_index = ml->nodeindices;
double* __restrict__ data = ml->data;
const double* __restrict__ voltage = nt->_actual_v;
Datum* __restrict__ indexes = ml->pdata;
ThreadDatum* __restrict__ thread = ml->_thread;
GluSynapse_Instance* __restrict__ inst = (GluSynapse_Instance*) ml->instance;
int start = 0;
int end = nodecount;
#pragma ivdep
for (int id = start; id < end; id++) {
int node_id = node_index[id];
double v = voltage[node_id];
double minf_VDCC, hinf_VDCC;
minf_VDCC = 1.0/(1.0+exp(((GluSynapse_global.vhm_VDCC-GluSynapse_global.ljp_VDCC)-v)/GluSynapse_global.km_VDCC));
hinf_VDCC = 1.0/(1.0+exp(((GluSynapse_global.vhh_VDCC-GluSynapse_global.ljp_VDCC)-v)/GluSynapse_global.kh_VDCC));
minf_VDCC = inst->Use_GB[id];
Eigen::Matrix<double, 10, 1> X;
X[0] = inst->A_AMPA[id];
X[1] = inst->B_AMPA[id];
X[2] = inst->A_NMDA[id];
X[3] = inst->B_NMDA[id];
X[4] = inst->m_VDCC[id];
X[5] = inst->h_VDCC[id];
X[6] = inst->cai_CR[id];
X[7] = inst->effcai_GB[id];
X[8] = inst->Rho_GB[id];
X[9] = inst->Use_GB[id];
struct functor {
NrnThread* nt;
GluSynapse_Instance* inst;
int id;
functor(NrnThread* nt, GluSynapse_Instance* inst, int id) : nt(nt), inst(inst), id(id) {}
void operator()(const Eigen::Matrix<double, 10, 1>& X, Eigen::Matrix<double, 10, 1>& F, Eigen::Matrix<double, 10, 10>& Jm) const{
double* J = Jm.data();
F[0] = -inst->A_AMPA[id]+X[0]*nt->_dt/GluSynapse_global.tau_r_AMPA+X[0];
F[1] = -inst->B_AMPA[id]+X[1]*nt->_dt/inst->tau_d_AMPA[id]+X[1];
F[2] = -inst->A_NMDA[id]+X[2]*nt->_dt/GluSynapse_global.tau_r_NMDA+X[2];
F[3] = -inst->B_NMDA[id]+X[3]*nt->_dt/GluSynapse_global.tau_d_NMDA+X[3];
F[4] = (nt->_dt*(X[4]-minf_VDCC)+GluSynapse_global.mtau_VDCC*(X[4]-inst->m_VDCC[id]))/GluSynapse_global.mtau_VDCC;
F[5] = (nt->_dt*(X[5]-hinf_VDCC)+GluSynapse_global.htau_VDCC*(X[5]-inst->h_VDCC[id]))/GluSynapse_global.htau_VDCC;
F[6] = (FARADAY*GluSynapse_global.tau_ca_CR*inst->volume_CR[id]*(X[6]-inst->cai_CR[id])+nt->_dt*(FARADAY*inst- >volume_CR[id]*(X[6]-GluSynapse_global.min_ca_CR)+499999.9999999999*GluSynapse_global.gamma_ca_CR*GluSynapse_global.tau_ca_CR*(inst- >ica_NMDA[id]+inst->ica_VDCC[id])))/(FARADAY*GluSynapse_global.tau_ca_CR*inst->volume_CR[id]);
F[7] = (nt->_dt*(X[7]+GluSynapse_global.tau_effca_GB*( -X[6]+GluSynapse_global.min_ca_CR))+GluSynapse_global. tau_effca_GB*(X[7]-inst->effcai_GB[id]))/GluSynapse_global.tau_effca_GB;
F[8] = (0.001*nt->_dt*(X[8]*inst->depress_GB[id]*GluSynapse_global.gamma_d_GB+X[8]*(X[8]-1.0)*(X[8]-GluSynapse_global. rho_star_GB)+GluSynapse_global.gamma_p_GB*inst->potentiate_GB[id]*(X[8]-1.0))+GluSynapse_global.tau_GB*( -inst->Rho_GB[id]+X[8]))/ GluSynapse_global.tau_GB;
F[9] = (0.001*nt->_dt*( -inst->Use_d_GB[id]+X[8]*(inst->Use_d_GB[id]-inst->Use_p_GB[id])+X[9])+GluSynapse_global. tau_Use_GB*( -inst->Use_GB[id]+X[9]))/GluSynapse_global.tau_Use_GB;
J[0] = nt->_dt/GluSynapse_global.tau_r_AMPA+1.0;
J[10] = 0.0;
J[20] = 0.0;
J[30] = 0.0;
In this case functor doesn't have access to :
double minf_VDCC, hinf_VDCC;
This bug prevents optimizing nmodl files.
Local variables are faster than assigned variables.
My new nmodl_preprocessor program automatically converts variables from assigned to local.
The best way to fix this bug is to just give up on trying to factor out the newton solver into a reusable peice of code. Instead you should generate a new instance of the newton solver algorithm for every model that uses it, and then inline the functor directly into the newton solver. For an example see the linear solver src/codegen/codegen_cpp_visitor.cpp CodegenCVisitor::print_eigen_linear_solver.
