leadfield.cc

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#include <unordered_set>
 
#include "leadfield.h"
#include "electrics.h"
#include "SF_init.h"
 
#ifdef WITH_CALIPER
#include "caliper/cali.h"
#else
#include "caliper_hooks.h"
#endif
 
namespace opencarp {
 
Leadfield::Leadfield()
{
  if(param_globals::lf_vmfile && std::string(param_globals::lf_vmfile).length() > 0)
    vm_file = param_globals::lf_vmfile;
  else
    vm_file = param_globals::vofile;
 
  ecg_timedt = param_globals::ecg_timedt;
 
  if(param_globals::lf_dir)
    lf_dir = param_globals::lf_dir;
}
 
Leadfield::~Leadfield()
{
  for(sf_vec* z : Z) delete z;
  if(Im) delete Im;
}
 
void Leadfield::collect_lead_names(const SF::vector<stimulus> & stimuli)
{
  lead_names.clear();
  std::unordered_set<std::string> seen;
  for(const stimulus & s : stimuli)
    if(s.phys.type == LF_I && seen.insert(s.name).second)
      lead_names.push_back(s.name);
}
 
void Leadfield::write_leadfield_igb(const std::string & lead_name, sf_vec * phie_i)
{
  int rank = get_rank();
  set_dir(POSTPROC);
  std::string filename = "LF_" + lead_name + ".igb";
 
  FILE* f_out = nullptr;
  int err = 0;
  if(rank == 0) {
    f_out = fopen(filename.c_str(), "w");
    if(!f_out) {
      log_msg(0, 5, 0, "Error: Cannot write to %s", filename.c_str());
      err = 1;
    } else {
      IGBheader head;
      head.x(phie_i->gsize());
      head.y(1);
      head.z(1);
      head.t(1);
      head.dim_t(0.0);
      head.type(IGB_FLOAT);
      head.systeme("little_endian");
      head.unites_x("um");
      head.unites("mV");
      head.fileptr(f_out);
      head.write();
    }
  }
  if(get_global(err, MPI_SUM)) EXIT(1);
  phie_i->write_binary<float>(f_out);
 
  if(rank == 0) {
    fclose(f_out);
    log_msg(0, 0, 0, "  - Wrote %s (nodes=%zu)", filename.c_str(), phie_i->gsize());
  }
}
 
void Leadfield::apply_unit_curr_lead(const std::string & lead, Electrics & elec)
{
  elec.ellip_solver.phiesrc->set(0.0);
 
  bool stim_found = false;
  for(const stimulus & s : elec.stimuli) {
    if(s.name == lead && s.phys.type == LF_I) {
      // Normalise each electrode patch by its own enclosed volume so that
      // multi-electrode leads (e.g. WCT-referenced precordial leads) keep the
      // correct relative weighting. Scale a local copy to avoid mutating the
      // shared stimulus state.
      stimulus scaled = s;
      SF_real vol = get_volume_from_nodes(*elec.ellip_solver.mass_e, scaled.electrode.vertices);
      if(vol > 0.0)
        scaled.pulse.strength *= 1.e12 / vol;
      apply_stim_to_vector(scaled, *elec.ellip_solver.phiesrc, true);
      stim_found = true;
    }
  }
 
  if(stim_found) {
    elec.ellip_solver.mass_e->mult(*elec.ellip_solver.phiesrc, *elec.ellip_solver.currtmp);
    elec.ellip_solver.phiesrc->deep_copy(*elec.ellip_solver.currtmp);
  }
}
 
void Leadfield::construct(Electrics & elec)
{
  CALI_CXX_MARK_FUNCTION;
 
  if(elec.ellip_solver.phie_mat == nullptr) {
    elec.ellip_solver.init();
 
    set_elec_tissue_properties(elec.mtype, Electrics::extra_grid, NULL);
    region_mask(extra_elec_msh, elec.mtype[Electrics::extra_grid].regions,
                elec.mtype[Electrics::extra_grid].regionIDs, true, "gregion_e", elec.logger);
 
    set_elec_tissue_properties(elec.mtype, Electrics::intra_grid, NULL);
    region_mask(intra_elec_msh, elec.mtype[Electrics::intra_grid].regions,
                elec.mtype[Electrics::intra_grid].regionIDs, true, "gregion_i", elec.logger);
 
    elec.ellip_solver.rebuild_stiffness(elec.mtype, elec.stimuli, elec.logger);
    elec.ellip_solver.rebuild_mass(elec.logger);
  }
 
  elec.ellip_solver.phie->set(0.0);
  elec.ellip_solver.phie_i->set(0.0);
 
  collect_lead_names(elec.stimuli);
 
  int num_leads = lead_names.size();
  if(num_leads == 0) {
    log_msg(0, 4, 0, "Warning: No Leadfield stimuli found. Skipping.");
    return;
  }
 
  sf_mesh & imesh = get_mesh(intra_elec_msh);
  if(Z.empty()) {
    Z.resize(num_leads, nullptr);
    for(int i = 0; i < num_leads; i++)
      SF::init_vector(&Z[i], imesh, 1, sf_vec::algebraic);
  }
 
  SF::scattering* p2c = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
  SF::scattering* i2e = get_scattering(intra_elec_msh, extra_elec_msh, 1);
 
  log_msg(0, 0, 0, "Computing %d leadfield matrices...", num_leads);
 
  for(int lead_idx = 0; lead_idx < num_leads; lead_idx++) {
    std::string current_lead = lead_names[lead_idx];
    log_msg(0, 1, 0, "  - Solving lead %d/%d: %s", lead_idx + 1, num_leads, current_lead.c_str());
 
    apply_unit_curr_lead(current_lead, elec);
    (*elec.ellip_solver.lin_solver)(*elec.ellip_solver.phie, *elec.ellip_solver.phiesrc);
 
    if(elec.ellip_solver.lin_solver->reason < 0) {
      log_msg(0, 5, 0, "Error: Leadfield solver failed to converge for lead %s after %d iterations with residual %g! Reason: %s.", current_lead.c_str(),
              static_cast<int>(elec.ellip_solver.lin_solver->niter), elec.ellip_solver.lin_solver->final_residual,
              elec.ellip_solver.lin_solver->reason_str());
      EXIT(1);
    }
 
    i2e->backward(*elec.ellip_solver.phie, *elec.ellip_solver.phie_i);
 
    // Z must stay in solver (algebraic) ordering to match Im in apply();
    // phie_i is permuted to canonical only for the IGB write below.
    Z[lead_idx]->deep_copy(*elec.ellip_solver.phie_i);
 
    if(p2c) (*p2c)(*elec.ellip_solver.phie_i, true);
    write_leadfield_igb(current_lead, elec.ellip_solver.phie_i);
  }
  log_msg(0, 0, 0, "Leadfield computation complete.");
}
 
void Leadfield::load(Electrics & elec)
{
  collect_lead_names(elec.stimuli);
 
  int num_leads = lead_names.size();
  if(num_leads == 0) {
    log_msg(0, 5, 0, "Error: 'load' called but no leads defined in config!");
    EXIT(1);
  }
 
  sf_mesh & imesh = get_mesh(intra_elec_msh);
  if(Z.empty()) {
    Z.resize(num_leads, nullptr);
    for(int i = 0; i < num_leads; i++)
      SF::init_vector(&Z[i], imesh, 1, sf_vec::algebraic);
  }
 
  std::string dir = lf_dir;
  if(!dir.empty() && dir.back() != '/') dir += "/";
 
  size_t global_nodes = Z[0]->gsize();
 
  SF::scattering* p2c = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
 
  int err = 0;
  for(int i = 0; i < num_leads; i++) {
    std::string name = lead_names[i];
    std::string fpath = dir + "LF_" + name + ".igb";
 
    FILE* f = nullptr;
    if(get_rank() == 0) {
      f = fopen(fpath.c_str(), "r");
      if(!f) {
        log_msg(0, 5, 0, "Error: Cannot open file %s", fpath.c_str());
        err = 1;
      } else {
        IGBheader igb(f, true);
        if(igb.type() == 0) {
          log_msg(0, 5, 0, "Error: Corrupt or unrecognized IGB header in %s", fpath.c_str());
          fclose(f); f = nullptr; err = 1;
        } else if((size_t)igb.x() != global_nodes) {
          log_msg(0, 5, 0, "Error: Node count mismatch in %s! Expected %zu, got %d",
                  fpath.c_str(), global_nodes, igb.x());
          fclose(f); f = nullptr; err = 1;
        }
      }
    }
 
    if(get_global(err, MPI_SUM)) EXIT(1);
    size_t nread = Z[i]->read_binary<float>(f);
    if(get_rank() == 0) fclose(f);
    if(nread != Z[i]->gsize()) {
      log_msg(0, 5, 0, "Error: Incomplete read from %s! Expected %zu elements, got %zu.",
              fpath.c_str(), Z[i]->gsize(), nread);
      EXIT(1);
    }
    if(p2c) (*p2c)(*Z[i], false);
 
    log_msg(0, 0, 0, "  - Loaded %s", fpath.c_str());
  }
}
 
void Leadfield::apply(sf_mat & Ki, sf_vec & vm)
{
  if(Z.empty()) return;
 
  if(Im == nullptr)
    SF::init_vector(&Im, get_mesh(intra_elec_msh), 1, sf_vec::algebraic);
 
  Ki.mult(vm, *Im);
 
  // V_L(t) = Z_L . (Ki vm); dot() reduces across ranks.
  for(size_t i = 0; i < Z.size(); i++) {
    double ecg = Z[i]->dot(*Im);
    if(get_rank() == 0)
      trace_data[i].push_back((float)ecg);
  }
}
 
int Leadfield::run(Electrics & elec)
{
  if(!get_physics(elec_phys) || !mesh_is_registered(intra_elec_msh) || !mesh_is_registered(extra_elec_msh)) {
    log_msg(0, 5, 0, "Error: Leadfield requires active EP physics! Aborting.");
    EXIT(1);
  }
 
  elec.output_manager.close_files_and_cleanup();
  set_dir(POSTPROC);
  MPI_Barrier(SF_COMM);
 
  if(lf_dir.length() > 0) {
    load(elec);
  } else {
    log_msg(0, 0, 0, "Computing leadfields...");
    construct(elec);
  }
 
  int num_leads = lead_names.size();
  if(num_leads == 0) {
    log_msg(0, 5, 0, "No leads detected. Exiting.");
    EXIT(1);
  }
 
  if(get_rank() == 0) {
    trace_data.resize(num_leads);
    for(auto & vec : trace_data) vec.clear();
  }
 
  sf_vec* vm_step = NULL;
  SF::init_vector(&vm_step, get_mesh(intra_elec_msh), 1, sf_vec::algebraic);
  SF::scattering* sc = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
 
  FILE* fd = nullptr;
  int total_frames = 0;
  float file_dt = 0.0;
  int stride = 1;
  int err = 0;
 
  set_dir(OUTPUT);
 
  if(get_rank() == 0) {
    FILE_SPEC fs = f_open(vm_file.c_str(), "r");
    if(fs) {
      IGBheader head;
      head.fileptr(fs->fd);
      head.read();
      total_frames = head.t();
      file_dt = head.inc_t();
 
      if((size_t)head.x() != vm_step->gsize()) {
        log_msg(0, 5, 0, "Error: Vm dimension in %s does not match mesh. Aborting.", vm_file.c_str());
        err = 1;
      }
 
      if(file_dt > 0.0)
        stride = (int)((ecg_timedt / file_dt) + 0.5);
      if(stride < 1) stride = 1;
 
      log_msg(0, 1, 0, "  VM file: dt=%.4f ms, target=%.4f ms, stride=%d",
              file_dt, ecg_timedt, stride);
      fd = fs->fd;
      delete(fs);
    } else {
      log_msg(0, 5, 0, "Error: Could not open %s", vm_file.c_str());
      err = 1;
    }
  }
  if(get_global(err, MPI_SUM)) EXIT(1);
 
  set_dir(POSTPROC);
  MPI_Bcast(&total_frames, 1, MPI_INT,   0, SF_COMM);
  MPI_Bcast(&file_dt,      1, MPI_FLOAT, 0, SF_COMM);
  MPI_Bcast(&stride,       1, MPI_INT,   0, SF_COMM);
 
  if(total_frames == 0) {
    delete vm_step;
    return 0;
  }
 
  log_msg(0, 0, 0, "Computing ECG for %d leads, %d frames...", num_leads, total_frames);
  sf_mat & Ki = *elec.parab_solver.rhs_parab;
 
  for(int t = 0; t < total_frames; t++) {
    size_t nread = vm_step->read_binary<float>(fd);
    if(nread != vm_step->gsize()) {
      log_msg(0, 5, 0, "Error: %s read incomplete data slice from %s!", __func__, vm_file.c_str());
      EXIT(1);
    }
 
    bool forward = false;
    if(sc) (*sc)(*vm_step, forward);
 
    if(t % stride == 0)
      apply(Ki, *vm_step);
  }
 
  if(get_rank() == 0 && fd) fclose(fd);
  delete vm_step;
 
  int write_err = 0;
  if(get_rank() == 0) {
    for(size_t i = 0; i < lead_names.size(); i++) {
      std::string fname = "ECG_" + lead_names[i] + ".dat";
      FILE* f = fopen(fname.c_str(), "w");
      if(f) {
        for(size_t t = 0; t < trace_data[i].size(); t++) {
          double time = t * file_dt * stride;
          fprintf(f, "%g %g\n", time, trace_data[i][t]);
        }
        fclose(f);
      } else {
        log_msg(0, 5, 0, "Error: Cannot write %s", fname.c_str());
        write_err = 1;
      }
    }
  }
 
  if(!write_err)
    log_msg(0, 0, 0, "ECG post-processing complete.");
  return write_err;
}
 
} // namespace opencarp