ION_IF.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
// ----------------------------------------------------------------------------
 
/*
 *  Specify an interface for an IMP
 *
 *  To use an IMP by itself outside of MULTI_IF.h, one must \n
 *  -# alloc_IIF()
 *  -# set the rdata pointers in the ::IMPDataStruct
 *  -# initialize_params()
 *  -# initialize_IIF()
 *  -# call the compute() method of the IMP
 *
 *  The IMP library tries to avoid undo computation by using lookup
 *  tables to avoid computation of costly mathematical functions.
 *  Tables indices may be computed based on vltage, calcium concentrations,
 *  or any other quantity.
 *
 *  One common scenario is commonly encountered with first-order differential
 *  equations:
 *  \f[ \frac{dz}{dt} = \frac{z_{\infty}(x) - z}{\tau_z(x)} \f]
 *
 *  Using an exponential solution assuming that
 *  \f$z_{\inf}(x)\f$ and \f$\tau_z(x)\f$ constant over a time step:
 *
 *  \f[ \begin{split}
 *  z_{i+1} & = z_{\infty}(x) - (z_{\infty}(x)-z_i)
 *              \exp\Large (-\Delta t/\tau_z(x) \Large ) \\
 *         & = A(x) + B(x) z_i
 *  \end{split}\f]
 *  where A(x) and B(x) are entered into a look up table which is
 *  indexed by \em x:
 *  \f[ \begin{split}
 *    A(x) &=  z_{\infty}(x)\Large ( 1 - \exp(-\Delta t/\tau_z(x)\Large ) \\
 *         &= -z_{\infty}(x) expm1\Large ( -\Delta t/\tau_z(x) \Large ) \\
 *    B(x) &=  \exp(-\Delta t/\tau_z(x))
 *  \end{split}\f]
 *
 */
#include "limpet_types.h"
#include "ION_IF.h"
#include "petsc_utils.h" // for EXIT
#include <stdarg.h>
 
#ifdef USE_CVODE
#include <nvector/nvector_serial.h>
#include <cvode/cvode.h>
#include <cvode/cvode_diag.h>
#endif
 
namespace limpet {
 
using ::opencarp::dupstr;
using ::opencarp::FILE_SPEC;
using ::opencarp::log_msg;
using ::opencarp::Salt_list;
 
#ifdef USE_HDF5
FILE_SPEC _nc_logf = (FILE_SPEC)calloc(1, sizeof(fileptr));
#else
FILE_SPEC _nc_logf = NULL;
#endif
 
 
char* tokstr_r(char *s1, const char *s2, char **lasts)
{
  char *ret;
 
  if (s1 == NULL)
    s1 = *lasts;
  while(*s1 && strchr(s2, *s1))
    ++s1;
  if(*s1 == '\0')
    return NULL;
  ret = s1;
  while(*s1 && !strchr(s2, *s1))
    ++s1;
  if(*s1)
    *s1++ = '\0';
  *lasts = s1;
  return ret;
}
 
 
#ifndef offsetof
#define offsetof(TYPE, MEMBER) ((size_t) &((TYPE *)0)->MEMBER)
#endif
 
// struct exception_context the_exception_context[1];
 
node_index_t* curr_node_list;   
 
IonIfBase::IonIfBase(const IonType& type, Target target, node_count_t num_node, const std::vector<std::reference_wrapper<IonType>>& plugins)
  : _type(type), _target(type.select_target(target)), _num_node(num_node) {
  this->_tstp = {};
  this->_reqdat = type.reqdat();
  this->_moddat = type.moddat();
  this->_plugins = {};
 
  for (const auto& plugin : plugins) {
    auto child = plugin.get().make_ion_if(target, num_node, {});
    child->_parent = this;
    this->_plugins.push_back(child);
 
    this->_reqdat |= child->_reqdat;
    this->_moddat |= child->_moddat;
  }
}
 
IonIfBase::~IonIfBase() {
  for (auto plugin : this->_plugins) {
    plugin->get_type().destroy_ion_if(plugin);
  }
 
  this->_type.destroy(*this);
  if (this->_tables_d != nullptr) {
    deallocate_on_target<LUT>(this->_target, this->_tables_d);
  }
}
 
const IonType& IonIfBase::get_type() const {
  return this->_type;
}
 
node_count_t IonIfBase::get_num_node() const {
  return this->_num_node;
}
 
std::size_t IonIfBase::get_num_threads() const {
  std::size_t num_threads = 0;
  switch (this->_target) {
    // CPU target threads number depends on wether OpenMP is enabled or not
    case Target::CPU:
    case Target::MLIR_CPU:
#ifdef _OPENMP
      num_threads = omp_get_max_threads();
#else
      num_threads = 1;
#endif
      break;
    // GPU targets threads numbers is the number of cells
    case Target::MLIR_CUDA:
    case Target::MLIR_ROCM:
      num_threads = this->get_num_node();
      break;
    default:
      throw std::logic_error("The IonIf execution target " + std::to_string(this->_target) + " is invalid");
      break;
  }
  return num_threads;
}
 
IonIfBase* IonIfBase::parent() const {
  return this->_parent;
}
 
void IonIfBase::set_parent(IonIfBase* parent) {
  this->_parent = parent;
}
 
std::vector<IonIfBase*>& IonIfBase::plugins() {
  return this->_plugins;
}
 
uint32_t IonIfBase::get_reqdat() const {
  return this->_reqdat;
}
 
uint32_t IonIfBase::get_moddat() const {
  return this->_moddat;
}
 
void IonIfBase::set_moddat(uint32_t data) {
  this->_moddat = data;
}
 
float IonIfBase::get_dt() const {
  return this->dt;
}
 
void IonIfBase::set_dt(float dt) {
  this->dt = dt;
}
 
ts& IonIfBase::get_tstp() {
  return this->_tstp;
}
 
std::vector<LUT>& IonIfBase::tables() {
  return this->_tables;
}
 
size_t IonIfBase::get_n_tables_d() const {
  return this->_n_tables_d;
}
 
void IonIfBase::set_target(Target target) {
  if (is_concrete(this->get_type().select_target(target))) {
      this->_target = target;
  }
  else {
      throw std::invalid_argument("new target set for IMP is unavailable or not concrete (AUTO, UNKNWOWN, ...)");
  }
}
 
void IonIfBase::initialize_params() {
  this->_type.initialize_params(*this);
 
  for (auto& plugin : this->_plugins) {
    plugin->initialize_params();
  }
}
 
void IonIfBase::initialize(double dt, GlobalData_t **impdat) {
  this->dt = dt;
 
  if (this->_num_node) {
    this->_type.construct_tables(*this);
    // If this is a GPU target, we need to copy the LUT definitions to the
    // device
    if (is_gpu(this->get_target()) && this->tables().size() > 0) {
      // Allocate on device and copy
      this->_tables_d = allocate_on_target<LUT>(this->get_target(), this->tables().size());
      // TODO replace this with some kind of memcpy (needs to be adjusted
      // depending on CUDA / HIP
      for (size_t i = 0; i < this->tables().size(); ++i) {
        this->_tables_d[i] = this->tables()[i];
      }
      this->_n_tables_d = this->tables().size();
    }
  }
 
  this->_type.initialize_sv(*this, impdat);
  this->ldata = impdat;
 
  for (auto& plugin : this->_plugins) {
    plugin->initialize(dt, impdat);
  }
}
 
void IonIfBase::compute(node_index_t start, node_index_t end, GlobalData_t **data) {
  this->_type.compute(this->_target, start, end, *this, data);
}
 
char* IonIfBase::fill_buf(char *buf, int* n, opencarp::Salt_list *l) const {
  int iif_sz = this->get_sv_size();
 
  for (auto& plugin : this->_plugins) {
    iif_sz += plugin->get_sv_size();
  }
 
  buf = (char *) realloc(buf, *n + l->nitems*iif_sz) + *n;
  *n += l->nitems * iif_sz;
 
  return buf - *n + l->nitems * iif_sz;
}
 
node_count_t IonIfBase::restore(FILE_SPEC in, node_count_t n, const node_index_t *pos, IIF_Mask_t *mask,
                    size_t *offset, IMPinfo *impinfo, const global_node_index_t* loc2canon) {
  if( !impinfo->compatible )
    return n;
 
  char *ptr  = (char *)(this->get_sv_address());
  long  base = ftell(in->fd);
  node_count_t mismatch = 0;
 
  // Restoration of the model data need to take care of the data layout
  // optimization. (if it is disabled, the vector size should be 1 anyways
  std::size_t vec_size = this->get_type().dlo_vector_size();
  for( node_index_t i=0; i<n; i+=vec_size ) {
    node_index_t index = i / vec_size;
    global_node_index_t canon = loc2canon[pos[i]];
    if(mask[canon] == this->miifIdx) {
      fseek(in->fd, base+offset[canon], SEEK_SET);
      fread(ptr+index*impinfo->sz, impinfo->sz, 1, in->fd);
 
      for(int j=0; j<impinfo->nplug; j++) {
        if(impinfo->plug[j].compatible) {
          fread((char*)(this->_plugins[impinfo->plug[j].map]->get_sv_address())+index*impinfo->plug[j].sz,
                impinfo->plug[j].sz, 1, in->fd);
        }
        else
          fseek(in->fd, impinfo->plug[j].sz, SEEK_CUR);
      }
    }
    else
      mismatch++;
  }
  return mismatch;
}
 
const std::vector<IonIfBase::sv_field>& IonIfBase::sv_layout() const {
  if (!_sv_layout_valid) {
    char** list = NULL;
    int nsv = this->get_type().get_sv_list(&list);
    _sv_layout.reserve(nsv);
    for (int i = 0; i < nsv; i++) {
      int off = 0, sz = 0, type = 0;
      char* type_name = NULL;
      this->get_type().get_sv_offset(list[i], &off, &sz);
      this->get_type().get_sv_type(list[i], &type, &type_name);
      free(type_name);
      _sv_layout.push_back({list[i], off, sz, type});
      free(list[i]);
    }
    free(list);
    _sv_layout_valid = true;
  }
  return _sv_layout;
}
 
void IonIfBase::get_sv_layout(std::vector<std::pair<int, int>>& fields) const {
  const std::vector<sv_field>& layout = this->sv_layout();
  fields.clear();
  fields.reserve(layout.size());
  for (const auto& f : layout)
    fields.emplace_back(f.off, f.sz);
}
 
size_t IonIfBase::get_sv_per_node_size() const {
  size_t total = 0;
  for (const auto& f : this->sv_layout())
    total += static_cast<size_t>(f.sz);
  return total;
}
 
uint64_t IonIfBase::sv_fingerprint() const {
  constexpr uint64_t kFNVOffsetBasis = 0xcbf29ce484222325ULL;  // 14695981039346656037
  constexpr uint64_t kFNVPrime       = 0x100000001b3ULL;       // 1099511628211
  uint64_t h = kFNVOffsetBasis;
  auto mix = [&h](unsigned char b) { h ^= b; h *= kFNVPrime; };
  for (const auto& f : this->sv_layout()) {
    for (unsigned char c : f.name) mix(c);
    mix(0);  // name terminator, so {"ab","c"} and {"a","bc"} differ
    for (int b = 0; b < 4; b++)
      mix((static_cast<uint32_t>(f.sz) >> (8 * b)) & 0xff);
    for (int b = 0; b < 4; b++)
      mix((static_cast<uint32_t>(f.type) >> (8 * b)) & 0xff);
  }
  return h;
}
 
node_count_t IonIfBase::restore_per_node(FILE_SPEC in, node_count_t n, const node_index_t *pos,
                    IIF_Mask_t *mask, size_t *offset, IMPinfo *impinfo,
                    const global_node_index_t* loc2canon) {
  if( !impinfo->compatible )
    return n;
 
  long base = ftell(in->fd);
  node_count_t mismatch = 0;
 
  // field layouts are model properties: gather them once, not per node
  std::vector<std::pair<int, int>> main_fields;
  this->get_sv_layout(main_fields);
 
  std::vector<std::vector<std::pair<int, int>>> plug_fields(impinfo->nplug);
  for (int j = 0; j < impinfo->nplug; j++)
    if (impinfo->plug[j].compatible)
      this->_plugins[impinfo->plug[j].map]->get_sv_layout(plug_fields[j]);
 
  for( node_index_t i=0; i<n; i++ ) {
    global_node_index_t canon = loc2canon[pos[i]];
    if(mask[canon] == this->miifIdx) {
      fseek(in->fd, base+offset[canon], SEEK_SET);
      // a short read means the checkpoint is truncated within the state payload,
      // which the header-level checks in restore_state cannot catch; abort rather
      // than scatter stale bytes into the state. A zero-width field (e.g. a
      // placeholder SV whose per-lane size rounds to 0 under DLO) reads nothing;
      // fread reports 0 items for size 0, which must not be read as truncation.
      auto scatter = [&](char* dst, std::size_t sz) {
        if (sz && fread(dst, sz, 1, in->fd) != 1) {
          log_msg(NULL, 5, 0, "Truncated ionic-state checkpoint during restore");
          EXIT(EXIT_FAILURE);
        }
      };
 
      this->for_each_sv_field(i, main_fields, scatter);
      for(int j=0; j<impinfo->nplug; j++) {
        if(impinfo->plug[j].compatible)
          this->_plugins[impinfo->plug[j].map]->for_each_sv_field(i, plug_fields[j], scatter);
        else
          fseek(in->fd, impinfo->plug[j].sz, SEEK_CUR);
      }
    }
    else
      mismatch++;
  }
  return mismatch;
}
 
int IonIfBase::dump_luts(bool zipped) {
  std::string name;
  std::string ext = zipped ? ".gz" : "";
 
  int dcnt = 0;
  for (int i=0; i < this->_tables.size(); i++) {
    // determine file name
    if (strcmp(this->_tables[i].name, ""))
      name = this->_type.get_name() + "_LUT_" + this->_tables[i].name + ".bin" + ext;
    else
      name = this->_type.get_name() + "_LUT_" + std::to_string(i) + ".bin" + ext;
 
    // dump table
    int err = LUT_dump(&this->_tables[i], name.c_str());
    if (!err) dcnt++;
  }
  return dcnt;
}
 
void IonIfBase::destroy_luts() {
  for (auto& lut : this->_tables) {
    destroy_lut(&lut, this->_target);
  }
 
  std::vector<LUT>().swap(this->_tables);
}
 
void IonIfBase::tune(const char *im_par, const char *plugs, const char *plug_par) {
  char *opar, *oplg, *plg, *nplg, *parlst, *nparlst;
 
  if( im_par && *im_par != '\0' ) {
    log_msg( _nc_logf, 0, 0, "Ionic model: %s", this->_type.get_name().c_str());
    this->_type.tune(*this, im_par);
  }
  opar = parlst = dupstr(plug_par);
  oplg = plg    = dupstr(plugs);
 
  while( plg!=NULL && *plg!='\0' ) {
    nparlst = get_next_list( parlst, ':' );
    nplg    = get_next_list( plg, ':' );
 
    log_msg( _nc_logf, 0, 0, "Plug-in: %s", plg );
 
    if( parlst==NULL || *parlst=='\0' ) {
      parlst = nparlst;
      plg    = nplg;
      continue;
    }
 
    // find apropriate plugin
    IonType* plugin = get_ion_type(plg);
    int j;
 
    for(j = 0; j < this->_plugins.size(); j++) {
      if(this->_plugins[j]->get_type() == *plugin) break;
    }
 
    if(j == this->_plugins.size()) {
      log_msg( _nc_logf, 2, 0, "Plugin %s not used with IM", plg );
      plg    = nplg;
      parlst = nparlst;
      continue;
    }
 
    plugin->tune(*this->_plugins[j], parlst);
    plg    = nplg;
    parlst = nparlst;
  }
 
  log_msg( _nc_logf, 0, 0, "" );
  free( opar );
  free( oplg );
}
 
int IonIfBase::read_svs(FILE* file) {
  if(this->get_num_node())
    return this->_type.read_svs(*this, file);
  else
    return 1;
}
 
int IonIfBase::write_svs(FILE* file, node_index_t node) {
  return this->_type.write_svs(*this, file, node);
}
 
void IonIfBase::copy_plugins_from(IonIfBase& other) {
  this->_plugins = {};
 
  for (auto& plugin : other.plugins()) {
    auto copy = plugin->get_type().make_ion_if(this->_target, plugin->get_num_node(), {});
    copy->_parent = this;
    this->_plugins.push_back(copy);
    copy->copy_SVs_from(*plugin, true);
  }
}
 
void IonIfBase::for_each(const std::function<void(IonIfBase&)>& consumer) {
  consumer(*this);
 
  for (auto& plugin : this->_plugins) {
    consumer(*plugin);
  }
}
 
// Is this ever used?
void initialize_ts(Target target, ts *tstp, int ng, int *skp, double dt)
{
  // initialize step counter that will run linearly as long the
  // simulation goes
  tstp->cnt = -1;
  tstp->ng  = ng;
  tstp->tcg = allocate_on_target<tc_grp>(target, ng);
 
  // initialize time constant grouping
  // fast variables, use dt=dt
  tstp->tcg[0].skp    = 1;
  tstp->tcg[0].dt     = dt;
  tstp->tcg[0].update = 1;
 
  for (int i=1;i<ng;i++)  {
    tstp->tcg[i].skp    = skp[i];
    tstp->tcg[i].dt     = tstp->tcg[i].skp*dt;
  }
}
 
 
void update_ts(ts *ptstp)
{
  ptstp->cnt++;
 
  tc_grp *ptcg = ptstp->tcg;
  for (int i=1;i<ptstp->ng;i++)
    ptcg[i].update = !(ptstp->cnt%ptcg[i].skp);
}
 
 
#ifndef _GNU_SOURCE
void *memmem( void *haystack, int sz_hay, void *needle, int sz_n )
{
  if( !sz_n ) return NULL;
 
  char *h = (char *)haystack;
 
  for( int i=0; i>sz_hay-sz_n+1; i++, h++ )
    if( !memcmp( h, needle, sz_n ) )
      return (void *)h;
     break;
 
  return NULL;
}
#endif
 
 
char *get_next_list( char *lst, char delimiter )
{
  if( lst == NULL || *lst == '\0' )
    return NULL;
 
  while( *lst != delimiter && *lst != '\0' )
    lst++;
 
  if( *lst != '\0' ) {
    *lst = '\0';
    return lst+1;
  } else
  return lst;
}
 
 
 bool verify_flags( const char *flags, const char* given )
 {
   char *flag, *ptr, *gvn_cp = dupstr( given );
 
   flag = tokstr_r( gvn_cp, "|", &ptr );
   while( flag ) {
     if( !flag_set( flags, flag ) )
       break;
     flag = tokstr_r( NULL, "|", &ptr );
   }
 
   free( gvn_cp );
   return flag ? false : true;
}
 
bool flag_set( const char *flags, const char *target )
{
  if( !flags ) return false;
 
  char *f  = dupstr( flags ), *last, *pos;
 
  pos = tokstr_r( f, "|", &last );
  while( pos && strcmp( pos, target ) ) {
    pos = tokstr_r( NULL, "|", &last );
  }
  free( f );
 
  return pos? true : false;
}
 
 
#ifdef USE_CVODE
 
void __bogus_function_for_cvode() {
 
    int flag;
    int N_CVODE = 1;
#if SUNDIALS_VERSION_MAJOR < 4
    void* cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
#elif SUNDIALS_VERSION_MAJOR < 6
    void* cvode_mem = CVodeCreate(CV_BDF);
#elif SUNDIALS_VERSION_MAJOR < 7
    SUNContext sunctx;
    MPI_Comm comm = PETSC_COMM_WORLD;
    if( SUNContext_Create(&comm, &sunctx )<0 ){
      assert(0);
    }
    void* cvode_mem = CVodeCreate(CV_BDF, sunctx);
#else  // version from 7.0.0
    SUNContext sunctx;
    SUNComm comm = PETSC_COMM_WORLD;
    if( SUNContext_Create( comm, &sunctx )<0 ){
      assert(0);
    }
    void* cvode_mem = CVodeCreate(CV_BDF, sunctx);
#endif
    assert(cvode_mem != NULL);
#if SUNDIALS_VERSION_MAJOR >= 6
    N_Vector cvode_y = N_VNew_Serial(N_CVODE, sunctx);
#else
    N_Vector cvode_y = N_VNew_Serial(N_CVODE);
#endif
    flag = CVodeInit(cvode_mem, NULL, 0, cvode_y);
    assert(flag == CV_SUCCESS);
    N_VDestroy(cvode_y);
 
    flag = CVodeSStolerances(cvode_mem, 1e-5, 1e-6);
    assert(flag == CV_SUCCESS);
    flag = CVodeSetMaxStep(cvode_mem, 1);
    assert(flag == CV_SUCCESS);
    flag = CVodeSetUserData(cvode_mem, NULL);
    assert(flag == CV_SUCCESS);
    flag = CVDiag(cvode_mem);
    assert(flag == CV_SUCCESS);
 
  NV_Ith_S(cvode_y,0) = 0;
 
  {
    int CVODE_flag;
    CVODE_flag = CVodeReInit(NULL, 1, NULL);
    assert(CVODE_flag == CV_SUCCESS);
    CVODE_flag = CVodeSetInitStep(NULL, 1);
    assert(CVODE_flag == CV_SUCCESS);
#if SUNDIALS_VERSION_MAJOR >= 7
    sunrealtype tret;
#else
    realtype tret;
#endif
    CVODE_flag = CVode(cvode_mem, 1, NULL, &tret, CV_NORMAL);
    assert(CVODE_flag == CV_SUCCESS);
  }
}
#endif
 
}  // namespace limpet