NBinning.cxx

#include <string>
#include <vector>
#include "TAxis.h"
#include <TAttAxis.h>
#include <TObjArray.h>
#include "NBinningDef.h"
#include "NBinningPoint.h"
#include "NDimensionalExecutor.h"
#include "NUtils.h"
#include "NLogger.h"
 
#include "NBinning.h"
 
ClassImp(Ndmspc::NBinning);
 
namespace Ndmspc {
NBinning::NBinning() : TObject()
{
}
NBinning::NBinning(std::vector<TAxis *> axes) : TObject(), fAxes(axes)
{
 
  Initialize();
}
NBinning::NBinning(TObjArray * axes) : TObject()
{
 
  if (axes == nullptr) {
    NLogError("NBinning(TObjArray * axes) : axes is nullptr");
    return;
  }
 
  for (int i = 0; i < axes->GetEntriesFast(); i++) {
    TAxis * axis = dynamic_cast<TAxis *>(axes->At(i));
    if (axis) {
      fAxes.push_back((TAxis *)axis->Clone());
    }
    else {
      NLogError("NBinning: Axis %d is not a TAxis", i);
    }
  }
 
  Initialize();
}
NBinning::~NBinning()
{
 
  // Cleanup fDefinitions
  for (auto & def : fDefinitions) {
    delete def.second;
  }
  fDefinitions.clear();
 
  // Cleanup fAxes
  for (auto & axis : fAxes) {
    if (axis) {
      delete axis; // Delete the axis if it was created by this class
    }
  }
 
  delete fMap;
  delete fContent;
}
void NBinning::Reset()
{
  fMap->Reset();
  fContent->Reset();
  if (fPoint != nullptr) {
    fPoint->Reset();
  }
}
 
void NBinning::Initialize()
{
 
  if (fAxes.size() == 0) {
    NLogError("No axes provided");
    return;
  }
 
  int dim = fAxes.size();
  // Calculate maximyn=m of nbins
  int     nbinsMax     = 0;
  int     nContentDims = 0;
  Binning binningType;
  for (int i = 0; i < dim; i++) {
    int nbins = fAxes[i]->GetNbins();
    if (nbins > nbinsMax) {
      nbinsMax = nbins;
    }
    // NLogDebug("Axis %d: %s nbins=%d", i, fAxes[i]->GetName(), nbins);
    // if (nbins == 1 || axes[i]->IsAlphanumeric()) {
    TString axisname(fAxes[i]->GetName());
    bool    isUser = axisname.Contains("/U");
    if (isUser) {
      axisname.ReplaceAll("/U", "");
      fAxes[i]->SetName(axisname.Data());
    }
 
    if (isUser) {
      binningType = Binning::kUser;
      nContentDims++;
    }
    else if (fAxes[i]->IsAlphanumeric()) {
      // TODO: Optimized beiing for alphanumeric axis (maybe no need to fix)
      // binningType = Binning::kSingle;
      // nContentDims++;
      binningType = Binning::kMultiple;
      nContentDims += 3;
    }
    else {
      binningType = Binning::kMultiple;
      nContentDims += 3;
    }
    fBinningTypes.push_back(binningType);
    fAxisTypes.push_back(AxisType::kFixed);
  }
  int        dimBinning   = 4;
  Int_t *    nbinsBinning = new Int_t[dimBinning];
  Double_t * xminBinning  = new Double_t[dimBinning];
  Double_t * xmaxBinning  = new Double_t[dimBinning];
  // set first item nbinsBinning to dim min to zero and max to dim
 
  nbinsBinning[0] = dim;
  xminBinning[0]  = 0;
  xmaxBinning[0]  = dim;
  for (int i = 1; i < dimBinning; i++) {
    nbinsBinning[i] = nbinsMax;
    xminBinning[i]  = 0;
    xmaxBinning[i]  = nbinsMax;
  }
 
  if (fMap) {
    delete fMap;
  }
  fMap = new THnSparseI("ngntBinningMap", "NGnTree binning map", dimBinning, nbinsBinning, xminBinning, xmaxBinning);
  if (fMap == nullptr) {
    NLogError("Cannot create binnings !!!");
    return;
  }
  fMap->GetAxis(0)->SetNameTitle("dim", "dimension");
  fMap->GetAxis(1)->SetNameTitle("rebin", "rebins");
  fMap->GetAxis(2)->SetNameTitle("start", "rebins start");
  fMap->GetAxis(3)->SetNameTitle("bin", "bin id");
 
  for (int i = 0; i < dim; i++) {
    fMap->GetAxis(0)->SetBinLabel(i + 1, fAxes[i]->GetName());
  }
 
  int        dimBinningContent   = nContentDims;
  Int_t *    nbinsBinningContent = new Int_t[dimBinningContent];
  Double_t * xminBinningContent  = new Double_t[dimBinningContent];
  Double_t * xmaxBinningContent  = new Double_t[dimBinningContent];
  int        iContentDim         = 0;
  for (int i = 0; i < dim; i++) {
    if (fBinningTypes[i] == Binning::kSingle) {
      // NLogDebug("Binning %d: %d", i, axes[idim]->GetNbins());
      nbinsBinningContent[iContentDim] = fAxes[i]->GetNbins();
      xminBinningContent[iContentDim]  = 0;
      xmaxBinningContent[iContentDim]  = fAxes[i]->GetNbins();
      iContentDim++;
      NLogTrace("[S] Binning %d: %d", i, fAxes[i]->GetNbins());
    }
    else if (fBinningTypes[i] == Binning::kMultiple) {
      for (int j = 0; j < 3; j++) {
        nbinsBinningContent[iContentDim + j] = fAxes[i]->GetNbins();
        xminBinningContent[iContentDim + j]  = 0;
        xmaxBinningContent[iContentDim + j]  = fAxes[i]->GetNbins();
      }
      iContentDim += 3;
      NLogTrace("[M] Binning %d: %d", i, fAxes[i]->GetNbins());
    }
    else if (fBinningTypes[i] == Binning::kUser) {
      // NLogDebug("Binning %d: %d", i, axes[idim]->GetNbins());
      nbinsBinningContent[iContentDim] = fAxes[i]->GetNbins();
      xminBinningContent[iContentDim]  = 0;
      xmaxBinningContent[iContentDim]  = fAxes[i]->GetNbins();
      iContentDim++;
      NLogTrace("[U] Binning %d: %d", i, fAxes[i]->GetNbins());
    }
    else {
      NLogError("Invalid binning type %d", fBinningTypes[i]);
      return;
    }
  }
  fContent = new THnSparseI("ngntBinningContent", "NGnTree binning content", dimBinningContent, nbinsBinningContent,
                            xminBinningContent, xmaxBinningContent);
 
  int                      iContentAxis = 0;
  std::vector<std::string> types        = {"rebin", "start", "bin"};
  for (int i = 0; i < dim; i++) {
    std::string name  = fAxes[i]->GetName();
    std::string title = fAxes[i]->GetName();
    // NLogTrace("Binning %d: %d %d", i, axes[i]->GetNbins(), fBinningTypes[i]);
    if (fBinningTypes[i] == Binning::kSingle) {
      fContent->GetAxis(iContentAxis)->SetNameTitle(name.c_str(), title.c_str());
      iContentAxis++;
    }
    else if (fBinningTypes[i] == Binning::kMultiple) {
      for (int j = 0; j < 3; j++) {
        int         imod = j % 3;
        std::string n    = name + "_" + types[imod];
        std::string t    = title + " (" + types[imod] + ")";
        fContent->GetAxis(iContentAxis)->SetNameTitle(n.c_str(), t.c_str());
        iContentAxis++;
      }
      // NLogDebug("Binning %d: %d", i, axes[idim]->GetNbins());
    }
    else if (fBinningTypes[i] == Binning::kUser) {
      fContent->GetAxis(iContentAxis)->SetNameTitle(name.c_str(), title.c_str());
      iContentAxis++;
    }
  }
  for (int i = 0; i < fContent->GetNdimensions(); i++) {
    NLogTrace("Axis[fContent] %d: %s nbins=%d", i, fContent->GetAxis(i)->GetName(), fContent->GetAxis(i)->GetNbins());
  }
 
  fPoint = new NBinningPoint(this);
 
  delete[] nbinsBinning;
  delete[] xminBinning;
  delete[] xmaxBinning;
}
 
void NBinning::Print(Option_t * option) const
{
 
  if (fMap == nullptr) {
    NLogError("THnSparse fMap is null");
    return;
  }
 
  TString opt(option);
  opt.ToLower();
 
  NLogInfo("NBinning base axes:");
  for (size_t i = 0; i < fAxes.size(); i++) {
    NLogInfo("  Axis %d: name='%s' title='%s' nbins=%d min=%.3f max=%.3f", i, fAxes[i]->GetName(), fAxes[i]->GetTitle(),
             fAxes[i]->GetNbins(), fAxes[i]->GetXmin(), fAxes[i]->GetXmax());
  }
 
  if (opt.Contains("A")) {
    NLogInfo("NBinning map name='%s' title='%s'", fMap->GetName(), fMap->GetTitle());
    NLogInfo("    dimensions: %d", fMap->GetNdimensions());
    // loop over all axes and print name title
    for (int i = 0; i < fMap->GetNdimensions(); i++) {
      NLogInfo("    [%d] name='%s' title='%s' nbins=%d min=%.3f max=%.3f", i, fMap->GetAxis(i)->GetName(),
               fMap->GetAxis(i)->GetTitle(), fMap->GetAxis(i)->GetNbins(), fMap->GetAxis(i)->GetXmin(),
               fMap->GetAxis(i)->GetXmax());
    }
    NLogInfo("    filled bins = %lld", fMap->GetNbins());
    NLogInfo("NBinning content name='%s' title='%s'", fContent->GetName(), fContent->GetTitle());
    NLogInfo("    dimensions: %d", fContent->GetNdimensions());
    // loop over all axes and print name title
    for (int i = 0; i < fContent->GetNdimensions(); i++) {
      NLogInfo("    [%d] name='%s' title='%s' nbins=%d min=%.3f max=%.3f", i, fContent->GetAxis(i)->GetName(),
               fContent->GetAxis(i)->GetTitle(), fContent->GetAxis(i)->GetNbins(), fContent->GetAxis(i)->GetXmin(),
               fContent->GetAxis(i)->GetXmax());
    }
  }
  NLogInfo("NBinning content:");
  NLogInfo("    filled bins = %lld", fContent->GetNbins());
  // print list of definitions names from fDefinitions
  NLogInfo("NBinning definitions:");
  for (const auto & kv : fDefinitions) {
    if (kv.first == fCurrentDefinitionName) {
      NLogInfo("  '%s' (current)", kv.first.c_str());
    }
    else {
      NLogInfo("  '%s'", kv.first.c_str());
    }
  }
 
  if (option && opt.Contains("all")) {
    NLogInfo("NBinning content:");
    PrintContent(option);
  }
 
  // // loop over definition and print
  // for (const auto & kv : fDefinition) {
  //   NLogInfo("Binning '%s':", kv.first.c_str());
  //   for (const auto & v : kv.second) {
  //     NLogInfo("  %s", NUtils::GetCoordsString(v, -1).c_str());
  //   }
  // }
 
  // fMap->Print(option);
 
  // // Using executor print all bins
  // NDimensionalExecutor executor(fMap, true);
  // // NDimensionalExecutor executor(fMap, false);
  // auto                          binning_task = [this, &all_bins](const std::vector<int> & coords) {
  //   fMap->Print();
  //   NUtils::PrintPointSafe(coords, -1);
  // };
  //
  // executor.Execute(binning_task);
 
  // size_t num_threads_to_use = 0;
  // // set num_threads_to_use to the number of available threads
  // if (num_threads_to_use < 1) num_threads_to_use = std::thread::hardware_concurrency();
  // NLogInfo("Using %zu threads ...", num_threads_to_use);
  // std::vector<NThreadData> thread_data_vector(num_threads_to_use);
  // for (size_t i = 0; i < thread_data_vector.size(); ++i) {
  //   thread_data_vector[i].SetAssignedIndex(i);
  // }
  //
  // // 2. Define the lambda (type in signature updated)
  // auto parallel_task = [](const std::vector<int> & coords, NThreadData & thread_obj) {
  //   NUtils::PrintPointSafe(coords, thread_obj.GetAssignedIndex());
  //   thread_obj.Process(coords);
  //   // thread_obj.Print();
  // };
  //
  // // 3. Call ExecuteParallel (template argument updated)
  // executor.ExecuteParallel<NThreadData>(parallel_task, thread_data_vector);
}
void NBinning::PrintContent(Option_t * option) const
{
 
  std::map<std::string, std::vector<int>> binningMap;
 
  TString opt(option);
  // loop over all selected bins via ROOT iterarot for THnSparse
  THnSparse *                                     cSparse = fContent;
  Int_t *                                         bins    = new Int_t[cSparse->GetNdimensions()];
  Long64_t                                        linBin  = 0;
  std::unique_ptr<ROOT::Internal::THnBaseBinIter> iter{cSparse->CreateIter(true /*use axis range*/)};
  while ((linBin = iter->Next()) >= 0) {
    Double_t    v         = cSparse->GetBinContent(linBin, bins);
    Long64_t    idx       = cSparse->GetBin(bins);
    std::string binCoords = NUtils::GetCoordsString(NUtils::ArrayToVector(bins, cSparse->GetNdimensions()), -1);
    NLogTrace("Bin %lld: %f %s idx=%lld content=%f", linBin, v, binCoords.c_str(), idx, v);
 
    size_t iAxis   = 0;
    bool   isValid = false;
    // for (int i = 0; i < cSparse->GetNdimensions(); i += 3) {
    int index = 0;
    for (iAxis = 0; iAxis < fAxes.size(); iAxis++) {
      int min;
      int max;
      // Print type of binning
      // NLogTrace("Axis %d: %s [%s]", iAxis, fAxes[iAxis]->GetName(),
      //                fBinningTypes[iAxis] == Binning::kSingle ? "S" : "M");
      if (fBinningTypes[iAxis] == Binning::kSingle) {
        isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], 1, 1, bins[index], min, max);
        binCoords += TString::Format(" | (S%c) %s %d %d %d [%d,%d] [%d,%d]", GetAxisTypeChar(iAxis),
                                     fAxes[iAxis]->GetName(), 1, 1, bins[index], min, max, 1, fAxes[iAxis]->GetNbins())
                         .Data();
        binningMap[fAxes[iAxis]->GetName()] = {1, 1, bins[index]};
        index++;
      }
      else if (fBinningTypes[iAxis] == Binning::kMultiple) {
        isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], bins[index], bins[index + 1], bins[index + 2], min, max);
        binCoords +=
            TString::Format(" | (M%c) %s %d %d %d [%d,%d] [%d,%d]", GetAxisTypeChar(iAxis), fAxes[iAxis]->GetName(),
                            bins[index], bins[index + 1], bins[index + 2], min, max, 1, fAxes[iAxis]->GetNbins())
                .Data();
        binningMap[fAxes[iAxis]->GetName()] = {bins[index], bins[index + 1], bins[index + 2]};
        index += 3;
      }
      else if (fBinningTypes[iAxis] == Binning::kUser) {
        isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], fAxes[iAxis]->GetNbins(), 1, 1, min, max);
        binCoords += TString::Format(" | (U) %s %d %d %d [%d,%d] [%d,%d]", fAxes[iAxis]->GetName(),
                                     fAxes[iAxis]->GetNbins(), 1, 1, min, max, 1, fAxes[iAxis]->GetNbins())
                         .Data();
        // binningMap[fAxes[iAxis]->GetName()] = {1, 1, bins[index]};
        binningMap[fAxes[iAxis]->GetName()] = {fAxes[iAxis]->GetNbins(), 1, 1};
        index++;
      }
      else {
        NLogError("[NBinning::PrintContent] Unknown binning type [%d]", fBinningTypes[iAxis]);
        continue;
      }
      if (!isValid) {
        NLogError("Cannot get axis range for axis %d", iAxis);
        continue;
      }
    }
    if (!isValid) {
      // NLogError("Cannot get axis range for axis %d", iAxis);
      NLogError("Bin %lld: %s  [not valid for axis=%s]", linBin, binCoords.c_str(), fAxes[iAxis]->GetName());
      continue;
    }
    NLogDebug("Bin %lld: %s", linBin, binCoords.c_str());
  }
  delete[] bins;
}
 
Long64_t NBinning::FillAll(NBinningDef * def)
{
  Long64_t nBinsFilled = 0;
 
  // fContent->Reset();
 
  std::vector<int>                           mins(fMap->GetAxis(0)->GetNbins(), 1);
  std::vector<int>                           maxs(fMap->GetAxis(0)->GetNbins(), 0);
  std::vector<std::vector<std::vector<int>>> content(fMap->GetAxis(0)->GetNbins());
 
  // loop over all selected bins via ROOT iterarot for THnSparse
  THnSparse *                                     cSparse = fMap;
  Int_t *                                         p       = new Int_t[cSparse->GetNdimensions()];
  Long64_t                                        linBin  = 0;
  std::unique_ptr<ROOT::Internal::THnBaseBinIter> iter{cSparse->CreateIter(true /*use axis range*/)};
  // Print number of filled bins
  // NLogInfo("Filled bins = %lld", cSparse->GetNbins());
  while ((linBin = iter->Next()) >= 0) {
    Double_t v = cSparse->GetBinContent(linBin, p);
    // continue;
    int idx = p[0] - 1;
    NLogTrace("NBinning::FillAll: Bin %lld: %d %d %d %d type=%d content=%f", linBin, p[0], p[1], p[2], p[3],
              fBinningTypes[idx], v);
    if (fBinningTypes[idx] == Binning::kSingle) {
 
      content[idx].push_back({p[0], p[3]});
    }
    else if (fBinningTypes[idx] == Binning::kMultiple) {
      // NLogDebug("Binning %d: %d", i, axes[idim]->GetNbins());
      content[idx].push_back({p[0], p[1], p[2], p[3]});
    }
    else if (fBinningTypes[idx] == Binning::kUser) {
      // TODO: Fix nbins for user binning (it should ok now)
      content[idx].push_back({p[0], 1});
    }
    else {
      NLogError("NBinning::FillAll: Unknown binning type %d", fBinningTypes[idx]);
      continue;
    }
    maxs[idx] = maxs[idx] + 1;
  }
  delete[] p;
 
  Long64_t nTotalBins = 1;
  // loop over content vector and set axis types
  for (size_t i = 0; i < content.size(); i++) {
    if (content[i].size() == 0) {
      NLogWarning("NBinning::FillAll: No content for binning %zu", i);
      continue;
    }
 
    if (content[i].size() > 1) {
      TAxis * axis  = fAxes[i];
      fAxisTypes[i] = AxisType::kVariable;
      NLogTrace("NBinning::FillAll: Axis id=%d name=%s bins=%zu", i, axis->GetName(), content[i].size());
      // check if i is persent in variable axes of def
      bool found = false;
      for (size_t j = 0; j < def->GetVariableAxes().size(); j++) {
        if (def->GetVariableAxes()[j] == static_cast<int>(i)) {
          NLogTrace("NBinning::FillAll: Axis id=%d name=%s already in variable axes of def", i, axis->GetName());
          found = true;
          break;
        }
      }
      NLogTrace("NBinning::FillAll: Axis id=%d name=%s set to variable found=%d", i, axis->GetName(), found);
      if (!found) def->AddVariableAxis(i);
    }
    nTotalBins *= content[i].size();
  }
 
  def->Print();
  NLogDebug("NBinning::FillAll: Filling total of %lld bins ...", nTotalBins);
 
  auto start_par = std::chrono::high_resolution_clock::now();
  // Loop over all binning combinations
  NDimensionalExecutor executor(mins, maxs);
  auto binning_task = [&content, &nBinsFilled, &nTotalBins, start_par, def, this](const std::vector<int> & coords) {
    std::vector<int> pointContentVector;
    NLogTrace("Binning task: %s", NUtils::GetCoordsString(coords, -1).c_str());
    for (size_t i = 0; i < coords.size(); i++) {
      NLogTrace("  Binning %zu: coord=%d content size=%zu", i, coords[i], content[i].size());
      if (content[i][coords[i] - 1].size() == 2) {
        pointContentVector.push_back(content[i][coords[i] - 1][1]);
      }
      else {
        pointContentVector.push_back(content[i][coords[i] - 1][1]);
        pointContentVector.push_back(content[i][coords[i] - 1][2]);
        pointContentVector.push_back(content[i][coords[i] - 1][3]);
      }
    }
 
    // NUtils::PrintPointSafe(pointContentVector, -1);
    Int_t nContentDims = fContent->GetNdimensions();
    NLogTrace("NBinning::FillAll: pointContentVector dims=%zu nContentDims=%d", pointContentVector.size(),
              nContentDims);
 
    auto pointContent = std::make_unique<Int_t[]>(nContentDims);
 
    // Int_t pointContent[nContentDims];
    NUtils::VectorToArray(pointContentVector, pointContent.get());
    Long64_t pointContentBin = fContent->GetBin(pointContent.get());
 
    // ids.push_back(pointContentBin);
    if (def) {
      NBinningPoint point(this);
      fContent->GetBinContent(pointContentBin, point.GetCoords());
      point.RecalculateStorageCoords();
      Long64_t linBin = def->GetContent()->GetBin(point.GetStorageCoords());
      def->GetContent()->SetBinContent(linBin, pointContentBin);
      def->GetIds().push_back(pointContentBin);
    }
 
    NLogTrace("NBinning::FillAll: Setting content bin %lld", pointContentBin);
 
    fContent->SetBinContent(pointContentBin, 1);
    nBinsFilled++;
    // NLogDebug("NBinning::FillAll: Filled bin %lld: %s", nBinsFilled,
    //                NUtils::GetCoordsString(pointContentVector, -1).c_str());
    NLogTrace("NBinning::FillAll: Filled bin %lld: %lld", nBinsFilled, nTotalBins);
    int refreshRate = nTotalBins / 100;
    if (refreshRate == 0) refreshRate = nTotalBins;
    if (nBinsFilled % (refreshRate) == 0) Ndmspc::NUtils::ProgressBar(nBinsFilled, nTotalBins, start_par, "I    ");
    // NLogDebug("NBinning::FillAll: [%3.2f%%] nBinsFilled=%lld", (double)nBinsFilled / nTotalBins * 100,
    //                nBinsFilled);
  };
  executor.Execute(binning_task);
 
  auto                                      end_par      = std::chrono::high_resolution_clock::now();
  std::chrono::duration<double, std::milli> par_duration = end_par - start_par;
 
  NLogInfo("NBinning::FillAll: Filled %lld bins in %s s", nTotalBins,
           NUtils::FormatTime(par_duration.count() / 1000).c_str());
 
  fMap->Reset();
 
  return nBinsFilled;
}
 
bool NBinning::AddBinning(size_t id, std::vector<int> binning, size_t n)
{
 
  if (fAxes.size() == 0) {
    NLogError("AddBinning: No axes defined !!!");
    return false;
  }
 
  if (id <= 0 || id > fAxes.size()) {
    NLogError("AddBinning: Invalid binning id %d", id);
    return false;
  }
 
  binning.insert(binning.begin(), id);
  Int_t * point = new Int_t[fMap->GetNdimensions()];
  NUtils::VectorToArray(binning, point);
  if (binning.size() == 2) {
    point[3] = point[1];
    point[2] = 1;
    point[1] = 1;
  }
  for (size_t i = 0; i < n; i++) {
    NLogTrace("Adding binning %d: %d %d %d %d", i, point[0], point[1], point[2], point[3]);
    fMap->SetBinContent(point, 1);
    point[3] += 1;
  }
  delete[] point;
 
  return true;
}
bool NBinning::AddBinningVariable(size_t id, std::vector<int> mins)
{
  if (id <= 0 || id > fAxes.size()) {
    NLogError("AddBinningVariable: Invalid binning id %d", id);
    return false;
  }
 
  // lopp over all mins and set binning
  for (size_t i = 1; i < mins.size(); i++) {
 
    if (mins[i] < 1 || mins[i] > fAxes[id - 1]->GetNbins() + 1) {
      // NLogWarning("AddBinningVariable: Invalid binning value mins=%d", mins[i]);
      return true;
    }
    // if (mins[i] == fAxes[id - 1]->GetNbins() || i == mins.size() - 1) mins[i]++;
    int rebin = mins[i] - mins[i - 1];
    if (rebin < 1) {
      NLogError("AddBinningVariable: Invalid binning value rebin=%d", rebin);
      return false;
    }
 
    int rebin_start = mins[i - 1] % rebin;
    int bin         = mins[i] / rebin;
 
    // int rebin_start = ((mins[i - 1]) % rebin) + 1;
    // int bin         = (mins[i] / rebin);
    if (rebin == 1) {
      rebin_start = 1;
      bin         = mins[i - 1];
    }
    // rebin_start + 1 is just that it looks better in the output
    // NLogDebug("Adding partial binning %d: %d %d %d %d [%d,%d]", i - 1, id, rebin, rebin_start + 1, bin,
    //                mins[i - 1], mins[i]);
    AddBinning(id, {rebin, rebin_start, bin}, 1);
  }
 
  return true;
}
bool NBinning::AddBinningViaBinWidths(size_t id, std::vector<std::vector<int>> widths)
{
  // NLogTrace("Adding binning via bin widths for id=%d", id);
  std::vector<int> mins;
  // Print vector of widths
  mins.push_back(1);
  for (auto & w : widths) {
    TAxis * axis    = fAxes[id - 1];
    int     width   = w[0];
    int     nWidths = w.size() > 1 ? w[1] : axis->GetNbins() / width;
 
    // NLogDebug("  width=%d nWidths=%d last min=%d", width, nWidths, mins[mins.size() - 1]);
    for (int iWidth = 0; iWidth < nWidths; iWidth++) {
      if (mins[mins.size() - 1] > axis->GetNbins()) {
        // NLogWarning("NBinning::AddBinningViaBinWidths: Exceeding axis nbins %d", axis->GetNbins());
        break;
      }
      mins.push_back(mins[mins.size() - 1] + width);
      // NLogDebug("  width=%d nWidths=%d iWidth=%d -> min=%d", width, nWidths, iWidth, mins[mins.size() - 1]);
    }
  }
 
  NLogTrace("NBinning::AddBinningViaBinWidths: Adding binning via bin widths: %s",
            NUtils::GetCoordsString(mins, -1).c_str());
 
  return AddBinningVariable(id, mins);
}
 
std::vector<std::vector<int>> NBinning::GetCoordsRange(std::vector<int> c) const
{
  std::vector<std::vector<int>> coordsRange;
  std::vector<int>              ids;
  std::vector<int>              mins;
  std::vector<int>              maxs;
  // int                           iAxis = 0;
  bool isValid;
  int  index = 0;
  for (size_t iAxis = 0; iAxis < fAxes.size(); iAxis++) {
    // for (int i = 0; i < fContent->GetNdimensions(); i += 3) {
    int min;
    int max;
    if (fBinningTypes[iAxis] == Binning::kSingle) {
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], 1, 1, c[index], min, max);
      index++;
    }
    else if (fBinningTypes[iAxis] == Binning::kMultiple) {
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], c[index], c[index + 1], c[index + 2], min, max);
      index += 3;
    }
    else if (fBinningTypes[iAxis] == Binning::kUser) {
      // isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], 1, 1, c[index], min, max);
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], fAxes[iAxis]->GetNbins(), 1, 1, min, max);
      index++;
    }
    else {
      NLogError("NBinning::GetCoordsRange: Unknown binning type %d", fBinningTypes[iAxis]);
      continue;
    }
    if (!isValid) {
      // NLogError("Cannot get axis range for axis %d", iAxis);
      return {};
    }
    // print min max
    NLogDebug("Axis %d: %s [%d,%d]", iAxis, fAxes[iAxis]->GetName(), min, max);
    ids.push_back(iAxis); // +1 because id starts from 1
    mins.push_back(min);
    maxs.push_back(max);
  }
 
  coordsRange.push_back(ids);
  coordsRange.push_back(mins);
  coordsRange.push_back(maxs);
  return coordsRange;
}
 
Binning NBinning::GetBinningType(size_t i) const
{
  if (i >= fBinningTypes.size()) {
    NLogError("Invalid binning type %d", i);
    // return Binning::kSingle;
    return Binning::kUndefined;
  }
  return fBinningTypes[i];
}
std::vector<std::vector<int>> NBinning::GetAxisRanges(std::vector<int> c) const
{
  std::vector<std::vector<int>> axisRanges;
 
  bool isValid;
  int  index = 0;
  for (size_t iAxis = 0; iAxis < fAxes.size(); iAxis++) {
    // for (int i = 0; i < fContent->GetNdimensions(); i += 3) {
    int min;
    int max;
    if (fBinningTypes[iAxis] == Binning::kSingle) {
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], 1, 1, c[index], min, max);
      index++;
    }
    else if (fBinningTypes[iAxis] == Binning::kMultiple) {
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], c[index], c[index + 1], c[index + 2], min, max);
      index += 3;
    }
    else if (fBinningTypes[iAxis] == Binning::kUser) {
      isValid = NUtils::GetAxisRangeInBase(fAxes[iAxis], fAxes[iAxis]->GetNbins(), 1, 1, min, max);
      index++;
    }
    else {
      NLogError("[NBinning::GetAxisRanges] Unknown binning type [%d]", fBinningTypes[iAxis]);
      continue;
    }
    if (!isValid) {
      // NLogError("Cannot get axis range for axis %d", iAxis);
      return {};
    }
    // print min max
    NLogTrace("Axis %d: %s [%d,%d]", iAxis, fAxes[iAxis]->GetName(), min, max);
    axisRanges.push_back({(int)iAxis, min, max});
  }
 
  return axisRanges;
}
std::vector<int> NBinning::GetAxisBinning(size_t axisId, const std::vector<int> & c) const
{
  std::vector<int> binning;
 
  int index = 0;
  for (size_t iAxis = 0; iAxis < fAxes.size(); iAxis++) {
    // for (int i = 0; i < fContent->GetNdimensions(); i += 3) {
    // int min;
    // int max;
    if (fBinningTypes[iAxis] == Binning::kSingle) {
      if (iAxis == axisId) {
        binning.push_back(1);
        binning.push_back(1);
        binning.push_back(c[index]);
        break;
      }
      index++;
    }
    else if (fBinningTypes[iAxis] == Binning::kMultiple) {
      if (iAxis == axisId) {
        binning.push_back(c[index]);
        binning.push_back(c[index + 1]);
        binning.push_back(c[index + 2]);
        break;
      }
      index += 3;
    }
    else if (fBinningTypes[iAxis] == Binning::kUser) {
      if (iAxis == axisId) {
        binning.push_back(c[index]);
        binning.push_back(1);
        binning.push_back(1);
        break;
      }
      index++;
    }
    else {
      NLogError("[NBinning::GetAxisBinning] Unknown binning type [%d]", fBinningTypes[iAxis]);
      continue;
    }
  }
 
  return binning;
}
 
bool NBinning::GetAxisRange(size_t axisId, double & min, double & max, std::vector<int> c) const
{
  int    minBin  = 0;
  int    maxBin  = 0;
  Bool_t isValid = GetAxisRangeInBase(axisId, minBin, maxBin, c);
 
  min = fAxes[axisId]->GetBinLowEdge(minBin);
  max = fAxes[axisId]->GetBinUpEdge(maxBin);
  return isValid;
}
 
bool NBinning::GetAxisRangeInBase(size_t axisId, int & min, int & max, std::vector<int> c) const
{
  if (axisId >= fAxes.size()) {
    NLogError("Invalid axis id %d", axisId);
    return false;
  }
 
  // Generate c if not provided
 
  int    minBin  = 0;
  int    maxBin  = 0;
  Bool_t isValid = false;
  if (fBinningTypes[axisId] == Binning::kSingle) {
    isValid = NUtils::GetAxisRangeInBase(fAxes[axisId], 1, 1, c[0], minBin, maxBin);
  }
  else if (fBinningTypes[axisId] == Binning::kMultiple) {
    isValid = NUtils::GetAxisRangeInBase(fAxes[axisId], c[0], c[1], c[2], minBin, maxBin);
  }
 
  min = minBin;
  max = maxBin;
  return isValid;
}
 
TObjArray * NBinning::GenerateListOfAxes()
{
 
  // TODO: Verify memery leak and its posible consequence for user
 
  TObjArray * axesArray = new TObjArray();
 
  // loop over all axes and create TObjArray
  for (size_t i = 0; i < fAxes.size(); i++) {
    TAxis * axis    = (TAxis *)fAxes[i];
    TAxis * axisNew = (TAxis *)axis->Clone();
 
    std::string name = axis->GetName();
    NLogTrace("NBinning::GetListOfAxes(): Binning '%s': %d", name.c_str(), axis->GetNbins());
 
    auto bins = std::make_unique<double[]>(axis->GetNbins() + 1);
    // double bins[axis->GetNbins() + 1];
    int count     = 0;
    int iBin      = 1; // start from bin 1
    bins[count++] = axis->GetBinLowEdge(1);
 
    NBinningDef *                                        def        = GetDefinition();
    std::map<std::string, std::vector<std::vector<int>>> definition = def->GetDefinition();
    for (auto & v : definition.at(name)) {
      NLogTrace("  %s", NUtils::GetCoordsString(v, -1).c_str());
 
      int n = v.size() > 1 ? v[1] : axis->GetNbins() / v[0];
      for (int i = 0; i < n; i++) {
        iBin += v[0];
        bins[count++] = axis->GetBinLowEdge(iBin);
      }
    }
    // loop over bins and print
    for (int i = 0; i < count; i++) {
      NLogTrace("  %s: %d %f", axis->GetName(), i + 1, bins[i]);
    }
    axisNew->Set(count - 1, bins.get());
    axesArray->Add(axisNew);
  }
 
  return axesArray;
}
 
bool NBinning::SetAxisType(size_t axisId, AxisType at)
{
  if (axisId >= fAxes.size()) {
    NLogError("Invalid axis id %d", axisId);
    return false;
  }
  fAxisTypes[axisId] = at;
  return true;
}
 
AxisType NBinning::GetAxisType(size_t i) const
{
  if (i >= fAxisTypes.size()) {
    NLogError("Invalid axis type %d", i);
    return AxisType::kUndefined; // Undefined
  }
  return fAxisTypes[i];
}
char NBinning::GetAxisTypeChar(size_t i) const
{
  if (i >= fAxisTypes.size()) {
    NLogError("Invalid axis type %d", i);
    return 'X'; // Undefined
  }
  switch (fAxisTypes[i]) {
  case AxisType::kFixed: return 'F';
  case AxisType::kVariable: return 'V';
  case AxisType::kUndefined:
  default: return 'X';
  }
}
 
std::vector<int> NBinning::GetAxisIndexes(AxisType type) const
{
  std::vector<int> ids;
  for (size_t i = 0; i < fAxes.size(); i++) {
    if (fAxisTypes[i] == type) {
      NLogTrace("Axis %d: %s type=%d", i, fAxes[i]->GetName(), fAxisTypes[i]);
      ids.push_back(i);
    }
    else {
      NLogTrace("Axis %d: %s type=%d [not selected]", i, fAxes[i]->GetName(), fAxisTypes[i]);
    }
  }
  return ids;
}
 
std::vector<TAxis *> NBinning::GetAxesByType(AxisType type) const
{
  std::vector<TAxis *> axes;
  for (size_t i = 0; i < fAxes.size(); i++) {
    if (fAxisTypes[i] == type) {
      NLogTrace("Axis %d: %s type=%d", i, fAxes[i]->GetName(), fAxisTypes[i]);
      axes.push_back(fAxes[i]);
    }
    else {
      NLogTrace("Axis %d: %s type=%d [not selected]", i, fAxes[i]->GetName(), fAxisTypes[i]);
    }
  }
  return axes;
}
 
std::vector<int> NBinning::GetAxisIndexesByNames(std::vector<std::string> names) const
{
  std::vector<int> ids;
  for (auto & name : names) {
    for (size_t i = 0; i < fAxes.size(); i++) {
      if (fAxes[i]->GetName() == name) {
        NLogTrace("Axis %d: %s type=%d", i, fAxes[i]->GetName(), fAxisTypes[i]);
        ids.push_back(i);
        break; // break inner loop if found
      }
    }
  }
  return ids;
}
 
std::vector<std::string> NBinning::GetAxisNamesByType(AxisType type) const
{
  std::vector<std::string> names;
  for (size_t i = 0; i < fAxes.size(); i++) {
    if (fAxisTypes[i] == type) {
      NLogTrace("Axis %d: %s type=%d", i, fAxes[i]->GetName(), fAxisTypes[i]);
      names.push_back(fAxes[i]->GetName());
      continue;
    }
    else {
      NLogTrace("Axis %d: %s type=%d [not selected]", i, fAxes[i]->GetName(), fAxisTypes[i]);
    }
  }
  return names;
}
std::vector<std::string> NBinning::GetAxisNamesByIndexes(std::vector<int> ids) const
{
  std::vector<std::string> names;
  for (auto & id : ids) {
    // FIXME: change to size_t
    if (id < 0 || id >= (int)fAxes.size()) {
      NLogError("Invalid axis id %d", id);
      continue;
    }
    NLogTrace("Axis %d: %s type=%d", id, fAxes[id]->GetName(), fAxisTypes[id]);
    names.push_back(fAxes[id]->GetName());
  }
  return names;
}
 
NBinningDef * NBinning::GetDefinition(const std::string & name)
{
 
 
  if (name.empty()) {
    if (fCurrentDefinitionName.empty()) {
      return nullptr;
    }
    // check if current definition exists
    if (fDefinitions.find(fCurrentDefinitionName) == fDefinitions.end()) {
      NLogError("NBinning::GetDefinition: Current definition '%s' not found", fCurrentDefinitionName.c_str());
      return nullptr;
    }
    NLogTrace("NBinning::GetDefinition: Using current definition '%s'", fCurrentDefinitionName.c_str());
    return fDefinitions.at(fCurrentDefinitionName);
  }
 
  if (fDefinitions.find(name) == fDefinitions.end()) {
    NLogError("NBinning::GetDefinition: Definition '%s' not found", name.c_str());
    return nullptr;
  }
  NLogTrace("NBinning::GetDefinition: Using definition '%s'", name.c_str());
  fCurrentDefinitionName = name;
  return fDefinitions.at(name);
}
 
void NBinning::AddBinningDefinition(std::string name, std::map<std::string, std::vector<std::vector<int>>> binning,
                                    bool forceDefault)
{
  if (fDefinitions.find(name) != fDefinitions.end()) {
    NLogError("Binning definition '%s' already exists", name.c_str());
    return;
  }
 
  NLogInfo("NBinning::AddBinningDefinition: Adding binning definition '%s'", name.c_str());
  NBinningDef * def = new NBinningDef(name, binning, this);
  fDefinitionNames.push_back(name);
  fDefinitions[name]     = def;
  fCurrentDefinitionName = name;
  fMap->Reset();
  // loop over binning axes
  for (size_t i = 0; i < fAxes.size(); i++) {
    TAxis *     axis     = fAxes[i];
    std::string axisName = axis->GetName();
    // NLogDebug("Axis %zu: name='%s' title='%s' nbins=%d min=%.3f max=%.3f", i, axisName.c_str(),
    //                        axis->GetTitle(), axis->GetNbins(), axis->GetXmin(), axis->GetXmax());
    if (!binning[axisName].empty()) {
      AddBinningViaBinWidths(i + 1, binning[axisName]);
      SetAxisType(i, AxisType::kVariable);
      // GetDefinition()[axisName] = binning[axisName];
      // def->SetAxisDefinition(axisName, binning[axisName]);
    }
    else {
      AddBinning(i + 1, {axis->GetNbins(), 1, 1}, 1);
      // GetDefinition()[axisName] = {{axis->GetNbins()}};
      // def->SetAxisDefinition(axisName, {{axis->GetNbins()}});
    }
  }
 
  // Long64_t nStart = fContent->GetNbins();
  // def->SetStartId(nStart);
  // std::vector<Long64_t> entries;
  Long64_t nFilled = FillAll(def);
  // def->Print();
  NLogTrace("NBinning::AddBinningDefinition: Filled %lld bins for definition '%s'", nFilled,
            NUtils::GetCoordsString(def->GetIds(), -1).c_str());
  // def->SetNEntries(nFilled);
 
  if (forceDefault || fCurrentDefinitionName.empty()) {
    fCurrentDefinitionName = name;
    NLogTrace("Set default binning definition '%s'", fCurrentDefinitionName.c_str());
  }
  else {
    NLogTrace("Added binning definition '%s' [not default]", name.c_str());
  }
}
 
NBinningPoint * NBinning::GetPoint()
{
 
  if (fPoint == nullptr) {
    fPoint = new NBinningPoint(this);
  }
  return fPoint;
}
 
NBinningPoint * NBinning::GetPoint(size_t id, std::string binning)
{
 
  // TODO: Optimize speed here
 
  NBinningDef * def = GetDefinition(binning);
  if (def == nullptr) {
    NLogError("NBinning::FillPoint: Definition '%s' not found", binning.c_str());
    return nullptr;
  }
 
  if (fPoint == nullptr) {
    fPoint = new NBinningPoint(this);
  }
 
  Long64_t linBin = def->GetId(id);
  if (linBin < 0) {
    NLogError("NBinning::GetPoint: Invalid linear bin %lld for id=%d in definition '%s'", linBin, id, binning.c_str());
    return nullptr;
  }
  if (linBin >= fContent->GetNbins()) {
    NLogError("NBinning::GetPoint: Linear bin %lld out of range [0,%lld) for id=%d in definition '%s'", linBin,
              fContent->GetNbins(), id, binning.c_str());
    return nullptr;
  }
  fContent->GetBinContent(linBin, fPoint->GetCoords());
  fPoint->RecalculateStorageCoords(linBin);
 
  return fPoint;
}
 
bool NBinning::SetCfg(const json & cfg)
{
  GetPoint()->SetCfg(cfg);
 
  return true;
}
 
void NBinning::SetCurrentDefinitionName(const std::string & name)
{
 
  if (name.empty()) {
    if (!fDefinitionNames.empty()) {
      fCurrentDefinitionName = fDefinitionNames[0];
      NLogInfo("Binning definition name='' set to first available definition '%s'", fCurrentDefinitionName.c_str());
      return;
    }
    else {
      NLogError("Binning definition name='' and no definitions are available !!! Ignorred ...");
      return;
    }
  }
 
  if (fDefinitions.find(name) == fDefinitions.end()) {
    NLogError("Binning definition '%s' not found", name.c_str());
    return;
  }
 
  // check if name is already current
  if (fCurrentDefinitionName.compare(name) == 0) {
    NLogTrace("Binning definition '%s' is already current", name.c_str());
    return;
  }
 
  // Set current definition name
  fCurrentDefinitionName = name;
 
  // Reset Point
  if (fPoint != nullptr) {
    fPoint = new NBinningPoint(this);
  }
  fPoint->Reset();
}
 
} // namespace Ndmspc