ACDCGen.h

// -*- C++ -*-
//
// ACDCGen.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef ACDCGen_H
#define ACDCGen_H

#include <algorithm>
#include "ACDCGenConfig.h"
#include "ACDCTraits.h"
#include "ACDCGenCell.h"
#include "ThePEG/Utilities/Exception.h"

namespace ACDCGenerator {

template <typename Rnd, typename FncPtr>
class ACDCGen {

public:

  typedef Rnd RndType;
  typedef ACDCRandomTraits<RndType> RndTraits;
  typedef FncPtr FncPtrType;
  typedef vector<ACDCGenCell*> CellVector;
  typedef vector<FncPtrType> FncVector;
  typedef vector<DimType> DimVector;
  typedef DimVector::size_type size_type;
  typedef ACDCFncTraits<FncPtrType> FncTraits;

public:

  inline ACDCGen(Rnd * r);

  inline ACDCGen();

  inline ~ACDCGen();

  inline bool addFunction(DimType dim, FncPtrType f, double maxrat = -1.0);

  inline void clear();

public:

  inline FncPtrType generate();

  inline void reject();

  inline const DVector & lastPoint() const;

  inline double lastF() const;

  inline FncPtrType lastFunction() const;

  inline size_type last() const;

public:

  inline void eps(double newEps);

  inline void margin(double newMargin);

  inline void nTry(size_type newNTry);

  inline void maxTry(long);

public:

  inline double integral(FncPtrType f = FncPtrType()) const;

  inline double integralErr(FncPtrType f = FncPtrType()) const;

  inline long n() const;

  inline long N() const;

  inline double efficiency() const;

  inline int nBins() const;

  inline int depth() const;

  inline double maxInt() const;

  inline double eps() const;

  inline double margin() const;

  inline size_type nTry() const;

  inline long maxTry() const;

  inline bool cheapRandom() const;

  inline size_type size() const;

  inline bool compensating();

  inline long compleft() const;

  vector<ACDCGenCellInfo> extractCellInfo() const;

public:

  inline void cheapRandom(bool b);

  inline void setRnd(Rnd * r);

  inline double rnd() const;

  inline double rnd(double lo, double up) const;

  inline void rnd(const DVector & lo, const DVector & up, DVector & r)const;

  inline void rnd(DimType D, DVector & r) const;

  inline long rndInt(long x) const;

public:

  template <typename POStream>
  void output(POStream &) const;

  template <typename PIStream>
  void input(PIStream &);

private:

  inline double doMaxInt();

  inline const FncVector & functions() const;

  inline FncPtrType function(size_type i) const;

  inline const DimVector & dimensions() const;

  inline DimType dimension(size_type i) const;

  inline DimType lastDimension() const;

  inline const CellVector & cells() const;

  inline ACDCGenCell * cell(size_type i) const;

  inline ACDCGenCell * lastPrimary() const;

  inline const DVector & sumMaxInts() const;

  inline ACDCGenCell * lastCell() const;


  inline void chooseCell(DVector & lo, DVector & up);

  inline void  compensate(const DVector & lo, const DVector & up);

private:

  RndType * theRnd;
  
  long theNAcc;

  long theN;

  vector<long> theNI;

  DVector theSumW;

  DVector theSumW2;

  double theEps;

  double theMargin;

  size_type theNTry;

  long theMaxTry;

  bool useCheapRandom;

  FncVector theFunctions;

  DimVector theDimensions;

  CellVector thePrimaryCells;

  DVector theSumMaxInts;

  size_type theLast;

  ACDCGenCell * theLastCell;

  DVector theLastPoint;

  double theLastF;

  struct Level {

    long lastN;

    double g;

    ACDCGenCell * cell;

    size_type index;

    DVector up;
    DVector lo;

  };

  typedef vector<Level> LevelVector;

  LevelVector levels;


  struct Slicer {

    Slicer(DimType, ACDCGen &, const DVector &, const DVector &);

    Slicer(DimType Din, const Slicer & s, ACDCGenCell * cellin,
           const DVector & loin, const DVector & xselin, const DVector & upin,
           double fselin);

    ~Slicer();

    void divideandconquer();

    void init();

    void slice();

    double shiftmaxmin();

    void dohalf(DimType);

    void checkdiag(ACDCGenCell * cell, DimType d, double lod, double upd);

    DimType D;

    DVector lo;
    DVector up;

    DVector xcl;
    DVector xcu;

    DVector xhl;
    DVector xhu;

    DVector fhl;

    DVector fhu;

    DVector xsel;

    double fsel;

    ACDCGenCell * current;

    ACDCGenCell * first;

    DVector firstlo;
    DVector firstup;

    FncPtr f;

    double epsilon;

    double margin;

    multimap<double,DimType> rateslice;

    double minf;

    bool wholecomp;

  };

public:

  static size_type maxsize;

private:

  ACDCGen(const ACDCGen &);

  ACDCGen & operator=(const ACDCGen &) = delete;

};

}

#include "ACDCGen.icc"

#endif