#include "bcbg2.hpp"
#include "constants.hpp"
void BCBG2::SingleChannelNucleus::initialize_single_channel_nucleus(float Smax, float Sini, float vh, float k, int damping, float* dists, float* diams, int compartment_nb, const char* id)
{
int i,j;
// initialization of single_channel_nucleus parameters
this->Smax = Smax; // in Hz
//this->vh = vh;
this->vh = vh*1e-3; // convert vh from mV to V
//this->k = k;
this->k = k*1e3; // convert vh from mV⁻¹ to V⁻¹
this->kvh = k*vh; // kvh is dimensionless
this->id = (char*)malloc((1+strlen(id))*sizeof(char));
for (i=0;i<strlen(id);i++) {
this->id[i]=id[i];
}
this->id[strlen(id)] = '\0';
this->is_damped = damping;
// Computation of distance attenuation factor <
this->Rm = 20000.*1e-4; // conversion Ohm.cm² -> Ohm.m²
this->Ri = 200.*1e-2; // conversion Ohm.cm -> Ohm.m
this->membrane_area = 0.;
this->distances.resize(0); this->diameters.resize(0);
for (i=0; i<compartment_nb; i++) {this->distances.push_back(dists[i]); this->diameters.push_back(diams[i]); this->membrane_area += 2.*3.141592654*diams[i]*dists[i];}
// Computation of distance attenuation factor >
//nb_afferents = afferents_number;
// initialization of counter for afferents
n = 0;
// allocation of arrays
S.resize(bg_.max_tau);
N_in.resize(0);
A_in.resize(0);
D_in.resize(0);
C_in.resize(0);
ADC_in.resize(0);
DD_in.resize(0);
D2_in.resize(0);
nu_in.resize(0);
dtADC_in.resize(0);
dtDD_in.resize(0);
dtD2_in.resize(0);
Sign_in.resize(0);
T_in.resize(0);
H_in_RK4.resize(0);
Hp_in_RK4.resize(0);
// initialization to Sini for the initial discharge rate
for(i=0;i<bg_.max_tau;i++) {
S[i]=Sini;
}
}
void BCBG2::SingleChannelNucleus::set_afferent(float A, float D, int Sign, float C, int T, float to_be_ignored, float distance, SingleChannelNucleus* N)
{
SingleChannelNucleus::set_afferent(A, D, Sign, C, T, distance, N); // Dummy function, calls other "set_afferent"
}
void BCBG2::SingleChannelNucleus::set_afferent(float A, float D, int Sign, float C, int T, float distance, SingleChannelNucleus* N)
{
this->initialize_new_afferent();
A_in[n] = A*compute_distance_factor(distance)*1e-3; // convert A in V (from mV)
D_in[n] = D*1e3; // convert D in s⁻¹ (from ms⁻¹)
C_in[n] = C;
nu_in[n] = 0.;
ADC_in[n] = A_in[n]*D_in[n]*C;
DD_in[n] = D_in[n]*D_in[n];
D2_in[n] = 2.*D_in[n];
T_in[n] = T - bg_.max_tau;
Sign_in[n] = Sign;
N_in[n] = N;
dtADC_in[n] = ADC_in[n]*bg_.dt;
dtDD_in[n] = DD_in[n]*bg_.dt;
dtD2_in[n] = D2_in[n]*bg_.dt;
n++;
}
void BCBG2::SingleChannelNucleus::set_afferent(float nu, int T, SingleChannelNucleus* N)
{
this->initialize_new_afferent();
A_in[n] = 0;
D_in[n] = 0;
C_in[n] = 0;
nu_in[n] = nu;
ADC_in[n] = 0;
DD_in[n] = 0;
D2_in[n] = 0;
Sign_in[n] = 0;
T_in[n] = T - bg_.max_tau;
N_in[n] = N;
n++;
}
float BCBG2::SingleChannelNucleus::compute_distance_factor(float distance)
{
/* calculations of the distance factor */
float end_of_pseudocompart, total_length, lambda, L, X, pos, distance_factor;
total_length = 0.;
for (int j=0; j<distances.size(); j++) {
total_length += distances[j];
}
distance_factor = 1.;
pos = distance*total_length; //conversion ratio -> m
for (int i=distances.size()-1; i>=0; i--) {
end_of_pseudocompart = 0.;
for (int j=0; j<i; j++) {
end_of_pseudocompart += distances[j];
}
if (pos > end_of_pseudocompart) {
lambda = sqrt((diameters[i]*Rm)/(4.*Ri));
L = (distances[i]) / lambda;
X = (pos - end_of_pseudocompart) / lambda;
distance_factor *= cosh(L-X) / cosh(L);
pos = end_of_pseudocompart;
}
}
return distance_factor;
}
void BCBG2::SingleChannelNucleus::initialize_new_afferent()
{
N_in.resize(n+1);
A_in.resize(n+1);
D_in.resize(n+1);
C_in.resize(n+1);
ADC_in.resize(n+1);
DD_in.resize(n+1);
D2_in.resize(n+1);
nu_in.resize(n+1);
dtADC_in.resize(n+1);
dtDD_in.resize(n+1);
dtD2_in.resize(n+1);
Sign_in.resize(n+1);
T_in.resize(n+1);
H_in_RK4.resize(n+1); H_in_RK4[n].assign(bg_.max_tau,0);
Hp_in_RK4.resize(n+1); Hp_in_RK4[n].assign(bg_.max_tau,0);
}
void BCBG2::SingleChannelNucleus::set_dt(float value)
{
dt=value;
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus()
{
int i;
float sum_v = 0;
for (i=0; i<n; i++) {
Hp_in[i] = Hp_in_old[i] + dtADC_in[i] * N_in[i]->get_S((bg_.tmod - T_in[i]) % bg_.max_tau) - dtD2_in[i] * Hp_in_old[i] - dtDD_in[i] * H_in_old[i];
H_in[i] = H_in_old[i] + Hp_in[i] * dt;
H_in_old[i] = H_in[i];
Hp_in_old[i] = Hp_in[i];
sum_v += H_in[i] * Sign_in[i];
}
S[(bg_.tmod+1) % bg_.max_tau] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_stabilize(int steps)
{
int i,s;
int t = 1;
int t_1 = 0;
std::vector <float> virtual_S;
float H_previous = 0.;
float H = 0.;
float Hp_previous = 0.;
float Hp = 0.;
float input_S;
for (i=0; i<n; i++) {
input_S = dtADC_in[i] * N_in[i]->get_S(0);
for (s=steps; s; s--) {
// fastest version
Hp = Hp_previous + input_S - dtD2_in[i] * Hp_previous - dtDD_in[i] * H_previous;
H_previous = H_previous + dt * Hp_previous;
Hp_previous = Hp;
////// explicit Euler
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] +
// (
// virtual_S[i]
// - dtD2_in[i] * Hp_in_RK4[i][t_1]
// - dtDD_in[i] * H_in_RK4[i][t_1]
// );
////// explicit Euler
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1];
//std::cout << "just added " << dt * Hp_in_RK4[i][t_1] << " but Hp was " << Hp_in_RK4[i][t_1] << std::endl;
//t_1++;
//t++;
//if (t_1 == bg_.max_tau) {
// t_1 = 0;
//}
//if (t == bg_.max_tau) {
// t = 0;
//}
}
Hp_in_RK4[i].assign(bg_.max_tau,Hp_previous); H_in_RK4[i].assign(bg_.max_tau,H_previous);
}
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_evo_Hp()
{
int i;
float sum_v = 0;
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
int t_2 = (t - 2 + bg_.max_tau) % bg_.max_tau;
int t_3 = (t - 3 + bg_.max_tau) % bg_.max_tau;
int t_4 = (t - 4 + bg_.max_tau) % bg_.max_tau;
for (i=0; i<n; i++) {
k1 = dt * (
ADC_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4])
- DD_in[i] * H_in_RK4[i][t_4]
);
k2 = dt * (
ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+2.*k1)
- DD_in[i] * H_in_RK4[i][t_2]
);
k3 = dt * (
ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+3.*k2)
- DD_in[i] * H_in_RK4[i][t_1]
);
Hp_in_RK4[i][t] = Hp_in_RK4[i][t_4] + (8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3;
H_in_RK4[i][t] = H_in_RK4[i][t_1];
sum_v += H_in_RK4[i][t] * Sign_in[i];
}
S[t] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_euler()
{
int i;
float sum_v = 0;
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
for (i=0; i<n; i++) {
// This gives the same discharge rates as Tsirogiannis et al 2010
////// explicit Euler
Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] +
(
dtADC_in[i] * N_in[i]->get_S((t - T_in[i] - 1 + bg_.max_tau ) % bg_.max_tau)
- dtD2_in[i] * Hp_in_RK4[i][t_1]
- dtDD_in[i] * H_in_RK4[i][t_1]
);
////// explicit Euler
H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1];
sum_v += H_in_RK4[i][t] * Sign_in[i];
}
S[t] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_rk3()
{
int i;
float sum_v = 0;
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
int t_2 = (t - 2 + bg_.max_tau) % bg_.max_tau;
int t_4 = (t - 4 + bg_.max_tau) % bg_.max_tau;
for (i=0; i<n; i++) {
k1 = dt * (
ADC_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4])
- DD_in[i] * H_in_RK4[i][t_4]
);
k1bis = dt * (Hp_in_RK4[i][t_4]);
k2 = dt * (
ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+2.*k1)
- DD_in[i] * (H_in_RK4[i][t_4]+2.*k1bis)
);
k2bis = dt * (Hp_in_RK4[i][t_4]+2.*k1bis);
k3 = dt * (
ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+3.*k2)
- DD_in[i] * (H_in_RK4[i][t_4]+3.*k2bis)
);
k3bis = dt * (Hp_in_RK4[i][t_4]+3.*k2bis);
Hp_in_RK4[i][t] = Hp_in_RK4[i][t_4] + (8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3;
H_in_RK4[i][t] = H_in_RK4[i][t_4] + (8./9.) * k1bis + (4./3.) * k2bis + (16./9.) * k3bis;
sum_v += H_in_RK4[i][t] * Sign_in[i];
}
S[t] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_evo()
{
int i;
float sum_v = 0;
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
int t_2 = (t - 2 + bg_.max_tau) % bg_.max_tau;
int t_4 = (t - 4 + bg_.max_tau) % bg_.max_tau;
for (i=0; i<n; i++) {
#ifdef INTEGRATIONEULER
// This gives the same discharge rates as Tsirogiannis et al 2010
////// explicit Euler
Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] +
(
dtADC_in[i] * N_in[i]->get_S((t - T_in[i] - 1 + bg_.max_tau ) % bg_.max_tau)
- dtD2_in[i] * Hp_in_RK4[i][t_1]
- dtDD_in[i] * H_in_RK4[i][t_1]
);
////// explicit Euler
H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1];
#endif
//// This gives the same discharge rates as Tsirogiannis et al 2010
//////// explicit Euler
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1];
//////// explicit (?) Euler
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] +
// (
// dtADC_in[i] * N_in[i]->get_S(t_1)
// - dtD2_in[i] * Hp_in_RK4[i][t_1]
// - dtDD_in[i] * H_in_RK4[i][t_1]
// );
//// midpoint rule
//k1 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_2])
// - DD_in[i] * H_in_RK4[i][t_2]
// );
//k2 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_2]+k1)
// - DD_in[i] * H_in_RK4[i][t_1]
// );
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_2] + 2. * k2;
//k1 = dt * (Hp_in_RK4[i][t_2]);
//k2 = dt * (Hp_in_RK4[i][t_1]);
//H_in_RK4[i][t] = H_in_RK4[i][t_2] + 2. * k2;
//// Bogacki-Shampine
//// work on a simple example, not on the whole circuit
//k1 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_4])
// - DD_in[i] * H_in_RK4[i][t_4]
// );
//k2 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_4]+2.*k1)
// - DD_in[i] * H_in_RK4[i][t_2]
// );
//k3 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_4]+3.*k2)
// - DD_in[i] * H_in_RK4[i][t_1]
// );
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_4] + (8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3;
//k1 = dt * (Hp_in_RK4[i][t_4]);
//k2 = dt * (Hp_in_RK4[i][t_2]);
//k3 = dt * (Hp_in_RK4[i][t_1]);
//H_in_RK4[i][t] = H_in_RK4[i][t_4] + (8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3;
#ifdef INTEGRATIONRK3
// Bogacki-Shampine
k1 = dt * (
ADC_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4])
- DD_in[i] * H_in_RK4[i][t_4]
);
k1bis = dt * (Hp_in_RK4[i][t_4]);
k2 = dt * (
ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+2.*k1)
- DD_in[i] * (H_in_RK4[i][t_4]+2.*k1bis)
);
k2bis = dt * (Hp_in_RK4[i][t_4]+2.*k1bis);
k3 = dt * (
ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- D2_in[i] * (Hp_in_RK4[i][t_4]+3.*k2)
- DD_in[i] * (H_in_RK4[i][t_4]+3.*k2bis)
);
k3bis = dt * (Hp_in_RK4[i][t_4]+3.*k2bis);
Hp_in_RK4[i][t] = Hp_in_RK4[i][t_4] + (8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3;
H_in_RK4[i][t] = H_in_RK4[i][t_4] + (8./9.) * k1bis + (4./3.) * k2bis + (16./9.) * k3bis;
#endif
//// Newmark (tested only on a subproblem)
//float gamma=0.5;
//float beta=0.25;
//Hpp_in_RK4[i][t] =
// (
// dtADC_in[i] * N_in[i]->get_S((t - T_in[i] - 1 + bg_.max_tau ) % bg_.max_tau)
// - dtD2_in[i] * Hp_in_RK4[i][t_1]
// - dtDD_in[i] * H_in_RK4[i][t_1]
// );
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1] + ((dt * dt)/2) * ((1-2*beta)*Hpp_in_RK4[i][t_1] + 2*beta*Hpp_in_RK4[i][t]);
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] + dt * ((1-gamma)*Hpp_in_RK4[i][t_1] + gamma*Hpp_in_RK4[i][t]);
// Below are other integration schemes of high-order:
//k1 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-3 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_3])
// - DD_in[i] * H_in_RK4[i][t_3]
// );
//k2 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_3]+k1)
// - DD_in[i] * H_in_RK4[i][t_2]
// );
//k3 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_3]+k2)
// - DD_in[i] * H_in_RK4[i][t_2]
// );
//k4 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_3]+2.*k3)
// - DD_in[i] * H_in_RK4[i][t_1]
// );
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_3] + (1./3.) * (k1 + 2. * k2 + 2. * k3 + k4);
//k1 = dt * (Hp_in_RK4[i][t_3]);
//k2 = dt * (Hp_in_RK4[i][t_2]);
//k3 = dt * (Hp_in_RK4[i][t_2]);
//k4 = dt * (Hp_in_RK4[i][t_1]);
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + (1./6.) * (k1 + 2. * k2 + 2. * k3 + k4);
//k1 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_4])
// - DD_in[i] * H_in_RK4[i][t_4]
// );
//k2 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_2])
// - DD_in[i] * H_in_RK4[i][t_2]
// );
//k3 = dt * (
// ADC_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_1])
// - DD_in[i] * H_in_RK4[i][t_1]
// );
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_4] + ((8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3);
//k1 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t-90 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90])
// - DD_in[i] * H_in_RK4[i][t_90]
// );
//k2 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t-72 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90] + 1./5. * k1)
// - DD_in[i] * H_in_RK4[i][t_72]
// );
//k3 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t-63 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90] + 3./40. * k1 + 9./40. * k2)
// - DD_in[i] * H_in_RK4[i][t_63]
// );
//k4 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t-18 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90] + 44./45. * k1 - 56./15. * k2 + 32./9. * k3)
// - DD_in[i] * H_in_RK4[i][t_18]
// );
//k5 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t-10 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90] + 19372./6561. * k1 - 25360./2187. * k2 + 64448./6561. * k3 - 212./729. * k4)
// - DD_in[i] * H_in_RK4[i][t_10]
// );
//k6 = dt90 * (
// ADC_in[i] * N_in[i]->get_S((t - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - D2_in[i] * (Hp_in_RK4[i][t_90] + 9017./3168. * k1 - 355./33. * k2 + 46732./5247. * k3 + 49./176. * k4 - 5103./18656. * k5)
// - DD_in[i] * (H_in_RK4[i][t])
// );
//Hp_in_RK4[i][t] =
// Hp_in_RK4[i][t_90]
// + 35./384. * k1 + 500./1113. * k3 + 125./192. * k4 - 2187./6784. * k5 + 11./84. * k6;
sum_v += H_in_RK4[i][t] * Sign_in[i];
}
S[t] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_damping_vana()
{
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
int t_2 = (t - 2 + bg_.max_tau) % bg_.max_tau;
int t_4 = (t - 4 + bg_.max_tau) % bg_.max_tau;
if (is_damped) {
//// Bogacki-Shampine (ode23 in matlab) for damping
//k1 = dt * (Phip_in_RK4[t_4]);
//k2 = dt * (Phip_in_RK4[t_2] + 2. * dt * k1);
//k3 = dt * (Phip_in_RK4[t_1] + 3. * dt * k2);
//Phi_in_RK4[t] = Phi_in_RK4[t_4] + ((8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3);
// explicit Euler for damping Phi
Phi_in_RK4[t] = Phi_in_RK4[t_1] + dt * Phip_in_RK4[t_1];
// explicit Euler for damping Phip
Phip_in_RK4[t] = Phip_in_RK4[t_1] + dt *
(
15625 * this->get_S((t-1 +bg_.max_tau) % bg_.max_tau)
- 250 * Phip_in_RK4[t_1]
- 15625 * Phi_in_RK4[t_1]
);
}
}
void BCBG2::SingleChannelNucleus::update_single_channel_nucleus_vana()
{
// should be optimised by setting the dt factor inside parameters
int i;
float sum_v = 0;
float gamma=0.5;
float beta=0.25;
int t = bg_.tmod;
int t_1 = (t - 1 + bg_.max_tau) % bg_.max_tau;
int t_2 = (t - 2 + bg_.max_tau) % bg_.max_tau;
int t_3 = (t - 3 + bg_.max_tau) % bg_.max_tau;
int t_4 = (t - 4 + bg_.max_tau) % bg_.max_tau;
int t_5 = (t - 5 + bg_.max_tau) % bg_.max_tau;
int t_6 = (t - 6 + bg_.max_tau) % bg_.max_tau;
int t_7 = (t - 7 + bg_.max_tau) % bg_.max_tau;
int c=-1;
int t_90c = (t - 90 + c + bg_.max_tau) % bg_.max_tau;
int t_72c = (t - 72 + c + bg_.max_tau) % bg_.max_tau;
int t_63c = (t - 63 + c + bg_.max_tau) % bg_.max_tau;
int t_18c = (t - 18 + c + bg_.max_tau) % bg_.max_tau;
int t_10c = (t - 10 + c + bg_.max_tau) % bg_.max_tau;
int tc = (t - 1 + c + bg_.max_tau) % bg_.max_tau;
int t_90 = (t - 90 + bg_.max_tau) % bg_.max_tau;
int t_72 = (t - 72 + bg_.max_tau) % bg_.max_tau;
int t_63 = (t - 63 + bg_.max_tau) % bg_.max_tau;
int t_18 = (t - 18 + bg_.max_tau) % bg_.max_tau;
int t_10 = (t - 10 + bg_.max_tau) % bg_.max_tau;
float dt90 = 90*dt;
for (i=0; i<n; i++) {
//// explicit Euler
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1];
//// explicit trapezoidal (Heun)
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + 0.5 * dt *
// (
// Hp_in_RK4[i][t_1]
// +
// (H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1])
// );
// Bogacki-Shampine (ode23)
k1 = dt * (Hp_in_RK4[i][t_4]);
k2 = dt * (Hp_in_RK4[i][t_2] + 2. * dt * k1);
k3 = dt * (Hp_in_RK4[i][t_1] + 3. * dt * k2);
H_in_RK4[i][t] = H_in_RK4[i][t_4] + ((8./9.) * k1 + (4./3.) * k2 + (16./9.) * k3);
//// 2-steps Adams-Bashforth
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * ((3./2.)*Hp_in_RK4[i][t_1] - (1./2.)*Hp_in_RK4[i][t_2]);
//// 3-steps Adams-Bashforth ALT
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * ((23./12.)*Hp_in_RK4[i][t_1] - (4./3.)*Hp_in_RK4[i][t_2] + (5./12.)*Hp_in_RK4[i][t_3]);
//// 4-steps Adams-Bashforth
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * ((55./24.)*Hp_in_RK4[i][t_1] - (59./24.)*Hp_in_RK4[i][t_2] + (37./24.)*Hp_in_RK4[i][t_3] - (3./8.)*Hp_in_RK4[i][t_4]);
//// 5-steps Adams-Bashforth
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * ((1901./720.)*Hp_in_RK4[i][t_1] - (1387./360.)*Hp_in_RK4[i][t_2] + (109./30.)*Hp_in_RK4[i][t_3] - (637./360.)*Hp_in_RK4[i][t_4] + (251./720.)*Hp_in_RK4[i][t_5]);
//// Dormand Prince (approximated)
//k1 = dt90 * (Hp_in_RK4[i][t_90]);
//k2 = dt90 * (Hp_in_RK4[i][t_72] + 1./5. * k1);
//k3 = dt90 * (Hp_in_RK4[i][t_63] + 3./40. * k1 + 9./40. * k2);
//k4 = dt90 * (Hp_in_RK4[i][t_18] + 44./45. * k1 - 56./15. * k2 + 32./9. * k3);
//k5 = dt90 * (Hp_in_RK4[i][t_10] + 19372./6561. * k1 - 25360./2187. * k2 + 64448./6561. * k3 - 212./729. * k4);
//k6 = dt90 * (Hp_in_RK4[i][t_1] + 9017./3168. * k1 - 355./33. * k2 + 46732./5247. * k3 + 49./176. * k4 - 5103./18656. * k5);
//H_in_RK4[i][t] = H_in_RK4[i][t_90] + 35./384. * k1 + 500./1113. * k3 + 125./192. * k4 - 2187./6784. * k5 + 11./84. * k6;
//// Dormand Prince
//k1 = dt90 * (Hp_in_RK4[i][t_90c]);
//k2 = dt90 * (Hp_in_RK4[i][t_72c] + 1./5. * k1);
//k3 = dt90 * (Hp_in_RK4[i][t_63c] + 3./40. * k1 + 9./40. * k2);
//k4 = dt90 * (Hp_in_RK4[i][t_18c] + 44./45. * k1 - 56./15. * k2 + 32./9. * k3);
//k5 = dt90 * (Hp_in_RK4[i][t_10c] + 19372./6561. * k1 - 25360./2187. * k2 + 64448./6561. * k3 - 212./729. * k4);
//k6 = dt90 * (Hp_in_RK4[i][tc] + 9017./3168. * k1 - 355./33. * k2 + 46732./5247. * k3 + 49./176. * k4 - 5103./18656. * k5);
//H_in_RK4[i][t] = H_in_RK4[i][t_90c] + 35./384. * k1 + 500./1113. * k3 + 125./192. * k4 - 2187./6784. * k5 + 11./84. * k6;
//// explicit Euler
//Hp_in_RK4[i][bg_.tmod] = Hp_in_RK4[i][t_1] + dt *
// (
// 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau)
// - 800 * Hp_in_RK4[i][t_1]
// - 102400 * H_in_RK4[i][t_1]
// );
//// explicit Euler with better estimation of H (needs prior computation of H)
//Hp_in_RK4[i][bg_.tmod] = Hp_in_RK4[i][t_1] + dt *
// (
// 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau)
// - 800 * Hp_in_RK4[i][t_1]
// - 102400 * H_in_RK4[i][t]
// );
//// "improved" explicit Euler (mix current input with current H and former Hp)
//Hp_in_RK4[i][bg_.tmod] = Hp_in_RK4[i][t_1] + dt *
// (
// 102400 * nu_in[i] * N_in[i]->get_S((t - T_in[i]) % bg_.max_tau)
// - 800 * Hp_in_RK4[i][t_1]
// - 102400 * H_in_RK4[i][t]
// );
//// Heun (explicit trapezoidal)
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] + 0.5 * dt *
// (
// (102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau) - 800 *
// Hp_in_RK4[i][t_1]
// - 102400 *
// H_in_RK4[i][t_1])
// +
// (102400 * nu_in[i] * N_in[i]->get_S((t - T_in[i]) % bg_.max_tau) - 800 *
// (
// Hp_in_RK4[i][t_1] + dt *
// (102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau) - 800 * Hp_in_RK4[i][t_1] - 102400 * H_in_RK4[i][t_1])
// )
// - 102400 *
// H_in_RK4[i][t_1])
// );
////bsf Heun (explicit trapezoidal) with better estimation of H (needs prior computation of H)
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] + 0.5 * dt *
// (
// (102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau) - 800 *
// Hp_in_RK4[i][t_1]
// - 102400 *
// H_in_RK4[i][t_1])
// +
// (102400 * nu_in[i] * N_in[i]->get_S((t - T_in[i]) % bg_.max_tau) - 800 *
// (
// Hp_in_RK4[i][t_1] + dt *
// (102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i]) % bg_.max_tau) - 800 * Hp_in_RK4[i][t_1] - 102400 * H_in_RK4[i][t_1])
// )
// - 102400 *
// H_in_RK4[i][t])
// );
//// RK4 with better estimation of Hp (needs prior computation of Hp)
//k1 = 102400 * nu_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2])
// - 102400 * H_in_RK4[i][t_2];
//k2 = 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2] + (dt * k1))
// - 102400 * H_in_RK4[i][t_1];
//k3 = 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2] + (dt * k2))
// - 102400 * H_in_RK4[i][t_1];
//k4 = 102400 * nu_in[i] * N_in[i]->get_S((t - T_in[i]) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2] + 2 * dt * k3)
// - 102400 * H_in_RK4[i][t];
//Hp_in_RK4[i][bg_.tmod] = Hp_in_RK4[i][t_2] + dt * (k1 + 2*k2 + 2*k3 + k4)/3;
//// 5-steps Adams-Bashforth
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] + dt * (
// (1901./720.)*
// (102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_1])
// - 102400 * H_in_RK4[i][t_1])
// - (1387./360.)*
// (102400 * nu_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2])
// - 102400 * H_in_RK4[i][t_2])
// + (109./30.)*
// (102400 * nu_in[i] * N_in[i]->get_S((t-3 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_3])
// - 102400 * H_in_RK4[i][t_3])
// - (637./360.)*
// (102400 * nu_in[i] * N_in[i]->get_S((t-4 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_4])
// - 102400 * H_in_RK4[i][t_4])
// + (251./720.)*
// (102400 * nu_in[i] * N_in[i]->get_S((t-5 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_5])
// - 102400 * H_in_RK4[i][t_5])
// );
//// Dormand Prince (approximated)
//k1 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-90 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90])
// - 102400 * H_in_RK4[i][t_90]
// );
//k2 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-72 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90] + 1./5. * k1)
// - 102400 * H_in_RK4[i][t_72]
// );
//k3 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-63 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90] + 3./40. * k1 + 9./40. * k2)
// - 102400 * H_in_RK4[i][t_63]
// );
//k4 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-18 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90] + 44./45. * k1 - 56./15. * k2 + 32./9. * k3)
// - 102400 * H_in_RK4[i][t_18]
// );
//k5 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-10 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90] + 19372./6561. * k1 - 25360./2187. * k2 + 64448./6561. * k3 - 212./729. * k4)
// - 102400 * H_in_RK4[i][t_10]
// );
//k6 = dt90 * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_90] + 9017./3168. * k1 - 355./33. * k2 + 46732./5247. * k3 + 49./176. * k4 - 5103./18656. * k5)
// - 102400 * (H_in_RK4[i][t-1])
// );
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_90] + 35./384. * k1 + 500./1113. * k3 + 125./192. * k4 - 2187./6784. * k5 + 11./84. * k6;
// Dormand Prince
k1 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t-90 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90])
- 102400 * H_in_RK4[i][t_90]
);
k2 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t-72 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90] + 1./5. * k1)
- 102400 * H_in_RK4[i][t_72]
);
k3 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t-63 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90] + 3./40. * k1 + 9./40. * k2)
- 102400 * H_in_RK4[i][t_63]
);
k4 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t-18 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90] + 44./45. * k1 - 56./15. * k2 + 32./9. * k3)
- 102400 * H_in_RK4[i][t_18]
);
k5 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t-10 - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90] + 19372./6561. * k1 - 25360./2187. * k2 + 64448./6561. * k3 - 212./729. * k4)
- 102400 * H_in_RK4[i][t_10]
);
k6 = dt90 * (
102400 * nu_in[i] * N_in[i]->get_S((t - T_in[i] + bg_.max_tau) % bg_.max_tau)
- 800 * (Hp_in_RK4[i][t_90] + 9017./3168. * k1 - 355./33. * k2 + 46732./5247. * k3 + 49./176. * k4 - 5103./18656. * k5)
- 102400 * (H_in_RK4[i][t])
);
Hp_in_RK4[i][t] =
Hp_in_RK4[i][t_90]
+ 35./384. * k1 + 500./1113. * k3 + 125./192. * k4 - 2187./6784. * k5 + 11./84. * k6;
//// Runge-Kutta integration scheme (as defined in http://mymathlib.webtrellis.net/c_source/diffeq/second_order/runge_kutta_2nd_order.c)
////k1 = h * f(x[i], y[i], yp[i])
//k1 = dt * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-2 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2])
// - 102400 * H_in_RK4[i][t_2]
// );
////k2 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i], yp[i]+k1/2)
//k2 = dt * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2] + k1/2.)
// - 102400 * (H_in_RK4[i][t_2] + dt * Hp_in_RK4[i][t_2])
// );
////k3 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i]+(h/4)*k1, yp[i]+k2/2)
//k3 = dt * (
// 102400 * nu_in[i] * N_in[i]->get_S((t-1 - T_in[i] + bg_.max_tau) % bg_.max_tau)
// - 800 * (Hp_in_RK4[i][t_2] + k2/2.)
// - 102400 * (H_in_RK4[i][t_2] + dt * Hp_in_RK4[i][t_2] + (dt/2.)*k1)
// );
////y[i+1] = y[i] + h * ( yp[i] + 1/6 (k1 + k2 + k3 )
//H_in_RK4[i][t] = H_in_RK4[i][t_2] + dt * (2. * Hp_in_RK4[i][t_2] + k1 + k2 + k3)/3.;
////yp[i+1] = yp[i] + 1/6 * ( k1 + 2 k2 + 2 k3 + k4 )
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_2] + (k1 + 2. * k2 + 2. * k3 + k4)/6.;
//// Newmark integration scheme
//Hpp_in_RK4[i][t] = (102400 * nu_in[i] * N_in[i]->get_S(t) - 800 * Hp_in_RK4[i][t_1] - 102400 * H_in_RK4[i][t_1]);
//H_in_RK4[i][t] = H_in_RK4[i][t_1] + dt * Hp_in_RK4[i][t_1] + ((dt * dt)/2) * ((1-2*beta)*Hpp_in_RK4[i][t_1] + 2*beta*Hpp_in_RK4[i][t]);
//Hp_in_RK4[i][t] = Hp_in_RK4[i][t_1] + dt * ((1-gamma)*Hpp_in_RK4[i][t_1] + gamma*Hpp_in_RK4[i][t]);
sum_v += H_in_RK4[i][t];
}
S[(bg_.tmod) % bg_.max_tau] = Smax / (1 + exp(kvh - k*sum_v));
}
void BCBG2::SingleChannelNucleus::display(std::ostream& os)
{
os << "self-description of " << get_name() << std::endl;
os << " Smax = " << Smax << std::endl;
os << " vh = " << vh << std::endl;
os << " k = " << k << std::endl;
for (size_t i=0;i<n;i++) {
os << " afferent n° " << i << " is from " << N_in[i]->get_name() << std::endl;
os << " A = " << A_in[i] << std::endl;
os << " D = " << D_in[i] << std::endl;
os << " T = " << T_in[i] << std::endl;
os << " C = " << C_in[i] << std::endl;
os << " nu = " << nu_in[i] << std::endl;
os << " dtADC = " << dtADC_in[i] << std::endl;
}
}
void BCBG2::SingleChannelNucleus::save() {
int i,j;
S_backup.resize(S.size());
for (i=0;i<S.size(); i++) {
S_backup[i] = S[i];
}
H_in_RK4_backup.resize(H_in_RK4.size());
for (i=0;i<H_in_RK4.size(); i++) {
H_in_RK4_backup[i].resize(H_in_RK4[i].size());
for (j=0;j<H_in_RK4[i].size(); j++) {
H_in_RK4_backup[i][j] = H_in_RK4[i][j];
}
}
Hp_in_RK4_backup.resize(Hp_in_RK4.size());
for (i=0;i<Hp_in_RK4.size(); i++) {
Hp_in_RK4_backup[i].resize(Hp_in_RK4[i].size());
for (j=0;j<Hp_in_RK4[i].size(); j++) {
Hp_in_RK4_backup[i][j] = Hp_in_RK4[i][j];
}
}
}
void BCBG2::SingleChannelNucleus::save(MemorySCN& mem) {
int i,j;
mem.S_backup.resize(S.size());
for (i=0;i<S.size(); i++) {
mem.S_backup[i] = S[i];
}
mem.H_in_RK4_backup.resize(H_in_RK4.size());
for (i=0;i<H_in_RK4.size(); i++) {
mem.H_in_RK4_backup[i].resize(H_in_RK4[i].size());
for (j=0;j<H_in_RK4[i].size(); j++) {
mem.H_in_RK4_backup[i][j] = H_in_RK4[i][j];
}
}
mem.Hp_in_RK4_backup.resize(Hp_in_RK4.size());
for (i=0;i<Hp_in_RK4.size(); i++) {
mem.Hp_in_RK4_backup[i].resize(Hp_in_RK4[i].size());
for (j=0;j<Hp_in_RK4[i].size(); j++) {
mem.Hp_in_RK4_backup[i][j] = Hp_in_RK4[i][j];
}
}
}
void BCBG2::SingleChannelNucleus::load() {
int i,j;
S.resize(S_backup.size());
for (i=0;i<S.size(); i++) {
S[i] = S_backup[i];
}
H_in_RK4.resize(H_in_RK4_backup.size());
for (i=0;i<H_in_RK4_backup.size(); i++) {
H_in_RK4[i].resize(H_in_RK4_backup[i].size());
for (j=0;j<H_in_RK4_backup[i].size(); j++) {
H_in_RK4[i][j] = H_in_RK4_backup[i][j];
}
}
Hp_in_RK4.resize(Hp_in_RK4_backup.size());
for (i=0;i<Hp_in_RK4_backup.size(); i++) {
Hp_in_RK4[i].resize(Hp_in_RK4_backup[i].size());
for (j=0;j<Hp_in_RK4_backup[i].size(); j++) {
Hp_in_RK4[i][j] = Hp_in_RK4_backup[i][j];
}
}
}
void BCBG2::SingleChannelNucleus::load(MemorySCN& mem) {
int i,j;
S.resize(mem.S_backup.size());
for (i=0;i<S.size(); i++) {
S[i] = mem.S_backup[i];
}
H_in_RK4.resize(mem.H_in_RK4_backup.size());
for (i=0;i<mem.H_in_RK4_backup.size(); i++) {
H_in_RK4[i].resize(mem.H_in_RK4_backup[i].size());
for (j=0;j<mem.H_in_RK4_backup[i].size(); j++) {
H_in_RK4[i][j] = mem.H_in_RK4_backup[i][j];
}
}
Hp_in_RK4.resize(mem.Hp_in_RK4_backup.size());
for (i=0;i<mem.Hp_in_RK4_backup.size(); i++) {
Hp_in_RK4[i].resize(mem.Hp_in_RK4_backup[i].size());
for (j=0;j<mem.Hp_in_RK4_backup[i].size(); j++) {
Hp_in_RK4[i][j] = mem.Hp_in_RK4_backup[i][j];
}
}
}