clamp.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 <stdbool.h>
#include <assert.h>

#include "MULTI_ION_IF.h"
#include "clamp.h"

#include "petsc_utils.h" // TODO: for EXIT

namespace limpet {

using ::opencarp::dupstr;
using ::opencarp::log_msg;
using ::opencarp::sf_vec;
using ::opencarp::timer_manager;

bool
initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f,
                 int trans, float *duration)
{
  bool fmode   = (bool) strcmp(f, "");
  bool isClamp = dur || fmode;
  memset(cl, 0, sizeof(Clamp));

  if (fmode) {
    cl->file      = dupstr(f);
    *duration     = trace_duration(&cl->tr, cl->file);
  }
  else  {
    cl->dur       = dur;
    cl->val_pre   = ini_val;
    cl->val       = cl_val;
    cl->val_post  = ini_val;
    cl->start     = start;
    cl->tau       = 0.1;  // clamp time constant in ms
    cl->transient = trans;
  }

  return isClamp;
}

void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
{
  if(cl->transient && ! tm->trigger(cl->transient)) return;

  pMIIF->getRealData();

  SF_real     *v   = pMIIF->procdata[Vm];
  static double  v0;
  static bool    fst = true;

  if (fst) {
    v0  = v[0];
    fst = false;
    if (cl->file) {
      read_trace(&cl->tr, cl->file);
      resample_trace(&cl->tr, tm->time_step);
    }
    if(cl->transient)
      cl->val_pre = v0;
  }

  double val;
  if (!cl->file) {
    double time   = tm->time;

    if (time <= cl->start)
      val = cl->val_pre + (v0-cl->val_pre)*exp(-time / cl->tau);
    else if (time > cl->start && time < cl->start+cl->dur)
      val = cl->val + (cl->val_pre-cl->val)*exp(-(time - cl->start)/cl->tau);
    else
      val = cl->val_post + (cl->val-cl->val_post)*exp(-(time-(cl->start+cl->dur))/cl->tau);
  } else
    val = cl->tr.s[tm->d_time];

  for (int i=0; i<pMIIF->numNode;i++)
    v[i] = val;

  pMIIF->releaseRealData();
}


void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
{
  memset(cl, 0, sizeof(Clamp));

  read_trace(&cl->tr, cl->file=file);
  resample_trace(&cl->tr, dt);
  cl->sv = strdup(sv);

  cl->impDataID = IMPdataLabel2Index(sv);
  cl->transient = -1000000;       // last one was a real long time ago
}


void
sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
{
  int d_time = tm->d_time;

  if(trigger)
    cl->transient = d_time;

  if(d_time - cl->transient < cl->tr.N) {
    if(cl->impDataID>=0)  {
      if(miif->gdata[cl->impDataID] == NULL)  {
        log_msg(NULL, 5, 0, "Cannot clamp %s as the global vector %s is not used!!\n\n", cl->sv, cl->sv);
        EXIT(-1);
      }
      miif->gdata[cl->impDataID]->set(cl->tr.s[d_time - cl->transient]);
    }
    else  {
      int offset;
      int temparray[1] = {0};

      SVgetfcn svget = miif->IIF[0]->get_type().get_sv_offset(cl->sv, &offset, temparray);
      if(!svget) {
        log_msg(NULL, 5, 0, "\nCannot clamp %s as there is no %s!!\n\n", cl->sv, cl->sv);
        EXIT(1);
      }
      SVputfcn svput= getPutSV(svget);
      for(int i=0; i<miif->N_Nodes[0]; i++)
        svput( *miif->IIF[0], i, offset, cl->tr.s[d_time - cl->transient]);
    }
  }
}


void
AP_clamp(Clamp *cl, timer_manager *tm, sf_vec* v, bool trigger)
{
  if(trigger)
    cl->transient = tm->d_time;

  if(tm->d_time-cl->transient < cl->tr.N)
    v->set(cl->tr.s[tm->d_time-cl->transient]);
}


int
process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
{
  *clamps = NULL;
  int clno=0;

  if(!strlen(SVs) || !strlen(files))
    return 0;

  char *svlist = strdup(SVs);
  char *flist  = strdup(files);
  char *sptr, *fptr;

  // count the number of :'s
  int   ncol = 0;
  char *p    = SVs;
  while((p = strchr(p, ':'))!=NULL)  {
    p++;
    ncol++;
  }

  *clamps = (Clamp*)calloc(ncol+1, sizeof(Clamp));

  char* sv   = strtok_r(svlist, ":", &sptr);
  char* file = strtok_r(flist, ":", &fptr);
  while(sv != NULL  && file!= NULL) {
    initialize_sv_clamp((*clamps)+clno++, sv, file, dt);
    sv   = strtok_r(NULL, ":", &sptr);
    file = strtok_r(NULL, ":", &fptr);
  }

  if(sv != file) {
    log_msg(NULL, 5, 0, "Mismatch between SV list and file list");
    exit(1);
  }

  free(svlist);
  free(flist );

  return clno;
}

}  // namespace limpet