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 <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];

int* curr_node_list;   

IonIfBase::IonIfBase(const IonType& type, Target target, int 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;
}

int 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(int start, int 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;
}

int IonIfBase::restore(FILE_SPEC in, int n, const int *pos, IIF_Mask_t *mask,
                    size_t *offset, IMPinfo *impinfo, const int* loc2canon) {
  if( !impinfo->compatible )
    return n;

  char *ptr  = (char *)(this->get_sv_address());
  long  base = ftell(in->fd);
  int   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( int i=0; i<n; i+=vec_size ) {
    int index = i / vec_size;
    int 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;
}

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, int 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