BmnMSCMatrix.h

Go to the documentation of this file.

#ifndef BMNMSCMATRIX_H
#define BMNMSCMATRIX_H
 
#include <chrono>
#include <map>
#include <vector>
// Boost
#include <boost/dynamic_bitset.hpp>
// #include <boost/numeric/ublas/matrix_sparse.hpp>
// ROOT
// #include <TMatrixTSparse.h>
#include <TObject.h>
// FairRoot
#include <FairLogger.h>
// BmnRoot
#include "BmnAliases.h"
#include "BmnFunctionSet.h"
 
using SysClock = std::chrono::system_clock;
using SysPoint = std::chrono::time_point<SysClock>;
 
using std::string;
using std::vector;
 
// namespace ublas = boost::numeric::ublas;
 
// template class TMatrixTSparse<uint8_t>;
// template class ublas::compressed_matrix<uint8_t>;
 
/**
 * Storage for all MSC counters of the spill
 * @param times_count
 * @param scaler_count
 * it is assumed that time points entered by SetTime(...) are chronologically ordered
 */
template<typename IntType = uint8_t>
class BmnMSCMatrix : public TObject
{
    typedef std::pair<SysPoint, SysPoint> TimeRange;   // [begin,end) range
    typedef vector<TimeRange> TimeRangeList;
    //    typedef map<SysPoint, IntType> TimeLine;
    typedef std::pair<SysPoint, IntType> TimeLineEntry;   // <timeslice start, count>
    typedef vector<TimeLineEntry> TimeLine;
    template<typename TimeLocationPair>   // for TimeRange, TimeLineEntry
    struct LessThanForFirst
    {
        bool operator()(const TimeLocationPair& entry, const SysPoint& key) const { return entry.first < key; }
        bool operator()(const SysPoint& key, const TimeLocationPair& entry) const { return key < entry.first; }
    };
    template<typename TimeLocationPair>
    struct LessThanForSecond
    {
        bool operator()(const TimeLocationPair& entry, const SysPoint& key) const { return entry.second < key; }
        bool operator()(const SysPoint& key, const TimeLocationPair& entry) const { return key < entry.second; }
    };
 
    typedef std::pair<SysPoint, size_t> TimeMapEntry;
    //    typedef vector<TimeMapEntry> Time2Pos;
    typedef map<SysPoint, size_t> Time2Pos;
    typedef map<string, TimeLine> Scaler2TimeDist;
 
    typedef std::map<uint8_t, TimeLine> ScalerIdx2TimeDist;
 
  public:
    BmnMSCMatrix(size_t times_count = 0,
                 //            uint16_t scaler_count = 0,
                 NameChMap scalerName2Index = NameChMap{})
        : fNRecords{times_count}
        //        , fNScalers{scalerName2Index.size()}   //        , fVector(fNRecords * fNScalers, 0)
        //                , fMatrix{fNRecords, fNScalers, false}
        //                , fMatrix{(int32_t)fNRecords, fNScalers}
        , fScalerName2index{scalerName2Index}
        , fMapVector(GetNScalers())
        , fTimes(fNRecords, SysPoint::min())
    //        , fBitVecSpill(fNRecords, false)
    //        , fBitVecBusy(fNRecords, false)
    {}
 
    //    IntType& Counts(size_t i_record, uint16_t i_scaler) { return fVector[i_record * fNScalers + i_scaler]; }
    //        auto& Counts() { return fMatrix; }
 
    //    /**
    //     * map (time point -> vector position)
    //     * @return
    //     */
    //    Time2Pos& TimeMap() { return fTimeMap; }
 
    const IntType GetCount(size_t record_id, const string& scaler_name)
    {
        if (record_id >= GetNRecords())
            return 0;
        TimeLine& row = ScalerTimeLine(scaler_name);
        SysPoint& time = fTimes[record_id];
        return GetCntByTime(row, time);
    }
    const IntType GetCount(size_t record_id, uint16_t scaler_id)
    {
        if (record_id >= GetNRecords())
            return 0;
        TimeLine& row = ScalerTimeLine(scaler_id);
        SysPoint& time = fTimes[record_id];
        return GetCntByTime(row, time);
        //            return fMatrix(record_id, scaler_id);
    }
    const IntType GetCount(const SysPoint& time, const string& scaler_name)
    {
        auto it = fScalerName2index.find(scaler_name);
        if (it == fScalerName2index.end()) {
            throw std::out_of_range("Scaler not found");
        }
        return GetCount(time, it->second);
    }
    const IntType GetCount(const SysPoint& time, const uint16_t scaler_id)
    {
        TimeLine& row = ScalerTimeLine(scaler_id);
        return GetCntByTime(row, time);
    }
 
    void SetCount(size_t record_id, uint16_t scaler_id, IntType val)
    {
        TimeLine& row = ScalerTimeLine(scaler_id);
        SetCount(record_id, row, val);
        //        row[time] = val;
        //            fMatrix(record_id, scaler_id) = val;
    }
 
    void SetCount(size_t record_id, const string& scaler_name, IntType val)
    {
        TimeLine& row = ScalerTimeLine(scaler_name);
        SetCount(record_id, row, val);
        //        row[time] = val;
    }
 
    SysPoint& GetTime(const size_t& record_id) { return fTimes[record_id]; }
    TTimeStamp GetTS(const size_t& record_id) { return BmnFunctionSet::TimePoint2TS(fTimes[record_id]); }
 
    /** we assume the user is filling the timestamps in the chronological order
     * on the basis of preliminary created time map */
    void SetTime(const size_t& record_id, const SysPoint& time)
    {
        fTimes[record_id] = time;
        //        fTimeMap[time] = record_id;
    }
    void SetTS(const size_t& record_id, const TTimeStamp& ts) { SetTime(record_id, BmnFunctionSet::TimeStamp2TP(ts)); }
 
    NameChMap& GetScalerName2indexMap() { return fScalerName2index; }
    //    NameChMap& ScalerName2indexMap() { return fScalerName2index; }
 
    TimeLine& ScalerTimeLine(const string& scaler_name)
    {
        auto it = fScalerName2index.find(scaler_name);
        if (it == fScalerName2index.end()) {
            throw std::out_of_range("Scaler name not found");
        }
        return ScalerTimeLine(it->second);
    }
    TimeLine& ScalerTimeLine(const uint16_t& scaler_id)
    {
        if (scaler_id >= fScalerName2index.size())
            throw std::out_of_range("Scaler index not found");
        return fMapVector[scaler_id];
    }
    auto GetScalerRange(const string& scaler_name,
                        const SysPoint& timeBegin = SysPoint::min(),
                        const SysPoint& timeEnd = SysPoint::max())
    {
        TimeLine& row = ScalerTimeLine(scaler_name);
        auto itStart = std::lower_bound(row.begin(), row.end(), timeBegin, LessThanForFirst<TimeLineEntry>{});
        auto itFinish = std::upper_bound(row.begin(), row.end(), timeEnd, LessThanForFirst<TimeLineEntry>{});
        return make_pair(itStart, itFinish);
    }
    /**
     * Returns the [begin,end) iterators for range of time point pairs representing
     * begin and end of a spill time interval overlapping with the given interval
     * @param timeBegin start time point
     * @param timeEnd end time point
     * @return
     */
    auto GetSpillRanges(const SysPoint& timeBegin = SysPoint::min(), const SysPoint& timeEnd = SysPoint::max())
    {
        return GetRanges(fSpillRangeVec, timeBegin, timeEnd);
    }
    /**
     * Returns the [begin,end) iterators for range of time point pairs representing
     * begin and end of a busy time interval overlapping with the given interval
     * @param timeBegin start time point
     * @param timeEnd end time point
     * @return
     */
    auto GetBusyRanges(const SysPoint& timeBegin = SysPoint::min(), const SysPoint& timeEnd = SysPoint::max())
    {
        return GetRanges(fBusyRangeVec, timeBegin, timeEnd);
    }
    // Scaler count between time points
    uint64_t GetCntInRange(const string& scaler_name,
                           const SysPoint& timeBegin = SysPoint::min(),
                           const SysPoint& timeEnd = SysPoint::max())
    {
        return GetCntInRange(ScalerTimeLine(scaler_name), timeBegin, timeEnd);
    }
 
    size_t GetRecordIdByTime(const SysPoint& timePoint)
    {
        auto it = std::lower_bound(fTimes.begin(), fTimes.end(), timePoint);
        if (it == fTimes.end())
            throw std::out_of_range("Time point not found");
        return std::distance(it, fTimes.begin());
        //        auto it = fTimeMap.find(timePoint);
        //        if (it == fTimeMap.end())
        //            throw std::out_of_range("Time point not found");
        //        return it->second;
    }
 
    bool GetIsSpill(const SysPoint& time)
    {
        return IsInRange(fSpillRangeVec, time);
        //        size_t& record_id = GetRecordIdByTime(timePoint);
        //        return fBitVecSpill[record_id];
    }
 
    bool GetIsBusy(const SysPoint& time)
    {
        return IsInRange(fBusyRangeVec, time);
        //        size_t& record_id = GetRecordIdByTime(timePoint);
        //        return fBitVecBusy[record_id];
    }
 
    auto& SpillRangeVec() { return fSpillRangeVec; }
 
    auto& BusyRangeVec() { return fBusyRangeVec; }
 
    //    auto& IsSpillBitVec() { return fBitVecSpill; }
    //
    //    auto& IsBusyBitVec() { return fBitVecBusy; }
 
    const size_t GetNRecords() { return fNRecords; }
 
    const uint16_t GetNScalers() { return fScalerName2index.size(); }
    //    const uint16_t GetNScalers() { return fNScalers; }
 
    uint32_t GetMissingHits() { return fMissingHits; }
 
    void SetMissingHits(uint32_t val) { fMissingHits = val; }
 
    /** Destructor **/
    virtual ~BmnMSCMatrix() {}
 
  private:
    static bool IsInRange(const TimeRangeList& rangeList, const SysPoint& time)
    {
        auto it = std::lower_bound(
            rangeList.begin(), rangeList.end(), time,
            [](const TimeRange& rangePair, const SysPoint& timePoint) { return rangePair.second < timePoint; });
        if (it == rangeList.end())
            return false;
        //                printf("time \n\t%s \n range\n\t%s\n\t%s\n\t%d\n",
        //                        BmnFunctionSet::TimePoint2String(time).data(),
        //                        BmnFunctionSet::TimePoint2String(it->first).data(),
        //                        BmnFunctionSet::TimePoint2String(it->second).data(),
        //                        ((time >= it->first)&&(time < it->second)));
        return ((time >= it->first) && (time < it->second));
    }
    /**
     * Returns the [begin,end) iterators for range of time point pairs representing
     * begin and end of a specific time interval overlapping with the given interval
     * @param rangeList the vector of time point pairs to work with
     * @param timeBegin start time point
     * @param timeEnd end time point
     * @return
     */
    static auto GetRanges(const TimeRangeList& rangeList,
                          const SysPoint& timeBegin = SysPoint::min(),
                          const SysPoint& timeEnd = SysPoint::max())
    {
        auto itBeginByFirsts =
            std::lower_bound(rangeList.begin(), rangeList.end(), timeBegin, LessThanForFirst<TimeRange>{});
        auto itBeginBySeconds =
            std::lower_bound(rangeList.begin(), rangeList.end(), timeBegin, LessThanForSecond<TimeRange>{});
        auto& itBegin = std::distance(itBeginByFirsts, itBeginBySeconds) > 0 ? itBeginByFirsts : itBeginBySeconds;
 
        auto itEndByFirsts =
            std::upper_bound(rangeList.begin(), rangeList.end(), timeEnd, LessThanForFirst<TimeRange>{});
        auto itEndBySeconds =
            std::upper_bound(rangeList.begin(), rangeList.end(), timeEnd, LessThanForSecond<TimeRange>{});
        auto& itEnd = std::distance(itEndByFirsts, itEndBySeconds) > 0 ? itEndBySeconds : itEndByFirsts;
        return make_pair(itBegin, itEnd);
    }
    static IntType GetCntByTime(const TimeLine& row, const SysPoint& time)
    {
        auto it = std::lower_bound(row.begin(), row.end(), time, LessThanForFirst<TimeLineEntry>{});
        if (it->first != time) {
            LOGF(error, "Time point not found: %s", BmnFunctionSet::TimePoint2String(time).data());
            return 0;
        }
        //        auto it = row.find(time);
        if (it == row.end())
            return 0;
        return it->second;
    }
    static uint64_t GetCntInRange(const TimeLine& row,
                                  const SysPoint& timeBegin = SysPoint::min(),
                                  const SysPoint& timeEnd = SysPoint::max())
    {
        if (timeEnd < timeBegin)
            return 0;
        auto it = row.begin()->first > timeBegin
                      ? row.begin()
                      : std::lower_bound(row.begin(), row.end(), timeBegin, LessThanForFirst<TimeLineEntry>{});
        //            row.lower_bound(timeBegin);
        if (it == row.end())
            return 0;
        auto itFinish = std::upper_bound(row.begin(), row.end(), timeEnd, LessThanForFirst<TimeLineEntry>{});
        //        auto itFinish = row.upper_bound(timeEnd);
        uint64_t sum{0};
        while ((it != itFinish))
            sum += (it++)->second;
        return sum;
    }
 
    void SetCount(const size_t& record_id, TimeLine& row, const IntType& val)
    {
        SysPoint& time = fTimes[record_id];
        auto it = std::lower_bound(row.begin(), row.end(), time, LessThanForFirst<TimeLineEntry>{});
        if (it == row.end())
            row.push_back(make_pair(time, val));
        else
            it->second = val;
    }
 
    const size_t fNRecords;
    //    const uint16_t fNScalers;
    NameChMap fScalerName2index;
    vector<TimeLine> fMapVector;
    //    std::map<IntType> fVector;
    //    std::vector<IntType> fVector;
    //    Time2Pos fTimeMap;   // (time -> position in the fTimes vector) map
    //    ublas::compressed_matrix<IntType> fMatrix; // WTF ROOT cannot work with it in a macro
    //    TMatrixTSparse<uint8_t> fMatrix; // "template" only supports a couple of prebuilt types
    vector<SysPoint> fTimes;
 
    TimeRangeList fSpillRangeVec;
    TimeRangeList fBusyRangeVec;
    //    boost::dynamic_bitset<> fBitVecSpill;
    //    boost::dynamic_bitset<> fBitVecBusy;
    uint32_t fMissingHits;
 
    //    Scaler2TimeDist scalers;
    //    TimeDist ext_cond;
 
    ClassDef(BmnMSCMatrix, 4);
};
 
#endif /* BMNMSCMATRIX_H */