bench.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#define PrMGLOBAL
#include <stdlib.h>
#include <stdio.h>
#include "cmdline.h"
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <limits.h>
#include <math.h>
#include <stdbool.h>
#include <vector>
#include <sstream>
#include <fstream>
#include <string>
#include <cmath>
#include "basics.h"
 
#include "MULTI_ION_IF.h"
 
#include "trace.h"
#include "clamp.h"
#include "ap_analyzer.h"
#include "restitute.h"
#include "stretch.h"
#include "sv_init.h"
#include "bench_utils.h"
#include "target.h"
#include "Lapicque.h"
 
#include "libgen.h"
#include "build_info.h"
 
#include "petsc_utils.h" //TODO for EXIT, FPRINT*
#include "SF_init.h" // for SF::init_xxx()
 
 
// globals
namespace opencarp {
namespace user_globals {
  SF::scatter_registry scatter_reg;
  std::map<mesh_t, sf_mesh> mesh_reg;
  std::map<SF::quadruple<int>, SF::index_mapping<mesh_int_t> > map_reg;
  // /// the physics
  // std::map<physic_t, Basic_physic*> physics_reg;
  timer_manager* tm_manager;
  // /// important solution vectors from different physics
  // std::map<datavec_t, sf_vec*> datavec_reg;
}  // namespace user_globals
}  // namespace opencarp
 
using namespace opencarp;
using namespace limpet;
 
#undef __FUNCT__
#define __FUNCT__ "main"
 
int main(int argc, char *argv[]) {
 
  // specify the number of regions and their names
  const int num_region = 1;
  enum { R1 }
  Region;
 
  // specifiy the ionic model for each region
  IonType* RegionDef;
 
  // specify the plugins for each region
  int num_plugins;
  IonTypeList RegionPlug = {};
 
  // specify the region for each node
  char *Regions;
 
 
  struct gengetopt_args_info params;
 
  // we allocate the vectors on the heap since they will be freed by the MIIF
  sf_vec* Vmv;
  sf_vec* I_ion;
  
  FILE *fhdls[NUM_IMP_DATA_TYPES+1];
  GVEC_DUMP gvd;
  IOCtrl io;
 
  char *PetscDBfile = COMPAT_PETSC_NULLPTR;  // could be used to pass in ~/.petscrc
  char *help_msg    = COMPAT_PETSC_NULLPTR;
  initialize_PETSc(&argc, argv, PetscDBfile, help_msg);
  COMPAT_PetscOptionsInsertString("-options_left no");
 
  if (cmdline_parser(argc, argv, &params) != 0)
    exit(1);
 
  // performance measurements
  double t1, t2, t3, t4, t5, t6, t7, t8;
  struct timeval tt1;
 
  event_timing timings[N_TIMINGS];
  initialize_timings(timings);
 
  // start measurementstrace_IA
  get_time(t1);
 
  // determine mode
  bool   RRC          = params.RRC_mode_counter;
  bool   Lapicque_crv = params.StimThresh_mode_counter;
  bool   restitute    = params.restitute_given;
 
  // assign command line params and dependent variables
  double dt           = params.dt_arg;
  double bcl          = Lapicque_crv ? params.Lap_bcl_arg : params.bcl_arg;
  double stim_start   = Lapicque_crv ? 0. : params.stim_start_arg;
  int    nstim        = Lapicque_crv ? params.Lap_numstim_arg : params.bcl_arg;
  double stim_dur     = params.stim_dur_arg;
  double stim_charge  = params.stim_curr_arg*dt;
  double dt_out       = params.dt_out_arg;
  double start_out    = params.start_out_arg;
  int    numNode      = params.num_arg;
  int    stim_assign  = params.stim_assign_flag;
  bool   analyze_AP   = params.APstatistics_given || params.restitute_given;
  bool   validate     = params.validate_given;
  bool   APclamp      = params.AP_clamp_file_given;
  bool   do_trace     = restitute ? params.res_trace_flag : (!params.no_trace_flag);
  TrgList stim_lst;
  float   duration = determine_duration(&params, &stim_lst);
 
  nstim = nstim ? nstim : -1;
  io.w2file = validate ? 1 : params.fout_given;
  io.wbin   = validate ? 1 : params.bin_flag;
 
  if(get_size() > numNode) {
    numNode = get_size();
 
    log_msg(NULL, 3, 0, "Warning: Number of nodes is less than number of processes used!");
    log_msg(NULL, 3, 0, "Setting number of nodes to %d.", numNode);
  }
 
  // stimulate with trace from file
  trace stim_trace;
  bool tr_stim = false;
  if (params.stim_file_given) {
    char *fname = dupstr(params.stim_file_arg);
    read_trace(&stim_trace, fname);
    resample_trace(&stim_trace, dt);
    tr_stim  = true;
    stim_dur = stim_trace.dur;
    if (bcl < stim_trace.dur)
      log_msg(NULL, 4, 0, "BCL less than stimulus duration!");
  }
#define RRC_MAX_N 50     
  if( RRC ) {
    mk_RRC_trace(&stim_trace, params.stim_curr_arg,params.RRC_delay_arg,
            -1., params.RRC_clamp_dur_arg );
    resample_trace(&stim_trace, dt);
    tr_stim  = true;
    stim_dur = stim_trace.dur;
    nstim    = RRC_MAX_N;
    bcl      = stim_dur+10.;
    duration = nstim*bcl;
  }
 
  // illuminate with trace from file
  trace light_trace;
  bool tr_light = false;
  if(params.light_file_given) {
    char* fname = dupstr(params.light_file_arg);
    read_trace(&light_trace, fname);
    resample_trace(&light_trace, dt);
    tr_light = true;
 
    int n_clipped_to_0 = 0;
    for(int i=0; i<light_trace.N; i++) {
      if( light_trace.s[i] < 0.0 ) {
        light_trace.s[i] = 0.0;
        n_clipped_to_0++;
      }
    }
    if( n_clipped_to_0 ) {
      FPRINTF(WORLD stderr, "WARNING: %d irradiance values in %s were negative (clipped to 0.0)\n", n_clipped_to_0, fname);
    }
  }
 
  // state variable clamp
  Clamp *sv_cl;
  int SVclamp = process_sv_clamps(params.clamp_SVs_arg, params.SV_clamp_files_arg, &sv_cl, dt);
 
  // action potential clamp
  Clamp ap_cl;
  if (APclamp){
    const char *Vm_clamp = "Vm";
    initialize_sv_clamp(&ap_cl, Vm_clamp, params.AP_clamp_file_arg, dt);
  }
 
  // voltage clamp experiment
  Clamp cl;
  bool IsVclamp = initialize_clamp(&cl, params.clamp_arg, params.clamp_ini_arg,
                                    params.clamp_start_arg, params.clamp_dur_arg,
                                    params.clamp_file_arg,
                                    !params.clamp_ini_given ? CLAMP_TM_IDX : 0, &duration);
  if (IsVclamp) {
    if (APclamp) {
      log_msg(NULL, 0, 0,  "info: cannot clamp voltage and have AP clamp at same time");
      exit(1);
    }
    if (cl.transient) {
      log_msg(NULL, 0, 0,  "info: Vm will be clamped to %f mV over the interval [%f, %f]",
              params.clamp_arg, params.clamp_start_arg, params.clamp_start_arg + params.clamp_dur_arg);
    } else {
      log_msg(NULL, 0, 0,  "info: Vm will be clamped from %f to %f mV over the interval [%f, %f]",
              params.clamp_ini_arg, params.clamp_arg, params.clamp_start_arg,
              params.clamp_start_arg + params.clamp_dur_arg);
    }
  }
 
  action_potential AP;
  if (analyze_AP) initialize_AP_analysis(&AP);
 
  Regions = static_cast<char *>(calloc(numNode, sizeof(char) ));
 
  SF::init_vector(&Vmv, numNode, -1);
  SF::init_vector(&I_ion, Vmv);
 
  int *stim_list = static_cast<int *>(malloc(numNode*sizeof(int) ));
  for (int i = 0; i < numNode; i++) {
    Regions[i]   = R1;
    stim_list[i] = i;
  }
 
  for (unsigned int i = 0; i < params.load_module_given; i++) {
#ifdef HAVE_DLOPEN
    load_ionic_module(params.load_module_arg[i]);
    if (!params.imp_given) {
 
      char *imp = strdup(basename(params.load_module_arg[i]) );
      char *dot = strrchr(imp, '.');
      if (dot)
        *dot = 0;
      params.imp_arg = imp;
    }
#else // ifdef HAVE_DLOPEN
    log_msg(NULL, 5, 0, "Compile with USE_DLOPEN to support dynamic module loading.");
    EXIT(-1);
#endif // ifdef HAVE_DLOPEN
  }
 
  if (get_ion_type(std::string(params.imp_arg)) == NULL) {
    fprintf(stderr, "Illegal IMP specified: %s\n", params.imp_arg);
    exit(1);
  } else {
    RegionDef = get_ion_type(std::string(params.imp_arg));
  }
  if (!get_plug_flag(params.plug_in_arg, &num_plugins, RegionPlug)) {
    fprintf(stderr, "Illegal plugin specified: %s\n", params.plug_in_arg);
    exit(1);
  }
  // Fetch target enum from command line string
  Target target = get_target_from_string(std::string(params.target_arg));
  // Automatic target selection
  if (target == Target::AUTO) {
    target = RegionDef->select_target(target);
  }
  // Check whether target exists and was generated for the given model and plugin
  if (target == Target::UNKNOWN) {
    fprintf(stderr, "Unkown target: %s\n", params.target_arg);
    fprintf(stderr, "Available targets are: %s\n", get_target_list_string().c_str());
    exit(1);
  }
  // Check if the main imp was generated for the requested target
  if (RegionDef->select_target(target) == Target::UNKNOWN) {
    fprintf(stderr, "Model %s was not generated for target %s\n",
        params.imp_arg, get_string_from_target(target).c_str());
    exit(1);
  }
  // Check the same thing if a plugin is specified
  if (!RegionPlug.empty()) {
    if (RegionPlug[0].get().select_target(target) == Target::UNKNOWN) {
      fprintf(stderr, "Plugin %s was not generated for target %s\n",
          params.plug_in_arg, get_string_from_target(target).c_str());
      exit(1);
    }
  }
 
  // data structure to manage multiple IIF's
  MULTI_IF MIIF;
  MULTI_IF doppel, *cMIIF = &MIIF;
  bool doppel_created = false;
 
  MIIF.N_IIF       = num_region;
  MIIF.iontypes    = {*RegionDef};
  MIIF.targets     = {target};
  MIIF.IIFmask     = Regions;
  MIIF.numplugs    = &num_plugins;
  MIIF.plugtypes   = {RegionPlug};
  MIIF.numNode     = Vmv->lsize();
  MIIF.gdata[Vm]   = Vmv;
  MIIF.gdata[Iion] = I_ion;
 
  // A manually imposed voltage protocol (--clamp / --clamp-file) owns Vm: it is updated
  // externally, so the ionic model must not integrate the membrane equation for Vm.
  MIIF.extUpdateVm = IsVclamp;
 
  double setup_time = timing(t2, t1);
  update_timing(timings+SETUP_IDX, setup_time);
 
  MIIF.initialize_MIIF();
 
  if (params.list_imps_flag || params.plugin_outputs_flag) {
    print_models(params.plugin_outputs_flag);
    exit(0);
  }
 
  if (params.imp_info_flag) {
    print_param_help(RegionDef, RegionPlug);
    exit(0);
  }
 
  fprintf(stderr, "\n*** GIT tag:              %s\n", GIT_COMMIT_TAG);
  fprintf(stderr, "*** GIT hash:             %s\n", GIT_COMMIT_HASH);
  fprintf(stderr, "*** GIT repo:             %s\n", GIT_PATH);
  fprintf(stderr, "*** dependency commits:   %s\n\n", SUBREPO_COMMITS);
 
  if (params.buildinfo_flag) {
    fprintf(stderr, "\n*** -buildinfo flag is deprecated and will be removed.\n");
    fprintf( stderr, "\n*** The build information is printed by default. \n\n");
    exit(0);
  }
 
  MIIF.IIF[0]->tune(params.imp_par_arg, params.plug_in_arg, params.plug_par_arg);//FIXME
 
  // now compute the tables
  MIIF.initialize_currents(dt, 1);
 
  if (analyze_AP) {
    // initialize AP analyzer
    double Vm0 = getCellVal(MIIF.gdata[Vm], 0);
    for (int i = 0; i < VM_HIST_LEN; i++)
      AP.vm_trc[i] = Vm0;
 
    char ap_stats_fname[1024];
    snprintf(ap_stats_fname, sizeof ap_stats_fname, "%s_AP_stats.dat", params.fout_arg);
    AP.stats = fopen(ap_stats_fname, "wt");
    print_AP_stats_header(&AP, AP.stats);
  }
 
  restitution r;
  int n_dopple;
  double *t_dopple;
  if (restitute) {
    char *res_file = dupstr(params.res_file_arg);
    restitution_trigger_list(res_file, &r, params.restitute_arg, &n_dopple, &t_dopple);
    if (r.dur > duration) duration = r.dur;
 
    // output restitution data
    char ap_rstats_fname[1024];
    snprintf(ap_rstats_fname, sizeof ap_rstats_fname, "%s_APD_restitution.dat", params.fout_arg);
    AP.rstats = fopen(ap_rstats_fname, "wt");
    print_AP_stats_header(&AP, AP.rstats);
  }
 
  // read initial state vector from file and/or assignment
  double t = 0.;
  // if (params.restore_given)
  //   t = MIIF.restore_state( params.restore_arg, &pl, true);
 
  // illumination stuff
  double light_irrad = params.light_irrad_arg;
  double light_start = params.light_start_arg;
  double light_end = params.duration_arg;
  int    light_nstim = params.light_numstim_arg;
  double light_bcl = params.light_bcl_arg;
  double light_dur = params.light_dur_arg;
 
  if(light_irrad < 0) {
    FPRINTF(WORLD stderr, "WARNING: irradiance value %.3f makes no sense (clipped to 0.0)\n", light_irrad);
    light_irrad = 0.0;
  }
  TrgList light_lst;
  if(params.light_times_given)  //user specified list
    determine_stim_list(params.light_times_arg, &light_lst, false);
 
  if (params.SV_init_given)
    initial_SVs(&MIIF, params.SV_init_arg, params.imp_arg, params.plug_in_arg, numNode);
   
  timer_manager tmo(dt, t, duration);
  tmo.timers.resize(N_TIMERS);
  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., CON_TM_IDX, "CON_TIMER");
  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0., SVD_TM_IDX, "SVD_TIMER");
 
  if (params.save_time_arg)
    tmo.initialize_singlestep_timer(params.save_time_arg, 0., STA_TM_IDX, "STA_TIMER");
 
  if (params.save_ini_time_given)
    tmo.initialize_singlestep_timer(params.save_ini_time_arg, 0., SSV_TM_IDX, "SSV_TIMER");
 
  if( RRC )
    tmo.initialize_eq_timer( 0, duration-1, RRC_MAX_N, bcl, 1., DOPPLE_TM_IDX, "DOPPLE_TIMER");
  if( Lapicque_crv )
    tmo.initialize_singlestep_timer( nstim*bcl, 0, DOPPLE_TM_IDX, "DOPPLE_TIMER");
  
  // use a trigger for stimulation
  if (params.restitute_given) {
    std::vector<double> trg;
    trg.assign(r.trigs.lst, r.trigs.lst + r.trigs.n);
    tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");
 
    if (params.res_state_vector_given) {
      trg.assign(r.saveState.lst, r.saveState.lst + r.saveState.n);
      tmo.initialize_neq_timer(trg, 0, RES_SAVE_TM_IDX, "RES_SAVE_TIMER");
    }
 
    if (!strcmp(params.restitute_arg, "S1S2f") ) {
      trg.assign(t_dopple, t_dopple + n_dopple);
      tmo.initialize_neq_timer(trg, 1, DOPPLE_TM_IDX, "DOPPLE_TIMER");
    }
  }
  else if (params.stim_times_given) {
    std::vector<double> trg;
    trg.assign(stim_lst.lst, stim_lst.lst + stim_lst.n);
    tmo.initialize_neq_timer(trg, stim_dur, STM_TM_IDX, "STIM_TIMER");
  }
  else {
    tmo.initialize_eq_timer(stim_start, stim_start+duration, nstim, bcl, stim_dur, STM_TM_IDX, "STIM_TIMER");
  }
 
  // illumination timers
  if (MIIF.gdata[illum] != NULL && !tr_light) {
    if (params.light_times_given) {
      std::vector<double> trg;
      trg.assign(light_lst.lst, light_lst.lst + light_lst.n);
      tmo.initialize_neq_timer(trg, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
    } else {
      tmo.initialize_eq_timer(light_start, light_end, light_nstim,
                              light_bcl, light_dur, LIGHT_TM_IDX, "LIGHT_TIMER", NULL);
    }
  }
 
  tmo.initialize_eq_timer(start_out, t+duration, 0, dt_out, 0, TRACE_TM_IDX, "TRACE_TIMER");
 
  if (IsVclamp && cl.transient)
    tmo.initialize_singlestep_timer(params.clamp_start_arg, params.clamp_dur_arg, CLAMP_TM_IDX, "CLAMP_TIMER");
  if (!restitute && params.doppel_on_given) {
    tmo.initialize_singlestep_timer(params.doppel_on_arg, params.doppel_dur_arg, DOPPLE_TM_IDX, "DOPPLE_TIMER");
  }
 
  // read initial state from single cell file and spread it out to all cells
  if (params.read_ini_file_given) {
    if (read_sv(&MIIF, R1, params.read_ini_file_arg) != 0) {
      return 1;  // exit with error code.
    }
  }
 
  // stretch experiment
  stretch s;
  memset(&s, 0, sizeof(stretch));
  if (params.strain_given && MIIF.use_stretch()) {
    initializePulseStretch(params.strain_arg, params.strain_time_arg,
                           params.strain_dur_arg, params.strain_rate_arg,
                           params.strain_rate_arg, &s);
 
    sf_vec* lambdavec; SF::init_vector(&lambdavec, MIIF.gdata[Vm]);
    MIIF.gdata[Lambda] = lambdavec;
    SF_real *l = lambdavec->ptr();
 
    for (int i = 0; i < lambdavec->lsize(); i++)
      l[i] = s.pulse.sr;
 
    lambdavec->release_ptr(l);
  }
 
  // Lapique curve
  Lapicque *lapicque = NULL;
  if( Lapicque_crv ) lapicque = new Lapicque(&tmo,params.Lap_min_arg,params.Lap_max_arg,
                        params.Lapicque_arg, params.thresh_tol_arg, params.thresh_tm_arg );
 
  // set up assigning stimulus current to ion concentrations
  std::stringstream stim_species(params.stim_species_arg);
  std::stringstream stim_ratios(params.stim_ratios_arg);
  std::string specie, ratio;
  std::vector<std::string> species_list;
  std::vector<float> ratios_list;
  float ratio_sum = 0;      
  if (stim_assign) {
    // split stim_species by ":" and then iterate over them
    while(std::getline(stim_species, specie, ':'))
    {
       species_list.push_back(specie);
    }
    while(std::getline(stim_ratios, ratio, ':'))
    {
       ratios_list.push_back(std::stof(ratio));
    }
    if (species_list.size() !=  ratios_list.size()) {
      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);
      exit(1);    
    }
    for (auto i : ratios_list) 
      ratio_sum += i;
    if (std::abs(ratio_sum - 1) >= 0.001)
      log_msg(NULL, 3, 0, "Warning: Sum of ion species ratios for stimulus assignement != 1. Normalizing to sum %.1f.", ratio_sum);
  }
 
  // dumping of state variables goes here
  char sv_dump_names[1024];
  if (!params.fout_given && !params.validate_flag)
    snprintf(sv_dump_names, sizeof sv_dump_names, "%s_%s", params.fout_arg, params.imp_arg);
  else
    strcpy(sv_dump_names, params.fout_arg);
 
  if (params.validate_flag) {
    dump_all(&MIIF, R1, params.imp_arg, params.plug_in_arg, t, dt, sv_dump_names);
  } else {
    MIIF.sv_dump_add_by_name_list(R1, params.imp_arg, sv_dump_names,
                             params.imp_sv_dump_arg, params.plug_in_arg,
                             params.plug_sv_dump_arg, t, dt_out);
  }
 
  // dumping of LUTs
  if (params.dump_lut_flag)
    MIIF.dump_luts_MIIF(true);
 
  // initialize tracing
  int trace_nodes[] = {0};
  if (do_trace)
    open_trace(&MIIF, sizeof(trace_nodes)/sizeof(trace_nodes[0]), trace_nodes, &params.trace_no_arg, NULL);
 
  // dumping global vectors to file?
  open_globalvec_dump(fhdls, &gvd, &MIIF, sv_dump_names, &io);
 
  // initialization time
  double ini_time = timing(t3, t2);
  update_timing(timings + INIT_IDX, ini_time);
  if (tmo.timers[DOPPLE_TM_IDX] != nullptr) {
    doppel_MIIF(&MIIF, &doppel);
    doppel_created = true;
  }
 
  int gpu_Vm_update = 1;
  double counter = 0.0;
  double epsilon = 0.0000001f;
  Bisect RRC_find( 0, INFINITY, params.RRC_tol_arg );
  float RRC_est = params.stim_curr_arg/10.;
  bool  RRC_repol;
 
  /***********************************************/
  /***************** MAIN LOOP *******************/
  /***********************************************/
 
  // Print simulation target
  fprintf(stderr, "Running simulation on target %s\n",      
                get_string_from_target(target).c_str());
  while (!tmo.elapsed())
  {
    get_time(t4);
    t = tmo.time;
 
    
    if (tmo.triggered_now(DOPPLE_TM_IDX) ) {
      printf("[%f] Switching to new ion structure\n", tmo.time);
      if( RRC ) {
        doppel_update(&doppel, &MIIF);
        if( tmo.timers[DOPPLE_TM_IDX]->trigger_count>1 ) {
          if( RRC_find.update( RRC_est, !RRC_repol ) ){
            std::cout << "\nRRC found: " << RRC_est << std::endl;
            break;
          }
        }
        set_trace_amp( &stim_trace, RRC_est, params.RRC_delay_arg );
        std::cout << RRC_est << " " << RRC_repol << std::endl;
        RRC_repol = false;
      } else if( Lapicque_crv ) {
        if( lapicque->init(DOPPLE_TM_IDX, STM_TM_IDX) ){
          bool fire = cMIIF->gdata[Vm]->get(0) >= params.vthreshold_arg;
          float stimcurr = stim_charge/dt;
          if( lapicque->check( fire, stimcurr, DOPPLE_TM_IDX, STM_TM_IDX ) ) 
            break;
          stim_charge = stimcurr*dt;
          doppel_update(&doppel, &MIIF);
        } else {
          doppel_update(&MIIF, &doppel);
        }
      } else {
        doppel_update(&MIIF, &doppel);
        cMIIF = &doppel;
      }
    }
 
    cMIIF->extUpdateVm = IsVclamp;   // external Vm control whenever a manual protocol is active
 
    if (IsVclamp)
      clamp_signal(cMIIF, &cl, &tmo);
 
    if (APclamp)
      AP_clamp(&ap_cl, &tmo, cMIIF->gdata[Vm], tmo.triggered_now(STM_TM_IDX) );
 
    for (int i = 0; i < SVclamp; i++) {
      bool trigger;
      if (params.SV_I_trigger_flag)
        trigger =  tmo.triggered_now(STM_TM_IDX);
      else
        trigger = !tmo.d_time;
 
      sv_clamp(sv_cl+i, &tmo, cMIIF, trigger);
    }
 
    // output selected state variables
    cMIIF->dump_svs(tmo.timers[SVD_TM_IDX]);
    globalvec_dump(fhdls, &gvd, cMIIF, &tmo, &io, numNode);
 
    // save state of all cells to binary file
    if (tmo.trigger(STA_TM_IDX))
      cMIIF->dump_state(params.save_file_arg, tmo.time, mesh_t::unset_msh, false, 0);
 
    // save state of cell 0 to text file
    if (tmo.trigger(SSV_TM_IDX))
      save_sv(cMIIF, R1, params.save_ini_file_arg);
 
    // save state within an restitution protocol
    if (tmo.trigger(RES_SAVE_TM_IDX))
      restitution_save_sv(cMIIF, R1, &r, &AP);
 
    // add stim current
    if (tmo.trigger(STM_TM_IDX)) {
      gpu_Vm_update = 0;
      if (tr_stim) {
        stim_charge = stim_trace.s[tmo.trigger_elapse(STM_TM_IDX)]*dt;
      }
      else if (params.stim_volt_given) {
        SF_real *d = cMIIF->gdata[Vm]->ptr();
        // peak_vm = last_vm = d[0];
        stim_charge = (params.stim_volt_arg-d[0])/params.resistance_arg*dt;
        cMIIF->gdata[Vm]->release_ptr(d);
      }
 
#ifndef ALGEBRAIC
      *(cMIIF->gdata[Vm]) += stim_charge;
#endif // ifndef ALGEBRAIC
 
      if (stim_assign) {
        int iIon = 0;
        for (auto ion : species_list)
        {
          cMIIF->transmem_stim_species(-stim_charge*(ratios_list[iIon]/ratio_sum), ion.c_str(), params.surface_to_volume_arg, stim_list, numNode);
          FPRINTF(WORLD stderr, "[stim_assign] @ t = %.3f ms; charge: %f pA/pF*ms; type: %s\n", tmo.time, stim_charge, ion.c_str());
          iIon++;
        }
      }
    }
 
    if (MIIF.gdata[illum] != NULL) {
      MIIF.gdata[illum]->set(tmo.trigger(LIGHT_TM_IDX) ? light_irrad : 0.0);
 
      if(tmo.triggered_now(LIGHT_TM_IDX) && light_irrad) {
        FPRINTF(WORLD stderr, "[stim] illum ON  @ t = %.3f ms; E_e = %.3f mW/mm^2\n", tmo.time, light_irrad);
      }
      if(tmo.trigger_end(LIGHT_TM_IDX) && light_irrad) {
        FPRINTF(WORLD stderr, "[stim] illum OFF @ t = %.3f ms\n", tmo.time);
      }
    }
 
    // strain
    if (params.strain_given && cMIIF->use_stretch())
      apply_stretch(cMIIF, &s, &tmo);
 
    if (analyze_AP) {
      if (restitute) AP.pmat = r.trigs.pmat[AP.beat];
      check_events(getCellVal(cMIIF->gdata[Vm], 0), &AP, &tmo);
    }
 
    get_time(t5);
 
 
    if (is_gpu(target)) {
      if (gpu_Vm_update == 0) {
        cMIIF->compute_ionic_current(1, 1);
      } else if (gpu_Vm_update == 1) {
        cMIIF->compute_ionic_current(fabs(counter - 0.0) < epsilon || fabs(counter - 1.0) < epsilon  , fabs((counter+dt) - 1.0) < epsilon);
      }
    }
    else {
      cMIIF->compute_ionic_current();
    }
 
    double ode_time = timing(t6, t5);
    update_timing(timings+ODE_IDX, ode_time);
 
    if (do_trace && tmo.trigger(TRACE_TM_IDX))
      dump_trace(cMIIF, t);
 
    if (!restitute && !RRC && !Lapicque_crv && tmo.trigger_end(DOPPLE_TM_IDX) ) {
      log_msg(NULL, 0, 0, "[%f] Switching back to original ion structure\n", tmo.time);
      cMIIF = &MIIF;
    }
 
    tmo.update_timers();
 
    double loop_time = timing(t7, t4);
    update_timing(timings+LOOP_IDX, loop_time);
 
    gpu_Vm_update = 1;
    counter = counter + dt;
    if ( fabs(counter - 1.0) < epsilon ) {
      counter = 0.0;
    }
 
    if( RRC && tmo.trigger(STM_TM_IDX) &&  
        tmo.trigger_elapse(STM_TM_IDX)*tmo.time_step>params.RRC_delay_arg ){
      auto data = cMIIF->gdata[Vm]->ptr();
      if( *data <= params.RRC_repol_arg )
        RRC_repol = true;
      cMIIF->gdata[Vm]->release_ptr(data);
    }
  }//loop end
 
  // loop_current time
  double main_time = timing(t8, t3);
 
  // write header so that we know which files belong together
  const char *ExpID = params.fout_arg;
  write_dump_header(&gvd, &MIIF.svd, ExpID);
 
  // clean up everything
  if( Lapicque_crv ) {
    if( params.Lap_file_given )  {
      std::ofstream ofs(params.Lap_file_arg);
      lapicque->print(ofs);
    } else
      lapicque->print(std::cout);
  }
  if (tr_stim) free_trace(&stim_trace);
  if (analyze_AP) cleanup_AP_analysis(&AP);
  if (do_trace) close_trace(&MIIF);
  MIIF.close_svs_dumps();
  if (doppel_created)
    free_doppel(&doppel);
  MIIF.free_MIIF();
 
  close_globalvec_dump(fhdls, &gvd, &io);
 
  log_msg(NULL, 0, 0, "\n\n\nAll done!\n\n");
 
  // print timing data
  event_timing *t_setup = timings+SETUP_IDX;
  log_msg(NULL, 0, 0, "setup time          %.6f s",   (double)t_setup->tot);
  event_timing *t_init = timings+INIT_IDX;
  log_msg(NULL, 0, 0, "initialization time %.6f s",   (double)t_init->tot);
  event_timing *t_loop = timings+LOOP_IDX;
  log_msg(NULL, 0, 0, "main loop time      %.6f s",   (double)t_loop->tot);
  event_timing *t_ode = timings+ODE_IDX;
  log_msg(NULL, 0, 0, "total ode time      %.6f s\n", (double)t_ode->tot);
 
  log_msg(NULL, 0, 0, "mn/avg/mx loop time %.6f %.6f %.6f ms",
          (double)t_loop->mn*1e3, (double)t_loop->avg*1e3, (double)t_loop->mx*1e3);
  log_msg(NULL, 0, 0, "mn/avg/mx ODE time  %.6f %.6f %.6f ms",
          (double)t_ode->mn*1e3, (double)t_ode->avg*1e3, (double)t_ode->mx*1e3);
 
  log_msg(NULL, 0, 0, "real time factor    %.6f\n\n", duration/(t_ode->tot*1e3) );
 
 
  CHKERRQ(PetscFinalize());
  return 0;
} // main