CellGrid.h

// -*- C++ -*-
//
// CellGrid.hpp is a part of ExSample
// Copyright (C) 2012-2019 Simon Platzer, The Herwig Collaboration
//
// ExSample is licenced under version 3 of the GPL, see COPYING for details.
//
 
#ifndef EXSAMPLE_CellGrid_hpp_included
#define EXSAMPLE_CellGrid_hpp_included
 
#include <vector>
#include <utility>
#include <set>
#include <map>
#include <cstdlib>
#include <cassert>
#include <cmath>
 
#define units
 
#ifdef units
#include <iostream>
#include <string>
#endif
 
#include "Herwig/Utilities/XML/Element.h"
 
namespace ExSample {
 
  inline double sqr(double x) {
    return x*x;
  }
 
  class CellGrid {
 
  public:
 
    CellGrid()
      : theVolumeOrIntegral(0.0), theWeight(0.0) {}
 
    CellGrid(const std::vector<double>& newLowerLeft,
             const std::vector<double>& newUpperRight,
             double newWeight = 0.0);
 
    virtual CellGrid* makeInstance() const;
 
    virtual CellGrid* makeInstance(const std::vector<double>& newLowerLeft,
                                   const std::vector<double>& newUpperRight,
                                   double newWeight = 0.0) const;
 
    virtual ~CellGrid();
 
  public:
 
    void boundaries(const std::vector<double>& newLowerLeft,
                    const std::vector<double>& newUpperRight);
 
    const std::vector<double>& lowerLeft() const { return theLowerLeft; }
 
    const std::vector<double>& upperRight() const { return theUpperRight; }
 
    const std::vector<bool>& upperBoundInclusive() const { return theUpperBoundInclusive; }
 
    double volume(const std::vector<double>& lowerLeft,
                  const std::vector<double>& upperRight) const;
 
    bool isLeaf() const { return theChildren.empty(); }
 
    std::size_t depth() const;
 
    std::size_t size() const;
 
    virtual void split(std::size_t newSplitDimension, double newSplitCoordinate);
 
    std::pair<std::size_t,double> splitPoint() const;
 
    const CellGrid& firstChild() const;
 
    CellGrid& firstChild();
 
    const CellGrid& secondChild() const;
 
    CellGrid& secondChild();
 
    void splitCoordinates(std::size_t, std::set<double>&) const;
 
  public:
 
    bool active() const;
 
    double volume() const;
 
    double integral() const;
 
    double weight() const { return theWeight; }
 
    void weight(double w);
 
    void updateIntegral();
 
    void minimumSelection(double p = 0.1);
 
    bool contains(const std::vector<double>& point,
                  const std::vector<bool>& parameterFlags) const;
 
    double nonParametricVolume(const std::vector<double>& point,
                               const std::vector<bool>& parameterFlags) const;
 
    void updateIntegral(const std::vector<double>& point,
                        const std::vector<bool>& parameterFlags,
                        std::vector<bool>::iterator hashPosition);
 
    double projectInterval(const std::pair<double,double>& interval,
                           std::size_t dimension) const;
 
    std::map<std::pair<double,double>,double> 
    project(std::pair<double,double> interval,
            std::size_t dimension) const;
 
  public:
 
    virtual void fromXML(const XML::Element&);
 
    virtual XML::Element toXML() const;
 
  private:
 
    void doMinimumSelection(double r,
                            double ref);
 
    std::vector<double> theLowerLeft;
 
    std::vector<double> theUpperRight;
 
    std::vector<bool> theUpperBoundInclusive;
 
    double theVolumeOrIntegral;
 
    double theWeight;
 
    std::vector<CellGrid*> theChildren;
 
#ifdef units
 
  public:
 
    template<class RndGenerator>
    void randomGrid(RndGenerator& rnd, 
                    std::size_t maxDepth,
                    bool forceSplit = true) {
      if ( maxDepth > 0 ) {
        if ( !forceSplit )
          if ( rnd.rnd() < 0.5 )
            return;
        std::size_t dimension = (std::size_t)(std::floor(rnd.rnd()*lowerLeft().size()));
        double point = 0.5*(lowerLeft()[dimension]+upperRight()[dimension]);
        split(dimension,point);
        firstChild().weight(rnd.rnd());
        secondChild().weight(rnd.rnd());
        firstChild().randomGrid(rnd,maxDepth-1,false);
        secondChild().randomGrid(rnd,maxDepth-1,false);
      }
    }
 
    void dumpToC(std::ostream& os,
                 const std::string& name) const;
 
    void dumpPartToC(std::ostream& os,
                     std::string prefix = "") const;
 
#endif
 
  };
 
}
 
#endif // EXSAMPLE_CellGrid_hpp_included