restitute.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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "restitute.h"
#include "ap_analyzer.h"
#include "sv_init.h"

namespace limpet {

using ::opencarp::log_msg;

void get_protocol_definition(char *r_file, restitution *r, bool useS1S2);
int  read_restitution_protocol_def(char *r_file, restitution *r);

void restitution_trigger_list(char *r_file, restitution *r, char *protocol, int *n_dop, double **t_dop)      {
  // define protocol
  get_protocol_definition(r_file, r, !strncmp("S1S2", protocol, 4));
  if (!strcmp(protocol, "S1S2f") ) r->prtcl = S1S2_fast;
  *n_dop = 0;

  // compute dependent parameters based on protocol definition
  if (r->prtcl == S1S2) {
    // S1-S2 protocol
    int num_decs = (r->rtype.S1S2.S2_start-r->rtype.S1S2.S2_end)/r->rtype.S1S2.S2_dec;
    int num_stim = r->numppBeats+num_decs*(r->rtype.S1S2.beats_per_S2+1);
    r->trigs.lst = (double *)calloc(num_stim, sizeof(double));
    if (r->trigs.lst == NULL) {
      /* Memory could not be allocated, so print an error and exit. */
      fprintf(stderr, "Couldn't allocate memory\n");
      exit(EXIT_FAILURE);
    }
    r->trigs.pmat = (bool *)calloc(num_stim, sizeof(bool));

    // build trigger list for saving state vectors
    r->saveState.lst  = (double *)calloc(num_decs, sizeof(double));
    r->saveState.pmat = NULL;

    int trgIdx     = 0;
    double trgTime = 1.0;
    r->trigs.pmat[trgIdx]  = false;
    r->trigs.lst[trgIdx++] = trgTime;
    for (int i = 1; i < r->numppBeats; i++) {
      trgTime               += r->rtype.S1S2.bcl;
      r->trigs.pmat[trgIdx]  = false;
      r->trigs.lst[trgIdx++] = trgTime;
    }

    for (int i = 0; i < num_decs; i++) {
      for (int j = 0; j < r->rtype.S1S2.beats_per_S2; j++) {
        trgTime               += r->rtype.S1S2.bcl;
        r->trigs.pmat[trgIdx]  = false;
        r->trigs.lst[trgIdx++] = trgTime;
      }
      trgTime               += r->rtype.S1S2.S2_start-(i+1)*r->rtype.S1S2.S2_dec;
      r->trigs.pmat[trgIdx]  = true;
      r->trigs.lst[trgIdx++] = trgTime;
      r->saveState.lst[i]    = trgTime;
    }
    r->trigs.n     = trgIdx;
    r->dur         = trgTime + r->rtype.S1S2.bcl;
    r->saveState.n = num_decs;
  } else if (r->prtcl == S1S2_fast) {
    // fast S1-S2 protocol, we save the state after S1's and only apply S2's
    int num_decs = (r->rtype.S1S2.S2_start-r->rtype.S1S2.S2_end)/r->rtype.S1S2.S2_dec;
    int num_stim = r->numppBeats-1+2*num_decs;
    r->trigs.lst  = (double *)calloc(num_stim, sizeof(double));
    r->trigs.pmat = (bool *)calloc(num_stim, sizeof(bool)  );

    // build trigger list for saving state vectors
    r->saveState.lst  = (double *)calloc(num_decs, sizeof(double));
    r->saveState.pmat = NULL;

    int trgIdx     = 0;
    double trgTime = 1.0;
    r->trigs.pmat[trgIdx]  = false;
    r->trigs.lst[trgIdx++] = trgTime;
    for (int i = 1; i < r->numppBeats-1; i++) {
      trgTime               += r->rtype.S1S2.bcl;
      r->trigs.pmat[trgIdx]  = false;
      r->trigs.lst[trgIdx++] = trgTime;
    }

    for (int i = 0; i < num_decs; i++) {
      trgTime               += r->rtype.S1S2.bcl;
      r->trigs.pmat[trgIdx]  = false;
      r->trigs.lst[trgIdx++] = trgTime;
      trgTime               += r->rtype.S1S2.S2_start-(i+1)*r->rtype.S1S2.S2_dec;
      r->trigs.pmat[trgIdx]  = true;
      r->trigs.lst[trgIdx++] = trgTime;
      r->saveState.lst[i]    = trgTime;
    }
    r->trigs.n     = trgIdx;
    r->dur         = trgTime + r->rtype.S1S2.bcl;
    r->saveState.n = num_decs;

    *n_dop = num_decs;
    *t_dop = static_cast<double *>(malloc(num_decs*sizeof(double)));
    trgIdx = r->numppBeats-1;
    for (int i = 0; i < num_decs; i++, trgIdx += 2)
      (*t_dop)[i] = r->trigs.lst[trgIdx] - 1;
  } else {
    // dynamic protocol
    int num_decs = (r->rtype.dyn.bcl_start-r->rtype.dyn.bcl_end)/r->rtype.dyn.bcl_dec;
    int num_stim = r->numppBeats+num_decs*r->rtype.dyn.beats_per_bcl;
    r->trigs.lst  = (double *)calloc(num_stim, sizeof(double));
    r->trigs.pmat = (bool *)calloc(num_stim, sizeof(bool));

    // build trigger list for saving state vectors
    r->saveState.lst  = (double *)calloc(num_decs, sizeof(double));
    r->saveState.pmat = NULL;

    int trgIdx     = 0;
    double trgTime = 1.0;
    r->trigs.pmat[trgIdx]  = false;
    r->trigs.lst[trgIdx++] = trgTime;
    for (int i = 1; i < r->numppBeats; i++) {
      trgTime               += r->rtype.dyn.bcl_start;
      r->trigs.pmat[trgIdx]  = false;
      r->trigs.lst[trgIdx++] = trgTime;
    }

    for (int i = 0; i < num_decs; i++) {
      for (int j = 0; j < r->rtype.dyn.beats_per_bcl; j++) {
        trgTime += r->rtype.dyn.bcl_start-i*r->rtype.dyn.bcl_dec;
        if ((i > 0) && (j == 0))
          r->trigs.pmat[trgIdx] = true;
        else
          r->trigs.pmat[trgIdx] = false;
        r->trigs.lst[trgIdx++] = trgTime;
      }
      r->saveState.lst[i] = trgTime;
    }
    r->trigs.n     = trgIdx;
    r->dur         = trgTime + r->rtype.dyn.bcl_end;
    r->saveState.n = num_decs;
  }
} // restitution_trigger_list

void get_protocol_definition(char *r_file, restitution *r, bool useS1S2) {
  // file based
  if (strcmp(r_file, "")) {
    if (read_restitution_protocol_def(r_file, r)) {
      log_msg(NULL, 0, 0, "\n\nError reading protocol definition from %s",r_file);
    }
    else {
      // if specified protocol matches protocol in file we are done
      if( (useS1S2&&r->prtcl==S1S2) || (!useS1S2&&r->prtcl==DYNAMIC) )
        return;

      // file does not match, we don't use the protocol file and revert to defaults
      log_msg(NULL, 3, 0, "Restitution protocol definition file does not match.\n");
      log_msg(NULL, 3, 0, "Using default protocol definitions.\n");

      r->numppBeats =  20;
      if (useS1S2) {
        r->prtcl                   = S1S2;
        r->rtype.S1S2.bcl          = 1000;
        r->rtype.S1S2.S2_start     = 500;
        r->rtype.S1S2.S2_end       = 200;
        r->rtype.S1S2.beats_per_S2 =  12;
        r->rtype.S1S2.S2_dec       =  10;
      } else   {
        r->prtcl                   = DYNAMIC;
        r->rtype.dyn.bcl_start     = 400;
        r->rtype.dyn.bcl_end       = 200;
        r->rtype.dyn.beats_per_bcl =  10;
        r->rtype.dyn.bcl_dec       =  10;
      }
    }
  }
} // get_protocol_definition

int read_restitution_protocol_def(char *r_file, restitution *r)      {
  const int BUFSIZE = 256;
  char buf[256];

  FILE *pdef = fopen(r_file, "rt");

  if (!pdef) return -1;

  // skip possible empty or comment lines
  do {
    fgets(buf, BUFSIZE, pdef);
  } while ((*buf == '\n') || (*buf == '#'));

  int useS1S2;
  sscanf(buf, "%d", &useS1S2);

  if (useS1S2) {
    r->prtcl = S1S2;
    sscanf(fgets(buf, BUFSIZE, pdef), "%d", &r->numppBeats);
    if (r->numppBeats < CALIBRATION_BEAT) {
      log_msg(NULL, 2, 0, "#prepacing beats %d < #calibration beats %d.",
              r->numppBeats, CALIBRATION_BEAT);
      log_msg(NULL, 2, 0, "This may bias the calibration of AP statistics parameters.");
    }

    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.S1S2.bcl);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.S1S2.S2_start);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.S1S2.S2_end);
    sscanf(fgets(buf, BUFSIZE, pdef), "%d", &r->rtype.S1S2.beats_per_S2);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.S1S2.S2_dec);
  } else   {
    r->prtcl = DYNAMIC;
    sscanf(fgets(buf, BUFSIZE, pdef), "%d", &r->numppBeats);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.dyn.bcl_start);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.dyn.bcl_end);
    sscanf(fgets(buf, BUFSIZE, pdef), "%d", &r->rtype.dyn.beats_per_bcl);
    sscanf(fgets(buf, BUFSIZE, pdef), "%f", &r->rtype.dyn.bcl_dec);
  }

  fclose(pdef);

  return 0;
} // read_restitution_protocol_def

void restitution_save_sv(MULTI_IF *miif, int R1, restitution *r, action_potential *AP)      {
  char save_name[1024];
  static int cnt = 0;

  snprintf(save_name, sizeof save_name, "test_%d.sv", cnt++);
  save_sv(miif, R1, save_name);
}

}  // namespace limpet