signals.h

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
// ----------------------------------------------------------------------------

#ifndef SIGNAL
#define SIGNAL

#include <string>
#include <assert.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <sstream>

#include "SF_vector.h"
namespace opencarp {
namespace sig {

typedef double sReal;

// trailing underscore added to avoid clash with bidomain.h
// Stimulation.c needs to be modified to make use of traces as used in here
typedef enum ip_method { _LINEAR_IP, _NEAREST_IP } IpMeth_t;

template<typename V>
V mod(const V& a, const V& b)
{
    return (a % b + b) % b;
}

class time_trace_ffmt{
  public:
    bool num_samples_on;      
    bool header_on;           
    std::string header;       
    std::string col_sep;      
    int  fmt_id;              
    SF::vector<int> colw;    

    time_trace_ffmt() : time_trace_ffmt( true, false, "", " ", 1)
    {}

    time_trace_ffmt(bool num_samples_on_, bool header_on_, std::string header_, std::string col_sep_)
        : time_trace_ffmt(num_samples_on_, header_on_, header_, col_sep_, 1 )
    {}

    time_trace_ffmt(bool num_samples_on_, bool header_on_, std::string header_, std::string col_sep_, int fmt_id_)
        : num_samples_on(num_samples_on_), header_on(header_on_), header(header_), col_sep(col_sep_), fmt_id(fmt_id_)
    {
      // default column width is 8
      colw = {8, 8};
    }

    int parse_header(std::ifstream& inpf);
};

inline int time_trace_ffmt::parse_header(std::ifstream& inpf)
{
  int err = 0;

  // assume no comment, data section starts at the beginning
  std::streampos start_of_data = inpf.beg;

  // read in first line
  std::string line{};
  std::getline(inpf,line);

  // comment line?
  std::size_t found_hash = line.find_first_of("#");
  if(found_hash != std::string::npos)
  {
    // keep reading comment lines
    start_of_data = inpf.tellg();
  }
  else
  {
    // no header, check for number of samples
    sReal t, f;
    size_t items = sscanf(line.c_str(),"%lf %lf", &t, &f);

    if(items != 2)
    {
      // line holds number of samples
      int num_samples{0};
      size_t item = sscanf(line.c_str(),"%d", &num_samples);
      if(!item)
        err = -1;
      else  {
        num_samples_on = true;
        // put back
        //for(size_t i=0; i<line.size(); i++)
          //inpf.unget();
      }
    }
    else
    {
      // we found a time value pair
      num_samples_on = false;
    }

    // set stream pointer back to start of data section
    inpf.seekg(start_of_data,inpf.beg);
  }
  return err;
}

class time_trace{
  public:
    std::string label;      
    std::string t_unit;     
    std::string f_unit;     
    SF::vector<sReal> f;    
    SF::vector<sReal> t;    
    sReal dt;                
    sReal rtOff;             
    bool fade;              

    // constructors

    time_trace() : time_trace( 1000., 10.0e-3, 0.0)
    {}

    time_trace(sReal _dt) : time_trace(1000., _dt, 0.)
    {}

    time_trace(sReal _duration, sReal _dt) : time_trace(_duration, _dt, 0.)
    {}

    time_trace(sReal _duration, sReal _dt, sReal _rtOff, std::string _label) :
    time_trace(_duration, _dt, _rtOff, _label, "", "")
    {}

    time_trace(sReal _duration, sReal _dt, sReal _rtOff) : time_trace(_duration, _dt, _rtOff, "", "", "")
    {}

    time_trace(sReal _duration, sReal _dt, sReal _rtOff, std::string _label, std::string _t_unit, std::string _f_unit)
    {
      dt    = _dt;
      rtOff = _rtOff;
      int N = ceil(_duration/dt) + 1;

      f.assign(N, 0.0);
      t.resize(N);

      sReal ct = 0;
      for(size_t i=0; i<t.size(); i++) {
        t[i] = ct;
        ct += dt;
      }

      set_labels(_label, _t_unit, _f_unit);
    }

    time_trace& operator=(const time_trace& in)
    {
      label   =  in.label;
      t_unit  =  in.t_unit;
      f_unit  =  in.f_unit;
      f       =  in.f;
      t       =  in.t;
      dt      =  in.dt;
      rtOff   =  in.rtOff;
      fade    =  in.fade;

      return *this;
    }

    time_trace(const time_trace& a)
    {
      *this = a;
    }

    time_trace(const time_trace& a, sReal val) : time_trace(a)
    {
      // use standard copy constructor and initialize signal with constant value
      std::fill(f.begin(), f.end(), val);
    }

    std::string format_file_header(const time_trace_ffmt ffmt_spec);


    void operator *= (const time_trace& a);
    void operator *= (const sReal s);
    void operator += (const sReal a);
    void operator += (const time_trace& a);
    void operator <<= (int delta);
    void operator >>= (int delta);

    // member functions
    void set_labels(std::string _label)
    {
      label  = _label;
    }

    void set_labels(std::string _label, std::string _t_unit, std::string _f_unit)
    {
      label  = _label;
      setUnits(_t_unit, _f_unit);
    }

    void setUnits(std::string _t_unit, std::string _f_unit)
    {
      t_unit = _t_unit;
      f_unit = _f_unit;
    }

    void resize(int numSamples)
    {
      t.resize(numSamples);
      f.resize(numSamples);
    }

    sReal duration ()            
    {
      //return f.back();
      return t[t.size()-1];
    }

    size_t len() const          
    {
      return t.size();
    }

    sReal rta()                  
    {
      for(std::size_t i=0; i<len(); i++)
        t[i] += rtOff;

      rtOff = 0.;

      return t[0];
    }

    sReal zta()                  
    {
      rtOff = t[0];
      for(std::size_t i=0; i<len(); i++)
        t[i] -= rtOff;

      return -rtOff;
    }

    // signal processing

    SF::vector<sReal>& diff()
    {
      for(std::size_t i=0; i<len()-1; i++)
        f[i] = f[i+1] - f[i];

      // last sample has no difference
      f[len()-1] = 0.;

      return f;
    }

    SF::vector<sReal>& gradient()
    {
      auto dFdt = f;

      int N = len();
      int end = N-1;

      // take forward/backward difference at left/right edge
      dFdt[0] = (f[1] - f[0]) / (t[1] - t[0]);
      dFdt[end] = (f[end] - f[end-1]) / (t[end] - t[end-1]);

      // take centered differences on interior points
      for(int i=1; i<end-1; i++)
      {
        sReal fdF = (f[i  ] - f[i-1]) / (t[i  ] - t[i-1]);
        sReal bdF = (f[i+1] - f[i  ]) / (t[i+1] - t[i  ]);

        dFdt[i]  = (fdF + bdF) * 0.5;
      }

      // overwrite f
      f = dFdt;

      return f;
    }

    // time trace I/O
    int read_trace(const std::string fname);
    int read_trace(const std::string fname, bool unitize);

    int write_trace();
    int write_trace(const std::string fname);
    int write_trace(const std::string fname, time_trace_ffmt ffmt_spec);


    // shift operations
    size_t timeShift(sReal t_shift)
    {
      int dI = t_shift/dt;

      return dI;
    }


    // scaling functions

    void unitize()
    {
      assert(len()>0);

      sReal _mn = f[0];
      sReal _mx = f[0];

      for (std::size_t i=1;i<len();i++)  {
        if (f[i] > _mx)
          _mx = f[i];
        else
          if (f[i] < _mn)
            _mn = f[i];
      }

      //pair<SF::vector<int>::iterator, SF::vector<int>::iterator> mnmx;
      //eal __mx, __mn;
      //auto mnmx = std::minmax_element(f.begin(),f.end());
      //__mn = nmx.first;
      //__mx = mnmx.second;

      sReal unsc = fabs(_mx)>=fabs(_mn) ? 1./fabs(_mx):1./fabs(_mn);

      for (std::size_t i=0;i<len();i++)
        f[i] *= unsc;

    }


    // interpolation functions

    bool IsEquDistSampling()
    {
      sReal dt_o = t[1] - t[0];

      // differences in sampling should be less than 0.1%
      sReal r_err = 0.001*dt_o;

      for (std::size_t i=1; i<len()-1; i++)  {
        sReal dt_n = t[i+1] - t[i];
        if ( (dt_n - dt_o) > r_err)
          return false;
      }

      return true;
    }

    sReal fval_t(sReal _t)
    {
      if(_t < t[0])
        return fade ? 0. : f[0];

      if(_t > t.back())
        return fade ? 0. : f.back();

      int i = 0;
      while(t[i] <= _t) i++;

      sReal K = (f[i] - f[i-1]) / (t[i] - t[i-1]);
      sReal fval = f[i-1] + K * (_t - t[i-1]);

      return fval;
    }

    void resample(time_trace& trc)
    {
      for(std::size_t i=0; i<trc.len(); i++)
        trc.f[i] = fval_t(trc.t[i]);
    }

    void resample(time_trace& trc, IpMeth_t meth)
    {
      interp1(trc, meth);
    }

    void
    interp1(time_trace& trc, IpMeth_t meth)
    {

      // iterate over target array
      for (std::size_t i=0; i<trc.len(); i++)  {

        // data point to interpolate on
        sReal _t = trc.t[i];

        // are we within input x-range?
        if ( (_t < t[0]) || (_t  > t.back()) )
          continue;

        //binary search
        int i1 = 0;
        int i2 = len()-1;
        while(i2 > (i1+1)) {
          int mid = (i2+i1)/2;
          if (t[mid] <= _t) {
            i1 = mid;
          } else {
            i2 = mid;
          }
        }

        switch (meth)  {
          case _LINEAR_IP:
            trc.f[i] = f[i1] + (f[i2]-f[i1])/(t[i2]-t[i1])*(_t-t[i1]);
            break;
          case _NEAREST_IP:
            if ((trc.t[i]-t[i1]) >= (t[i2]-trc.t[i]))
              trc.f[i] = f[i1];
            else
              trc.f[i] = f[i2];
            break;
        }
      }
    }

};


inline void time_trace::operator *= (const time_trace& a)
{
  assert( len() == a.len() );

  for(size_t i=0; i<len(); i++)
    f[i] *= a.f[i];
}


inline void time_trace::operator *= (const sReal s)
{
  for(auto & v : f) v *= s;
}

inline void time_trace::time_trace::operator += (const sReal a)
{
  for(auto & v : f) v += a;
}


inline void time_trace::operator += (const time_trace& a)
{
  assert( len() == a.len() );

  for(size_t i=0; i<len(); i++)
    f[i] += a.f[i];
}

inline void time_trace::operator <<= (int delta)
{
  bool ring = delta < 0;  // ring shift?
  int  N    = len();

  if(ring)
    delta = -delta;

  // ring buffer shift
  SF::vector<sReal>tmp = f;
  for(int i=0; i<N; i++)
    f[mod((i-delta),N)] = tmp[i];

  // not ring, override
  if(!ring)
    for(std::size_t i=len()-delta;i<len();i++)
      f[i] = 0.;
}

inline void time_trace::operator >>= (int delta)
{
  bool ring = delta < 0;
  int  N    = len();

  if(ring)
    delta = -delta;

  // ring buffer shift
  SF::vector<sReal>tmp = f;
  for(int i=0; i<N; i++)
    f[mod((i+delta),N)] = tmp[i];

  // not ring, override
  if(!ring)
    for(int i=0;i<delta;i++)
      f[i] = 0.;
}

inline std::string time_trace::format_file_header(const time_trace_ffmt ffmt_spec)
{
    std::ostringstream sheader_;

    if(ffmt_spec.header_on)
    {
      // write trace header comment section, for now this is empty
      sheader_ << "# " << ffmt_spec.header << std::endl;

      // in this case we also print the format identifier
      sheader_ << "# version = " << ffmt_spec.fmt_id << std::endl;

      // write column table header
      sheader_ << "#"  << std::endl;
      sheader_ << "# TABLE_HEAD" << std::endl;
      sheader_ << "# A = time" << std::endl;
      sheader_ << "# B = " << label << std::endl;
      sheader_ << "# [" << std::setw(ffmt_spec.colw[0]) << "A" << "]";
      sheader_ <<   "[" << std::setw(ffmt_spec.colw[1]) << "B" << "]" << std::endl;
      sheader_ << "# [" << std::setw(ffmt_spec.colw[0]) << t_unit << "]";
      sheader_ <<   "[" << std::setw(ffmt_spec.colw[1]) << f_unit << "]" << std::endl;
      sheader_ << "#" << std::endl;
      sheader_ << "# DATA" << std::endl;
    }

    if(ffmt_spec.num_samples_on)
      // first line holds # of samples
      sheader_ << len() << std::endl;

    std::string sheader = sheader_.str();

    return sheader;
}


/*
sReal duration(const std::string fname)
{
  time_trace tmp;

  int err = tmp.read_trace(fname);

  return !err ? tmp.t.back() : -1;
}
*/

inline int time_trace::read_trace(std::string fname)
{
  int err = 0;

  std::ifstream inpf(fname);

  // If we couldn't open the output file stream for writing
  if (!inpf)
  {
      // Print an error and exit
      std::cerr << "Trace file " << fname << " could not be opened for reading!" << std::endl;
      err = -1;
  }
  else  {

    time_trace_ffmt tt_ffmt{};

    tt_ffmt.parse_header(inpf);

    // first line holds # of samples
    int numSamples{0};
    if(tt_ffmt.num_samples_on)
    {
      inpf >> numSamples;

      // adjust signal buffer size
      resize(numSamples);
    }

    // read time/sample pairs
    int i = 0;
    while( inpf >> t[i] >> f[i] )
      i++;

    if(tt_ffmt.num_samples_on)
    {
      if(i<numSamples || i<2)
        err = -2;
    }

    inpf.close();
  }

  if(err)
    return err;

  // fill trace data
  // get units from comments header, if available
  // set label from filename
  // resample file according to given dt

  return err;
}


inline int time_trace::read_trace(std::string fname, bool _unitize)
{
  int err = read_trace(fname);

  if(!err && _unitize)
    unitize();

  return err;
}

inline int time_trace::write_trace()
{
  return write_trace(label + ".trc");
}


inline int time_trace::write_trace(const std::string fname)
{
  time_trace_ffmt ffmt_spec;
  int err = write_trace(fname, ffmt_spec);

  return err;
}

inline int time_trace::write_trace(const std::string fname, time_trace_ffmt ffmt_spec)
{
  int err = 0;

  std::ofstream outf(fname);

  // If we couldn't open the output file stream for writing
  if (!outf) {
      // Print an error and exit
      std::cerr << "Trace file " << fname << " could not be opened for writing!" << std::endl;
      err = -1;
  }
  else {
    if(ffmt_spec.header_on)
    {
      // write trace file header
      outf << format_file_header(ffmt_spec); //ffmt_spec.header_string(len());
    }

    // write number of samples in file
    if(ffmt_spec.num_samples_on)
      outf << len() << std::endl;

    // data section
    // write time/sample pairs
    outf << std::fixed;
    for(std::size_t i=0; i<len(); i++)
    {
      outf << std::setw(ffmt_spec.colw[0]+3) << std::setprecision(3) << t[i] << ffmt_spec.col_sep;
      outf << std::setw(ffmt_spec.colw[1]+1) << std::setprecision(6) << f[i] << std::endl;
    }
    outf.close();
  }

  return err;
}

// define time_trace arithmetics
inline time_trace operator * (const time_trace& a, const time_trace& b)
{
  assert( a.len() == b.len() );

  time_trace _mult = a;

  for(std::size_t i=0; i<a.len(); i++)
    _mult.f[i] = a.f[i] * b.f[i];

  return _mult;
}

inline time_trace operator / (const time_trace& a, const time_trace& b)
{
  assert( a.len() == b.len() );

  time_trace _div = a;

  for(std::size_t i=0; i<a.len(); i++)
    _div.f[i] = a.f[i] / b.f[i];

  return _div;
}

inline time_trace operator + (const time_trace& a, const time_trace& b)
{
  assert( a.len() == b.len() );

  time_trace _add = a;
  for(std::size_t i=0; i<a.len(); i++)
    _add.f[i] = a.f[i] + b.f[i];

  return _add;
}

inline time_trace operator - (const time_trace& a, const time_trace& b)
{
  assert( a.len() == b.len() );

  time_trace _sub = a;
  for(std::size_t i=0; i<a.len(); i++)
    _sub.f[i] = a.f[i] - b.f[i];

  return _sub;
}

// truncated exponential pulse definition
class Pulse
{
  public:
  sReal   duration;            
  float  strength;            
  double d1;                  
  float  tau_edge;            
  float  tau_plat;            
  int    decay;               
  float  s2r;                 
  float  bias;                
  bool   constant;            

  // constructor
  Pulse()
  {
    duration =    1.0;
    strength =    1.0;
    d1       =    1.0;
    tau_edge =    0.1;
    tau_plat = 1000.0;
    decay    =    5.0;
    s2r      =    0.0;
    bias     =    0.0;
    constant =  false;
  }

  //waveForm defineWave(
};

// action potential foot definition
class APfoot
{
  public:
  sReal tau_f;       
  sReal A;           
  sReal B;           

  // constructor
  APfoot()
  {
    tau_f =   0.25; 
    A     =   1.00; 
    B     = -80.0;  
  }
};

// step function
class Step {
  public:
  sReal trig;
  bool rise;

  // constructor
  Step()
  {
    trig  = 0;
    rise = true;
  }
};

// sine waves
class sineWave {
  public:
  sReal   frq;         
  float  phase;       

  sineWave()
  {
    frq   = 1.;       
    phase = 0.;       
  }

  sineWave(sReal _frq, sReal _phase) : frq(_frq)
  {
    phase = _phase;
  }
};

// abstract trace sampling function
class sampleTraceFunc {
  public:
  virtual SF::vector<sReal> & sample (time_trace &trc) = 0;
};


class constFunc : public sampleTraceFunc
{
  public:
  sReal val;

  constFunc(sReal _val) : val(_val)
  {}

  virtual SF::vector<sReal> & sample (time_trace &trc)
  {
    for(std::size_t i=0; i<trc.len(); i++)
      trc.f[i] = val;

    return trc.f;
  }
};

class stepFunc : public sampleTraceFunc
{

  public:
  Step pars;

  stepFunc(Step _pars) : pars(_pars)
  {}

  virtual SF::vector<sReal> & sample (time_trace &trc)
  {
    sReal off = pars.rise ? 1. : 0.;
    sReal on  = pars.rise ? 0. : 1.;
    for(std::size_t i=0; i<trc.t.size(); i++)
      trc.f[i] = trc.t[i] < pars.trig ? on : off;

    return trc.f;
  }
};

class sineFunc : public sampleTraceFunc
{
  public:
  sineWave pars;

  sineFunc(sineWave _pars) : pars(_pars)
  {}

  virtual SF::vector<sReal> & sample (time_trace &trc)
  {
    sReal phase_shift = pars.phase/180*M_PI;
    for(std::size_t i=0; i<trc.len(); i++)
  trc.f[i] = sin(2*M_PI*pars.frq*trc.t[i] + phase_shift);

    return trc.f;
  }
};

class APfootFunc : public sampleTraceFunc
{
  public:
  APfoot pars;

  APfootFunc(APfoot _pars) : pars(_pars)
  {}

  virtual SF::vector<sReal> & sample (time_trace &trc)
  {
    for(std::size_t i=0; i<trc.len(); i++)
      trc.f[i] = pars.A*exp(trc.t[i]/pars.tau_f) + pars.B;

    return trc.f;
  }
};

class monophasicTruncExpFunc : public sampleTraceFunc
{
  public:
  Pulse pars;

  monophasicTruncExpFunc(Pulse _pars) : pars(_pars)
  {}

  virtual SF::vector<sReal> & sample (time_trace & trc)
  {
    double t0 = 0.;
    double t1 = t0 + pars.d1 - 5*pars.tau_edge;
    //double t2 = t0 + pars.d1;

    // first phase exponential rise and plateau
    time_trace zero(trc);
    //zero = 0.;
    time_trace phase_0 = zero;

    for(std::size_t i=0; i<trc.t.size(); i++)
    {
      sReal f0 = pars.tau_edge ? -expm1(-(trc.t[i]-t0)/pars.tau_edge):1.;
      sReal fp = exp(-(trc.t[i]-t0)/pars.tau_plat);

      phase_0.f[i]= f0*fp*pars.strength;
    }

    // end of plateau, exponential decay
    time_trace eps_0 = zero;
    time_trace eps_1 = zero;

    // step function
    Step stepPars;
    stepPars.trig = t1;
    stepPars.rise = false;

    // eps+, falling eps function 1 | 0
    stepFunc eps(stepPars);
    eps.sample(eps_0);

    // eps-, rising eps function 0 | 1
    eps.pars.rise = true;
    eps.sample(eps_1);

    // final falling phase
    time_trace phase_1 = zero;
    int  i1 = int(t1/trc.dt);
    sReal f0 = phase_0.f[i1];
    for(std::size_t i=i1; i<trc.len(); i++)
      phase_1.f[i] = f0*exp(-(trc.t[i]-t1)/pars.tau_edge);

    // assemble pulse components
    time_trace pulse  = phase_0 * eps_0 + phase_1 * eps_1;
    // this does not work, lvalue is time_trace &
    //trc = phase_0 * eps_0 + phase_1 * eps_1;
    trc = pulse;

    // add bias
    trc += pars.bias;

    return trc.f;
  }
};

class biphasicTruncExpFunc : public sampleTraceFunc
{
  public:
  Pulse pars;

  biphasicTruncExpFunc(Pulse _pars) : pars(_pars)
  {}

  virtual SF::vector<sReal> & sample (time_trace & trc)
  {
    // create first phase
    monophasicTruncExpFunc cdmFunc1(pars);
    cdmFunc1.sample(trc);

    // create second phase, determine tilt
    if(pars.s2r == 0.)
    {
      int ph1_i1 = int((pars.d1 - 5*pars.tau_edge)/trc.dt);
      pars.strength = -trc.f[ph1_i1];
    }
    else
      pars.strength = -pars.s2r*pars.strength;

    time_trace ph2(trc);
    //ph2 = 0.;
    monophasicTruncExpFunc cdmFunc2(pars);
    cdmFunc2.sample(ph2);

    // shift phase 2 relative to phase 1
    ph2 >>= (pars.d1 - 5*pars.tau_edge)/trc.dt;

    trc += ph2;

    return trc.f;
  }
};

}  // namespace sig
}  // namespace opencarp

#endif