clamp.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#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 (node_index_t 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(node_index_t 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