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 i_stim = 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;
255 #else // ifdef HAVE_DLOPEN 256 log_msg(NULL, 5, 0,
"Compile with USE_DLOPEN to support dynamic module loading.");
258 #endif // ifdef HAVE_DLOPEN 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 i_stim = (params.stim_volt_arg-d[0])/params.resistance_arg;
609 *(cMIIF->
gdata[Vm]) += i_stim;
610 #endif // ifndef ALGEBRAIC 614 for (
auto ion : species_list)
616 cMIIF->
transmem_stim_species( -i_stim*(ratios_list[iIon]/ratio_sum), ion.c_str(), params.surface_to_volume_arg, stim_list, numNode);
617 FPRINTF( WORLD stderr,
"[stim_assign] @ t = %.2f ms; magnitude: %fpA/pF; type: %s\n", tmo.
time, i_stim, ion.c_str() );
623 if ( MIIF.
gdata[illum] != NULL ) {
627 FPRINTF( WORLD stderr,
"[stim] illum ON @ t = %.2f ms; E_e = %.3f mW/mm^2\n", tmo.
time, light_irrad );
630 FPRINTF( WORLD stderr,
"[stim] illum OFF @ t = %.2f 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());
int trigger_elapse(int ID) const
int main(int argc, char *argv[])
int read_sv(MULTI_IF *, int, const char *)
void update_timing(event_timing *t, double event_duration)
void free_trace(trace *tr)
bool pmat
indiacte premature AP
void restitution_trigger_list(char *r_file, restitution *r, char *protocol, int *n_dop, double **t_dop)
std::string get_string_from_target(Target const target)
Get a string representation of a given target.
TrgList trigs
trigger list for defining stim sequence
virtual void release_ptr(S *&p)=0
void open_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, char *base_name, IOCtrl *io)
char * dupstr(const char *old_str)
virtual void add_scaled(const abstract_vector< T, S > &vec, S k)=0
void close_trace(MULTI_IF *MIIF)
void dump_all(MULTI_IF *MIIF, int reg, char *imp, char *plugs, double t, double ddt, char *fout)
double getCellVal(sf_vec *v, int ind)
int read_trace(trace *tr, const char *name)
SV_DUMP svd
state variable dump
void resample_trace(trace *tr, double dt)
std::vector< std::reference_wrapper< IonType > > IonTypeList
void sv_dump_add_by_name_list(int, char *, char *, char *, char *, char *, double, double)
void initialize_timings(event_timing *t)
FILE * rstats
output restitution statistics only
int * numplugs
number of plugins for each region
float determine_duration(struct gengetopt_args_info *p, TrgList *stim_lst)
determine time of last stimulus
double avg
average duration of event
void sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
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...
timing for initialization
void doppel_MIIF(MULTI_IF *orig, MULTI_IF *miif_doppel)
bool initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f, int trans, float *duration)
std::map< mesh_t, sf_mesh > mesh_reg
Registry for the different meshes used in a multi-physics simulation.
IIF_Mask_t * IIFmask
region for each node
void initialize_eq_timer(double istart, double iend, int ntrig, double iintv, double idur, int ID, const char *iname, const char *poolname=nullptr)
void dump_luts_MIIF(bool)
void print_AP_stats_header(action_potential *AP, FILE *outbuf)
virtual T lsize() const =0
timer_manager * tm_manager
a manager for the various physics timers
bool trigger_end(int ID) const
void save_sv(MULTI_IF *, int, const char *)
void print_param_help(IonType *im, IonTypeList &plugs)
timing for main loop (including IO and ODE solve)
IonTypeList iontypes
type for each region
char wbin
write to file in binary format
bool check_events(double vm, action_potential *AP, timer_manager *tm)
int get_plug_flag(char *plgstr, int *out_num_plugins, IonTypeList &out_plugins)
std::vector< IonTypeList > plugtypes
plugins types for each region
bool trigger(int ID) const
int N_IIF
how many different IIF's
void initialize_currents(double, int)
bool is_gpu(Target const target)
Checks if this is a GPU target.
TrgList saveState
instants at wich we save state vectors
void init_vector(SF::abstract_vector< T, S > **vec)
double mx
maximum duration of event
std::vector< IonIfBase * > IIF
array of IIF's
int n
number of pulses required for protocol
void apply_stretch(MULTI_IF *miif, stretch *s, timer_manager *tm)
virtual Target select_target(Target target) const =0
Gets a supported target from the given target.
int initialize_AP_analysis(action_potential *AP)
special value to handle unknown targets
double tot
total duration of all events
void doppel_update(MULTI_IF *orig, MULTI_IF *miif_doppel)
void dump_state(char *, float, opencarp::mesh_t gid, bool, unsigned int)
void initialize_singlestep_timer(double tg, double idur, int ID, const char *iname, const char *poolname=nullptr)
Define multiple ionic models to be used in different regions.
void initializePulseStretch(float strain, float onset, float duration, float rise, float fall, stretch *s)
IonType * get_ion_type(const std::string &name)
std::map< SF::quadruple< int >, SF::index_mapping< int > > map_reg
Registriy for the inter domain mappings.
double * lst
store instants of pulse delivery
void cleanup_AP_analysis(action_potential *AP)
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false)=0
double vm_trc[VM_HIST_LEN]
void free_doppel(MULTI_IF *m)
void determine_stim_list(char *stl, TrgList *trg, bool DIAs)
long d_time
current time instance index
std::string get_target_list_string()
Returns a string containing the list of available targets.
void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
void restitution_save_sv(MULTI_IF *miif, int R1, restitution *r, action_potential *AP)
void globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, timer_manager *tmo, IOCtrl *io, int numNode)
void get_time(double &tm)
FILE * stats
output statistics for each AP
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.
number of benchmark timings we use
V timing(V &t2, const V &t1)
Index mapping class. This is a bijective mapping.
int write_dump_header(GVEC_DUMP *gvd, SV_DUMP *svd, const char *ExpID)
V dot(const vec3< V > &p1, const vec3< V > &p2)
int dump_svs(opencarp::base_timer *)
void transmem_stim_species(float, const char *, float, int *, int)
opencarp::sf_vec * gdata[NUM_IMP_DATA_TYPES]
data used by all IMPs
Defines valid targets for an ionic model to run on and an allocator for allocating memory on a specif...
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Target get_target_from_string(std::string const str)
Returns a value from the Target enum from a given string.
void dump_trace(MULTI_IF *MIIF, limpet::Real time)
std::vector< Target > targets
target for each region
int numNode
local number of nodes
Target
enum that represents different targets to run ionic models on.
std::vector< base_timer * > timers
vector containing individual timers
double dur
total duration of protocol
bool triggered_now(int ID) const
void AP_clamp(Clamp *cl, timer_manager *tm, sf_vec *v, bool trigger)
void initialize_neq_timer(const std::vector< double > &itrig, double idur, int ID, const char *iname, const char *poolname=nullptr)
int load_ionic_module(const char *)
void close_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, IOCtrl *io)
double mn
minimum duration of event
Abstract class representing an ionic model type.
double SF_real
Use the general double as real type.
void initial_SVs(MULTI_IF *miif, char *SVs, char *imp, char *plgins, int num)
void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
bool extUpdateVm
flag indicating update function for Vm
Basic utility structs and functions, mostly IO related.
int process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
bool * pmat
store flag to indicate prematurity
int get_size(MPI_Comm comm=PETSC_COMM_WORLD)
manage input, output, resampling of traces
centralize time managment and output triggering
SF::scatter_registry scatter_reg
Registry for the different scatter objects.
The scatterer registry class.