Page MenuHomeHEPForge
Files F7877580

No OneTemporary
Actions

View Options

diff --git a/ACDC/ACDCGen.icc b/ACDC/ACDCGen.icc
--- a/ACDC/ACDCGen.icc
+++ b/ACDC/ACDCGen.icc
@@ -1,822 +1,822 @@
// -*- C++ -*-
//
// ACDCGen.icc is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2011 Leif Lonnblad
//
// ThePEG is licenced under version 2 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 )
if ( thePrimaryCells[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 map<double,DVector> PointMap;
+ 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[val] = x;
+ 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 - sumw*sumw/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];
os << thePrimaryCells[theLast]->getIndex(theLastCell);
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;
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-diff
Expires
Tue, Nov 19, 3:50 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805025
Default Alt Text
(25 KB)

Event Timeline