stimulate.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#include "stimulate.h"
#include "electrics.h"

#include "SF_init.h" // for SF::init_xxx()
#include "petsc_utils.h" // TODO for EXIT

namespace opencarp {

// stimulus convenience functions
// group stimulus types
namespace stim_info {
  std::set<int> voltage_stim {V_ex, GND_ex, V_ex_ol};
  std::set<int> current_stim {I_tm, I_ex, I_in, I_lat};
  std::set<int> dbc_stim     {V_ex, GND_ex, V_ex_ol};
  std::map<int,std::string> units;
}  // namespace stim_info

namespace user_globals {

  extern timer_manager* tm_manager;
}  // namespace user_globals

void init_stim_info(void)
{
  // start dealing with units, first try
  stim_info::units[I_tm]      = "uA/cm^2";   
  stim_info::units[I_ex]      = "uA/cm^3";   
  stim_info::units[V_ex]      = "mV";        
  stim_info::units[GND_ex]    = "mV";        
  stim_info::units[I_in]      = "uA/cm^3";   
  stim_info::units[V_ex_ol]   = "mV";        
  stim_info::units[I_tm_grad] = "";          
  stim_info::units[I_lat]     = "uA/cm^2";   
  stim_info::units[Vm_clmp]   = "mV";        
  stim_info::units[Phie_pre]  = "mV";        
}


bool is_voltage(stim_t type)        
{
  return stim_info::voltage_stim.count(type);
}

bool is_current(stim_t type)
{
  return stim_info::current_stim.count(type);
}

bool is_dbc(stim_t type)
{
  return stim_info::dbc_stim.count(type);
}

bool is_extra(stim_t type)
{
  switch(type) {
    default:
    case I_tm:
    case Vm_clmp:
      return false;

    case I_ex:
    case V_ex:
    case V_ex_ol:
    case GND_ex:
      return true;
  }
}


void stimulus::translate(int iidx)
{
  idx = iidx;

  // translate protocol
  Stim     & s       = param_globals::stim[idx];
  Stimulus & curstim = param_globals::stimulus[idx];

  s.ptcl.bcl      = curstim.bcl;
  s.ptcl.npls     = curstim.npls;
  s.ptcl.start    = curstim.start;
  s.ptcl.duration = curstim.duration;

  // translate pulse
  // create a name, old stimulus structure lacks this feature, no naming of pulse waveforms
  std::string pulse_name;
  pulse_name = "pulse_" + std::to_string(idx);
  s.pulse.name = (char*) calloc(pulse_name.length()+1, sizeof(char));
  strcpy(s.pulse.name, pulse_name.c_str());

  // old stimulus treats any pulse as truncated exponential shape
  s.pulse.shape = truncExpPulse;
  s.pulse.file  = strdup(curstim.pulse_file);
  s.pulse.bias  = curstim.bias;
  s.pulse.tau_plateau = curstim.tau_plateau;
  s.pulse.tau_edge    = curstim.tau_edge;
  s.pulse.s2          = curstim.s2;
  // transition to new tilt-based shape defintion
  s.pulse.tilt_ampl   = curstim.s2;
  s.pulse.tilt_time   = curstim.d1;
  s.pulse.strength    = curstim.strength;

  // strength is not a scalar, we prescribe a different strength at every node
  s.elec.vtx_fcn      = curstim.vtx_fcn;

  // translate circuit
  s.crct.type = curstim.stimtype;
  s.crct.balance = curstim.balance;
  s.crct.total_current = curstim.total_current;

  // translate electrode
  s.elec.geomID   = curstim.geometry;
  // s.elec.domain   = curstim.domain;
  s.elec.vtx_file = strdup(curstim.vtx_file);
  s.elec.dump_vtx_file = curstim.dump_vtx_file;

  // check whether we have a sound electrode geometry definition
  s.elec.geom_type = 2; // should be region_t block;
  s.elec.p0[0] = curstim.x0 - (curstim.ctr_def?curstim.xd/2.:0.);
  s.elec.p0[1] = curstim.y0 - (curstim.ctr_def?curstim.yd/2.:0.);
  s.elec.p0[2] = curstim.z0 - (curstim.ctr_def?curstim.zd/2.:0.);
  s.elec.p1[0] = s.elec.p0[0] + curstim.xd;
  s.elec.p1[1] = s.elec.p0[1] + curstim.yd;
  s.elec.p1[2] = s.elec.p0[2] + curstim.zd;

  // cp overall electrode name
  s.name = strdup(curstim.name);
}

void stimulus::setup(int iidx)
{
  idx = iidx;
  const Stim & cur_stim = param_globals::stim[idx];

  // label stimulus
  if(cur_stim.name != std::string(""))
    name = cur_stim.name;
  else
    name = "Stimulus_" + std::to_string(idx);

  // set up stimulus pulse
  pulse.setup(idx);

  // set up stimulation protocol
  ptcl.setup(idx,name);

  // set up physics of applied stimulus
  phys.setup(idx);

  // set up the spatial region of the stimulus
  electrode.setup(idx);
}

bool stimulus::is_active() const
{
  // extracellular potential stimuli are always active
  if(this->phys.type == V_ex)
    return true;

  return user_globals::tm_manager->trigger(ptcl.timer_id);
}


void stim_pulse::setup(int idx)
{
  const Stim & cur_stim = param_globals::stim[idx];
  const Pulse & pulse = cur_stim.pulse;

  // assign waveform type and parameters
  waveform_t wvt = truncExpPulse;  // default, we may expose this in prm as a choice

  // set waveform based on circuit type
  switch(cur_stim.crct.type) {
    case Vm_clmp:
    case GND_ex:
      wvt = constPulse; break;

    default:
      if(pulse.file && strlen(pulse.file)) wvt = arbPulse;
      else                                 wvt = truncExpPulse;
      break;
  }

  assign(pulse.strength, cur_stim.ptcl.duration, param_globals::dt, wvt);

  // sample pulse shape
  if(wvt == arbPulse) {
    sig::time_trace file_wave;
    int _err = 0;

    if(!get_rank()) {
      bool unitize = true;

      update_cwd();
      set_dir(INPUT);
      _err = file_wave.read_trace(pulse.file, unitize);
      set_dir(CURDIR);
    }

    int err = get_global(_err, MPI_MIN);
    if(err) {
      log_msg(user_globals::physics_reg[elec_phys]->logger, MAX_LOG_LEVEL, ECHO,
              "Failed reading pulse for stimulus[%d] from file %s.\n", idx, pulse.file);
    }
    else {
      get_global(file_wave.f);
      get_global(file_wave.t);

      file_wave.resample(this->wave);

      // check length of read in pulse and adjust stimulus settings accordingly
      if(wave.duration() != duration)
        duration = wave.duration();
    }
  }
  else {
    sample_wave_form(*this, idx);
  }

  // add unit labels
  //std::string sUnit = IsVoltageStim_(cur_stim.crct) ? "mV" : "uA/cm^3";
  //std::string sUnit = is_voltage((stim_t)cur_stim.crct.type) ? "mV" : "uA/cm^3";
  std::string sUnit = stim_info::units[(stim_t)cur_stim.crct.type];
  std::string tUnit = "ms";
  wave.setUnits(tUnit,sUnit);
}


void get_stim_list(const char* str_list, std::vector<double> & stlist)
{
  std::string tstim = str_list;
  std::vector<std::string> stims;
  split_string(tstim, ',', stims);

  stlist.resize(stims.size());
  size_t widx = 0, err = 0;

  for(size_t i=0; i<stims.size(); i++) {
    std::string & stim = stims[i];

    if(stim.size()) {
      if(stim[0] == '+') {
        double last = widx == 0 ? 0.0 : stlist[widx-1];
        stim.erase(0, 1); // erase first char

        // we user sscanf so that we can actually determine if the string was converted
        // into a double. atof does not allow for error-checking
        double cur = -1.0;
        size_t nread = sscanf(stim.c_str(), "%lf", &cur);
        if(nread)
          stlist[widx++] = cur + last;
      } else {
        double cur = -1.0;
        size_t nread = sscanf(stim.c_str(), "%lf", &cur);
        if(nread)
          stlist[widx++] = cur;
      }
    }
  }

  stlist.resize(widx);

  if(widx != stims.size())
    log_msg(0, 4, 0, "%s warning: Some values of %s could not be converted in stim times!",
            __func__, str_list);
}


void stim_protocol::setup(int idx, std::string name)
{
  const Stim & cur_stim = param_globals::stim[idx];
  const Protocol & ptcl = cur_stim.ptcl;

  if(strlen(ptcl.stimlist) == 0) {
    start = ptcl.start, npls = ptcl.npls, pcl = ptcl.bcl;

    // add timer
    timer_id = user_globals::tm_manager->add_eq_timer(ptcl.start,
               param_globals::tend, npls, pcl, ptcl.duration, name.c_str());
  } else {
    std::vector<double> stims;
    get_stim_list(ptcl.stimlist, stims);

    log_msg(0,0,ECHO | NONL, "stim %d: stimulating at times: ", idx);
    for(double t : stims)
      log_msg(0,0, ECHO | NONL, "%.2lf ", t);
    log_msg(0,0,0, "");
    log_msg(0,0,0, "stim %d: .bcl and .npls values will be ignored.", idx);

    timer_id = user_globals::tm_manager->add_neq_timer(stims, ptcl.duration, name.c_str());
  }
}


void stim_physics::setup(int idx)
{
  const Stim & cur_stim = param_globals::stim[idx];
  type = stim_t(cur_stim.crct.type);
  domain = stim_domain_t(cur_stim.elec.domain);
  unit = stim_info::units[cur_stim.crct.type];
  total_current = cur_stim.crct.total_current;

  const short dim = is_extra(stim_t(cur_stim.crct.type)) ? get_mesh_dim(extra_elec_msh) : get_mesh_dim(intra_elec_msh);

  switch(type) {
    case I_tm:
      scale = param_globals::operator_splitting ? user_globals::tm_manager->time_step : 1.0;
      break;

    case I_ex:
      scale = UM2_to_CM2 * (dim == 2 ? 1.0 : UM_to_CM);
      break;

    default: scale = 1.0; break;
  }
}


void sample_wave_form(stim_pulse& sp, int idx)
{
  const Stimulus & curstim = param_globals::stimulus[idx];

  if(sp.wform == truncExpPulse)
  {
    // monophasic pulse, first phase lasts for entire pulse duration
    sig::Pulse truncExp;
    truncExp.duration = sp.duration;
    truncExp.d1       = curstim.d1 * sp.duration;  // relative to absolute duration of d1
    truncExp.tau_edge = curstim.tau_edge;
    truncExp.tau_plat = curstim.tau_plateau;
    truncExp.decay    = 5.;

    // from relative duration of d1 we infer mono- or biphasic
    if(curstim.d1 == 1.)  {
      // monophasic pulse, first phase lasts for entire pulse duration
      sig::monophasicTruncExpFunc fnc(truncExp);
      fnc.sample(sp.wave);
    }
    else
    {
      // biphasic pulse, second phase lasts shorter than pulse duration
      sig::biphasicTruncExpFunc fnc(truncExp);
      fnc.sample(sp.wave);
    }
  }
  else if(sp.wform == constPulse)
  {
    sig::constFunc cnst_fnc(1);
    cnst_fnc.sample(sp.wave);
  }
  else
    assert(0);

  /*
  else if(sp.wform == squarePulse)
  {
    // generate step function using derived class
    Step stepPars;
    stepPars.trig = 0.25;
    stepPars.rise = true;

    stepFunc step_func(stepPars);
    s.set_labels("step_func");
    step_func.sample(s);
    s.write_trace();
  }
  else if(sp.wform == sinePulse)
  {
    // generate a sine wave
    sineWave  sinePars;
    sinePars.frq   = 1.;
    sinePars.phase = 0.;
    sineFunc sine(sinePars);
    s.set_labels("sine_func");
    sine.sample(s);
    s.write_trace();
  }
  else if(sp.wform == APfootPulse)
  {
    // generate AP foot
    APfoot footPars;
    footPars.tau_f = s.duration()/5.;   // use 5 time constants
    APfootFunc APft(footPars);
    s.set_labels("APfoot");
    APft.sample(s);
    s.write_trace();
  }
  else if(sp.wform == arbPulse)
  {
  }
  else
    log_msg(user_globals::  error);
  */
}

bool stimulus::value(double & v) const
{
  if(user_globals::tm_manager->trigger(ptcl.timer_id) == false) {
    // extracellular potential stimuli are always active, but return 0 when their
    // protocol is not active. Use V_ex_ol for extracellular potential stimuli that
    // get removed when they expire.
    if(this->phys.type == V_ex || this->phys.type == GND_ex) {
      v = 0.0;
      return true;
    }

    return false;
  }
  else {
    int i = user_globals::tm_manager->trigger_elapse(ptcl.timer_id);
    v = pulse.strength * phys.scale * pulse.wave.f[i];
    return true;
  }
}

void stim_electrode::setup(int idx)
{
  const Stim & curstim = param_globals::stim[idx];
  mesh_t        grid;

  // based on the stimulus domain type we select grid and physics
  switch(curstim.crct.type) {
    case I_ex:
    case V_ex:
    case V_ex_ol:
    case GND_ex:
      grid   = extra_elec_msh; break;

    case Illum:
    case Vm_clmp:
    case I_tm:
      grid   = intra_elec_msh; break;

    case Ignore_stim: return;

    default:
      log_msg(0,5,0, "stim_electrode::setup error: Can't determine domain from stim type! Aborting!", __func__);
      EXIT(1);
  }

  // get a logger from the physics
  FILE_SPEC logger = get_physics(elec_phys)->logger;

  // this are the criteria for choosing how we extract the stimulus electrode
  bool vertex_file_given = strlen(curstim.elec.vtx_file) > 0;
  bool tag_index_given   = curstim.elec.geomID > -1;
  // the mesh we need for computing the local vertex indices.
  const sf_mesh & mesh   = get_mesh(grid);

  if(vertex_file_given && tag_index_given)
    log_msg(0,3,0, "%s warning: More than one stimulus electrode definintions set in electrode %d", __func__, idx);

  if(vertex_file_given) {
    definition     = def_t::file_based;
    input_filename = curstim.elec.vtx_file;

    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from file %s", idx, input_filename.c_str());

    set_dir(INPUT);
    warn_when_passing_intra_vtx(input_filename);

    // we read the indices. they are being localized w.r.t. the provided numbering. In our
    // case this is always the reference numbering.
    if(curstim.elec.vtx_fcn)
      read_indices_with_data(vertices, scaling, input_filename, mesh, SF::NBR_REF, true, 1, PETSC_COMM_WORLD);
    else
      read_indices(vertices, input_filename, mesh, SF::NBR_REF, true, PETSC_COMM_WORLD);

    int gnum_idx = get_global(vertices.size(), MPI_SUM);
    if(gnum_idx == 0) {
      log_msg(0, 5, 0, "Stimulus %d: Specified vertices are not in stimulus domain! Aborting!", idx);
      EXIT(1);
    }
  }
  else if(tag_index_given) {
    definition = def_t::vol_based_tag;

    int tag = curstim.elec.geomID;
    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from tag %d", idx, tag);

    indices_from_region_tag(vertices, mesh, tag);
    // we restrict the indices to the algebraic subset
    SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());
  }
  else {
    definition = def_t::vol_based_shape;
    log_msg(logger, 0, 0, "Stimulus %d: Selecting vertices from shape.", idx);

    geom_shape shape;
    shape.type   = geom_shape::shape_t(curstim.elec.geom_type);
    shape.p0     = curstim.elec.p0;
    shape.p1     = curstim.elec.p1;
    shape.radius = curstim.elec.radius;

    bool nodal = true;
    indices_from_geom_shape(vertices, mesh, shape, nodal);
    // we restrict the indices to the algebraic subset
    SF::restrict_to_set(vertices, mesh.pl.algebraic_nodes());

    int gsize = vertices.size();
    if(get_global(gsize, MPI_SUM) == 0) {
      log_msg(0,5,0, "error: Empty stimulus[%d] electrode def! Aborting!", idx);
      EXIT(1);
    }
  }

  if(curstim.elec.dump_vtx_file) {
    SF::vector<mesh_int_t> glob_idx(vertices);
    mesh.pl.globalize(glob_idx);

    SF::vector<mesh_int_t> srt_idx;
    SF::sort_parallel(mesh.comm, glob_idx, srt_idx);

    size_t num_vtx = get_global(srt_idx.size(), MPI_SUM);
    int rank = get_rank();

    FILE_SPEC f = NULL;
    int err = 0;

    if(rank == 0) {
      char dump_name[1024];

      if (strlen(curstim.name)) {
        snprintf(dump_name, sizeof dump_name, "%s.vtx", curstim.name);
      } else {
        snprintf(dump_name, sizeof dump_name, "ELECTRODE_%d.vtx", idx);
      }

      f = f_open(dump_name, "w");
      if(!f) err++;
      else {
        fprintf(f->fd, "%zd\nextra\n", num_vtx);
      }
    }

    if(!get_global(err, MPI_SUM)) {
      print_vector(mesh.comm, srt_idx, 1, f ? f->fd : NULL);
    } else {
      log_msg(0, 4, 0, "error: stimulus[%d] cannot be dumped!");
    }

    // only root really does that
    f_close(f);
  }
}

void dbc_manager::recompute_dbcs()
{
  clear_active_dbc();

  for(const stimulus & s : stimuli) {
    if(is_dbc(s.phys.type) && s.is_active()) {

      // the local indices of the dbc
      const SF::vector<mesh_int_t> & dbc_idx = s.electrode.vertices;

      dbc_data* dbc_buff = new dbc_data();

      // the global petsc indices of the dbc
      dbc_buff->nod = new SF::vector<SF_int>(dbc_idx.size());

      const SF::vector<mesh_int_t> & petsc_nbr = mat.mesh_ptr()->get_numbering(SF::NBR_PETSC);

      size_t widx = 0;
      for(size_t i=0; i<dbc_idx.size(); i++)
        (*dbc_buff->nod)[widx++] = petsc_nbr[dbc_idx[i]];

      sf_vec* dbc; SF::init_vector(&dbc, *mat.mesh_ptr(), 1, sf_vec::algebraic);
      SF::init_vector(&dbc_buff->cntr, *mat.mesh_ptr(), 1, sf_vec::algebraic);

      if(s.electrode.scaling.size()) {
        // we have spatially heterogenous dirichlet values
        dbc->set(*dbc_buff->nod, s.electrode.scaling);
      } else {
        // we have spatially homogenous dirichlet values
        dbc->set(*dbc_buff->nod, 1.0);
      }

      mat.mult(*dbc, *dbc_buff->cntr);

      active_dbc[s.idx] = dbc_buff;
    }
  }
}

bool dbc_manager::dbc_update()
{
  // check if stim has come offline
  for(const auto & d : active_dbc) {
    if(stimuli[d.first].is_active() == false)
      return true;
  }

  // check if stim has come online
  for(const stimulus & s : stimuli) {
    if(is_dbc(s.phys.type) && s.is_active() && active_dbc.count(s.idx) == 0)
      return true;
  }

  return false;
}

void dbc_manager::enforce_dbc_lhs()
{
  sf_vec* dbc; 

  SF::init_vector(&dbc, *mat.mesh_ptr(), mat.dpn_row(), sf_vec::algebraic);

  for(auto it = active_dbc.begin(); it != active_dbc.end(); ++it) {
    const SF::vector<SF_int> & dbc_nod = *it->second->nod;
    // zero cols and rows in lhs mat
    dbc->set(1.0);
    dbc->set(dbc_nod, 0.0);
    // multiplicate dbc vector from left and right to matrix in order to zero cols and rows
    mat.mult_LR(*dbc, *dbc);
    // now add 1.0 to the diagonal at the dirichlet locations
    dbc->set(0.0);
    dbc->set(dbc_nod, 1.0);
    mat.diag_add(*dbc);
  }
}

void dbc_manager::enforce_dbc_rhs(sf_vec &rhs)
{
  for(auto it = active_dbc.begin(); it != active_dbc.end(); ++it) {
    const int & dbc_idx = it->first;
    const dbc_manager::dbc_data & dbc = *(it->second);
    double strength = 0.0;
    bool is_active = stimuli[dbc_idx].value(strength);

    // the whole idea of the dbc_manager and its active_dbc is
    // based on only storing the active dbcs. if the stimulus associated
    // to a dbc believed active is actually not active, we want to abort.
    assert(is_active);

    rhs.add_scaled(*dbc.cntr, -strength);

    if(stimuli[dbc_idx].electrode.scaling.size()) {
      SF::vector<SF_real> scaled_strength(stimuli[dbc_idx].electrode.scaling.size());

      size_t widx = 0;
      for(SF_real s : stimuli[dbc_idx].electrode.scaling)
        scaled_strength[widx++] = s * strength;

      rhs.set(*dbc.nod, scaled_strength);
    } else {
      rhs.set(*dbc.nod, strength);
    }
  }
}

}  // namespace opencarp