33 #include <sys/resource.h>
56 #include "build_info.h"
58 #include "petsc_utils.h"
64 namespace user_globals {
84 #define __FUNCT__ "main"
86 int main(
int argc,
char *argv[]) {
91 const int num_region = 1;
106 struct gengetopt_args_info params;
112 FILE *fhdls[NUM_IMP_DATA_TYPES+1];
116 char *PetscDBfile = COMPAT_PETSC_NULLPTR;
117 char *help_msg = COMPAT_PETSC_NULLPTR;
118 initialize_PETSc(&argc, argv, PetscDBfile, help_msg);
119 COMPAT_PetscOptionsInsertString(
"-options_left no");
121 if (cmdline_parser(argc, argv, ¶ms) != 0)
125 double t1, t2, t3, t4, t5, t6, t7, t8;
136 double stim_charge = params.stim_curr_arg*dt;
137 double bcl = params.bcl_arg;
138 double stim_start = params.stim_start_arg;
139 double stim_dur = params.stim_dur_arg;
140 int nstim = params.numstim_arg ? params.numstim_arg : -1;
141 double dt_out = params.dt_out_arg;
142 double start_out = params.start_out_arg;
143 int numNode = params.num_arg;
144 int stim_assign = params.stim_assign_flag;
145 bool extVmUpdate = params.ext_vm_update_given;
146 bool analyze_AP = params.APstatistics_given || params.restitute_given;
147 bool restitute = params.restitute_given;
148 bool validate = params.validate_given;
149 bool APclamp = params.AP_clamp_file_given;
150 bool do_trace = restitute ? params.res_trace_flag : (!params.no_trace_flag);
154 io.
w2file = validate ? 1 : params.fout_given;
155 io.
wbin = validate ? 1 : params.bin_flag;
160 log_msg(NULL, 3, 0,
"Warning: Number of nodes is less than number of processes used!");
161 log_msg(NULL, 3, 0,
"Setting number of nodes to %d.", numNode);
166 bool tr_stim =
false;
167 if (params.stim_file_given) {
168 char *fname =
dupstr(params.stim_file_arg);
172 stim_dur = stim_trace.
dur;
173 if (bcl < stim_trace.
dur)
174 log_msg(NULL, 4, 0,
"BCL less than stimulus duration!");
179 bool tr_light =
false;
180 if(params.light_file_given) {
181 char* fname =
dupstr(params.light_file_arg);
186 int n_clipped_to_0 = 0;
187 for(
int i=0; i<light_trace.
N; i++) {
188 if( light_trace.
s[i] < 0.0 ) {
189 light_trace.
s[i] = 0.0;
193 if( n_clipped_to_0 ) {
194 FPRINTF(WORLD stderr,
"WARNING: %d irradiance values in %s were negative (clipped to 0.0)\n", n_clipped_to_0, fname);
200 int SVclamp =
process_sv_clamps(params.clamp_SVs_arg, params.SV_clamp_files_arg, &sv_cl, dt);
211 bool IsVclamp =
initialize_clamp(&cl, params.clamp_arg, params.clamp_ini_arg,
212 params.clamp_start_arg, params.clamp_dur_arg,
213 params.clamp_file_arg,
217 log_msg(NULL, 0, 0,
"info: cannot clamp voltage and have AP clamp at same time");
221 log_msg(NULL, 0, 0,
"info: Vm will be clamped to %f mV over the interval [%f, %f]",
222 params.clamp_arg, params.clamp_start_arg, params.clamp_start_arg + params.clamp_dur_arg);
224 log_msg(NULL, 0, 0,
"info: Vm will be clamped from %f to %f mV over the interval [%f, %f]",
225 params.clamp_ini_arg, params.clamp_arg, params.clamp_start_arg,
226 params.clamp_start_arg + params.clamp_dur_arg);
233 Regions =
static_cast<char *
>(calloc(numNode,
sizeof(
char) ));
238 int *stim_list =
static_cast<int *
>(malloc(numNode*
sizeof(
int) ));
239 for (
int i = 0; i < numNode; i++) {
244 for (
unsigned int i = 0; i < params.load_module_given; i++) {
247 if (!params.imp_given) {
249 char *imp = strdup(basename(params.load_module_arg[i]) );
250 char *
dot = strrchr(imp,
'.');
253 params.imp_arg = imp;
256 log_msg(NULL, 5, 0,
"Compile with USE_DLOPEN to support dynamic module loading.");
261 if (
get_ion_type(std::string(params.imp_arg)) == NULL) {
262 fprintf(stderr,
"Illegal IMP specified: %s\n", params.imp_arg);
267 if (!
get_plug_flag(params.plug_in_arg, &num_plugins, RegionPlug)) {
268 fprintf(stderr,
"Illegal plugin specified: %s\n", params.plug_in_arg);
279 fprintf(stderr,
"Unkown target: %s\n", params.target_arg);
285 fprintf(stderr,
"Model %s was not generated for target %s\n",
290 if (!RegionPlug.empty()) {
292 fprintf(stderr,
"Plugin %s was not generated for target %s\n",
302 MIIF.
N_IIF = num_region;
309 MIIF.
gdata[Vm] = Vmv;
310 MIIF.
gdata[Iion] = I_ion;
315 double setup_time =
timing(t2, t1);
320 if (params.list_imps_flag || params.plugin_outputs_flag) {
325 if (params.imp_info_flag) {
330 fprintf(stderr,
"\n*** GIT tag: %s\n", GIT_COMMIT_TAG);
331 fprintf(stderr,
"*** GIT hash: %s\n", GIT_COMMIT_HASH);
332 fprintf(stderr,
"*** GIT repo: %s\n", GIT_PATH);
333 fprintf(stderr,
"*** dependency commits: %s\n\n", SUBREPO_COMMITS);
335 if (params.buildinfo_flag) {
336 fprintf(stderr,
"\n*** -buildinfo flag is deprecated and will be removed.\n");
337 fprintf( stderr,
"\n*** The build information is printed by default. \n\n");
341 MIIF.
IIF[0]->tune(params.imp_par_arg, params.plug_in_arg, params.plug_par_arg);
352 char ap_stats_fname[1024];
353 snprintf(ap_stats_fname,
sizeof ap_stats_fname,
"%s_AP_stats.dat", params.fout_arg);
354 AP.
stats = fopen(ap_stats_fname,
"wt");
362 char *res_file =
dupstr(params.res_file_arg);
364 if (r.
dur > duration) duration = r.
dur;
367 char ap_rstats_fname[1024];
368 snprintf(ap_rstats_fname,
sizeof ap_rstats_fname,
"%s_APD_restitution.dat", params.fout_arg);
369 AP.
rstats = fopen(ap_rstats_fname,
"wt");
379 double light_irrad = params.light_irrad_arg;
380 double light_start = params.light_start_arg;
381 double light_end = params.duration_arg;
382 int light_nstim = params.light_numstim_arg;
383 double light_bcl = params.light_bcl_arg;
384 double light_dur = params.light_dur_arg;
386 if(light_irrad < 0) {
387 FPRINTF(WORLD stderr,
"WARNING: irradiance value %.3f makes no sense (clipped to 0.0)\n", light_irrad);
391 if(params.light_times_given)
394 if (params.SV_init_given)
395 initial_SVs(&MIIF, params.SV_init_arg, params.imp_arg, params.plug_in_arg, numNode);
403 if (params.save_time_arg)
406 if (params.save_ini_time_given)
410 if (params.restitute_given) {
411 std::vector<double> trg;
415 if (params.res_state_vector_given) {
420 if (!strcmp(params.restitute_arg,
"S1S2f") ) {
421 trg.assign(t_dopple, t_dopple + n_dopple);
425 else if (params.stim_times_given) {
426 std::vector<double> trg;
427 trg.assign(stim_lst.
lst, stim_lst.
lst + stim_lst.
n);
435 if (MIIF.
gdata[illum] != NULL && !tr_light) {
436 if (params.light_times_given) {
437 std::vector<double> trg;
438 trg.assign(light_lst.
lst, light_lst.
lst + light_lst.
n);
442 light_bcl, light_dur,
LIGHT_TM_IDX,
"LIGHT_TIMER", NULL);
450 if (!restitute && params.doppel_on_given) {
455 if (params.read_ini_file_given) {
456 if (
read_sv(&MIIF, R1, params.read_ini_file_arg) != 0) {
463 memset(&s, 0,
sizeof(
stretch));
466 params.strain_dur_arg, params.strain_rate_arg,
467 params.strain_rate_arg, &s);
470 MIIF.
gdata[Lambda] = lambdavec;
473 for (
int i = 0; i < lambdavec->
lsize(); i++)
480 std::stringstream stim_species(params.stim_species_arg);
481 std::stringstream stim_ratios(params.stim_ratios_arg);
482 std::string specie, ratio;
483 std::vector<std::string> species_list;
484 std::vector<float> ratios_list;
488 while(std::getline(stim_species, specie,
':'))
490 species_list.push_back(specie);
492 while(std::getline(stim_ratios, ratio,
':'))
494 ratios_list.push_back(std::stof(ratio));
496 if (species_list.size() != ratios_list.size()) {
497 fprintf(stderr,
"Number of ion species for which stimulus should be assigned does not match number of ratio entries: %s | %s\n", params.stim_species_arg, params.stim_ratios_arg);
500 for (
auto i : ratios_list)
502 if (std::abs(ratio_sum - 1) >= 0.001)
503 log_msg(NULL, 3, 0,
"Warning: Sum of ion species ratios for stimulus assignement != 1. Normalizing to sum %.1f.", ratio_sum);
507 char sv_dump_names[1024];
508 if (!params.fout_given && !params.validate_flag)
509 snprintf(sv_dump_names,
sizeof sv_dump_names,
"%s_%s", params.fout_arg, params.imp_arg);
511 strcpy(sv_dump_names, params.fout_arg);
513 if (params.validate_flag) {
514 dump_all(&MIIF, R1, params.imp_arg, params.plug_in_arg, t, dt, sv_dump_names);
517 params.imp_sv_dump_arg, params.plug_in_arg,
518 params.plug_sv_dump_arg, t, dt_out);
522 if (params.dump_lut_flag)
526 int trace_nodes[] = {0};
528 open_trace(&MIIF,
sizeof(trace_nodes)/
sizeof(trace_nodes[0]), trace_nodes, ¶ms.trace_no_arg, NULL);
534 double ini_time =
timing(t3, t2);
538 int gpu_Vm_update = 1;
539 double counter = 0.0;
540 double epsilon = 0.0000001f;
553 printf(
"[%f] Switching to new ion structure\n", tmo.
time);
569 for (
int i = 0; i < SVclamp; i++) {
571 if (params.SV_I_trigger_flag)
576 sv_clamp(sv_cl+i, &tmo, cMIIF, trigger);
589 save_sv(cMIIF, R1, params.save_ini_file_arg);
601 else if (params.stim_volt_given) {
604 stim_charge = (params.stim_volt_arg-d[0])/params.resistance_arg*dt;
609 *(cMIIF->
gdata[Vm]) += stim_charge;
614 for (
auto ion : species_list)
616 cMIIF->
transmem_stim_species(-stim_charge*(ratios_list[iIon]/ratio_sum), ion.c_str(), params.surface_to_volume_arg, stim_list, numNode);
617 FPRINTF(WORLD stderr,
"[stim_assign] @ t = %.3f ms; magnitude: %f pA/pF; type: %s\n", tmo.
time, stim_charge, ion.c_str());
623 if (MIIF.
gdata[illum] != NULL) {
627 FPRINTF(WORLD stderr,
"[stim] illum ON @ t = %.3f ms; E_e = %.3f mW/mm^2\n", tmo.
time, light_irrad);
630 FPRINTF(WORLD stderr,
"[stim] illum OFF @ t = %.3f ms\n", tmo.
time);
647 if (gpu_Vm_update == 0) {
649 }
else if (gpu_Vm_update == 1) {
650 cMIIF->
compute_ionic_current(fabs(counter - 0.0) < epsilon || fabs(counter - 1.0) < epsilon , fabs((counter+dt) - 1.0) < epsilon);
661 double ode_time =
timing(t6, t5);
669 log_msg(NULL, 0, 0,
"[%f] Switching back to original ion structure\n", tmo.
time);
675 double loop_time =
timing(t7, t4);
679 counter = counter + dt;
680 if ( fabs(counter - 1.0) < epsilon ) {
686 double main_time =
timing(t8, t3);
689 const char *ExpID = params.fout_arg;
702 log_msg(NULL, 0, 0,
"\n\n\nAll done!\n\n");
706 log_msg(NULL, 0, 0,
"setup time %.6f s", (
double)t_setup->
tot);
708 log_msg(NULL, 0, 0,
"initialization time %.6f s", (
double)t_init->
tot);
710 log_msg(NULL, 0, 0,
"main loop time %.6f s", (
double)t_loop->
tot);
712 log_msg(NULL, 0, 0,
"total ode time %.6f s\n", (
double)t_ode->
tot);
714 log_msg(NULL, 0, 0,
"mn/avg/mx loop time %.6f %.6f %.6f ms",
715 (
double)t_loop->
mn*1e3, (
double)t_loop->
avg*1e3, (
double)t_loop->
mx*1e3);
716 log_msg(NULL, 0, 0,
"mn/avg/mx ODE time %.6f %.6f %.6f ms",
717 (
double)t_ode->
mn*1e3, (
double)t_ode->
avg*1e3, (
double)t_ode->
mx*1e3);
719 log_msg(NULL, 0, 0,
"real time factor %.6f\n\n", duration/(t_ode->
tot*1e3) );
722 CHKERRQ(PetscFinalize());
Define multiple ionic models to be used in different regions.
double SF_real
Use the general double as real type.
Basic utility structs and functions, mostly IO related.
int main(int argc, char *argv[])
virtual void release_ptr(S *&p)=0
virtual void add_scaled(const abstract_vector< T, S > &vec, S k)=0
virtual T lsize() const =0
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false, const bool local=false)=0
Index mapping class. This is a bijective mapping.
The scatterer registry class.
Abstract class representing an ionic model type.
virtual Target select_target(Target target) const =0
Gets a supported target from the given target.
int numNode
local number of nodes
bool extUpdateVm
flag indicating update function for Vm
int dump_svs(opencarp::base_timer *)
std::vector< IonIfBase * > IIF
array of IIF's
opencarp::sf_vec * gdata[NUM_IMP_DATA_TYPES]
data used by all IMPs
void sv_dump_add_by_name_list(int, char *, char *, char *, char *, char *, double, double)
int * numplugs
number of plugins for each region
std::vector< Target > targets
target for each region
std::vector< IonTypeList > plugtypes
plugins types for each region
IonTypeList iontypes
type for each region
void dump_state(char *, float, opencarp::mesh_t gid, bool, unsigned int)
SV_DUMP svd
state variable dump
void transmem_stim_species(float, const char *, float, int *, int)
void initialize_currents(double, int)
int N_IIF
how many different IIF's
void compute_ionic_current(bool flag_send=1, bool flag_receive=1)
GPU kernel to emulate the add_scaled call made to adjust the Vm values when the update to Vm is not m...
void dump_luts_MIIF(bool)
IIF_Mask_t * IIFmask
region for each node
centralize time managment and output triggering
void initialize_neq_timer(const std::vector< double > &itrig, double idur, int ID, const char *iname, const char *poolname=nullptr)
long d_time
current time instance index
bool trigger(int ID) const
void initialize_eq_timer(double istart, double iend, int ntrig, double iintv, double idur, int ID, const char *iname, const char *poolname=nullptr)
void initialize_singlestep_timer(double tg, double idur, int ID, const char *iname, const char *poolname=nullptr)
int trigger_elapse(int ID) const
bool trigger_end(int ID) const
std::vector< base_timer * > timers
vector containing individual timers
bool triggered_now(int ID) const
void init_vector(SF::abstract_vector< T, S > **vec)
void initial_SVs(MULTI_IF *miif, char *SVs, char *imp, char *plgins, int num)
void doppel_MIIF(MULTI_IF *orig, MULTI_IF *miif_doppel)
int write_dump_header(GVEC_DUMP *gvd, SV_DUMP *svd, const char *ExpID)
int get_plug_flag(char *plgstr, int *out_num_plugins, IonTypeList &out_plugins)
void doppel_update(MULTI_IF *orig, MULTI_IF *miif_doppel)
void globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, timer_manager *tmo, IOCtrl *io, int numNode)
void resample_trace(trace *tr, double dt)
int load_ionic_module(const char *)
void free_trace(trace *tr)
void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
Target
enum that represents different targets to run ionic models on.
@ UNKNOWN
special value to handle unknown targets
void sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
@ LOOP_IDX
timing for main loop (including IO and ODE solve)
@ ODE_IDX
timing for ODE solve
@ INIT_IDX
timing for initialization
@ N_TIMINGS
number of benchmark timings we use
@ SETUP_IDX
timing for setup phase
std::string get_string_from_target(Target const target)
Get a string representation of a given target.
void dump_all(MULTI_IF *MIIF, int reg, char *imp, char *plugs, double t, double ddt, char *fout)
void close_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, IOCtrl *io)
IonType * get_ion_type(const std::string &name)
int initialize_AP_analysis(action_potential *AP)
void restitution_trigger_list(char *r_file, restitution *r, char *protocol, int *n_dop, double **t_dop)
void save_sv(MULTI_IF *, int, const char *)
double getCellVal(sf_vec *v, int ind)
std::string get_target_list_string()
Returns a string containing the list of available targets.
void apply_stretch(MULTI_IF *miif, stretch *s, timer_manager *tm)
bool is_gpu(Target const target)
Checks if this is a GPU target.
void restitution_save_sv(MULTI_IF *miif, int R1, restitution *r, action_potential *AP)
void print_AP_stats_header(action_potential *AP, FILE *outbuf)
void initialize_timings(event_timing *t)
bool initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f, int trans, float *duration)
void initializePulseStretch(float strain, float onset, float duration, float rise, float fall, stretch *s)
void determine_stim_list(char *stl, TrgList *trg, bool DIAs)
void cleanup_AP_analysis(action_potential *AP)
std::vector< std::reference_wrapper< IonType > > IonTypeList
float determine_duration(struct gengetopt_args_info *p, TrgList *stim_lst)
determine time of last stimulus
int process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
int read_trace(trace *tr, const char *name)
void open_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, char *base_name, IOCtrl *io)
int read_sv(MULTI_IF *, int, const char *)
void close_trace(MULTI_IF *MIIF)
void dump_trace(MULTI_IF *MIIF, limpet::Real time)
void update_timing(event_timing *t, double event_duration)
void open_trace(MULTI_IF *MIIF, int n_traceNodes, int *traceNodes, int *label, opencarp::sf_mesh *imesh)
Set up ionic model traces at some global node numbers.
void free_doppel(MULTI_IF *m)
void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
bool check_events(double vm, action_potential *AP, timer_manager *tm)
void print_param_help(IonType *im, IonTypeList &plugs)
Target get_target_from_string(std::string const str)
Returns a value from the Target enum from a given string.
void AP_clamp(Clamp *cl, timer_manager *tm, sf_vec *v, bool trigger)
timer_manager * tm_manager
a manager for the various physics timers
std::map< SF::quadruple< int >, SF::index_mapping< int > > map_reg
Registriy for the inter domain mappings.
SF::scatter_registry scatter_reg
Registry for the different scatter objects.
std::map< mesh_t, sf_mesh > mesh_reg
Registry for the different meshes used in a multi-physics simulation.
V dot(const vec3< V > &p1, const vec3< V > &p2)
char * dupstr(const char *old_str)
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
void get_time(double &tm)
int get_size(MPI_Comm comm=PETSC_COMM_WORLD)
V timing(V &t2, const V &t1)
char wbin
write to file in binary format
int n
number of pulses required for protocol
double * lst
store instants of pulse delivery
bool * pmat
store flag to indicate prematurity
bool pmat
indiacte premature AP
FILE * rstats
output restitution statistics only
FILE * stats
output statistics for each AP
double vm_trc[VM_HIST_LEN]
double avg
average duration of event
double mx
maximum duration of event
double mn
minimum duration of event
double tot
total duration of all events
TrgList saveState
instants at wich we save state vectors
double dur
total duration of protocol
TrgList trigs
trigger list for defining stim sequence
manage input, output, resampling of traces
Defines valid targets for an ionic model to run on and an allocator for allocating memory on a specif...