ACDCGen.icc
No OneTemporary
Actions

Size

24 KB

Subscribers

None

ACDCGen.icc
View Options

	// -- C++ --
	//
	// ACDCGen.icc is a part of ThePEG - Toolkit for HEP Event Generation
	// Copyright (C) 1999-2017 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.
	//

	namespace ACDCGenerator {

	template <typename Rnd, typename FncPtr>
	inline ACDCGen<Rnd,FncPtr>::ACDCGen(Rnd * r)
	: theRnd(r), theNAcc(0), theN(0), theNI(1, 0),
	theSumW(1, 0.0), theSumW2(1, 0.0),
	theEps(100*std::numeric_limits<double>::epsilon()), theMargin(1.1),
	theNTry(100), theMaxTry(10000), useCheapRandom(false), theFunctions(1),
	theDimensions(1, 0), thePrimaryCells(1), theSumMaxInts(1, 0.0), theLast(0),
	theLastCell(0), theLastF(0.0) {
	maxsize = 0;
	}

	template <typename Rnd, typename FncPtr>
	inline ACDCGen<Rnd,FncPtr>::ACDCGen()
	: theRnd(0), theNAcc(0), theN(0), theNI(1, 0),
	theSumW(1, 0.0), theSumW2(1, 0.0),
	theEps(100*std::numeric_limits<double>::epsilon()), theMargin(1.1),
	theNTry(100), theMaxTry(10000), useCheapRandom(false), theFunctions(1),
	theDimensions(1, 0), thePrimaryCells(1), theSumMaxInts(1, 0.0), theLast(0),
	theLastCell(0), theLastF(0.0) {
	maxsize = 0;
	}

	template <typename Rnd, typename FncPtr>
	inline ACDCGen<Rnd,FncPtr>::~ACDCGen() {
	clear();
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::setRnd(Rnd * r) {
	theRnd = r;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::clear() {
	theNAcc = 0;
	theN = 0;
	theNI = vector<long>(1, 0);
	theSumW = DVector(1, 0.0);
	theSumW2 = DVector(1, 0.0);
	theFunctions = FncVector(1);
	theDimensions = DimVector(1, 0);
	for ( int i = 0, N = thePrimaryCells.size(); i < N; ++i )
	delete thePrimaryCells[i];
	thePrimaryCells = CellVector(1);
	theSumMaxInts = DVector(1, 0.0);
	theLast = 0;
	theLastCell = 0;
	theLastPoint.clear();
	theLastF = 0.0;
	levels.clear();
	}

	template <typename Rnd, typename FncPtr>
	inline bool ACDCGen<Rnd,FncPtr>::
	addFunction(DimType dim, FncPtrType fnc, double maxrat) {
	if ( maxrat < 0.0 ) maxrat = 1.0/nTry();
	typedef multimap<double,DVector> PointMap;
	theLast = theFunctions.size();
	theFunctions.push_back(fnc);
	theNI.push_back(0);
	theSumW.push_back(0.0);
	theSumW2.push_back(0.0);
	theDimensions.push_back(dim);

	// Generate nTry() points with non-zero function value
	DVector x(dim);
	PointMap pmap;
	long itry = 0;
	while ( pmap.size() < nTry() ) {
	if ( ++itry > maxTry() ) {
	thePrimaryCells.push_back(new ACDCGenCell(0.0));
	theSumMaxInts.push_back(theSumMaxInts.back() + cells().back()->doMaxInt());
	return false;
	}
	rnd(dim, x);
	double val = FncTraits::value(fnc, x);
	if ( val > 0.0 ) {
	pmap.insert(make_pair(val, x));
	itry = 0;
	}
	}

	// Create the root cell and set its overestimated function value to
	// the smallest non-zero value found
	double minf = pmap.begin()->first;
	double maxf = (--pmap.end())->first;
	minf = max(minf, maxrat*maxf);
	// thePrimaryCells.push_back(new ACDCGenCell(pmap.begin()->first));
	thePrimaryCells.push_back(new ACDCGenCell(minf));
	theLastF = pmap.begin()->first;
	pmap.erase(pmap.begin());
	theSumMaxInts.push_back(theSumMaxInts.back() + cells().back()->doMaxInt());

	// Start the divide-and-conquer procedure using the point with the
	// highest function value found.
	theLastF = (--pmap.end())->first;
	theLastPoint = (--pmap.end())->second;
	pmap.erase(--pmap.end());
	DVector up(dim, 1.0);
	DVector lo(dim, 0.0);
	theLastCell = cells().back()->getCell(lo, lastPoint(), up);
	if ( lastF() > lastCell()->g() ) {
	compensate(lo, up);
	levels.clear();
	}

	// For each other generated point check that it's below the
	// overestimated value of the corresponding cell. If not start the
	// divide-and-conquer using that point.
	while ( !pmap.empty() ) {
	theLastPoint = pmap.begin()->second;
	theLastF = pmap.begin()->first;
	pmap.erase(pmap.begin());
	DVector up(dim, 1.0);
	DVector lo(dim, 0.0);
	theLastCell = cells().back()->getCell(lo, lastPoint(), up);
	if ( lastF() > lastCell()->g() ) {
	compensate(lo, up);
	levels.clear();
	}
	}
	// cells().back()->smooth(1.0/nTry());
	return true;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::chooseCell(DVector & lo, DVector & up) {
	if ( compensating() ) {
	// If we are compensating, we must choose the cell to be compensated.
	up = levels.back().up;
	lo = levels.back().lo;
	theLastCell = levels.back().cell;
	theLast = levels.back().index;
	} else {
	// Otherwise, first choose the function to be used and choose the
	// corresponding root cell.
	theLast = upper_bound(sumMaxInts().begin(), sumMaxInts().end(),
	rnd()*sumMaxInts().back())
	- sumMaxInts().begin();
	if(theLast>=sumMaxInts().size()) {
	throw ThePEG::Exception() << "Selected a function outside the allowed range"
	<< " in ACDCGen::chooseCell(). This is usually due"
	<< " to a floating point error (nan or inf) in the"
	<< " calculation of the weight"
	<< ThePEG::Exception::abortnow;
	}
	up = DVector(lastDimension(), 1.0);
	lo = DVector(lastDimension(), 0.0);
	theLastCell = lastPrimary();
	}

	// Now select randomly a sub-cell of the chosen cell
	if ( cheapRandom() ) {
	theLastCell = lastCell()->generate(lo, up, theRnd);
	} else {
	DVector rndv(lastDimension());
	rnd(lastDimension(), rndv);
	theLastCell = lastCell()->generate(lo, up, rndv);
	}
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::FncPtrType
	ACDCGen<Rnd,FncPtr>::generate() {
	long itry = 0;
	while ( true ) {
	if ( ++itry > maxTry() ) return FncPtrType();
	++theN;

	// First choose a function and a cell to generate in.
	DVector up;
	DVector lo;
	chooseCell(lo, up);

	// Now choose a point in that cell according to a flat distribution.
	DimType D = lastDimension();
	theLastPoint.resize(D);
	for ( DimType d = 0; d < D; ++d ) theLastPoint[d] = rnd(lo[d], up[d]);

	// Calculate the function value in this point
	theLastF = FncTraits::value(lastFunction(), theLastPoint);
	if ( theLastF <= 0.0 ) continue;

	// If we are compensating we require the function value to be
	// above the previous overestimate of the function.
	if ( compensating() && lastF() < levels.back().g ) theLastF = 0.0;

	double w = lastF()/lastCell()->g();
	if ( w > 1.0 ) {
	// If the value was above the overestimate then we must start
	// the compensation procedure and the curren point is disregarded.
	--theN;
	compensate(lo, up);
	continue;
	}

	// Accept the point according to the ratio of the true and
	// overestimated function value.
	theSumW[last()] += w;
	theSumW2[last()] += w*w;
	++theNI[last()];
	if ( w > rnd() ) {
	++theNAcc;
	return lastFunction();
	}
	}
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::reject() {
	theSumW[last()] -= 1.0;
	theSumW2[last()] -= 1.0;
	--theNAcc;
	}

	template <typename Rnd, typename FncPtr>
	inline bool ACDCGen<Rnd,FncPtr>::compensating() {
	while ( levels.size() && levels.back().lastN < N() ) levels.pop_back();
	// Leave all levels which has reached there 'expiry date'.

	return !levels.empty();
	}

	template <typename Rnd, typename FncPtr>
	inline long ACDCGen<Rnd,FncPtr>::compleft() const {
	if ( levels.empty() ) return 0;
	long left = 1;
	for ( int i = 0, Ni = levels.size(); i < Ni; ++i )
	left = max(left, levels[i].lastN - N());

	return left;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::
	compensate(const DVector & lo, const DVector & up) {
	//Save the previous overestimated integral and create a new
	//compensation level.
	double i0 = maxInt();
	Level level;
	level.g = lastCell()->g();

	// Start the divide-and-conquer algorithm slicing up the selected
	// cell and specify it as the cell to compensate.
	Slicer slicer(lastDimension(), *this, lo, up);
	level.cell = slicer.first;
	level.index = last();
	level.up = slicer.firstup;
	level.lo = slicer.firstlo;

	// Now calculate the the new overestimated total integral and
	// calculate the number of attempted points needed to compensate.
	double rat = (doMaxInt())/i0;
	level.lastN = long(N()*rat);

	// If we are already compensating increase also the previous
	// compensation levels.
	for ( size_type i = 0; i < levels.size(); ++i )
	levels[i].lastN = long(levels[i].lastN*rat);
	levels.insert(levels.end(), level);
	maxsize = std::max(maxsize, levels.size());
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::eps(double newEps) {
	theEps = newEps;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::margin(double newMargin) {
	theMargin = newMargin;
	}

	template <typename Rnd, typename FncPtr>
	inline long ACDCGen<Rnd,FncPtr>::N() const {
	return theN;
	}

	template <typename Rnd, typename FncPtr>
	inline long ACDCGen<Rnd,FncPtr>::n() const {
	return theNAcc;
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::size_type
	ACDCGen<Rnd,FncPtr>::nTry() const {
	return theNTry;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::nTry(size_type newNTry) {
	theNTry = newNTry;
	}

	template <typename Rnd, typename FncPtr>
	inline long ACDCGen<Rnd,FncPtr>::maxTry() const {
	return theMaxTry;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::maxTry(long newMaxTry) {
	theMaxTry = newMaxTry;
	}

	template <typename Rnd, typename FncPtr>
	inline bool ACDCGen<Rnd,FncPtr>::cheapRandom() const {
	return useCheapRandom;
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::cheapRandom(bool b) {
	useCheapRandom = b;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::maxInt() const {
	return theSumMaxInts.back();
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::doMaxInt() {
	for ( size_type i = 1, imax = functions().size(); i < imax; ++i )
	theSumMaxInts[i] = sumMaxInts()[i - 1] + cells()[i]->doMaxInt();
	return maxInt();
	}

	template <typename Rnd, typename FncPtr>
	inline int ACDCGen<Rnd,FncPtr>::nBins() const {
	int sum = 0;
	for ( size_type i = 1; i < functions().size(); ++i )
	sum += cell(i)->nBins();
	return sum;
	}

	template <typename Rnd, typename FncPtr>
	inline int ACDCGen<Rnd,FncPtr>::depth() const {
	int mx = 0;
	for ( size_type i = 1; i < functions().size(); ++i )
	mx = max(mx, cell(i)->depth());
	return mx;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::efficiency() const {
	return N() > 0? double(n())/double(N()): 0.0;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::integral(FncPtrType f) const {
	if ( N() <= 0 ) return maxInt();
	double sumw = 0.0;
	for ( size_type i = 1; i < functions().size(); ++i )
	if ( functions()[i] == f \|\| !f ) sumw += theSumW[i];
	return maxInt()*sumw/N();
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::integralErr(FncPtrType f) const {
	if ( N() <= 0 ) return maxInt();
	double sumw2 = 0.0;
	double sumw = 0.0;
	for ( size_type i = 1; i < functions().size(); ++i )
	if ( functions()[i] == f \|\| !f ) {
	sumw2 += theSumW2[i];
	sumw += theSumW[i];
	}
	if ( f ) return maxInt()*sqrt(sumw2)/N();
	return maxInt()sqrt(max(0.,sumw2 - sumwsumw/N()))/N();
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::eps() const {
	return theEps;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::margin() const {
	return theMargin;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::rnd() const {
	return RndTraits::rnd(theRnd);
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::rnd(double lo, double up) const {
	return RndTraits::rnd(theRnd, lo, up);
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::
	rnd(const DVector & lo, const DVector & up, DVector & r) const {
	RndTraits::rnd(theRnd, lo.begin(), lo.end(), up.begin(), r.begin());
	}

	template <typename Rnd, typename FncPtr>
	inline void ACDCGen<Rnd,FncPtr>::
	rnd(DimType D, DVector & r) const {
	RndTraits::rnd(theRnd, D, r.begin());
	}

	template <typename Rnd, typename FncPtr>
	inline long ACDCGen<Rnd,FncPtr>::rndInt(long x) const {
	return RndTraits::rndInt(theRnd, x);
	}

	template <typename Rnd, typename FncPtr>
	inline const typename ACDCGen<Rnd,FncPtr>::FncVector &
	ACDCGen<Rnd,FncPtr>::functions() const {
	return theFunctions;
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::FncPtrType
	ACDCGen<Rnd,FncPtr>::function(size_type i) const {
	return functions()[i];
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::FncPtrType
	ACDCGen<Rnd,FncPtr>::lastFunction() const {
	return function(last());
	}

	template <typename Rnd, typename FncPtr>
	inline const typename ACDCGen<Rnd,FncPtr>::DimVector &
	ACDCGen<Rnd,FncPtr>::dimensions() const {
	return theDimensions;
	}

	template <typename Rnd, typename FncPtr>
	inline DimType ACDCGen<Rnd,FncPtr>::dimension(size_type i) const {
	return dimensions()[i];
	}

	template <typename Rnd, typename FncPtr>
	inline DimType ACDCGen<Rnd,FncPtr>::lastDimension() const {
	return dimension(last());
	}

	template <typename Rnd, typename FncPtr>
	inline const typename ACDCGen<Rnd,FncPtr>::CellVector &
	ACDCGen<Rnd,FncPtr>::cells() const {
	return thePrimaryCells;
	}

	template <typename Rnd, typename FncPtr>
	inline ACDCGenCell * ACDCGen<Rnd,FncPtr>::cell(size_type i) const {
	return cells()[i];
	}

	template <typename Rnd, typename FncPtr>
	inline ACDCGenCell * ACDCGen<Rnd,FncPtr>::lastPrimary() const {
	return cell(last());
	}

	template <typename Rnd, typename FncPtr>
	inline const DVector & ACDCGen<Rnd,FncPtr>::sumMaxInts() const {
	return theSumMaxInts;
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::size_type
	ACDCGen<Rnd,FncPtr>::last() const {
	return theLast;
	}

	template <typename Rnd, typename FncPtr>
	inline typename ACDCGen<Rnd,FncPtr>::size_type
	ACDCGen<Rnd,FncPtr>::size() const {
	return cells().size() - 1;
	}

	template <typename Rnd, typename FncPtr>
	inline ACDCGenCell * ACDCGen<Rnd,FncPtr>::lastCell() const {
	return theLastCell;
	}

	template <typename Rnd, typename FncPtr>
	inline const DVector & ACDCGen<Rnd,FncPtr>::lastPoint() const {
	return theLastPoint;
	}

	template <typename Rnd, typename FncPtr>
	inline double ACDCGen<Rnd,FncPtr>::lastF() const {
	return theLastF;
	}

	template <typename Rnd, typename FncPtr>
	typename ACDCGen<Rnd,FncPtr>::size_type ACDCGen<Rnd,FncPtr>::maxsize = 0;

	template <typename Rnd, typename FncPtr>
	vector<ACDCGenCellInfo> ACDCGen<Rnd,FncPtr>::extractCellInfo() const {
	vector<ACDCGenCellInfo> ret;
	for ( size_type i = 1; i < cells().size(); ++i ) {
	DVector lo(dimension(i), 0.0);
	DVector up(dimension(i), 1.0);
	cell(i)->extract(lo, up, ret);
	}
	return ret;
	}

	template <typename Rnd, typename FncPtr>
	ACDCGen<Rnd,FncPtr>::Slicer::
	Slicer(DimType Din, ACDCGen & gen, const DVector & loin, const DVector & upin)
	: D(Din), lo(loin), up(upin), xcl(loin), xcu(upin), xhl(loin), xhu(upin),
	fhl(Din, 0.0), fhu(Din, 0.0), xsel(gen.lastPoint()), fsel(gen.lastF()),
	current(gen.lastCell()), first(gen.lastCell()),
	firstlo(loin), firstup(upin),f(gen.lastFunction()),
	epsilon(gen.eps()), margin(gen.margin()), minf(0.0), wholecomp(false) {
	divideandconquer();
	}

	template <typename Rnd, typename FncPtr>
	ACDCGen<Rnd,FncPtr>::Slicer::~Slicer() {
	// Added for debugging purposes.
	}

	template <typename Rnd, typename FncPtr>
	void ACDCGen<Rnd,FncPtr>::Slicer::divideandconquer() {
	// If the current function value was just a little above the
	// overestimate, just increase the overestimate of this cell and
	// we're done.
	if ( fsel < current->g()*margin ) {
	current->g(current->g()*margin);
	return;
	}

	// First initialize and slice up the current cell and save the
	// resulting for the compensation procedure.
	init();
	slice();
	if ( !wholecomp ) {
	first = current;
	firstlo = lo;
	firstup = up;
	}

	// Find the largest function value in the current cell and as long
	// as it is above the current overestimate, increase the
	// overestimate and repeat the slicing.
	while ( shiftmaxmin() > current->g() ) {
	current->g(minf*margin);
	if ( current->g() > fsel ) return;
	init();
	slice();
	}
	}

	template <typename Rnd, typename FncPtr>
	double ACDCGen<Rnd,FncPtr>::Slicer::shiftmaxmin() {
	// Find the largest diagonal
	DVector test = xsel;
	double scale = 0.0;
	for ( DimType d = 0; d < D; ++d )
	if ( fhl[d] > fsel \|\| fhu[d] > fsel ) scale += 1.0;
	scale = sqrt(scale);
	for ( DimType d = 0; d < D; ++d ) {
	if ( fhl[d] > fsel && fhl[d] > fhu[d] )
	test[d] = test[d] + (xhl[d] - test[d])/scale;
	if ( fhu[d] > fsel && fhu[d] > fhl[d] )
	test[d] = test[d] + (xhu[d] - test[d])/scale;
	}

	// Find the largest value above overestimate
	DimType dsel = -1;
	double x = 0;
	minf = fsel;
	for ( DimType d = 0; d < D; ++d ) {

	// Find the point with the function minimum value above the
	// current overestimate.
	minf = std::min(minf, fhl[d]);
	minf = std::min(minf, fhu[d]);

	// Find points with the maximum function value and shift the
	// current point to it.
	if ( fhu[d] > fsel ) {
	fsel = fhu[d];
	dsel = d;
	x = xhu[d];
	}
	if ( fhl[d] > fsel ) {
	fsel = fhl[d];
	dsel = d;
	x = xhl[d];
	}
	}

	// Check also the largest diagonal
	// double ftest = (*f)(test);
	// if ( ftest > fsel ) {
	// xsel = test;
	// fsel = ftest;
	// dsel = -1;
	// }

	if ( dsel >= 0 ) xsel[dsel] = x;
	minf = std::max(minf, current->g());
	return fsel;
	}

	template <typename Rnd, typename FncPtr>
	void ACDCGen<Rnd,FncPtr>::Slicer::dohalf(DimType d) {
	xcl[d] = lo[d];

	// Find the point closest below the current point in the given
	// dimension with a function value below the current overestimate,
	// also find the one furthest away from the current point with a
	// function value above the current overestimate. Use a crude
	// iterative mid-point selection.
	while ( true ) {
	xhl[d] = (xsel[d] + xcl[d])*0.5;

	std::swap(xsel[d], xhl[d]);
	fhl[d] = FncTraits::value(f, xsel);
	std::swap(xsel[d], xhl[d]);

	if ( fhl[d] > current->g() ) break;
	if ( xsel[d] - xcl[d] < epsilon ) break;

	xcl[d] = xhl[d];
	}

	// Check if the current slicing is worth doing...
	double cut = ( up[d] - xcl[d] )/( up[d] - lo[d] );
	if ( cut < 1.0 - current->g()/fsel && cut > 0.0 )
	rateslice.insert(std::make_pair(cut, -1-d));

	// Find the point closest above the current point in the given
	// dimension with a function value below the current overestimate,
	// also find the one furthest away from the current point with a
	// function value above the current overestimate. Use a crude
	// iterative mid-point selection.
	xcu[d] = up[d];
	while ( true ) {
	xhu[d] = (xsel[d] + xcu[d])*0.5;

	std::swap(xsel[d], xhu[d]);
	fhu[d] = FncTraits::value(f, xsel);
	std::swap(xsel[d], xhu[d]);

	if ( fhu[d] > current->g() ) break;
	if ( xcu[d] - xsel[d] < epsilon ) break;

	xcu[d] = xhu[d];
	}

	// Check if the current slicing is worth doing...
	cut = ( xcu[d] - lo[d] )/( up[d] - lo[d] );
	if ( cut < 1.0 - current->g()/fsel && cut > 0.0 )
	rateslice.insert(std::make_pair(cut, 1+d));

	}

	template <typename Rnd, typename FncPtr>
	void ACDCGen<Rnd,FncPtr>::Slicer::init() {
	for ( DimType d = 0; d < D; ++d ) dohalf(d);
	}

	template <typename Rnd, typename FncPtr>
	void ACDCGen<Rnd,FncPtr>::Slicer::slice() {
	while ( !rateslice.empty() ) {
	// Perform the slicing which reduces the volume of the cell the
	// most first.
	DimType d = rateslice.begin()->second;
	rateslice.erase(rateslice.begin());
	if ( d > 0 ) {
	// Slice from above.
	d = d - 1;
	current->splitme(lo[d], xcu[d], up[d], d);
	checkdiag(current->upper(), d, xcu[d], up[d]);
	current = current->lower();
	up[d] = xcu[d];
	} else {
	// Slice from below..
	d = -d - 1;
	current->splitme(lo[d], xcl[d], up[d], d);
	checkdiag(current->lower(), d, lo[d], xcl[d]);
	current = current->upper();
	lo[d] = xcl[d];
	}
	dohalf(d);
	}
	}

	template <typename Rnd, typename FncPtr>
	void ACDCGen<Rnd,FncPtr>::Slicer::
	checkdiag(ACDCGenCell * cell, DimType dc, double lod, double upd) {
	return;
	// Look at the midpoint in the dc direction in which a cell has been
	// chopped off.
	if ( upd - lod <= epsilon ) return;
	DVector newlo = lo;
	DVector newup = up;
	newlo[dc] = lod;
	newup[dc] = upd;
	DVector newsel = xsel;
	newsel[dc] = 0.5*(lod + upd);
	double newfsel = FncTraits::value(f, newsel);
	double newfh = newfsel;
	DVector newxsel = newsel;
	vector<int> dir(D, 0);

	// For each other direction look at the mid point between the point
	// chosen above and the borders of the cell. Save the point which
	// gives the highest function value.
	for ( DimType d = 0; d < D; ++d ) {
	if ( d == dc ) continue;
	double xdum = 0.5*(newlo[d] + newsel[d]);
	swap(xdum, newsel[d]);
	double fh1 = FncTraits::value(f, newsel);
	if ( fh1 > newfsel ) {
	newfsel = fh1;
	newxsel = newsel;
	}
	if ( fh1 > newfh ) dir[d] = -1;
	swap(xdum, newsel[d]);
	xdum = 0.5*(newsel[d] + newup[d]);
	swap(xdum, newsel[d]);
	double fh2 = FncTraits::value(f, newsel);
	if ( fh2 > newfsel ) {
	newfsel = fh2;
	newxsel = newsel;
	}
	if ( fh2 > newfh && fh2 > fh1 ) dir[d] = 1;
	swap(xdum, newsel[d]);
	}

	// Now check along the diagonal where we found the highest values.
	for ( DimType d = 0; d < D; ++d ) {
	if ( dir[d] == 0 ) continue;
	if ( dir[d] > 0 ) newsel[d] = 0.5*(newsel[d] + newup[d]);
	else newsel[d] = 0.5*(newlo[d] + newsel[d]);
	}
	newfh = FncTraits::value(f, newsel);
	if ( newfh > newfsel ) {
	newfsel = newfh;
	newxsel = newsel;
	}

	if ( newfsel < cell->g() ) return;

	// If this the highest value is above the overestimate, also this
	// cell needs to be divided up and conquered.
	wholecomp = true;
	Slicer dummy(D, *this, cell, newlo, newxsel, newup, newfsel);
	}

	template <typename Rnd, typename FncPtr>
	ACDCGen<Rnd,FncPtr>::Slicer::
	Slicer(DimType Din, const Slicer & s, ACDCGenCell * cellin,
	const DVector & loin, const DVector & xselin, const DVector & upin,
	double fselin)
	: D(Din), lo(loin), up(upin), xcl(loin), xcu(upin), xhl(loin), xhu(upin),
	fhl(Din, 0.0), fhu(Din, 0.0), xsel(xselin), fsel(fselin),
	current(cellin), first(cellin),
	firstlo(loin), firstup(upin),f(s.f),
	epsilon(s.epsilon), margin(s.margin), minf(0.0), wholecomp(false) {
	divideandconquer();
	}

	template <typename Rnd, typename FncPtr>
	template <typename POStream>
	void ACDCGen<Rnd,FncPtr>::output(POStream & os) const {
	os << theNAcc << theN << theEps << theMargin << theNTry << theMaxTry
	<< useCheapRandom << theLast << theLastPoint << theLastF
	<< theFunctions.size() << levels.size();
	for ( int i = 1, N = theFunctions.size(); i < N; ++i )
	os << theFunctions[i] << theDimensions [i] << theSumMaxInts[i]
	<< *thePrimaryCells[i] << theNI[i] << theSumW[i] << theSumW2[i];
	if ( theLast > 0 ) // first entry in thePrimaryCells always points at 0x0
	os << thePrimaryCells[theLast]->getIndex(theLastCell);
	else
	os << -1l;
	for ( int i = 0, N = levels.size(); i < N; ++i )
	os << levels[i].lastN << levels[i].g << levels[i].index
	<< levels[i].up << levels[i].lo
	<< thePrimaryCells[levels[i].index]->getIndex(levels[i].cell);
	}

	template <typename Rnd, typename FncPtr>
	template <typename PIStream>
	void ACDCGen<Rnd,FncPtr>::input(PIStream & is) {
	clear();
	long fsize = 0;
	long lsize = 0;
	is >> theNAcc >> theN >> theEps >> theMargin >> theNTry >> theMaxTry
	>> useCheapRandom >> theLast >> theLastPoint >> theLastF >> fsize >> lsize;
	while ( --fsize ) {
	theFunctions.push_back(FncPtrType());
	theDimensions.push_back(DimType());
	theSumMaxInts.push_back(0.0);
	theNI.push_back(0);
	theSumW.push_back(0.0);
	theSumW2.push_back(0.0);
	thePrimaryCells.push_back(new ACDCGenCell(0.0));
	is >> theFunctions.back() >> theDimensions.back() >> theSumMaxInts.back()
	>> *thePrimaryCells.back() >> theNI.back()
	>> theSumW.back() >> theSumW2.back();
	}
	long index = -1;
	is >> index;
	if ( index == -1 )
	theLastCell = 0x0;
	else
	theLastCell = thePrimaryCells[theLast]->getCell(index);
	while ( lsize-- ) {
	levels.push_back(Level());
	is >> levels.back().lastN >> levels.back().g >> levels.back().index
	>> levels.back().up >> levels.back().lo >> index;
	levels.back().cell = thePrimaryCells[levels.back().index]->getCell(index);
	}
	}

	}

File Metadata

Mime Type: text/x-c++
Expires: Sat, Dec 21, 1:54 PM (15 h, 54 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4022629
Default Alt Text: ACDCGen.icc (24 KB)

ACDCGen.iccNo OneTemporaryActions

ACDCGen.iccView Options

File Metadata

Event Timeline

ACDCGen.icc
No OneTemporary
Actions

ACDCGen.icc
View Options