ActivationTracker.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "cardiac/common/ActivationTracker.hpp"
#include "cardiac/problem/CardiacProblem.hpp"
#include "cardiac/solver/CardiacSolver.hpp"
#include "pde/assemblers/FE/FEIntegrator.hpp"
 
ActivationTracker::ActivationTracker(CardiacSolver* solver) :
  m_problem             (nullptr),
  m_geom                (nullptr),
  m_timeWriter          (nullptr),
  m_nbIterPostProcess   (1),
  m_activationTimeData  ({}),
  m_activationNotYetSeen({}),
  m_activationSeen      ({}),
  m_foundPeak           ({}),
  m_foundAPD50          ({}),
  m_foundAPD            ({}),
  m_activatedVolume     ({}),
  m_tissueVolume        (0.),
  m_nToBeSeen           (0)
{
 
  CEPS_ABORT_IF(ceps::isNullPtr(solver),
    "Cannot create post processing unit from nullptr solver link"
  );
 
  m_problem    = solver->getCardiacProblem();
  m_geom       = m_problem->getGeometry();
  m_timeWriter = solver->getTimeWriter();
 
  if (ceps::isValidPtr(m_problem->getParameters()))
    setupWithParameters(m_problem->getParameters());
 
  initializeActivationMap();
 
}
 
ActivationTracker::~ActivationTracker()
{
  for (CepsUInt i=0; i<m_activationTimeData.size(); i++)
  CEPS_SAYS("Fraction of activated tissue for activation data set #" << i << ": " 
     << m_activatedVolume[i]/m_tissueVolume);
}
 
void
ActivationTracker::update(CepsInt  iter,
                          CepsReal t,
                          DHVecPtr solution,
                          DHVecPtr prevSolution
)
{
 
  if (ceps::isProfiling())
    m_profiler.start("actTrackUpdate","updating activation data");
 
  if (m_activationTimeData.empty() or iter%m_nbIterPostProcess!=0)
    return;
 
  solution         ->getLocalData();
  prevSolution     ->getLocalData();
  m_activationTimes->getLocalData();
  m_peakValues     ->getLocalData();
  m_peakTimes      ->getLocalData();
  m_APD50          ->getLocalData();
  m_APD            ->getLocalData();
  CepsUInt count = 0;
 
  // Loop through parameter set in definition of activation times
  for (auto it: m_activationTimeData)
  {
    // Access data
 
    CepsSet<CepsUInt> tmpSetOfSeens;
    CepsReal          threshold = m_activationTimeData[it.first][_threshold];
    CepsReal          apdpctg   = m_activationTimeData[it.first][_apdpctg  ];
    CepsReal          minV4Peak = m_activationTimeData[it.first][_minV4Peak];
 
    // A vector to compute currently activated volume
    m_activated->getLocalData();
 
    // Activation time: check if threshold has been exceeded
    // beware i is already in the form PROBLEM_DIM * j  + s
    for (auto dofId : m_activationNotYetSeen[it.first])
      if ((*solution)[dofId] > threshold)
      {
        tmpSetOfSeens.insert(dofId);
        (*m_activationTimes)[dofId] = t;
        (*m_activated      )[dofId] = 1.;
      }
 
    // we have to seek peak and APD  just among
    // for dof which had pass activation
    for (auto dofId : m_activationSeen[it.first])
    {
      if (not m_foundPeak[it.first].contains(dofId))
      {
        // v+ > v_step and dv / dt <0
        if (((*solution)[dofId]>minV4Peak) and (*solution)[dofId] < (*prevSolution)[dofId])
        {
          m_foundPeak[it.first].insert(dofId);
          (*m_peakValues)[dofId] = (*solution)[dofId];
          (*m_peakTimes) [dofId] = t;
        }
      }
      if ((*solution)[dofId]<threshold)
      {
        (*m_activated)[dofId] = 0.;
      }
    }
 
    // we have to seek peak and APD among dof which had pass activation
    for (auto dofId : m_foundPeak[it.first])
    {
      auto testPercent = [&] (double x, double perc) {
        double thres = ((*m_peakValues)[dofId]-threshold)*(1-perc)+threshold;
        return (x < thres);
      };
      // we assume that m_apdpctg belongs to [0.5,1]
      if (not m_foundAPD50[it.first].contains(dofId) 
          and testPercent((*solution)[dofId],0.5))
      {
        (*m_APD50)[dofId] = t;
        m_foundAPD50[it.first].insert(dofId);
      }
 
      if (m_foundAPD50[it.first].contains(dofId) 
      and not m_foundAPD[it.first].contains(dofId) 
      and testPercent((*solution)[dofId],apdpctg))
      {
        (*m_APD)[dofId] = t;
        m_foundAPD[it.first].insert(dofId);
      }
    }
 
    m_activated->releaseLocalData();
    {
      auto fe = ceps::runtimeCast<FiniteElements*>(m_problem->getSpatialDiscretization());
      FEIntegrator<FEMassAssembler> itg(fe,false,m_problem->getTissueAttributes(),{m_problem->getUnknown(it.first)});
      m_activatedVolume[count] = itg.integrate(m_activated);
    }
    count++;
 
    // update of seens
    for (auto dofId: tmpSetOfSeens)
    {
      m_activationNotYetSeen[it.first].erase (dofId);
      m_activationSeen      [it.first].insert(dofId);
    }
  }
 
  solution         ->releaseLocalData();
  prevSolution     ->releaseLocalData();
  m_activationTimes->releaseLocalData();
  m_peakValues     ->releaseLocalData();
  m_peakTimes      ->releaseLocalData();
  m_APD50          ->releaseLocalData();
  m_APD            ->releaseLocalData();
 
  if (ceps::isProfiling())
    m_profiler.stop("actTrackUpdate");
}
 
void
ActivationTracker::writeActivationMap()
{
  if (m_activationTimeData.empty())
    return;
 
  CEPS_SAYS("Now saving activation maps");
 
  CepsString fileName = m_problem->getOutputFileBase()+"_APD";
 
  VtkWriter* w = ceps::getNew<VtkWriter>(m_problem->getSpatialDiscretization(),fileName,1,
                               m_problem->getOutputFormat(),
                               m_problem->writesGlobalIndices());
 
  CepsVector<Unknown*> us;
  CepsVector<CepsString> namesAct,namesPV,namesPT,namesAPD50,namesAPD;
  for (auto it: m_activationTimeData)
  {
    Unknown* u = m_problem->getUnknown(it.first);
    us        .push_back(u);
    namesAct  .push_back("Activation time (" + u->getName() + ")");
    namesPV   .push_back("Peak value ("      + u->getName() + ")");
    namesPT   .push_back("Peak time ("       + u->getName() + ")");
    namesAPD50.push_back("APD50 ("           + u->getName() + ")");
    namesAPD  .push_back("APD ("             + u->getName() + ")");
  }
 
  w->addScalarData(m_activationTimes,namesAct  ,us);
  w->addScalarData(m_peakValues     ,namesPV   ,us);
  w->addScalarData(m_peakTimes      ,namesPT   ,us);
  w->addScalarData(m_APD50          ,namesAPD50,us);
  w->addScalarData(m_APD            ,namesAPD  ,us);
  w->write();
 
  ceps::destroyObject(w);
 
}
 
CepsBool
ActivationTracker::allPointsHaveBeenActivated() const
{
  if (m_activationTimeData.empty())
    return false;
  return allSeen(m_activationSeen);
}
 
CepsBool
ActivationTracker::allPointsHaveBeenRepolarized() const
{
  if (m_activationTimeData.empty())
    return false;
  return allSeen(m_foundAPD);
}
 
CepsString
ActivationTracker::getByeByeMessage() const
{
  return "Activation data can be found in file " + m_problem->getOutputFileBase() + "_APD.*";
}
 
CepsVector<CepsReal>&
ActivationTracker::getActivatedVolume()
{
  return m_activatedVolume;
}
 
void
ActivationTracker::setupWithParameters(InputParameters* p)
{
  if (p->hasKey("post processing perdiod"))
    m_nbIterPostProcess = CepsUInt(p->getReal("post processing period")/m_problem->getTimeStep());
 
  if (p->hasKey("activation time data"))
  {
    CepsVector<CepsString> options = ceps::split(p->getString("activation time data"),",");
    for (CepsString opts: options)
    {
      std::stringstream ss(opts);
      CepsUnknownIndex uId = ceps::readInt(ss,
        "activation time data: unable to read unknown id from string " + opts
      );
      CepsReal3D adt = ceps::readVertex(ss,
        "unable to read activation time data (3 real numbers) from input string " + opts
      );
      CEPS_ABORT_IF(ceps::isNullPtr(m_problem->getUnknown(uId)),
        "activation time data : cannot find unknown with ID " + ceps::toString(uId)
      );
 
      m_activationTimeData.insert(std::make_pair(uId,adt));
    }
  }
}
 
void
ActivationTracker::initializeActivationMap()
{
 
  if (m_activationTimeData.empty())
    return;
 
  // allocation for saving action potential markers
  auto sd = m_problem->getSpatialDiscretization();
  m_activationTimes = sd->newDofHaloVector();
  m_activated       = sd->newDofHaloVector();
  m_peakValues      = sd->newDofHaloVector();
  m_peakTimes       = sd->newDofHaloVector();
  m_APD50           = sd->newDofHaloVector();
  m_APD             = sd->newDofHaloVector();
 
  // Fill with nans in order to indicate non computed values
  m_activationTimes->fill(std::nan(""));
  m_activated      ->fill(          0.);
  m_peakValues     ->fill(std::nan(""));
  m_peakTimes      ->fill(std::nan(""));
  m_APD50          ->fill(std::nan(""));
  m_APD            ->fill(std::nan(""));
  m_activatedVolume.resize(m_activationTimeData.size(),0.);
 
  // Build sets of nodes with seen/unseen activation
  // (all unseen for now)
  
  auto dofs = m_problem->getSpatialDiscretization()->getDistributedDofs()->getOwned();
  CepsUInt i=0;
  for (auto it: m_activationTimeData)
  {
    Unknown* u = m_problem->getUnknown(it.first);
    CepsSet<CepsDofGlobalIndex> indices;
    for (DegreeOfFreedom* dof: dofs)
      if (dof->getUnknown()==u)
        indices.insert(dof->getGlobalIndex());
    m_activationNotYetSeen[it.first] = indices;
    m_activationSeen      [it.first] = {};
    m_foundPeak           [it.first] = {};
    m_foundAPD50          [it.first] = {};
    m_foundAPD            [it.first] = {};
    m_nToBeSeen += indices.size();
    std::stringstream ss;
    ss << "relative activated volume for activation time data set #" << i;
    m_timeWriter->addCustomData(ss.str(),&(m_activatedVolume[i]));
    i++;
  }
 
 
  // Compute the tissue volume once and for all
  auto     fe        = ceps::runtimeCast<FiniteElements*>(sd);
  DHVecPtr tissue    = sd->newDofHaloVector();
  auto     hdofs     = m_problem->getSpatialDiscretization()->getDistributedDofs()->getHalo();
  auto     tAttrs    = m_problem->getTissueAttributes();
  CepsBool allTissue = tAttrs.empty() or (tAttrs.size()==1 and tAttrs.contains(-1));
 
  tissue->fill(0.);
  tissue->getLocalData();
  for (auto toParse: {dofs,hdofs})
    for (auto dof: toParse)
      if (dof->hasOneOfAttributes(tAttrs) or allTissue)
        (*tissue)[dof->getGlobalIndex()] = 1.;
  tissue->releaseLocalData();
  {
    FEIntegrator<FEMassAssembler> itg(fe,false,tAttrs,{m_problem->getTMVUnknowns()[0]});
    m_tissueVolume = itg.integrate(tissue);
  }
 
}
 
CepsBool
ActivationTracker::allSeen(const CepsMap<CepsUnknownIndex,CepsSet<CepsDofGlobalIndex>>& arrays) const
{
  CepsUInt nbSeen = 0;
  for (auto aset : arrays)
    nbSeen += aset.second.size();
 
  CepsInt myAllSeen = (nbSeen==m_nToBeSeen) ? 1 : 0;
  CepsInt gAllSeen  = 0;
  MPI_Allreduce(&myAllSeen,&gAllSeen,1,MPI_INT,MPI_SUM,ceps::getCommunicator());
  return (CepsUInt)gAllSeen == ceps::getGridSize();
}