openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
clamp.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
19 #include <stdbool.h>
20 #include <assert.h>
21 
22 #include "MULTI_ION_IF.h"
23 #include "clamp.h"
24 
25 #include "petsc_utils.h" // TODO: for EXIT
26 
27 namespace limpet {
28 
32 using ::opencarp::timer_manager;
33 
44 bool
45 initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f,
46  int trans, float *duration)
47 {
48  bool fmode = (bool) strcmp(f, "");
49  bool isClamp = dur || fmode;
50  memset(cl, 0, sizeof(Clamp));
51 
52  if (fmode) {
53  cl->file = dupstr(f);
54  *duration = trace_duration(&cl->tr, cl->file);
55  }
56  else {
57  cl->dur = dur;
58  cl->val_pre = ini_val;
59  cl->val = cl_val;
60  cl->val_post = ini_val;
61  cl->start = start;
62  cl->tau = 0.1; // clamp time constant in ms
63  cl->transient = trans;
64  }
65 
66  return isClamp;
67 }
68 
69 void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
70 {
71  if(cl->transient && ! tm->trigger(cl->transient)) return;
72 
73  pMIIF->getRealData();
74 
75  SF_real *v = pMIIF->procdata[Vm];
76  static double v0;
77  static bool fst = true;
78 
79  if (fst) {
80  v0 = v[0];
81  fst = false;
82  if (cl->file) {
83  read_trace(&cl->tr, cl->file);
84  resample_trace(&cl->tr, tm->time_step);
85  }
86  if(cl->transient)
87  cl->val_pre = v0;
88  }
89 
90  double val;
91  if (!cl->file) {
92  double time = tm->time;
93 
94  if (time <= cl->start)
95  val = cl->val_pre + (v0-cl->val_pre)*exp(-time / cl->tau);
96  else if (time > cl->start && time < cl->start+cl->dur)
97  val = cl->val + (cl->val_pre-cl->val)*exp(-(time - cl->start)/cl->tau);
98  else
99  val = cl->val_post + (cl->val-cl->val_post)*exp(-(time-(cl->start+cl->dur))/cl->tau);
100  } else
101  val = cl->tr.s[tm->d_time];
102 
103  for (int i=0; i<pMIIF->numNode;i++)
104  v[i] = val;
105 
106  pMIIF->releaseRealData();
107 }
108 
109 
117 void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
118 {
119  memset(cl, 0, sizeof(Clamp));
120 
121  read_trace(&cl->tr, cl->file=file);
122  resample_trace(&cl->tr, dt);
123  cl->sv = strdup(sv);
124 
125  cl->impDataID = IMPdataLabel2Index(sv);
126  cl->transient = -1000000; // last one was a real long time ago
127 }
128 
129 
137 void
138 sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
139 {
140  int d_time = tm->d_time;
141 
142  if(trigger)
143  cl->transient = d_time;
144 
145  if(d_time - cl->transient < cl->tr.N) {
146  if(cl->impDataID>=0) {
147  if(miif->gdata[cl->impDataID] == NULL) {
148  log_msg(NULL, 5, 0, "Cannot clamp %s as the global vector %s is not used!!\n\n", cl->sv, cl->sv);
149  EXIT(-1);
150  }
151  miif->gdata[cl->impDataID]->set(cl->tr.s[d_time - cl->transient]);
152  }
153  else {
154  int offset;
155  int temparray[1] = {0};
156 
157  SVgetfcn svget = miif->IIF[0]->get_type().get_sv_offset(cl->sv, &offset, temparray);
158  if(!svget) {
159  log_msg(NULL, 5, 0, "\nCannot clamp %s as there is no %s!!\n\n", cl->sv, cl->sv);
160  EXIT(1);
161  }
162  SVputfcn svput= getPutSV(svget);
163  for(int i=0; i<miif->N_Nodes[0]; i++)
164  svput( *miif->IIF[0], i, offset, cl->tr.s[d_time - cl->transient]);
165  }
166  }
167 }
168 
169 
177 void
178 AP_clamp(Clamp *cl, timer_manager *tm, sf_vec* v, bool trigger)
179 {
180  if(trigger)
181  cl->transient = tm->d_time;
182 
183  if(tm->d_time-cl->transient < cl->tr.N)
184  v->set(cl->tr.s[tm->d_time-cl->transient]);
185 }
186 
187 
196 int
197 process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
198 {
199  *clamps = NULL;
200  int clno=0;
201 
202  if(!strlen(SVs) || !strlen(files))
203  return 0;
204 
205  char *svlist = strdup(SVs);
206  char *flist = strdup(files);
207  char *sptr, *fptr;
208 
209  // count the number of :'s
210  int ncol = 0;
211  char *p = SVs;
212  while((p = strchr(p, ':'))!=NULL) {
213  p++;
214  ncol++;
215  }
216 
217  *clamps = (Clamp*)calloc(ncol+1, sizeof(Clamp));
218 
219  char* sv = strtok_r(svlist, ":", &sptr);
220  char* file = strtok_r(flist, ":", &fptr);
221  while(sv != NULL && file!= NULL) {
222  initialize_sv_clamp((*clamps)+clno++, sv, file, dt);
223  sv = strtok_r(NULL, ":", &sptr);
224  file = strtok_r(NULL, ":", &fptr);
225  }
226 
227  if(sv != file) {
228  log_msg(NULL, 5, 0, "Mismatch between SV list and file list");
229  exit(1);
230  }
231 
232  free(svlist);
233  free(flist );
234 
235  return clno;
236 }
237 
238 } // namespace limpet
239 
GlobalData_t(* SVgetfcn)(IonIfBase &, int, int)
Definition: ion_type.h:48
int transient
Definition: clamp.h:37
char * dupstr(const char *old_str)
Definition: basics.cc:44
double dur
Definition: clamp.h:33
trace tr
Definition: clamp.h:36
int read_trace(trace *tr, const char *name)
Definition: trace.cc:43
GlobalData_t * procdata[NUM_IMP_DATA_TYPES]
data for this processor
Definition: MULTI_ION_IF.h:204
void resample_trace(trace *tr, double dt)
Definition: trace.cc:115
double val_pre
Definition: clamp.h:29
void sv_clamp(Clamp *cl, timer_manager *tm, MULTI_IF *miif, bool trigger)
Definition: clamp.cc:138
bool initialize_clamp(Clamp *cl, double cl_val, double ini_val, double start, double dur, const char *f, int trans, float *duration)
Definition: clamp.cc:45
double val_post
Definition: clamp.h:31
void(* SVputfcn)(IonIfBase &, int, int, GlobalData_t)
Definition: ion_type.h:49
SVputfcn getPutSV(SVgetfcn)
double start
Definition: clamp.h:32
std::vector< IonIfBase * > IIF
array of IIF&#39;s
Definition: MULTI_ION_IF.h:202
char * sv
Definition: clamp.h:38
Define multiple ionic models to be used in different regions.
double tau
Definition: clamp.h:34
double trace_duration(trace *tr, const char *f)
Definition: trace.cc:147
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false)=0
float * s
samples
Definition: trace.h:28
void clamp_signal(MULTI_IF *pMIIF, Clamp *cl, timer_manager *tm)
Definition: clamp.cc:69
int IMPdataLabel2Index(const char *sv)
SF::abstract_vector< SF_int, SF_real > sf_vec
Definition: sf_interface.h:49
int impDataID
Definition: clamp.h:39
int * N_Nodes
#nodes for each IMP
Definition: MULTI_ION_IF.h:200
opencarp::sf_vec * gdata[NUM_IMP_DATA_TYPES]
data used by all IMPs
Definition: MULTI_ION_IF.h:216
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Definition: basics.cc:72
double val
Definition: clamp.h:30
int N
number of samples
Definition: trace.h:26
int numNode
local number of nodes
Definition: MULTI_ION_IF.h:210
char * file
Definition: clamp.h:35
void AP_clamp(Clamp *cl, timer_manager *tm, sf_vec *v, bool trigger)
Definition: clamp.cc:178
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
void initialize_sv_clamp(Clamp *cl, const char *sv, char *file, double dt)
Definition: clamp.cc:117
int process_sv_clamps(char *SVs, char *files, Clamp **clamps, double dt)
Definition: clamp.cc:197