CaParameters.cxx

Go to the documentation of this file.

/* Copyright (C) 2022-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergey Gorbunov, Sergei Zharko [committer] */

 
#include "CaParameters.h"
 
#include "AlgoFairloggerCompat.h"
 
#include <boost/filesystem.hpp>
 
#include <iomanip>
 
#include <fmt/format.h>
 
using cbm::algo::ca::Config;
using cbm::algo::ca::EDetectorID;
using cbm::algo::ca::Iteration;
using cbm::algo::ca::Parameters;
using cbm::algo::ca::kfutils::CheckSimdVectorEquality;
 
namespace fs = boost::filesystem;
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
Parameters<DataT>::Parameters(ca::Setup<DataT>&& setup, ca::Config&& config,
                              ca::SearchWindowMapContainer&& searchWindows)
  : fSetup(std::move(setup))
  , fConfig(std::move(config))
  , fSwMaps(std::move(searchWindows))
{
}
 
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
void Parameters<DataT>::CheckConsistency() const
{
  LOG(info) << "Consistency test for L1 parameters object... ";
  /*
   * Array of active station IDs
   *
   * In the array of active station IDs, 
   *   (i)   sum of elements, which are not equal -1, must be equal the number of stations,
   *   (ii)  subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers starting with 0
   */
 
  // TODO: rewrite or remove
  //auto filterInactiveStationIDs = [](int x) { return x != -1; };
  //int nStationsCheck = std::count_if(fvGeoToActiveMap.cbegin(), fvGeoToActiveMap.cend(), filterInactiveStationIDs);
  //if (nStationsCheck != fNstationsActiveTotal) {
  //  std::stringstream msg;
  //  msg << "ca::Parameters: invalid object condition: array of active station IDs is not consistent "
  //      << "with the total number of stations (" << nStationsCheck << " vs. " << fNstationsActiveTotal << ' '
  //      << "expected)";
  //  throw std::logic_error(msg.str());
  //}
 
  // Check, if the subsequence of all the elements, which are not equal -1, must be a gapless subset of integer numbers
  // starting from 0. If it is, the testValue in the end must be equal the number of non -1 elements
 
  // TODO: rewrite or remove
  //std::vector<int> idsCheck(
  //  nStationsCheck);  // temporary vector containing a sequence of active station IDs without -1 elements
  //std::copy_if(fvGeoToActiveMap.cbegin(), fvGeoToActiveMap.cend(), idsCheck.begin(), filterInactiveStationIDs);
  //bool isStationIDsOk = true;
  //for (int id = 0; id < nStationsCheck; ++id) {
  //  isStationIDsOk = isStationIDsOk && idsCheck[id] == id;
  //}
  //if (!isStationIDsOk) {
  //  std::stringstream msg;
  //  msg << "ca::Parameters: invalid object condition: array of active station IDs is not a gapless subset "
  //      << "of integer numbers starting from 0:\n\t";
  //  for (auto id : fvGeoToActiveMap) {
  //    msg << std::setw(3) << std::setfill(' ') << id << ' ';
  //  }
  //  throw std::logic_error(msg.str());
  //}
 
  /*
   * Check magnetic field flags of the stations
   *
   * In a current version of tracking there are three configurations possible to be proceeded:
   *  A. All the stations are inside magnetic field
   *  B. There is no magnetic field in a setup
   *  C. All the first stations are inside magnetic field, all the last stations are outside the field
   * In all the cases the fieldStatus flags should be sorted containing all non-zero elements in the beginning
   * (representing stations placed into magnetic field) and all zero elements in the end of z-axis.
   */
 
  // TODO: rewrite for the active setup
  //bool ifFieldStatusFlagsOk = std::is_sorted(fStations.cbegin(), fStations.cbegin() + fNstationsActiveTotal,
  //                                           [&](const ca::Station<DataT>& lhs, const ca::Station<DataT>& rhs) {
  //                                             return bool(lhs.fieldStatus) > bool(rhs.fieldStatus);
  //                                           });
 
  //if (!ifFieldStatusFlagsOk) {
  //  std::stringstream msg;
  //  msg << "ca::Parameters: invalid object condition: L1 tracking is impossible for a given field configuration:\n";
  //  for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
  //    msg << "- station ID:  " << iSt << ",  field status: " << fStations[iSt].fieldStatus << '\n';
  //  }
  //  throw std::logic_error(msg.str());
  //}
 
 
  /*
   * Check if each station object is consistent itself, and if all of them are placed after the target
   * NOTE: If a station was not set up, it is accounted inconsistent (uninitialized). In the array of stations there are uninitialized
   *       stations possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over all the stations in array
   *       but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
   * TODO: Probably, we should introduce methods, which check the consistency of fully initialized objects such as ca::Station,
   *       L1MaterialInfo, etc. (S.Zharko)
   */
 
  // TODO: Rewrite for active setup
  //{
  //  auto trgZ = kfutils::simd::Cast<DataT, float>(fActiveSetup.GetTarget().GetZ(), 0);
  //  for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
  //    fStations[iSt].CheckConsistency();
  //    auto zSta = kfutils::simd::Cast<DataT, float>(fStations[iSt].fZ, 0);
  //    if (zSta < trgZ) {
  //      std::stringstream msg;
  //      msg << "ca::Parameters: station with global ID = " << iSt << " is placed before target "
  //          << "(z_st = " << zSta << " [cm] < z_targ = " << trgZ << " [cm])";
  //      throw std::logic_error(msg.str());
  //    }
  //  }
  //}
 
  /*
   * Check thickness maps
   * NOTE: If a ca::MaterialMap map was not set up, it is accounted inconsistent (uninitialized). In the array of thickness maps for each 
   *       there are uninitialized elements possible (with id > NstationsActiveTotal), thus one should NOT run the loop above over 
   *       all the stations in array but only until *(fNstationsActive.cend() - 1) (== NstationsActiveTotal).
   */
  // TODO: Provide these checks in the kf::Setup
  //for (int iSt = 0; iSt < fNstationsActiveTotal; ++iSt) {
  //  fThickMap[iSt].CheckConsistency();
  //}
 
  /*
   *  Check equality of station indices in setups:
   */
  // TODO: rewrite
  //{
  //  if (fActiveSetup.GetNofLayers() != this->GetNstationsActive()) {
  //    std::stringstream msg;
  //    msg << "ca::Parameters: number of stations in active kf setup (" << fActiveSetup.GetNofLayers()
  //        << ") differes from the actual number of active tracking stations (" << GetNstationsActive() << ')';
  //    throw std::logic_error(msg.str());
  //  }
 
  //  if (fGeometrySetup.GetNofLayers() != this->GetNstationsGeometry()) {
  //    std::stringstream msg;
  //    msg << "ca::Parameters: number of stations in geometry kf setup (" << fGeometrySetup.GetNofLayers()
  //        << ") differes from the actual number of geometry tracking stations (" << GetNstationsGeometry() << ')';
  //    throw std::logic_error(msg.str());
  //  }
 
  //  {
  //    std::stringstream msg;
  //    msg << "ca::Parameters: inconsistency in geometry station indexing in kf-setup and old init:";
  //    bool bConsistent = true;
  //    for (int iStGeo = 0; iStGeo < GetNstationsGeometry(); ++iStGeo) {
  //      auto [detId, locId]           = GetStationIndexLocal(iStGeo);
  //      auto [detIdSetup, locIdSetup] = fGeometrySetup.GetIndexMap().template GlobalToLocal<EDetectorID>(iStGeo);
  //      if (detId != detIdSetup || locId != locIdSetup) {
  //        bConsistent = false;
  //      }
  //      msg << "\n- (old) detId = " << static_cast<int>(detId) << ", locId = " << locId
  //          << " ---- (kf) detId = " << static_cast<int>(detIdSetup) << ", locId = " << locIdSetup;
  //    }
  //    if (!bConsistent) {
  //      throw std::logic_error(msg.str());
  //    }
  //  }
 
  //  {
  //    std::stringstream msg;
  //    msg << "ca::Parameters: inconsistency in active station indexing in kf-setup and old init:";
  //    bool bConsistent = true;
  //    for (int iStGeo = 0; iStGeo < GetNstationsGeometry(); ++iStGeo) {
  //      auto [detId, locId] = GetStationIndexLocal(iStGeo);
  //      int iStActive       = GetStationIndexActive(locId, detId);
  //      if (iStActive < 0) {
  //        continue;
  //      }
  //      auto [detIdSetup, locIdSetup] = fActiveSetup.GetIndexMap().template GlobalToLocal<EDetectorID>(iStActive);
  //      if (detId != detIdSetup || locId != locIdSetup) {
  //        bConsistent = false;
  //      }
  //      msg << "\n- (old) detId = " << static_cast<int>(detId) << ", locId = " << locId
  //          << " ---- (kf) detId = " << static_cast<int>(detIdSetup) << ", locId = " << locIdSetup;
  //    }
  //    if (!bConsistent) {
  //      throw std::logic_error(msg.str());
  //    }
  //  }
  //}
 
 
  /*
   *  Check iterations sequence
   *  1. Number of iterations should be larger then zero
   *  2. Each iteration should contain values within predefined limits
   *  3. Number of iterations with TrackFromTriplets flag turned on no more then 1
   *  4. If the TrackFromTriplets iteration exists, it should be the last one in the sequence
   */
  // TODO: move this to the ca::Config validation
  {
    const auto& iterations = GetCAIterations();
    if (iterations.empty()) {
      std::stringstream msg;
      msg << "L1Parameters: 0 track finder iterations were found. Please, define at least one iteration";
      throw std::logic_error(msg.str());
    }
 
    std::string names = "";
    for (const auto& iter : iterations) {
      if (!iter.Check()) {
        names += iter.GetName() + " ";
      }
    }
    if (names.size()) {
      std::stringstream msg;
      msg << "L1Parameters: some parameters are out of range for the following iterations: " << names;
      throw std::logic_error(msg.str());
    }
 
    int nTrackFromTripletIterations = std::count_if(iterations.begin(), iterations.end(),
                                                    [=](const auto& it) { return it.GetTrackFromTripletsFlag(); });
    if (nTrackFromTripletIterations > 1) {
      std::stringstream msg;
      msg << "L1Parameters: found " << nTrackFromTripletIterations
          << " iterations with GetTrackFromTripletsFlag() == true:\n";
      for (const auto& iter : iterations) {
        if (iter.GetTrackFromTripletsFlag()) {
          msg << '\t' << "- " << iter.GetName() << '\n';
        }
      }
      msg << "Only the one iteration can have GetTrackFromTripletsFlag() == true, and this iteration should be ";
      msg << "the last one";
      throw std::logic_error(msg.str());
    }
 
    if (nTrackFromTripletIterations == 1 && !(iterations.end() - 1)->GetTrackFromTripletsFlag()) {
      std::stringstream msg;
      msg << "L1Parameters: iteration with GetTrackFromTripletsFlag() == true is not the last in a sequence. ";
      msg << "The GetTrackFromTripletsFlag() value in the iterations sequence: \n";
      for (const auto& iter : iterations) {
        msg << '\t' << "- " << std::setw(15) << std::setfill(' ') << iter.GetName() << ' ';
        msg << std::setw(6) << std::setfill(' ') << iter.GetTrackFromTripletsFlag() << '\n';
      }
      throw std::logic_error(msg.str());
    }
  }
 
  LOG(info) << "Consistency test for L1 parameters object... \033[1;32mpassed\033[0m";
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
void Parameters<DataT>::Print(int /*verbosityLevel*/) const
{
  LOG(info) << ToString();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
template<typename DataT>
std::string Parameters<DataT>::ToString(int verbosity, int indentLevel) const
{
  using namespace constants;
  using std::setfill;
  using std::setw;
  std::stringstream msg{};
  if (verbosity < 1) {
    return msg.str();
  }
 
  constexpr char indentCh = '\t';
  std::string indent(indentLevel, indentCh);

  auto ClrFlg = [&](bool flag) { return (flag ? clrs::GN : clrs::RD); };
 
  msg << " ----- CA parameters list -----\n";
  msg << indent << clrs::CLb << "COMPILE TIME CONSTANTS:\n" << clrs::CL;
  msg << indent << indentCh << "Bits to code one station:           " << constants::size::StationBits << '\n';
  msg << indent << indentCh << "Bits to code one triplet:           " << constants::size::TripletBits << '\n';
  msg << indent << indentCh << "Max number of detectors:            " << constants::size::MaxNdetectors << '\n';
  msg << indent << indentCh << "Max number of stations:             " << constants::size::MaxNstations << '\n';
  msg << indent << indentCh << "Max number of triplets:             " << constants::size::MaxNtriplets << '\n';
  msg << indent << clrs::CLb << "GEOMETRY:\n" << clrs::CL;
  msg << indent << indentCh << clrs::CLb << "NUMBER OF STATIONS:\n" << clrs::CL;
  msg << indent << indentCh << indentCh << "Number of stations (Geometry): ";
  for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
    msg << setw(2) << this->GetNstationsGeometry(static_cast<EDetectorID>(iDet)) << ' ';
  }
  msg << " | total = " << setw(2) << this->GetNstationsGeometry();
  msg << '\n';
  msg << indent << indentCh << indentCh << "Number of stations (Active):   ";
  for (int iDet = 0; iDet < constants::size::MaxNdetectors; ++iDet) {
    msg << setw(2) << setfill(' ') << this->GetNstationsActive(static_cast<EDetectorID>(iDet)) << ' ';
  }
  msg << " | total = " << setw(2) << this->GetNstationsActive();
  msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Geometry station indices: ";
  for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
    bool isActive = fSetup.GeoToActStationId(iStGeo) != -1;
    msg << ClrFlg(isActive) << setw(3) << iStGeo << ' ';
  }
  msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Local station indices:    ";
  for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
    bool isActive       = fSetup.GeoToActStationId(iStGeo) != -1;
    auto [detId, locId] = fSetup.GeoToLocStationId(iStGeo);
    msg << ClrFlg(isActive) << setw(3) << locId << ' ';
  }
  msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Detector indices:         ";
  for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
    bool isActive       = fSetup.GeoToActStationId(iStGeo) != -1;
    auto [detId, locId] = fSetup.GeoToLocStationId(iStGeo);
    msg << ClrFlg(isActive) << setw(3) << static_cast<int>(detId) << ' ';
  }
  msg << '\n' << indent << indentCh << indentCh << clrs::CL << "Active station indices:   ";
  for (int iStGeo = 0; iStGeo < this->GetNstationsGeometry(); ++iStGeo) {
    int iStAct = fSetup.GeoToActStationId(iStGeo);
    msg << ClrFlg(iStAct != -1) << setw(3) << iStAct << ' ';
  }
  msg << clrs::CL << '\n';
 
  if (fConfig.DevUseParSearchWindows()) {
    msg << indent << "DOUBLET SEARCH WINDOWS:\n";
    for (int iSt = 1; iSt < fSetup.GetNofActStations(); ++iSt) {
      for (int iIter = 0; iIter < GetNcaIterations(); ++iIter) {
        const auto& iter = fConfig.GetIteration(iIter);
        for (int gap = 0; gap <= iter.GetMaxStationGap(); gap++) {
          msg << indent << "- station: " << iSt << ", iteration: " << iIter << ", gap: " << gap << '\n';
          msg << fSwMaps.DoubletSw(iIter, iSt, gap).ToString() << '\n';
        }
      }
    }
    msg << indent << "TRIPLET SEARCH WINDOWS:\n";
    for (int iSt = 1; iSt < fSetup.GetNofActStations(); ++iSt) {
      for (int iIter = 0; iIter < GetNcaIterations(); ++iIter) {
        const auto& iter = fConfig.GetIteration(iIter);
        for (int gapL = 0; gapL <= iter.GetMaxStationGap(); gapL++) {
          for (int gapM = 0; gapM <= iter.GetMaxStationGap(); gapM++) {
            msg << indent << "- station: " << iSt << ", iteration: " << iIter << ", gapL: " << gapL
                << ", gapM: " << gapM << '\n';
            msg << fSwMaps.TripletSw(iIter, iSt, gapL, gapM).ToString() << '\n';
          }
        }
      }
    }
  }
 
  msg << '\n';
  msg << fConfig.ToString(indentLevel) << '\n';
  return msg.str();
}
 
// ---------------------------------------------------------------------------------------------------------------------
//
void cbm::algo::ca::ParametersIO::Store(const Parameters<double>& parameters, const std::string& fileName)
{
  auto path = fs::absolute(fs::weakly_canonical(fileName));
  fs::create_directories(path.parent_path());
 
  std::ofstream ofs(path.string(), std::ios::binary);
  if (!ofs) {
    throw std::runtime_error(
      fmt::format("ca::ParametersIO::Store: failed opening file \"{}\" to store the setup", path.string()));
  }
  LOG(info) << "ca::ParametersIO: writing parameters into file \"" << path.string() << '\"';
  boost::archive::binary_oarchive oa(ofs);
  oa << parameters;
}
 
 
namespace cbm::algo::ca
{
  template class Parameters<fvec>;
  template class Parameters<float>;
  template class Parameters<double>;
}  // namespace cbm::algo::ca