text:
/*
 * PatternGenerator.cc
 *
 *  Created on: Nov 8, 2019
 *      Author: kbunkow
 */

#include "L1Trigger/L1TMuonOverlapPhase1/interface/Tools/PatternGenerator.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <boost/range/adaptor/reversed.hpp>

#include "TFile.h"
#include "TDirectory.h"

PatternGenerator::PatternGenerator(const edm::ParameterSet& edmCfg,
                                   const OMTFConfiguration* omtfConfig,
                                   GoldenPatternVec<GoldenPatternWithStat>& gps)
    : PatternOptimizerBase(edmCfg, omtfConfig, gps), eventCntPerGp(gps.size(), 0) {
  edm::LogImportant("l1tOmtfEventPrint") << "constructing PatternGenerator, type: "
                                         << edmCfg.getParameter<string>("patternGenerator") << std::endl;

  if (edmCfg.getParameter<string>("patternGenerator") == "patternGen")
    initPatternGen();
}

PatternGenerator::~PatternGenerator() {}

void PatternGenerator::initPatternGen() {
  //reseting the golden patterns
  unsigned int i = 0;
  for (auto& gp : goldenPatterns) {
    gp->setKeyNumber(i++);  //needed  if patterns were added

    if (gp->key().thePt == 0)
      continue;

    gp->reset();

    int statBinsCnt = 1024;  //gp->getPdf()[0][0].size() * 8; //TODO should be big enough to comprise the pdf tails
    gp->iniStatisitics(statBinsCnt, 1);  //TODO
  }

  edm::LogImportant("l1tOmtfEventPrint") << "PatternGenerator::initPatternGen():" << __LINE__
                                         << " goldenPatterns.size() " << goldenPatterns.size() << std::endl;

  //GoldenPatternResult::setFinalizeFunction(3); TODO why it was this one????
  // edm::LogImportant("l1tOmtfEventPrint") << "reseting golden pattern !!!!!" << std::endl;

  //setting all pdf to 1, this will cause that  when the OmtfProcessor process the input, the result will be based only on the number of fired layers,
  //and then the omtfCand will come from the processor that has the biggest number of fired layers
  //however, if the GoldenPatternResult::finalise3() is used - which just count the number of muonStubs (but do not check if it is valid, i.e. fired the pdf)
  // - the below does not matter
  for (auto& gp : goldenPatterns) {
    for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
      for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
        for (unsigned int iBin = 0; iBin < gp->getPdf()[iLayer][iRefLayer].size(); iBin++) {
          gp->pdfAllRef[iLayer][iRefLayer][iBin] = 1;
        }
      }
    }
  }

  //TODO uncomment if filling the ptDeltaPhiHist is needed
  /*  ptDeltaPhiHists.resize(2);
  for(unsigned int iCharge = 0; iCharge <= 1; iCharge++) {
    for(unsigned int iLayer = 0; iLayer < omtfConfig->nLayers(); ++iLayer) { //for the moment filing only ref layer, remove whe
      if(iLayer == 0 || iLayer == 2 || iLayer == 4 || iLayer == 6 || iLayer == 7 || iLayer == 10 || iLayer == 11  || iLayer == 16 || //refLayars
         iLayer == 1 || iLayer == 3 || iLayer == 5  ) //bending layers
      {
        ostringstream name;
        name<<"ptDeltaPhiHist_ch_"<<iCharge<<"_Layer_"<<iLayer;
        int phiFrom = -10;
        int phiTo   = 300; //TODO
        int phiBins = phiTo - phiFrom;

        if(iCharge == 1) {
          phiFrom = -300; //TODO
          phiTo = 10;
        }

        TH2I* ptDeltaPhiHist = new TH2I(name.str().c_str(), name.str().c_str(), 400, 0, 200, phiBins, phiFrom -0.5, phiTo -0.5);
        //cout<<"BinLowEdge "<<ptDeltaPhiHist->GetYaxis()->GetBinLowEdge(100)<<" BinUpEdge "<<ptDeltaPhiHist->GetYaxis()->GetBinUpEdge(100);
        ptDeltaPhiHists[iCharge].push_back(ptDeltaPhiHist);
      }
      else
        ptDeltaPhiHists[iCharge].push_back(nullptr);
    }
  }*/

  /* cannot be called  here, will cause crash
  edm::LogImportant("OMTFReconstruction")<<" PatternGenerator constructor - patterns after modification "<<std::endl;
  for(auto& gp : goldenPatterns) {
    edm::LogImportant("OMTFReconstruction")<<gp->key()<<" "
        <<omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom
        <<" - "<<omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo<<" GeV"<<std::endl;
  }*/
}

void PatternGenerator::updateStat() {
  //cout<<__FUNCTION__<<":"<<__LINE__<<" omtfCand "<<*omtfCand<<std::endl;;
  AlgoMuon* algoMuon = omtfCand.get();
  if (!algoMuon) {
    edm::LogImportant("l1tOmtfEventPrint") << ":" << __LINE__ << " algoMuon is null" << std::endl;
    throw runtime_error("algoMuon is null");
  }

  double ptSim = simMuon->momentum().pt();
  int chargeSim = (abs(simMuon->type()) == 13) ? simMuon->type() / -13 : 0;

  unsigned int exptPatNum = omtfConfig->getPatternNum(ptSim, chargeSim);
  GoldenPatternWithStat* exptCandGp = goldenPatterns.at(exptPatNum).get();  // expected pattern

  eventCntPerGp[exptPatNum]++;

  //edm::LogImportant("l1tOmtfEventPrint")<<"\n" <<__FUNCTION__<<": "<<__LINE__<<" exptCandGp "<<exptCandGp->key()<<" candProcIndx "<<candProcIndx<<" ptSim "<<ptSim<<" chargeSim "<<chargeSim<<std::endl;

  int pdfMiddle = 1 << (omtfConfig->nPdfAddrBits() - 1);

  //iRefHit is the index of the hit
  for (unsigned int iRefHit = 0; iRefHit < exptCandGp->getResults()[candProcIndx].size(); ++iRefHit) {
    auto& gpResult = exptCandGp->getResults()[candProcIndx][iRefHit];

    unsigned int refLayer = gpResult.getRefLayer();

    if (gpResult.getFiredLayerCnt() >= 3) {
      for (unsigned int iLayer = 0; iLayer < gpResult.getStubResults().size(); iLayer++) {
        //updating statistic for the gp which should have fired

        bool fired = false;
        if (gpResult.getStubResults()[iLayer].getMuonStub()) {
          if (omtfConfig->isBendingLayer(iLayer)) {
            if (gpResult.getStubResults()[iLayer].getMuonStub()->qualityHw >= 4)  //TODO change quality cut if needed
              fired = true;
          } else
            fired = true;
        }

        if (fired) {  //the result is not empty
          int phiDist = gpResult.getStubResults()[iLayer].getPdfBin();
          phiDist += exptCandGp->meanDistPhiValue(iLayer, refLayer) - pdfMiddle;
          //removing the shift applied in the GoldenPatternBase::process1Layer1RefLayer

          //TODO uncomment if filling ptDeltaPhiHists is needed
          /*
          unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefHit];
          if(ptDeltaPhiHists[iCharge][iLayer] != nullptr &&
              (iLayer == refLayerLogicNum || omtfConfig->getLogicToLogic().at(iLayer) == (int)refLayerLogicNum) )
            ptDeltaPhiHists[iCharge][iLayer]->Fill(ttAlgoMuon->getPt(), phiDist); //TODO correct
           */

          phiDist += exptCandGp->getStatistics()[iLayer][refLayer].size() / 2;

          //edm::LogImportant("l1tOmtfEventPrint")<<__FUNCTION__<<":"<<__LINE__<<" refLayer "<<refLayer<<" iLayer "<<iLayer<<" phiDist "<<phiDist<<" getPdfBin "<<gpResult.getStubResults()[iLayer].getPdfBin()<<std::endl;
          if (phiDist > 0 && phiDist < (int)(exptCandGp->getStatistics()[iLayer][refLayer].size())) {
            //updating statistic for the gp which found the candidate
            //edm::LogImportant("l1tOmtfEventPrint")<<__FUNCTION__<<":"<<__LINE__<<" updating statistic "<<std::endl;
            exptCandGp->updateStat(iLayer, refLayer, phiDist, 0, 1);
          }
        } else {  //if there is no hit at all in a given layer, the bin = 0 is filled
          int phiDist = 0;
          exptCandGp->updateStat(iLayer, refLayer, phiDist, 0, 1);
        }
      }
    }
  }
}

void PatternGenerator::observeEventEnd(const edm::Event& iEvent,
                                       std::unique_ptr<l1t::RegionalMuonCandBxCollection>& finalCandidates) {
  if (simMuon == nullptr || omtfCand->getGoldenPatern() == nullptr)  //no sim muon or empty candidate
    return;

  if (abs(simMuon->momentum().eta()) < 0.8 || abs(simMuon->momentum().eta()) > 1.24)
    return;

  PatternOptimizerBase::observeEventEnd(iEvent, finalCandidates);

  updateStat();
}

void PatternGenerator::endJob() {
  if (edmCfg.getParameter<string>("patternGenerator") == "modifyClassProb")
    modifyClassProb(1);
  else if (edmCfg.getParameter<string>("patternGenerator") == "groupPatterns")
    groupPatterns();
  else if (edmCfg.getParameter<string>("patternGenerator") == "patternGen") {
    upadatePdfs();
    writeLayerStat = true;
  } else if (edmCfg.getParameter<string>("patternGenerator") == "patternGenFromStat") {
    std::string rootFileName = edmCfg.getParameter<edm::FileInPath>("patternsROOTFile").fullPath();
    edm::LogImportant("l1tOmtfEventPrint") << "PatternGenerator::endJob() rootFileName " << rootFileName << std::endl;
    TFile inFile(rootFileName.c_str());
    TDirectory* curDir = (TDirectory*)inFile.Get("layerStats");

    ostringstream ostrName;
    for (auto& gp : goldenPatterns) {
      if (gp->key().thePt == 0)
        continue;

      int statBinsCnt = 1024;  //= gp->getPdf()[0][0].size() * 8; //TODO should be big enough to comprise the pdf tails
      gp->iniStatisitics(statBinsCnt, 1);  //TODO

      for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
        for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
          ostrName.str("");
          ostrName << "histLayerStat_PatNum_" << gp->key().theNumber << "_refLayer_" << iRefLayer << "_Layer_"
                   << iLayer;

          TH1I* histLayerStat = (TH1I*)curDir->Get(ostrName.str().c_str());

          if (histLayerStat) {
            for (int iBin = 0; iBin < statBinsCnt; iBin++) {
              gp->updateStat(iLayer, iRefLayer, iBin, 0, histLayerStat->GetBinContent(iBin + 1));
            }
          } else {
            edm::LogImportant("l1tOmtfEventPrint")
                << "PatternGenerator::endJob() - reading histLayerStat: histogram not found " << ostrName.str()
                << std::endl;
          }
        }
      }
    }

    TH1* simMuFoundByOmtfPt_fromFile = (TH1*)inFile.Get("simMuFoundByOmtfPt");
    for (unsigned int iGp = 0; iGp < eventCntPerGp.size(); iGp++) {
      eventCntPerGp[iGp] = simMuFoundByOmtfPt_fromFile->GetBinContent(simMuFoundByOmtfPt_fromFile->FindBin(iGp));
      edm::LogImportant("l1tOmtfEventPrint")
          << "PatternGenerator::endJob() - eventCntPerGp: iGp" << iGp << " - " << eventCntPerGp[iGp] << std::endl;
    }

    //TODO chose the desired grouping ///////////////
    int group = 0;
    int indexInGroup = 0;
    for (auto& gp : goldenPatterns) {
      indexInGroup++;
      gp->key().setGroup(group);
      gp->key().setIndexInGroup(indexInGroup);
      //indexInGroup is counted from 1

      edm::LogImportant("l1tOmtfEventPrint")
          << "setGroup(group): group " << group << " indexInGroup " << indexInGroup << std::endl;

      if (gp->key().thePt <= 10 && indexInGroup == 2) {  //TODO
        indexInGroup = 0;
        group++;
      }

      if (gp->key().thePt > 10 && indexInGroup == 4) {  //TODO
        indexInGroup = 0;
        group++;
      }
    }  /////////////////////////////////////////////

    upadatePdfs();

    modifyClassProb(1);

    //groupPatterns(); IMPORTANT don't call grouping here, just set the groups above!!!!

    reCalibratePt();
    this->writeLayerStat = true;
  }

  PatternOptimizerBase::endJob();
}

void PatternGenerator::upadatePdfs() {
  //TODO setting the DistPhiBitShift i.e. grouping of the pdfBins
  for (auto& gp : goldenPatterns) {
    if (gp->key().thePt == 0)
      continue;
    for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
      for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
        if (gp->getDistPhiBitShift(iLayer, iRefLayer)) {
          throw runtime_error(
              string(__FUNCTION__) + ":" + to_string(__LINE__) +
              "gp->getDistPhiBitShift(iLayer, iRefLayer) != 0 -  cannot change DistPhiBitShift then!!!!");
        }

        if ((gp->key().thePt <= 10) && (iLayer == 1 || iLayer == 3 || iLayer == 5)) {
          gp->setDistPhiBitShift(1, iLayer, iRefLayer);
        } else
          gp->setDistPhiBitShift(0, iLayer, iRefLayer);

        //watch out: the shift in a given layer must be the same for patterns in one group
        //todo  make the setting of the shift on the group base
        //TODO set the DistPhiBitShift
        /*if( (gp->key().thePt <= 10) && (iLayer == 3 || iLayer == 5 ) && (iRefLayer == 0 || iRefLayer == 2 || iRefLayer == 6 || iRefLayer == 7)) {
          gp->setDistPhiBitShift(3, iLayer, iRefLayer);
        }
        else if( (gp->key().thePt <= 10) && ( iLayer == 1 || iLayer == 3 || iLayer == 5 ) ) {
          gp->setDistPhiBitShift(2, iLayer, iRefLayer);
        }
        else if( ( (gp->key().thePt <= 10) && (iLayer == 7 ||iLayer == 8 || iLayer == 17 ) ) ) {
          gp->setDistPhiBitShift(1, iLayer, iRefLayer);
        }
        else if( (gp->key().thePt <= 10) && (iLayer == 10 || iLayer == 11 || iLayer == 12 || iLayer == 13) && (iRefLayer == 1)) {
          gp->setDistPhiBitShift(1, iLayer, iRefLayer);
        }*/
      }
    }
  }

  double minHitCntThresh = 0.001;
  //Calculating meanDistPhi
  for (auto& gp : goldenPatterns) {
    if (gp->key().thePt == 0)
      continue;

    int minHitCnt = minHitCntThresh * eventCntPerGp[gp->key().number()];  // //TODO tune threshold <<<<<<<<<<<<<<<<<<
    edm::LogImportant("l1tOmtfEventPrint")
        << "PatternGenerator::upadatePdfs() Calculating meanDistPhi " << gp->key() << " eventCnt "
        << eventCntPerGp[gp->key().number()] << " minHitCnt " << minHitCnt << std::endl;

    for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
      for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
        //calculate meanDistPhi
        double meanDistPhi = 0;
        double count = 0;
        for (unsigned int iBin = 1; iBin < gp->getStatistics()[iLayer][iRefLayer].size(); iBin++) {
          //iBin = 0 is reserved for the no hit
          meanDistPhi += iBin * gp->getStatistics()[iLayer][iRefLayer][iBin][0];
          count += gp->getStatistics()[iLayer][iRefLayer][iBin][0];
        }

        if (count != 0) {
          meanDistPhi /= count;

          meanDistPhi -= (gp->getStatistics()[iLayer][iRefLayer].size() / 2);

          if (count < minHitCnt)
            meanDistPhi = 0;
          else
            edm::LogImportant("l1tOmtfEventPrint")
                << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << " iLayer " << iLayer << " iRefLayer "
                << iRefLayer << " count " << count << " meanDistPhi " << meanDistPhi << endl;
        }
        gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
      }
    }
  }

  OMTFConfiguration::vector2D patternGroups = omtfConfig->getPatternGroups(goldenPatterns);
  edm::LogImportant("l1tOmtfEventPrint") << "patternGroups:" << std::endl;
  for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
    edm::LogImportant("l1tOmtfEventPrint") << "patternGroup " << std::setw(2) << iGroup << " ";
    for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
      edm::LogImportant("l1tOmtfEventPrint") << i << " patNum " << patternGroups[iGroup][i] << " ";
    }
    edm::LogImportant("l1tOmtfEventPrint") << std::endl;
  }

  //averaging the meanDistPhi for the gp belonging to the same group
  for (unsigned int iLayer = 0; iLayer < goldenPatterns.at(0)->getPdf().size(); ++iLayer) {
    for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns.at(0)->getPdf()[iLayer].size(); ++iRefLayer) {
      //unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
      //if(refLayerLogicNum == iLayer)
      {
        for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
          double meanDistPhi = 0;
          int mergedCnt = 0;
          for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
            auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
            meanDistPhi += gp->meanDistPhiValue(iLayer, iRefLayer);
            if (gp->meanDistPhiValue(iLayer, iRefLayer) != 0)
              mergedCnt++;
          }

          if (mergedCnt) {
            meanDistPhi /= mergedCnt;
            //because for some gps the statistics can be too low, and then the meanDistPhiValue is 0, so it should not contribute
            for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
              auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
              gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
              edm::LogImportant("l1tOmtfEventPrint")
                  << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " " << gp->key()
                  << " iLayer " << iLayer << " iRefLayer " << iRefLayer << " meanDistPhi after averaging "
                  << meanDistPhi << endl;
            }
          }
        }
      }
    }
  }

  //calculating the pdfs
  for (auto& gp : goldenPatterns) {
    if (gp->key().thePt == 0)
      continue;

    //TODO tune threshold <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
    int minHitCnt = 2 * minHitCntThresh * eventCntPerGp[gp->key().number()];

    for (unsigned int iLayer = 0; iLayer < gp->getPdf().size(); ++iLayer) {
      for (unsigned int iRefLayer = 0; iRefLayer < gp->getPdf()[iLayer].size(); ++iRefLayer) {
        {
          double norm = 0;
          for (unsigned int iBin = 0; iBin < gp->getStatistics()[iLayer][iRefLayer].size();
               iBin++) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
            norm += gp->getStatistics()[iLayer][iRefLayer][iBin][0];
          }

          int pdfMiddle = gp->getPdf()[iLayer][iRefLayer].size() / 2;
          int statBinGroupSize = 1 << gp->getDistPhiBitShift(iLayer, iRefLayer);
          for (unsigned int iBinPdf = 0; iBinPdf < gp->getPdf()[iLayer][iRefLayer].size(); iBinPdf++) {
            double pdfVal = 0;
            if (iBinPdf > 0) {
              int groupedBins = 0;
              for (int i = 0; i < statBinGroupSize; i++) {
                int iBinStat =
                    statBinGroupSize * ((int)(iBinPdf)-pdfMiddle) + i + gp->meanDistPhiValue(iLayer, iRefLayer);

                iBinStat += (gp->getStatistics()[iLayer][iRefLayer].size() / 2);

                if (iBinStat >= 0 && iBinStat < (int)gp->getStatistics()[iLayer][iRefLayer].size()) {
                  pdfVal += gp->getStatistics()[iLayer][iRefLayer][iBinStat][0];
                  groupedBins++;
                  //cout<<__FUNCTION__<<": "<<__LINE__<<" "<<gp->key()<<" iLayer "<<iLayer<<" iBinStat "<<iBinStat<<" iBinPdf "<<iBinPdf<<" statVal "<<gp->getStatistics()[iLayer][iRefLayer][iBinStat][0]<<endl;
                }
              }
              if (norm > minHitCnt) {
                pdfVal /= (norm * statBinGroupSize);
              } else
                pdfVal = 0;
              /*edm::LogImportant("l1tOmtfEventPrint")
                                << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << "calculating pdf: iLayer " << iLayer
                                << " iRefLayer " << iRefLayer //<< " norm " << std::setw(5) << norm
                                << " pdfVal " << pdfVal << endl;*/
            } else {  //iBinPdf == 0 i.e. no hit
              int iBinStat = 0;
              if (norm > 0) {
                pdfVal = gp->getStatistics()[iLayer][iRefLayer][iBinStat][0] / norm;
              }
              edm::LogImportant("l1tOmtfEventPrint")
                  << __FUNCTION__ << ": " << __LINE__ << " " << gp->key() << "calculating pdf: iLayer " << iLayer
                  << " iRefLayer " << iRefLayer << " norm " << std::setw(5) << norm << " no hits cnt " << std::setw(5)
                  << gp->getStatistics()[iLayer][iRefLayer][iBinStat][0] << " pdfVal " << pdfVal << endl;
            }

            double minPdfValFactor = 1;
            const double minPlog = log(omtfConfig->minPdfVal() * minPdfValFactor);
            const double pdfMaxVal = omtfConfig->pdfMaxValue();

            int digitisedVal = 0;
            if (pdfVal >= omtfConfig->minPdfVal() * minPdfValFactor) {
              digitisedVal = rint(pdfMaxVal - log(pdfVal) / minPlog * pdfMaxVal);
            }

            gp->setPdfValue(digitisedVal, iLayer, iRefLayer, iBinPdf);
            //cout<<__FUNCTION__<<": "<<__LINE__<<" "<<gp->key()<<" iLayer "<<iLayer<<" iBinPdf "<<iBinPdf<<" pdfVal "<<pdfVal<<" digitisedVal "<<digitisedVal<<endl;
          }
        }
      }
    }
  }
}

void PatternGenerator::saveHists(TFile& outfile) {
  outfile.mkdir("ptDeltaPhiHists")->cd();
  //TODO uncomment if ptDeltaPhiHists are needed
  /*  for(unsigned int iCharge = 0; iCharge <= 1; iCharge++) {
    for(unsigned int iLayer = 0; iLayer < omtfConfig->nLayers(); ++iLayer) { //for the moment filing only ref layer, remove whe
      if(ptDeltaPhiHists[iCharge][iLayer]) {
        ptDeltaPhiHists[iCharge][iLayer]->Write();
      }
    }
  }*/
}

void PatternGenerator::modifyClassProb(double step) {
  edm::LogImportant("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " Correcting P(C_k) " << std::endl;
  unsigned int iPdf = omtfConfig->nPdfBins() / 2;  // <<(omtfConfig->nPdfAddrBits()-1);
  for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns[0]->getPdf()[0].size(); ++iRefLayer) {
    unsigned int refLayerLogicNumber = omtfConfig->getRefToLogicNumber()[iRefLayer];
    if (iRefLayer == 0 || iRefLayer == 2)  //DT
      step = 1.5;
    else if (iRefLayer == 5)  //DT
      step = 1.5;
    else if (iRefLayer == 1)  //CSC
      step = 1.5;
    else if (iRefLayer == 3)  //CSC
      step = 1.5;
    else if (iRefLayer == 5)  //RE2/3
      step = 1.5;
    else if (iRefLayer == 6 || iRefLayer == 7)  //bRPC
      step = 1.5;

    edm::LogImportant("l1tOmtfEventPrint")
        << __FUNCTION__ << ":" << __LINE__ << " RefLayer " << iRefLayer << " step " << step << std::endl;
    for (int sign = -1; sign <= 1; sign++) {
      for (auto& gp : boost::adaptors::reverse(goldenPatterns)) {
        if (gp->key().thePt == 0 || gp->key().theCharge != sign)
          continue;

        double ptFrom = omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom;
        double ptTo = omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo;

        double ptRange = ptTo - ptFrom;

        double minPdfValFactor = 0.1;
        double minPlog = log(omtfConfig->minPdfVal() * minPdfValFactor);
        double pdfMaxVal = omtfConfig->pdfMaxValue();

        pdfMaxVal /= 3.;
        minPlog *= 2;

        //last bin of the ptRange goes to 10000, so here we change it to 1000
        if (ptRange > 800)
          ptRange = 800;

        double norm = 0.001;
        double classProb = vxIntegMuRate(ptFrom, ptRange, 0.82, 1.24) * norm;

        int digitisedVal = rint(pdfMaxVal - log(classProb) / minPlog * pdfMaxVal);

        int newPdfVal = digitisedVal;  //gp->getPdf()[refLayerLogicNumber][iRefLayer][iPdf]

        if (ptFrom == 0)
          newPdfVal += 15;
        if (ptFrom == 3.5)
          newPdfVal += 15;
        if (ptFrom == 4)
          newPdfVal += 12;
        if (ptFrom == 4.5)
          newPdfVal += 9;
        if (ptFrom == 5)
          newPdfVal += 7;
        if (ptFrom == 6)
          newPdfVal += 4;
        if (ptFrom == 7)
          newPdfVal += 2;

        if (ptFrom == 100)
          newPdfVal = 16;
        if (ptFrom == 200)
          newPdfVal = 22;

        gp->setPdfValue(newPdfVal, refLayerLogicNumber, iRefLayer, iPdf);

        edm::LogImportant("l1tOmtfEventPrint")
            << gp->key() << " " << omtfConfig->getPatternPtRange(gp->key().theNumber).ptFrom << " - "
            << omtfConfig->getPatternPtRange(gp->key().theNumber).ptTo << " GeV"
            << " ptRange " << ptRange << " RefLayer " << iRefLayer << " newPdfVal " << newPdfVal << std::endl;
      }
    }
  }
}

void PatternGenerator::reCalibratePt() {
  edm::LogImportant("l1tOmtfEventPrint") << __FUNCTION__ << ": " << __LINE__ << " reCalibratePt" << std::endl;
  std::map<int, float> ptMap;
  //for Patterns_0x0009_oldSample_3_10Files_classProb2.xml
  ptMap[7] = 4.0;
  ptMap[8] = 4.5;
  ptMap[9] = 5.0;
  ptMap[10] = 5.5;
  ptMap[11] = 6.0;
  ptMap[13] = 7.0;
  ptMap[15] = 8.5;
  ptMap[17] = 10.0;
  ptMap[21] = 12.0;
  ptMap[25] = 14.0;
  ptMap[29] = 16.0;
  ptMap[33] = 18.5;
  ptMap[37] = 21.0;
  ptMap[41] = 23.0;
  ptMap[45] = 26.0;
  ptMap[49] = 28.0;
  ptMap[53] = 30.0;
  ptMap[57] = 32.0;
  ptMap[61] = 36.0;
  ptMap[71] = 40.0;
  ptMap[81] = 48.0;
  ptMap[91] = 54.0;
  ptMap[101] = 60.0;
  ptMap[121] = 70.0;
  ptMap[141] = 82.0;
  ptMap[161] = 96.0;
  ptMap[201] = 114.0;
  ptMap[401] = 200.0;

  for (auto& gp : goldenPatterns) {
    if (gp->key().thePt == 0)
      continue;

    int newPt = omtfConfig->ptGevToHw(ptMap[gp->key().thePt]);
    edm::LogImportant("l1tOmtfEventPrint") << gp->key().thePt << " -> " << newPt << std::endl;

    gp->key().setPt(newPt);
  }
}

void PatternGenerator::groupPatterns() {
  int group = 0;
  int indexInGroup = 0;
  for (auto& gp : goldenPatterns) {
    indexInGroup++;
    gp->key().setGroup(group);
    gp->key().setIndexInGroup(indexInGroup);
    //indexInGroup is counted from 1

    edm::LogImportant("l1tOmtfEventPrint")
        << "setGroup(group): group " << group << " indexInGroup " << indexInGroup << std::endl;

    if (gp->key().thePt <= 12 && indexInGroup == 2) {  //TODO
      indexInGroup = 0;
      group++;
    }

    if (gp->key().thePt > 12 && indexInGroup == 4) {  //TODO
      indexInGroup = 0;
      group++;
    }
  }

  OMTFConfiguration::vector2D patternGroups = omtfConfig->getPatternGroups(goldenPatterns);
  edm::LogImportant("l1tOmtfEventPrint") << "patternGroups:" << std::endl;
  for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
    edm::LogImportant("l1tOmtfEventPrint") << "patternGroup " << std::setw(2) << iGroup << " ";
    for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
      edm::LogImportant("l1tOmtfEventPrint") << i << " patNum " << patternGroups[iGroup][i] << " ";
    }
    edm::LogImportant("l1tOmtfEventPrint") << std::endl;
  }

  int pdfBins = exp2(omtfConfig->nPdfAddrBits());

  for (unsigned int iLayer = 0; iLayer < goldenPatterns.at(0)->getPdf().size(); ++iLayer) {
    for (unsigned int iRefLayer = 0; iRefLayer < goldenPatterns.at(0)->getPdf()[iLayer].size(); ++iRefLayer) {
      //unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
      //if(refLayerLogicNum == iLayer)
      {
        //averaging the meanDistPhi for the gp belonging to the same group
        for (unsigned int iGroup = 0; iGroup < patternGroups.size(); iGroup++) {
          double meanDistPhi = 0;
          int mergedCnt = 0;
          for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
            auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
            meanDistPhi += gp->meanDistPhiValue(iLayer, iRefLayer);
            if (gp->meanDistPhiValue(iLayer, iRefLayer) != 0)
              mergedCnt++;
            edm::LogImportant("l1tOmtfEventPrint")
                << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " " << gp->key()
                << " iLayer " << iLayer << " iRefLayer " << iRefLayer << " old meanDistPhiValue "
                << gp->meanDistPhiValue(iLayer, iRefLayer) << endl;
          }

          if (mergedCnt) {
            meanDistPhi /= mergedCnt;
            meanDistPhi = (int)meanDistPhi;

            //because for some gps the statistics can be too low, and then the meanDistPhiValue is 0, so it should not contribute
            for (unsigned int i = 0; i < patternGroups[iGroup].size(); i++) {
              auto gp = goldenPatterns.at(patternGroups[iGroup][i]).get();
              unsigned int refLayerLogicNum = omtfConfig->getRefToLogicNumber()[iRefLayer];
              if (refLayerLogicNum != iLayer) {
                int shift = meanDistPhi - gp->meanDistPhiValue(iLayer, iRefLayer);
                edm::LogImportant("l1tOmtfEventPrint")
                    << __FUNCTION__ << ": " << __LINE__ << " iGroup " << iGroup << " numInGroup " << i << " "
                    << gp->key() << " iLayer " << iLayer << " iRefLayer " << iRefLayer
                    << " new meanDistPhi after averaging " << meanDistPhi << " old meanDistPhiValue "
                    << gp->meanDistPhiValue(iLayer, iRefLayer) << " shift " << shift << endl;

                if (shift < 0) {
                  for (int iBin = 1 - shift; iBin < pdfBins;
                       iBin++) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
                    auto pdfVal = gp->pdfValue(iLayer, iRefLayer, iBin);
                    gp->setPdfValue(pdfVal, iLayer, iRefLayer, iBin + shift);
                    edm::LogImportant("l1tOmtfEventPrint")
                        << "     iBin " << iBin << "  iBin + shift " << iBin + shift << " pdfVal " << pdfVal << endl;
                  }
                  for (int iBin = pdfBins + shift; iBin < pdfBins; iBin++) {
                    gp->setPdfValue(0, iLayer, iRefLayer, iBin);
                  }
                } else if (shift > 0) {
                  for (int iBin = pdfBins - 1 - shift; iBin > 0;
                       iBin--) {  //iBin = 0 i.e. no hit is included here, to have the proper norm
                    auto pdfVal = gp->pdfValue(iLayer, iRefLayer, iBin);
                    gp->setPdfValue(pdfVal, iLayer, iRefLayer, iBin + shift);
                    edm::LogImportant("l1tOmtfEventPrint")
                        << "     iBin " << iBin << "  iBin + shift " << iBin + shift << " pdfVal " << pdfVal << endl;
                  }
                  for (int iBin = shift; iBin > 0; iBin--) {
                    gp->setPdfValue(0, iLayer, iRefLayer, iBin);
                  }
                }
              }

              gp->setMeanDistPhiValue(round(meanDistPhi), iLayer, iRefLayer);
            }
          }
        }
      }
    }
  }
}

ext: cc

max_stars_repo_name: malbouis/cmssw

max_issues_repo_name: gartung/cmssw

max_forks_repo_name: gartung/cmssw


text:
/* Siconos is a program dedicated to modeling, simulation and control
 * of non smooth dynamical systems.
 *
 * Copyright 2020 INRIA.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
*/

#include "BlockMatrixIterators.hpp"
#include <iostream>                                 // for operator<<, ostream
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include "SiconosAlgebra.hpp"
#include "SimpleMatrixFriends.hpp"
#include "BlockMatrix.hpp"
#include "SimpleMatrix.hpp"
#include "SiconosVector.hpp"
#include "SiconosMatrixSetBlock.hpp" // for setBlock
#include "SiconosAlgebraTools.hpp" // for isComparableTo
#include "Tools.hpp"
#include "SiconosException.hpp"

using  std::cout;
using std::endl;
using Siconos::Algebra::isComparableTo;

// =================================================
//                CONSTRUCTORS
// =================================================

BlockMatrix::BlockMatrix(const SiconosMatrix &m): SiconosMatrix(Siconos::BLOCK), _dimRow(0), _dimCol(0)
{
  _tabRow.reset(new Index());
  _tabCol.reset(new Index());
  if(m.isBlock())
  {
    const BlockMatrix& mB = static_cast<const BlockMatrix&>(m);
    unsigned int nbRows = m.numberOfBlocks(0);
    unsigned int nbCols = m.numberOfBlocks(1);
    _tabRow->reserve(nbRows);
    _tabCol->reserve(nbCols);

    // mat construction
    _mat.reset(new BlocksMat(nbRows, nbCols, nbRows * nbCols));

    unsigned int i, j;
    ConstBlocksIterator1 it1;
    ConstBlocksIterator2 it2;
    bool firstLoop = true;
    // We scan all the blocks of m ...
    for(it1 = mB._mat->begin1(); it1 != mB._mat->end1(); ++it1)
    {
      _dimRow += (*(it1.begin()))->size(0);
      _tabRow->push_back(_dimRow);
      for(it2 = it1.begin(); it2 != it1.end(); ++it2)
      {
        i = it2.index1();
        j = it2.index2();
        if((*it2)->isBlock())   // if the current matrix is a blockMatrix
          _mat->insert_element(i, j, std::shared_ptr<SiconosMatrix>(new BlockMatrix(**it2)));
        else
          _mat->insert_element(i, j, std::shared_ptr<SiconosMatrix>(new SimpleMatrix(**it2)));
        // _dimCol must be incremented only at first "column-loop"
        if(firstLoop)
        {
          _dimCol += (*it2)->size(1);
          _tabCol->push_back(_dimCol);
        }
      }
      firstLoop = false;
    }
  }
  else // if m is a SimpleMatrix
  {
    _tabRow->reserve(1);
    _tabCol->reserve(1);
    // _mat construction
    _mat.reset(new BlocksMat(1, 1, 1));
    _mat->insert_element(0, 0, std::shared_ptr<SiconosMatrix>(new SimpleMatrix(m)));

    _dimRow = m.size(0);
    _dimCol = m.size(1);
    _tabRow->push_back(_dimRow);
    _tabCol->push_back(_dimCol);
  }
}

BlockMatrix::BlockMatrix(const BlockMatrix &m): SiconosMatrix(Siconos::BLOCK), _dimRow(0), _dimCol(0)
{
  unsigned int nbRows = m.numberOfBlocks(0);
  unsigned int nbCols = m.numberOfBlocks(1);
  _tabRow.reset(new Index());
  _tabCol.reset(new Index());
  _tabRow->reserve(nbRows);
  _tabCol->reserve(nbCols);

  // _mat construction
  _mat.reset(new BlocksMat(nbRows, nbCols, nbRows * nbCols));

  unsigned int i, j;
  // We scan all the blocks of m ...
  ConstBlocksIterator1 it1;
  ConstBlocksIterator2 it2;
  bool firstLoop = true;
  // We scan all the blocks of m ...
  for(it1 = m._mat->begin1(); it1 != m._mat->end1(); ++it1)
  {
    _dimRow += (*(it1.begin()))->size(0);
    _tabRow->push_back(_dimRow);
    for(it2 = it1.begin(); it2 != it1.end(); ++it2)
    {
      i = it2.index1();
      j = it2.index2();
      if((*it2)->isBlock())   // if the current _matrix is a blockMatrix
        _mat->insert_element(i, j, std::shared_ptr<SiconosMatrix>(new BlockMatrix(**it2)));
      else
        _mat->insert_element(i, j, std::shared_ptr<SiconosMatrix>(new SimpleMatrix(**it2)));

      // _dimCol must be incremented only at first "column-loop"
      if(firstLoop)
      {
        _dimCol += (*it2)->size(1);
        _tabCol->push_back(_dimCol);
      }
    }
    firstLoop = false;
  }
}

BlockMatrix::BlockMatrix(const std::vector<SP::SiconosMatrix >& m, unsigned int row, unsigned int col):
  SiconosMatrix(Siconos::BLOCK), _dimRow(0), _dimCol(0)
{
  if(m.size() != (row * col))
    THROW_EXCEPTION("number of blocks inconsistent with provided dimensions.");

  _tabRow.reset(new Index());
  _tabCol.reset(new Index());
  _tabRow->reserve(row);
  _tabCol->reserve(col);

  // _mat construction
  _mat.reset(new BlocksMat(row, col, row * col));

  unsigned int k = 0;
  bool firstRowLoop = true;
  bool firstColLoop = true;

  for(unsigned int i = 0; i < row; ++i)
  {
    for(unsigned int j = 0; j < col; ++j)
    {
      (*_mat)(i, j) = m[k++];

      // _dimCol must be incremented only at first "column-loop"
      if(firstColLoop)
      {
        _dimCol += m[k - 1]->size(1);
        _tabCol->push_back(_dimCol);
      }
      if(firstRowLoop)
      {
        _dimRow += m[k - 1]->size(0);
        _tabRow->push_back(_dimRow);
        firstRowLoop = false;
      }
    }
    firstColLoop = false;
    firstRowLoop = true;
  }
}

BlockMatrix::BlockMatrix(SP::SiconosMatrix A, SP::SiconosMatrix B, SP::SiconosMatrix C, SP::SiconosMatrix D):
  SiconosMatrix(Siconos::BLOCK), _dimRow(0), _dimCol(0)
{
  if(A->size(0) != B->size(0) || C->size(0) != D->size(0) ||  A->size(1) != C->size(1) ||  B->size(1) != D->size(1))
    THROW_EXCEPTION("inconsistent sizes between A, B, C or D SiconosMatrices.");

  // _mat = [ A B ]
  //       [ C D ]

  // _mat construction
  _mat.reset(new BlocksMat(2, 2, 4));

  _tabRow.reset(new Index());
  _tabCol.reset(new Index());
  _tabRow->reserve(2);
  _tabCol->reserve(2);

  (*_mat)(0, 0) = A;
  (*_mat)(0, 1) = B;
  (*_mat)(1, 0) = C;
  (*_mat)(1, 1) = D;
  _dimRow = A->size(0);
  _tabRow->push_back(_dimRow);
  _dimRow += C->size(0);
  _tabRow->push_back(_dimRow);
  _dimCol = A->size(1);
  _tabCol->push_back(_dimCol);
  _dimCol += B->size(1);
  _tabCol->push_back(_dimCol);

}

BlockMatrix::~BlockMatrix()
{

  _mat->clear();

  _tabRow->clear();
  _tabCol->clear();
}

// =================================================
//    get number of blocks
// =================================================

unsigned int BlockMatrix::numberOfBlocks(unsigned int dim) const
{
  if(dim == 0)
    return _tabRow->size();
  else
    return _tabCol->size();
}

// =================================================
//        get Ublas component (dense ...)
// =================================================

// return the boost dense _matrix of the block (i, j)
const DenseMat  BlockMatrix::getDense(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp;
  tmp = (*_mat)(row, col);

  if(tmp->num() != Siconos::DENSE)
    THROW_EXCEPTION("the matrix at (row, col) is not a Dense matrix");

  return (tmp->getDense());
}

// return the boost triangular matrix of the block (i, j)
const TriangMat BlockMatrix::getTriang(unsigned int row, unsigned int col) const
{
  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::TRIANGULAR)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Triangular matrix");
  }
  return (tmp->getTriang());
}

// return the boost symmetric matrix of the block (i, j)
const SymMat BlockMatrix::getSym(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SYMMETRIC)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Symmmetric matrix");
  }
  return (tmp->getSym());
}

// return the boost sparse matrix of the block (i, j)
const SparseMat  BlockMatrix::getSparse(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SPARSE)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Sparse matrix");
  }
  return (tmp->getSparse());
}

// return the boost sparse matrix of the block (i, j)
const SparseCoordinateMat  BlockMatrix::getSparseCoordinate(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SPARSE_COORDINATE)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Sparse matrix");
  }
  return (tmp->getSparseCoordinate());
}
// return the boost banded matrix of the block (i, j)
const BandedMat  BlockMatrix::getBanded(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::BANDED)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Banded matrix");
  }
  return (tmp->getBanded());
}

// return the boost zero matrix of the block (i, j)
const ZeroMat  BlockMatrix::getZero(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::ZERO)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Zero matrix");
  }
  return (tmp->getZero());
}

// return the boost identity matrix of the block (i, j)
const IdentityMat  BlockMatrix::getIdentity(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::IDENTITY)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Identity matrix");
  }
  return (tmp->getIdentity());
}

// The following functions return the corresponding pointers
DenseMat*  BlockMatrix::dense(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::DENSE)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Dense matrix");
  }

  return (tmp->dense());
}

TriangMat* BlockMatrix::triang(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::TRIANGULAR)
  {
    THROW_EXCEPTION("The matrix at (row, col) is not a Triangular matrix");
  }
  return (tmp->triang());
}
SymMat* BlockMatrix::sym(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SYMMETRIC)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Symmmetric matrix");
  }
  return (tmp->sym());
}

SparseMat*  BlockMatrix::sparse(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SPARSE)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Sparse matrix");
  }
  return (tmp->sparse());
}
SparseCoordinateMat*  BlockMatrix::sparseCoordinate(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::SPARSE_COORDINATE)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Sparse coordinate matrix");
  }
  return (tmp->sparseCoordinate());
}

BandedMat*  BlockMatrix::banded(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::BANDED)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Banded matrix");
  }
  return (tmp->banded());
}

ZeroMat*  BlockMatrix::zero_mat(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::ZERO)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Zero matrix");
  }
  return (tmp->zero_mat(row, col));
}

IdentityMat*  BlockMatrix::identity(unsigned int row, unsigned int col) const
{

  SP::SiconosMatrix tmp = (*_mat)(row, col);
  if(tmp->num() != Siconos::IDENTITY)
  {
    THROW_EXCEPTION("the matrix at (row, col) is not a Identity matrix");
  }
  return (tmp->identity());
}

double* BlockMatrix::getArray(unsigned int i, unsigned int j) const
{
  SP::SiconosMatrix tmp = (*_mat)(i, j);
  return tmp->getArray();
}

// ===========================
//       fill matrix
// ===========================

void BlockMatrix::zero()
{
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      (*it2)->zero();
    }
  }
}

void BlockMatrix::randomize()
{
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      (*it2)->randomize();
    }
  }
}

void BlockMatrix::randomize_sym()
{
  THROW_EXCEPTION("not yet implemented for block matrices.");
}

void BlockMatrix::eye()
{
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;

  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      if(it2.index1() == it2.index2())
        (*it2)->eye();
      else
        (*it2)->zero();
    }
  }
}

unsigned int BlockMatrix::size(unsigned int index) const
{
  if(index == 0) return _dimRow;
  else return _dimCol;
};

//=======================
// set matrix dimensions
//=======================

void BlockMatrix::resize(unsigned int, unsigned int, unsigned int, unsigned int, bool)
{
  THROW_EXCEPTION("forbidden for block matrices.");
}

//=======================
//       get norm
//=======================

double BlockMatrix::normInf()const
{
  double sum = 0, norm = 0;
  for(unsigned int i = 0; i < size(0); i++)
  {
    for(unsigned int j = 0; j < size(1); j++)
    {
      sum += (*this)(i, j);
    }
    if(fabs(sum) > norm) norm = fabs(sum);
    sum = 0;
  }
  return norm;
}

//=====================
// screen display
//=====================

void BlockMatrix::display(void)const
{
  std::cout << "==========> BlockMatrix (" << numberOfBlocks(0) << " X " << numberOfBlocks(1) << " blocks): " << std::endl;
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      (*it2)->display();
    }
  }
  std::cout << "===========================================================================================" << std::endl;
}
void BlockMatrix::displayExpert(bool brief)const
{
  std::cout << "==========> BlockMatrix (" << numberOfBlocks(0) << " X " << numberOfBlocks(1) << " blocks): " << std::endl;
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      (*it2)->displayExpert(brief);
    }
  }
  std::cout << "===========================================================================================" << std::endl;
}

//=====================
// convert to a string
//=====================

std::string BlockMatrix::toString() const
{
  return ::toString(*this);
}

//=====================
// convert to an ostream
//=====================

std::ostream& operator<<(std::ostream& os, const BlockMatrix& bm)
{
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  os << "[" << bm.numberOfBlocks(0) << "," << bm.numberOfBlocks(1) << "](";
  for(it = bm._mat->begin1(); it != bm._mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
    {
      if(it2 != it.begin()) os << ",";
      if(*it2) os << **it2;
      else os << "(nil)";
    }
  }
  os << ")";
  return os;
}

//=============================
// Elements access (get or set)
//=============================

double& BlockMatrix::operator()(unsigned int row, unsigned int col)
{
  unsigned int nbRow = 0;
  unsigned int nbCol = 0;

  while(row >= (*_tabRow)[nbRow] && nbRow < _tabRow->size())
    nbRow ++;

  while(col >= (*_tabCol)[nbCol] && nbCol < _tabCol->size())
    nbCol ++;

  unsigned int posRow = row;
  unsigned int posCol = col;

  if(nbRow != 0)
    posRow -= (*_tabRow)[nbRow - 1];
  if(nbCol != 0)
    posCol -= (*_tabCol)[nbCol - 1];

  SP::SiconosMatrix tmp = (*_mat)(nbRow, nbCol);
  return (*tmp)(posRow, posCol);
}

double BlockMatrix::operator()(unsigned int row, unsigned int col) const
{

  unsigned int nbRow = 0;
  unsigned int nbCol = 0;

  while(row >= (*_tabRow)[nbRow] && nbRow < _tabRow->size())
    nbRow ++;

  while(col >= (*_tabCol)[nbCol] && nbCol < _tabCol->size())
    nbCol ++;

  unsigned int posRow = row;
  unsigned int posCol = col;

  if(nbRow != 0)
    posRow -= (*_tabRow)[nbRow - 1];
  if(nbCol != 0)
    posCol -= (*_tabCol)[nbCol - 1];

  SP::SiconosMatrix tmp = (*_mat)(nbRow, nbCol);
  return (*tmp)(posRow, posCol);
}

double BlockMatrix::getValue(unsigned int row, unsigned int col) const
{
  unsigned int nbRow = 0;
  unsigned int nbCol = 0;

  while(row >= (*_tabRow)[nbRow] && nbRow < _tabRow->size())
    nbRow ++;

  while(col >= (*_tabCol)[nbCol] && nbCol < _tabCol->size())
    nbCol ++;

  unsigned int posRow = row;
  unsigned int posCol = col;

  if(nbRow != 0)
    posRow -= (*_tabRow)[nbRow - 1];
  if(nbCol != 0)
    posCol -= (*_tabCol)[nbCol - 1];

  SP::SiconosMatrix tmp = (*_mat)(nbRow, nbCol);
  return (*tmp)(posRow, posCol);
}

void BlockMatrix::setValue(unsigned int row, unsigned int col, double value)
{
  unsigned int nbRow = 0;
  unsigned int nbCol = 0;

  while(row >= (*_tabRow)[nbRow] && nbRow < _tabRow->size())
    nbRow ++;

  while(col >= (*_tabCol)[nbCol] && nbCol < _tabCol->size())
    nbCol ++;

  unsigned int posRow = row;
  unsigned int posCol = col;

  if(nbRow != 0)
    posRow -= (*_tabRow)[nbRow - 1];
  if(nbCol != 0)
    posCol -= (*_tabCol)[nbCol - 1];

  SP::SiconosMatrix tmp = (*_mat)(nbRow, nbCol);
  (*tmp)(posRow, posCol) = value;
}

//============================================
// Access (get or set) to blocks of elements
//============================================

void BlockMatrix::getRow(unsigned int r, SiconosVector &v) const
{
  unsigned int numRow = 0, posRow = r, start = 0, stop = 0;

  if(r > _dimRow)
    THROW_EXCEPTION("row number is out of range");

  // Verification of the size of the result vector
  if(v.size() != _dimCol)
    THROW_EXCEPTION("inconsistent sizes");

  // Find the row-block number where "r" is
  while(r >= (*_tabRow)[numRow] && numRow < _tabRow->size())
    numRow ++;

  // Computation of the value of the index row into this block
  if(numRow != 0)
    posRow -= (*_tabRow)[numRow - 1];

  for(unsigned int j = 0; j < _tabCol->size(); j++)
  {
    start = stop;
    SP::SiconosMatrix tmp = (*_mat)(numRow, j);
    stop += tmp->size(1);
    ublas::subrange(*(v.dense()), start, stop) = ublas::row(*(tmp->dense()), posRow);
  }
}

void BlockMatrix::getCol(unsigned int c, SiconosVector &v) const
{
  unsigned int numCol = 0, posCol = c, start = 0, stop = 0;

  if(c > _dimCol)
    THROW_EXCEPTION("column number is out of range");

  // Verification of the size of the result vector
  if(v.size() != _dimRow)
    THROW_EXCEPTION("inconsistent sizes");

  // Find the column-block number where "c" is
  while(c >= (*_tabCol)[numCol] && numCol < _tabCol->size())
    numCol ++;

  // Computation of the value of the index column into this block
  if(numCol != 0)
    posCol -= (*_tabCol)[numCol - 1];

  for(unsigned int i = 0; i < _tabRow->size(); i++)
  {
    start = stop;
    SP::SiconosMatrix tmp = (*_mat)(i, numCol);
    stop += tmp->size(0);
    ublas::subrange(*(v.dense()), start, stop) = ublas::column(tmp->getDense(), posCol);
  }
}

void BlockMatrix::setRow(unsigned int r, const SiconosVector &v)
{

  unsigned int numRow = 0, posRow = r, start = 0, stop = 0;

  if(v.size() != _dimCol)
    THROW_EXCEPTION("inconsistent sizes");

  while(r >= (*_tabRow)[numRow] && numRow < _tabRow->size())
    numRow ++;

  if(numRow != 0)
    posRow -= (*_tabRow)[numRow - 1];

  for(unsigned int j = 0; j < _tabCol->size(); j++)
  {
    start = stop;
    SP::SiconosMatrix tmp = (*_mat)(numRow, j);
    stop += tmp->size(1);
    ublas::row(*(tmp->dense()), posRow) = ublas::subrange(*(v.dense()), start, stop);
  }
}

void BlockMatrix::setCol(unsigned int col, const SiconosVector &v)
{

  unsigned int numCol = 0, posCol = col, start = 0, stop = 0;

  if(v.size() != _dimRow)
    THROW_EXCEPTION("inconsistent sizes");

  while(col >= (*_tabCol)[numCol] && numCol < _tabCol->size())
    numCol ++;

  if(numCol != 0)
    posCol -= (*_tabCol)[numCol - 1];

  for(unsigned int i = 0; i < _tabRow->size(); i++)
  {
    start = stop;
    SP::SiconosMatrix tmp = (*_mat)(i, numCol);
    stop += tmp->size(0);
    ublas::column(*(tmp->dense()), posCol) = ublas::subrange(*(v.dense()), start, stop);
  }
}

void BlockMatrix::addSimple(unsigned int& indRow, unsigned int& indCol, const SiconosMatrix& m)
{
  // Add a part of m (starting from (indRow,indCol) to the current matrix.
  // m must be a SimpleMatrix.

  // At the end of the present function, indRow (resp. indCol) is equal to indRow + the corresponding dimension of the added sub-matrix.

  unsigned int row = m.size(0) - indRow; // number of rows of the block to be added.
  unsigned int col = m.size(1) - indCol; // number of columns of the block to be added.
  unsigned int initCol = indCol;

  if(row > _dimRow || col > _dimCol) THROW_EXCEPTION("invalid ranges");

  Siconos::UBLAS_TYPE numM = m.num();

  // iterators through this
  BlocksMat::iterator1 it1;
  BlocksMat::iterator2 it2;
  unsigned int currentRow = 0, currentCol = 0, currentNum;
  for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
  {
    for(it2 = it1.begin(); it2 != it1.end(); ++it2)
    {
      if((*it2)->isBlock())   // if the sub-block is also a BlockMatrix ...
        (std::static_pointer_cast<BlockMatrix>(*it2))->addSimple(indRow, indCol, m);

      else
      {
        currentCol = (*it2)->size(1);
        currentRow = (*it2)->size(0);
        currentNum = (*it2)->num();
        if(numM != currentNum) THROW_EXCEPTION("inconsistent types.");

        if(numM == Siconos::DENSE)
          noalias(*(*it2)->dense()) += ublas::subrange(*m.dense(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::TRIANGULAR)
          noalias(*(*it2)->triang()) += ublas::subrange(*m.triang(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::SYMMETRIC)
          noalias(*(*it2)->sym()) += ublas::subrange(*m.sym(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::SPARSE)
          noalias(*(*it2)->sparse()) += ublas::subrange(*m.sparse(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::BANDED)
          noalias(*(*it2)->banded()) += ublas::subrange(*m.banded(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::ZERO) {}
        else
          THROW_EXCEPTION("inconsistent types.");
      }
      indCol += currentCol;
    }
    indRow += currentRow;
    indCol = initCol;
  }
}

void BlockMatrix::subSimple(unsigned int& indRow, unsigned int& indCol, const SiconosMatrix& m)
{
  // subtract a part of m (starting from (indRow,indCol) to the current matrix.
  // m must be a SimpleMatrix.

  // At the end of the present function, indRow (resp. indCol) is equal to indRow + the corresponding dimension of the subtracted sub-matrix.

  unsigned int row = m.size(0) - indRow; // number of rows of the block to be added.
  unsigned int col = m.size(1) - indCol; // number of columns of the block to be added.
  unsigned int initCol = indCol;
  if(row > _dimRow || col > _dimCol) THROW_EXCEPTION("invalid ranges");

  Siconos::UBLAS_TYPE numM = m.num();

  // iterators through this
  BlocksMat::iterator1 it1;
  BlocksMat::iterator2 it2;
  unsigned int currentRow = 0, currentCol = 0;
  Siconos::UBLAS_TYPE currentNum;
  for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
  {
    for(it2 = it1.begin(); it2 != it1.end(); ++it2)
    {
      if((*it2)->isBlock())   // if the sub-block is also a BlockMatrix ...
        (std::static_pointer_cast<BlockMatrix>(*it2))->subSimple(indRow, indCol, m);

      else
      {
        currentCol = (*it2)->size(1);
        currentRow = (*it2)->size(0);
        currentNum = (*it2)->num();
        if(numM != currentNum) THROW_EXCEPTION("inconsistent types.");

        if(numM == Siconos::DENSE)
          noalias(*(*it2)->dense()) -= ublas::subrange(*m.dense(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::TRIANGULAR)
          noalias(*(*it2)->triang()) -= ublas::subrange(*m.triang(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::SYMMETRIC)
          noalias(*(*it2)->sym()) -= ublas::subrange(*m.sym(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::SPARSE)
          noalias(*(*it2)->sparse()) -= ublas::subrange(*m.sparse(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::BANDED)
          noalias(*(*it2)->banded()) -= ublas::subrange(*m.banded(), indRow, indRow + currentRow, indCol, indCol + currentCol);
        else if(numM == Siconos::ZERO) {}
        else
          THROW_EXCEPTION("inconsistent types.");
      }
      indCol += currentCol;
    }
    indRow += currentRow;
    indCol = initCol;
  }
}

//===============
//  Assignment
//===============

BlockMatrix& BlockMatrix::operator = (const SiconosMatrix &m)
{
  if(&m == this) return *this;  // auto-assignment.

  if(m.size(0) != _dimRow || m.size(1) != _dimCol)
    THROW_EXCEPTION("Left and Right values have inconsistent sizes.");

  // Warning: we do not reallocate the blocks, but only copy the values. This means that
  // all blocks are already allocated and that dim of m and mat are to be consistent.
  // Thus, _tabRow and _tabCol remains unchanged.
  // If m and mat are not "block-consistent", we use the () operator for a component-wise copy.

  if(m.isBlock())
  {
    if(isComparableTo(*this, m))
    {
      const BlockMatrix& mB = static_cast<const BlockMatrix&>(m);
      // iterators through this
      BlocksIterator1 it1;
      BlocksIterator2 it2;
      // iterators through m
      ConstBlocksIterator1 itM1 = mB._mat->begin1();
      ConstBlocksIterator2 itM2;

      for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
      {
        itM2 = itM1.begin();
        for(it2 = it1.begin(); it2 != it1.end(); ++it2)
        {
          (**it2) = (**itM2);
          itM2++; // increment column pos. in m.
        }
        itM1++; // increment row pos. in m.
      }
    }
    else
    {
      for(unsigned int i = 0; i < _dimRow; ++i)
        for(unsigned int j = 0; j < _dimCol; ++j)
          (*this)(i, j) = m(i, j);
    }
  }
  else // if m is a SimpleMatrix
  {
    BlocksIterator1 it;
    BlocksIterator2 it2;
    unsigned int posRow = 0;
    unsigned int posCol = 0;
    Index subDim(2);
    Index subPos(4);

    for(it = _mat->begin1(); it != _mat->end1(); ++it)
    {
      for(it2 = it.begin(); it2 != it.end(); ++it2)
      {
        // a sub-block of m is copied into this
        subDim[0] = (*it2)->size(0);
        subDim[1] = (*it2)->size(1);
        subPos[0] = posRow;
        subPos[1] = posCol;
        subPos[2] = 0;
        subPos[3] = 0;
        setBlock(createSPtrConstSiconosMatrix(m), *it2, subDim, subPos);
        posCol += subDim[1];
      }
      posRow += (*it)->size(0);
      posCol = 0;
    }
  }

  return *this;
}

BlockMatrix& BlockMatrix::operator = (const BlockMatrix &m)
{
  if(&m == this) return *this;  // auto-assignment.

  if(m.size(0) != _dimRow || m.size(1) != _dimCol)
    THROW_EXCEPTION("Left and Right values have inconsistent sizes.");

  // Warning: we do not reallocate the blocks, but only copy the values. This means that
  // all blocks are already allocated and that dim of m and mat are to be consistent.
  // Thus, _tabRow and _tabCol remains unchanged.
  // If m and mat are not "block-consistent", we use the () operator for a componet-wise copy.

  if(isComparableTo(*this, m))
  {
    const BlockMatrix& mB = static_cast<const BlockMatrix&>(m);
    // iterators through this
    BlocksIterator1 it1;
    BlocksIterator2 it2;
    // iterators through m
    ConstBlocksIterator1 itM1 = mB._mat->begin1();
    ConstBlocksIterator2 itM2;

    for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
    {
      itM2 = itM1.begin();
      for(it2 = it1.begin(); it2 != it1.end(); ++it2)
      {
        (**it2) = (**itM2);
        itM2++; // increment column pos. in m.
      }
      itM1++; // increment row pos. in m.
    }
  }
  else
  {
    for(unsigned int i = 0; i < _dimRow; ++i)
      for(unsigned int j = 0; j < _dimCol; ++j)
        (*this)(i, j) = m(i, j);
  }
  return *this;
}

BlockMatrix& BlockMatrix::operator = (const DenseMat &m)
{
  THROW_EXCEPTION("Not yet implemented.");
  return *this;
}

//=================================
// Op. and assignment (+=, -= ... )
//=================================

BlockMatrix& BlockMatrix::operator += (const SiconosMatrix &m)
{
  if(&m == this)
  {
    BlocksMat::iterator1 it1;
    BlocksMat::iterator2 it2;
    for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
    {
      for(it2 = it1.begin(); it2 != it1.end(); ++it2)
      {
        **it2 += **it2;
      }
    }
    return *this;
  }

  if(m.size(0) != _dimRow || m.size(1) != _dimCol)
    THROW_EXCEPTION("Left and Right values have inconsistent sizes.");

  if(m.isBlock())
  {
    if(isComparableTo(m, *this))
    {
      const BlockMatrix& mB = static_cast<const BlockMatrix&>(m);
      // iterators through this
      BlocksMat::iterator1 it1;
      BlocksMat::iterator2 it2;
      // iterators through m
      BlocksMat::const_iterator1 itM1 = mB._mat->begin1();
      BlocksMat::const_iterator2 itM2;

      for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
      {
        itM2 = itM1.begin();
        for(it2 = it1.begin(); it2 != it1.end(); ++it2)
        {
          (**it2) += (**itM2);
          itM2++; // increment column pos. in m.
        }
        itM1++; // increment row pos. in m.
      }
    }
    else
    {
      for(unsigned int i = 0; i < _dimRow; ++i)
        for(unsigned int j = 0; j < _dimCol; ++j)
          (*this)(i, j) += m(i, j);
    }
  }
  else // if m is a SimpleMatrix
  {
    unsigned int indRow = 0, indCol = 0;
    addSimple(indRow, indCol, m); // a sub-block of m is added to each block of this.
  }
  return *this;
}

BlockMatrix& BlockMatrix::operator -= (const SiconosMatrix &m)
{
  if(&m == this)
  {
    BlocksMat::iterator1 it1;
    BlocksMat::iterator2 it2;
    for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
    {
      for(it2 = it1.begin(); it2 != it1.end(); ++it2)
      {
        **it2 -= **it2;
      }
    }
    return *this;
  }

  if(m.size(0) != _dimRow || m.size(1) != _dimCol)
    THROW_EXCEPTION("Left and Right values have inconsistent sizes.");

  if(m.isBlock())
  {
    if(isComparableTo(m, *this))
    {
      const BlockMatrix& mB = static_cast<const BlockMatrix&>(m);
      // iterators through this
      BlocksMat::iterator1 it1;
      BlocksMat::iterator2 it2;
      // iterators through m
      BlocksMat::const_iterator1 itM1 = mB._mat->begin1();
      BlocksMat::const_iterator2 itM2;

      for(it1 = _mat->begin1(); it1 != _mat->end1(); ++it1)
      {
        itM2 = itM1.begin();
        for(it2 = it1.begin(); it2 != it1.end(); ++it2)
        {
          (**it2) -= (**itM2);
          itM2++; // increment column pos. in m.
        }
        itM1++; // increment row pos. in m.
      }
    }
    else
    {
      for(unsigned int i = 0; i < _dimRow; ++i)
        for(unsigned int j = 0; j < _dimCol; ++j)
          (*this)(i, j) -= m(i, j);
    }
  }
  else // if m is a SimpleMatrix
  {
    unsigned int indRow = 0, indCol = 0;
    subSimple(indRow, indCol, m); // a sub-block of m is subtracted to each block of this.
  }
  return *this;
}

void BlockMatrix::trans()
{
  THROW_EXCEPTION("not yet implemented.");
}

void BlockMatrix::trans(const SiconosMatrix &m)
{
  THROW_EXCEPTION("not yet implemented.");
}

void BlockMatrix::PLUFactorizationInPlace()
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}
void BlockMatrix::Factorize()
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}

void BlockMatrix::PLUInverseInPlace()
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}

void BlockMatrix::PLUForwardBackwardInPlace(SiconosMatrix &B)
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}
void BlockMatrix::Solve(SiconosMatrix &B)
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}

void BlockMatrix::PLUForwardBackwardInPlace(SiconosVector &B)
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}
void BlockMatrix::Solve(SiconosVector &B)
{
  THROW_EXCEPTION("not yet implemented for Block Matrices.");
}

SP::SiconosMatrix BlockMatrix::block(unsigned int row, unsigned int col)
{
  return (*_mat)(row, col);
}

SPC::SiconosMatrix BlockMatrix::block(unsigned int row, unsigned int col) const
{
  return std::shared_ptr<SiconosMatrix>((*_mat)(row, col));
}

size_t BlockMatrix::nnz(double tol)
{
  size_t nnz = 0;
  BlocksMat::iterator1 it;
  BlocksMat::iterator2 it2;
  for(it = _mat->begin1(); it != _mat->end1(); ++it)
  {
    for(it2 = it.begin(); it2 != it.end(); ++it2)
      nnz += (**it2).nnz();
  }
  return nnz;
}

ext: cpp

max_stars_repo_name: fperignon/sandbox

max_issues_repo_name: fperignon/sandbox

max_forks_repo_name: fperignon/sandbox


text:
#include <components/python/python.hpp>

#include <iostream>

#include <boost/uuid/uuid_io.hpp>
#include <boost/algorithm/string.hpp>
#include <pybind11/pybind11.h>
#include <nlohmann/json.hpp>

#include <detail/context.hpp>
#include <detail/context_manager.hpp>
#include <detail/data_set.hpp>
#include <detail/file_manager.hpp>
#include <detail/file_system.hpp>

#include <../log/log.hpp>
using json = nlohmann::json;

namespace components {

    constexpr static char init_script[] = R"__(
        import sys, os
        from importlib import import_module

        class Worker(pyrocketjoe.celery.apps.Worker):
            def __init__(self, app) -> None:
                super().__init__(app)

        pyrocketjoe.celery.apps.Worker = Worker

        sys.modules['pyrocketjoe'] = pyrocketjoe
    )__";

    constexpr static char load_script[] = R"__(
        import sys, os
        from importlib import import_module

        sys.path.insert(0, os.path.dirname(path))
        module_name, _ = os.path.splitext(path)
        import_module(os.path.basename(module_name))
    )__";

    namespace nl = nlohmann;
    using namespace py::literals;

    python_t::python_t(
        const components::python_sandbox_configuration& configuration,
        log_t& log)
        : python_{}
        , pyrocketjoe_{"pyrocketjoe"}
        , file_manager_{std::make_unique<python::detail::file_manager>()}
        , context_manager_{std::make_unique<python::detail::context_manager>(*file_manager_)}{
        log_ = log.clone();
        log_.info("processing env python start ");
        log_.info(fmt::format("Mode : {0}", configuration.mode_));
        script_path_ = configuration.script_path_;
        if (!configuration.jupyter_connection_path_.empty()) {
            log_.info(fmt::format("jupyter connection path : {0}", configuration.jupyter_connection_path_.string()));
        }

        if (mode_ == sandbox_mode_t::script) {
            init();
        }

        log_.info("processing env python finish ");
    }

    auto python_t::init() -> void {
        python::detail::add_file_system(pyrocketjoe_, file_manager_.get());
        ///python::detail::add_mapreduce(pyrocketjoe, context_manager_.get());

        log_.info("auto python_t::init() -> void {");
    }

    python_t::~python_t() = default;

    void python_t::run_script(const std::vector<std::string>& args) {
        /// TODO: alternative PySys_SetArgv https://stackoverflow.com/questions/18245140/how-do-you-use-the-python3-c-api-for-a-command-line-driven-app
        py::list tmp;

        auto it = args.begin();
        auto end = args.end();

        it = std::next(it);

        for (; it != end; ++it) {
            tmp.append(*it);
        }

        py::module::import("sys").add_object("argv", tmp);
        py::exec(load_script, py::globals(), py::dict("path"_a = script_path_.string()));
    }

    void python_t::call_task(const std::string& task, const std::string& body) {
        std::vector<std::string> parts;
        boost::algorithm::split(parts, task, boost::algorithm::is_any_of("."));
        if (parts.size() != 2) {
            throw std::runtime_error("Task ID is invalid");
        }

        auto jbody = json::parse(body);
        if (!jbody.is_array()) {
            throw std::runtime_error("Task body is not JSON array");
        }
        auto jargs = jbody[0];
        if (!jargs.is_array()) {
            throw std::runtime_error("Task body 1st item is not JSON array");
        }
        py::list args;
        for (auto& jarg : jargs) {
            if (jarg.is_number()) {
                args.append(jarg.get<int>());
            } else if (jarg.is_string()) {
                args.append(jarg.get<std::string>());
            }
        }

        py::module::import(parts[0].c_str()).attr(parts[1].c_str())(*args);
    }

} // namespace components

max_stars_repo_name: jinntechio/RocketJoe

max_issues_repo_name: jinntechio/RocketJoe

max_forks_repo_name: jinntechio/RocketJoe