bench_utils.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 "bench_utils.h"
#include "petsc_compat.h"

#include "libgen.h"

namespace limpet {

using ::opencarp::dupstr;
using ::opencarp::f_open;
using ::opencarp::get_rank;
using ::opencarp::is_big_endian;
using ::opencarp::log_msg;
using ::opencarp::sf_vec;
using ::opencarp::timer_manager;

namespace {

void write_global(sf_vec* v, GVEC_DUMP *gvd, FILE *out, MULTI_IF *pMIIF, IOCtrl *io, int numNode) {
  char stdsep[] = " ";
  char fsep[]   = "\n";
  static char sep[2];
  SF_real *buf;

  if (io->first) {
    // choose separator
    if (io->wsplt)
      strcpy(sep, fsep);
    else
      strcpy(sep, stdsep);

    io->first = 0;
  }

  buf = v->ptr();

#ifndef ENABLE_TESTING
  if (!get_rank()) {
    if (io->w2file) {
      if (io->wbin) {
        // fwrite(&val, 1, sizeof(float), out );
        fwrite(buf, 1, data_type_sizes[gvd->dtype[0]], out);
      } else {
        fprintf(out, "%+.8e%s", buf[0], sep);
      }
    }
  }

  // stdout only
  if (io->w2stdout)
    log_msg(NULL, 0, NONL, "%+.8e%s", buf[0], stdsep);
#else
  for(int i = 0; i< 8; i++)
  {
    if (io->w2stdout) {
      int cell_node = rand()%numNode;
      log_msg(NULL, 0, NONL, "\n\t(Cell: %d) ", cell_node);
      log_msg(NULL, 0, NONL, "%+.6e%s", buf[cell_node], stdsep);
    }
  }
#endif

  v->release_ptr(buf);
}

int double_cmp(const void *a, const void *b) {
  return *((double *)a)-*((double *)b);
}

}  // unnamed namespace

void determine_stim_list(char *stl, TrgList *trg, bool DIAs) {
  trg->n = 1;
  for (size_t i = 1; i < strlen(stl)-1; i++)
    trg->n += stl[i] == ',';
  trg->lst = static_cast<double *>(malloc(trg->n * sizeof(double) ));

  char *sp, *token;
  token = tokstr_r(stl, ",", &sp);
  int i = 0;
  while (token != NULL) {
    char *endp;
    trg->lst[i++] = strtod(token, &endp);
    if (endp == token) {
      fprintf(stderr, "Error in stimulus timing list\n");
      exit(1);
    }
    if (DIAs && (i > 1)) trg->lst[i-1] += trg->lst[i-2];
    token = tokstr_r(NULL, ",", &sp);
  }
  qsort(trg->lst, trg->n, sizeof(double), double_cmp);
}

/*  write a header file for state variable dumps
 *  This is to determine with postprocessing routines
 *  which files belong together to a particular experiment.
 *  Particularly, this will be used to gather validation output
 *  in a single hdf5 file.
 *
 * \param svd   pointer to state variable dumping structure
 * \param ExpID string to identify an experiment
 *
 * \return 0 if header written successfully, -1 otherwise
 *
 */
int write_dump_header(GVEC_DUMP *gvd, SV_DUMP *svd, const char *ExpID) {
  int retval = 0;

  // did we dump some state variables
  if (!svd->active)
    return retval;

  if (!get_rank()) {
  char hd_fname[512];
  snprintf(hd_fname, sizeof hd_fname, "%s_header.txt", ExpID);
  FILE *fh = fopen(hd_fname, "wt");
  if (fh) {
    fprintf(fh, "%d  # is bigendian\n", is_big_endian() );

    // write global vectors first
    for (int i = 0; i < NUM_IMP_DATA_TYPES+1; i++) {
      if (gvd->hdls[i]) {
        fprintf(fh, "%32s %10s %2d %10d\n", basename(gvd->fn[i]),
                data_type_names[gvd->dtype[i]],
                data_type_sizes[gvd->dtype[i]], gvd->n_dumps);
      }
    }

    // now state variables
    for (int i = 0; i < svd->n; i++) {
      fprintf(fh, "%32s %10s %2d %10d\n", basename(svd->fn[i]), data_type_names[svd->dtype[i]],
              data_type_sizes[svd->dtype[i]], svd->n_dumps);
    }
    fclose(fh);
    } else {
      fprintf(stderr, "Could not open %s.\n", hd_fname);
      retval = -1;
    }
  }

  PetscBarrier(COMPAT_PETSC_NULLPTR);
  return retval;
} // write_dump_header

void open_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, char *base_name, IOCtrl *io) {
  char nbuf[1024];
  char outname[] = {"GVECS"};
  char ext[5];
  char wmode[3];
  const char *nptr;


  // initialize
  io->first = 1;
  memset(gvd, 0, sizeof(GVEC_DUMP) );
  gvd->n_dumps = 0;

  // strategy
  if (io->w2file) {
    // writing global vector output to a file
    if (io->wbin) {
      io->wsplt    = 1;
      io->w2stdout = 0;
      strcpy(ext, ".bin");
      strcpy(wmode, "wb");
    } else {
      io->wsplt    = 0;
      io->w2stdout = 0;
      strcpy(ext, ".txt");
      strcpy(wmode, "wt");
    }
  } else {
    // output to stdout only
    io->wsplt    = 0;
    io->w2stdout = 1;
  }

  // open file handles
  if (io->w2file) {
    int p;
    for (p = 0; p < NUM_IMP_DATA_TYPES; p++) {
      if (pMIIF->gdata[p] == NULL) {
        fhdls[p]     = NULL;
        gvd->hdls[p] = NULL;
        gvd->fn[p]   = NULL;
      } else {
        if (io->wsplt) {
          nptr = imp_data_names[p];
          snprintf(nbuf, sizeof nbuf, "%s.%s%s", base_name, nptr, ext);
          gvd->fn[p]    = dupstr(nbuf);
          gvd->dtype[p] = dtype_Real;
        } else {
          nptr = &outname[0];
          snprintf(nbuf, sizeof nbuf, "%s%s", base_name, ext);
        }

        if (io->wsplt || !p) {
          fhdls[p]     = fopen(nbuf, wmode);
          gvd->hdls[p] = f_open(nbuf, wmode);
        } else {
          // not splitting, all valid hdls are the same
          // everything ends up in the same file
          // fhdls[0] is Vm which is always present
          fhdls[p] = fhdls[0];
        }
      }
    }

    // time is not a gdata vector -> repeat the same crap as above
    if (io->wsplt) {
      snprintf(nbuf, sizeof nbuf, "%s.t%s", base_name, ext);
      gvd->fn[p]    = dupstr(nbuf);
      gvd->dtype[p] = dtype_Real;
      fhdls[p]      = fopen(nbuf, wmode);
      gvd->hdls[p]  = f_open(nbuf, wmode);
    } else {
      fhdls[p] = fhdls[0];
    }
  }

  if (!get_rank()) {
    fprintf(stderr, "Outputting the following quantities at each time: \n");
    fprintf(stderr, "%10s\t", "Time");
    for (int p = 0; p < NUM_IMP_DATA_TYPES; p++)
      if (pMIIF->gdata[p] != NULL)
        fprintf(stderr, "%10s\t", imp_data_names[p]);
    fprintf(stderr, "\n\n");
  }

  PetscBarrier(COMPAT_PETSC_NULLPTR);
} // open_globalvec_dump

void globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, MULTI_IF *pMIIF, timer_manager *tmo, IOCtrl *io, int numNode) {
  int p;
  float ft = tmo->time;
  double t = tmo->time;

  if (!tmo->trigger(CON_TM_IDX))
    return;
  else
    gvd->n_dumps++;

  if (io->w2stdout) {
    // write time only if we are not writing to file
    // this is Ed's preference
    log_msg(NULL, 0, NONL, "%10.3f ", ft);
  }

  if (!get_rank()) {
    if (io->w2file) {
      if (io->wbin)
        fwrite(&tmo->time, 1, data_type_sizes[gvd->dtype[NUM_IMP_DATA_TYPES]], gvd->hdls[NUM_IMP_DATA_TYPES]->fd);
      // fwrite(&tmo->tm, 1, data_type_sizes[gvd->dtype[NUM_IMP_DATA_TYPES]], fhdls[NUM_IMP_DATA_TYPES]);
      else
        fprintf(fhdls[NUM_IMP_DATA_TYPES], "%10.3f ", t);
    }
  } // globalvec_dump

  for (p = 0; p < NUM_IMP_DATA_TYPES; p++)
    if (pMIIF->gdata[p])
      write_global(pMIIF->gdata[p], gvd, fhdls[p], pMIIF, io, numNode);

  if (io->w2stdout)
    log_msg(NULL, 0, 0, "");


  if (!get_rank()) {
    // we write everything into one file
    if (!io->wsplt && io->w2file)
      fprintf(fhdls[0], "\n");
  }
}

void close_globalvec_dump(FILE **fhdls, GVEC_DUMP *gvd, IOCtrl *io) {
  int p;

  if (!io->w2file)
    return;

  if (!io->wsplt) {
    if (fhdls[0])
      fclose(fhdls[0]);
  } else {
    for (p = 0; p < NUM_IMP_DATA_TYPES+1; p++) {
      if (fhdls[p])
        fclose(fhdls[p]);

      f_close(gvd->hdls[p]);
    }
  }
}


/* dump all state variables for IMP and all PLUGS
 *
 * \param MIIF
 * \param reg   region number
 * \param imp   name of IMP
 * \param plugs plug-ins
 * \param t     start of dump time
 * \param dt    dump interval
 */
void dump_all(MULTI_IF *MIIF, int reg, char *imp, char *plugs, double t, double ddt, char *fout) {
  MIIF->svd.active = 1;
  MIIF->svd.intv   = ddt;
  MIIF->svd.t_dump = t;

  char **sv;
  char   file[1024];
  char   IMPS[2048];

  snprintf(IMPS, sizeof IMPS, "%s:%s", imp, plugs);  // make one string with all IMPS
  char *plgs = dupstr(IMPS), *p;
  while ( (p = get_next_list(plgs, ':')) != NULL) {
    int nsv = get_ion_type(std::string(plgs))->get_sv_list(&sv);
    for (int i = 0; i < nsv; i++) {
      snprintf(file, sizeof file, "%s_%s.%s", fout, plgs, sv[i]);
      MIIF->sv_dump_add_by_name(reg, plgs, sv[i], plgs, file);
    }
    free(sv);
    plgs = p;
  }
}

void initialize_timings(event_timing *t) {
  for (int i = 0; i < N_TIMINGS; i++) {
    t[i].mn    = 10000000.;
    t[i].mx    =        0.;
    t[i].avg   =        0.;
    t[i].tot   =        0.;
    t[i].count =        0;
  }
}

void update_timing(event_timing *t, double event_duration) {
  if (event_duration < t->mn)
    t->mn = event_duration;
  if (event_duration > t->mx)
    t->mx = event_duration;
  t->tot += event_duration;
  t->avg  = t->tot/++t->count;
}

/* Read value from global vector for a particular cell
 *
 * \param v     global vector
 * \param ind   index of cell for which we read value
 *
 * \return value of cell ind in vector v
 */
double getCellVal(sf_vec *v, int ind) {
  return v->get(ind);
}

/* initalize state variables
 *
 * The SVs string has the following format
 *
 * extern_var1=a:imp_sv0=b,imp_sv1=c::plug2_sv=d
 *
 * where extern_var is a global variable, eg. Vm, imp_sv0 and imp_sv1 are state variables in the IMP,
 * and plug2_sv is a state variable in the second plugin. Any number of value pairs can be added,
 * each separated by a ",".
 *
 * \param MIIF   ionic interface
 * \param SVs    state variables to change
 * \param imp    IMP being used
 * \param plgins colon spearated plug-in list
 * \param num    number of cells
 */
void initial_SVs(MULTI_IF *miif, char *SVs, char *imp, char *plgins, int num) {
  SF::vector<SF_int>  ind (num);
  SF::vector<SF_real> vals(num);

  char    impsv[4096];
  char   *colptr;
  int     mode     = 0;
  char   *last_col = SVs;

  do {
    colptr = strchr(last_col, ':');
    char *tmp;
    if (colptr)
      tmp = strndup(last_col, colptr-last_col);
    else
      tmp = strdup(last_col);
    char *iptr = NULL;
    char *pair = strtok_r(tmp, ",", &iptr);
    while (pair) {
      char sv[2048];
      double value;
      if (sscanf(pair, "%[^= ]=%lf", sv, &value) != 2) {
        log_msg(_nc_logf, 5, 0, "bad SV initializer: %s", pair);
        exit(1);
      }

      for (int i = 0; i < num; i++) {
        ind[i]  = i;
        vals[i] = value;
      }

      if (!mode) {
        strcpy(impsv, sv);
      } else if (mode == 1) {
        snprintf(impsv, sizeof impsv, "%s.%s", imp, sv);
      } else if (mode > 1) {
        char *ptr = plgins;
        for (int j = 0; j < mode-2; j++)
          ptr = strchr(ptr, ':') + 1;
        char *ptr1 = strchr(ptr, ':');
        snprintf(impsv, sizeof impsv, "%.*s.%s", (int)(ptr1 ? ptr1-ptr : strlen(ptr)), ptr, sv);
      }

      miif->adjust_MIIF_variables(impsv, ind, vals);
      pair = strtok_r(NULL, ",", &iptr);
    }
    free(tmp);
    last_col = colptr+1;
    mode++;
  } while (colptr != NULL);
} // initial_SVs

/* print out the parameter help
 *
 * \param
 */
void print_param_help(IonType* im, IonTypeList& plugs) {
  printf("Ionic model:\n");
  printf("------------\n");
  im->print_metadata();
  im->print_params();
  char **sv;
  int n = im->get_sv_list(&sv);
  printf("\tState variables:\n");
  for (int i = 0; i < n; i++)
    printf("\t\t%20s\n", sv[i]);

  if (!plugs.empty()) {
    printf("\n\n");
    printf("Plugins:\n");  
    printf("--------\n");
  }
  for (auto& plug : plugs) {
    plug.get().print_metadata();
    plug.get().print_params();
    n = plug.get().get_sv_list(&sv);
    printf("\tState variables:\n");
    for (int j = 0; j < n; j++)
      printf("\t\t%20s\n", sv[j]);
    printf("\n");
  }
}

float determine_duration(struct gengetopt_args_info *p, TrgList *stim_lst) {
  if (p->stim_times_given)
    determine_stim_list(p->stim_times_arg, stim_lst, p->DIA_flag);

  if (p->duration_given)
    return p->duration_arg;
  else if (p->restitute_given)
    return p->duration_arg;
  else if (p->stim_times_given) // user specified list
    return stim_lst->lst[stim_lst->n-1]+p->past_stim_arg;
  else
    return p->stim_start_arg+p->bcl_arg*(p->numstim_arg-1)+p->past_stim_arg;
}

}  // namespace limpet