00001
00002
00003
00004
00005
00006
00007
00008
00009 #ifndef THEPEG_CompSelector_H
00010 #define THEPEG_CompSelector_H
00011
00012
00013
00014
00015 #include "ThePEG/Utilities/Selector.h"
00016
00017 namespace ThePEG {
00018
00046 template <typename T, typename WeightType = double>
00047 class CompSelector {
00048
00049 public:
00050
00058 CompSelector(double newMargin = 1.1, double newTolerance = 1.0e-6)
00059 : N(0), last(), theMargin(newMargin), theTolerance(newTolerance) {}
00061
00062 public:
00063
00075 WeightType insert(WeightType d, const T & t) {
00076 reset();
00077 return selector.insert(d, t);
00078 }
00079
00092 template <typename RNDGEN>
00093 T & select(RNDGEN & rnd) throw(range_error) {
00094 ++N;
00095 if ( !compensating() ) last = selector.select(rnd);
00096 return last;
00097 }
00098
00105 WeightType reweight(double & weight) {
00106 if ( abs(weight) > 1.0 + tolerance() ) {
00107
00108
00109 WeightType oldtot = selector.sum();
00110 WeightType oldmax = oldtot - selector.erase(last);
00111 WeightType newmax = oldmax*abs(weight)*margin();
00112 WeightType newtot = selector.insert(newmax, last);
00113 double rat = newmax/oldmax;
00114
00115
00116 Level level;
00117 level.weight = 1.0/rat;
00118 level.lastN = long(N*newtot/oldtot);
00119
00120
00121
00122 for ( int i = 0, M = levels.size(); i < M; ++i ) {
00123 levels[i].lastN = long(levels[i].lastN*newtot/oldtot);
00124 levels[i].weight /= rat;
00125 }
00126 levels.push_back(level);
00127 weight /= rat;
00128 return newmax;
00129 }
00130
00131
00132
00133 if ( compensating() ) if ( abs(weight) < levels.back().weight ) weight = 0.0;
00134
00135 return WeightType();
00136 }
00137
00142 void reset() {
00143 N = 0;
00144 levels.clear();
00145 last = T();
00146 }
00147
00151 void clear() {
00152 selector.clear();
00153 reset();
00154 }
00155
00160 void margin(double m) { theMargin = m; }
00161
00166 void tolerance(double t) { theTolerance = t; }
00168
00169
00175 bool compensating() {
00176
00177 while ( levels.size() && levels.back().lastN < N ) levels.pop_back();
00178 return !levels.empty();
00179 }
00180
00185 long compleft() const { return levels.empty()? 0: levels.back().lastN - N; }
00186
00193 WeightType sum() const { return selector.sum(); }
00194
00199 double margin() const { return theMargin; }
00200
00205 double tolerance() const { return theTolerance; }
00207
00213 template <typename OStream>
00214 void output(OStream & os) const {
00215 os << selector << N << last << theMargin << theTolerance << levels.size();
00216 for ( int i = 0, M = levels.size(); i < M; ++i )
00217 os << levels[i].lastN << levels[i].weight;
00218 }
00219
00223 template <typename IStream>
00224 void input(IStream & is) {
00225 long M;
00226 is >> selector >> N >> last >> theMargin >> theTolerance >> M;
00227 levels.resize(M);
00228 for ( int i = 0; i < M; ++i ) is >> levels[i].lastN >> levels[i].weight;
00229 }
00231
00232 private:
00233
00237 struct Level {
00238
00243 long lastN;
00244
00248 double weight;
00249
00250 };
00251
00252 private:
00253
00257 Selector<T,WeightType> selector;
00258
00262 long N;
00263
00267 T last;
00268
00273 double theMargin;
00274
00279 double theTolerance;
00280
00284 vector<Level> levels;
00285
00286 };
00287
00291 template <typename OStream, typename T, typename WeightType>
00292 inline OStream & operator<<(OStream & os,
00293 const CompSelector<T,WeightType> & s) {
00294 s.output(os);
00295 return os;
00296 }
00297
00301 template <typename IStream, typename T, typename WeightType>
00302 inline IStream & operator>>(IStream & is,
00303 CompSelector<T,WeightType> & s) {
00304 s.input(is);
00305 return is;
00306 }
00307
00308 }
00309
00310 #endif