File indexing completed on 2026-08-06 09:24:22
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 namespace exsample {
0012
0013 template<class Function, class Random>
0014 void exponential_generator<Function,Random>::initialize() {
0015 adaption_info_.dimension = function_->dimension();
0016 adaption_info_.lower_left = function_->support().first;
0017 adaption_info_.upper_right = function_->support().second;
0018 if (adaption_info_.adapt.empty())
0019 adaption_info_.adapt = std::vector<bool>(adaption_info_.dimension,true);
0020 evolution_variable_ = function_->evolution_variable();
0021 evolution_cutoff_ = function_->evolution_cutoff();
0022 sample_variables_ = function_->variable_flags();
0023 sample_other_variables_ = sample_variables_;
0024 sample_other_variables_[evolution_variable_] = false;
0025 last_point_.resize(adaption_info_.dimension);
0026 parametric_selector_ = parametric_selector(&last_point_,sample_other_variables_);
0027 exponent_selector_ = parametric_selector(&last_point_,sample_variables_);
0028 missing_accessor_ = parametric_missing_accessor(&last_parameter_bin_);
0029 parametric_sampler_ = parametric_sampling_selector<rnd_generator<Random> >
0030 (&last_point_,&last_parameter_bin_,sample_other_variables_,rnd_gen_);
0031 if (initialized_) return;
0032 splits_ = 0;
0033 for ( std::size_t k = 0; k < adaption_info_.dimension; ++k ) {
0034 if ( sample_other_variables_[k] )
0035 continue;
0036 parameter_splits_[k].push_back(adaption_info_.lower_left[k]);
0037 parameter_splits_[k].push_back(adaption_info_.upper_right[k]);
0038 }
0039 root_cell_ =
0040 binary_tree<cell>(cell(adaption_info_.lower_left,
0041 adaption_info_.upper_right,
0042 sample_other_variables_,
0043 adaption_info_));
0044 root_cell_.value().info().explore(rnd_gen_,adaption_info_,function_);
0045 root_cell_.value().integral(root_cell_.value().info().volume() * root_cell_.value().info().overestimate());
0046 last_exponent_integrand_.resize(1);
0047 check_events_ = adaption_info_.presampling_points;
0048 initialized_ = true;
0049 }
0050
0051 template<class Function, class Random>
0052 bool exponential_generator<Function,Random>::split () {
0053 if (adaption_info_.freeze_grid <= accepts_)
0054 return false;
0055 if (compensating_)
0056 return false;
0057 if (!(*last_cell_).info().bad(adaption_info_)) return false;
0058 bool dosplit = false;
0059 std::pair<std::size_t,double> sp =
0060 (*last_cell_).info().get_split(adaption_info_,dosplit);
0061 if (!dosplit) return false;
0062 if (!adaption_info_.adapt[sp.first]) return false;
0063 if (splits_ == parameter_hash_bits/2)
0064 return false;
0065 ++splits_;
0066 last_cell_.node().split((*last_cell_).split(sp,rnd_gen_,function_,adaption_info_,
0067 sample_other_variables_));
0068 if ( !sample_other_variables_[sp.first] ) {
0069 if ( std::find(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end(),sp.second)
0070 == parameter_splits_[sp.first].end() ) {
0071 parameter_splits_[sp.first].push_back(sp.second);
0072 std::sort(parameter_splits_[sp.first].begin(),parameter_splits_[sp.first].end());
0073 if ( sp.first == evolution_variable_ ) {
0074 last_exponent_integrand_.push_back(0.);
0075 }
0076 }
0077 }
0078 did_split_ = true;
0079 last_point_ = function_->parameter_point();
0080 root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
0081 exponents_.clear();
0082 get_exponent();
0083 return true;
0084 }
0085
0086 template<class Function, class Random>
0087 void exponential_generator<Function,Random>::get_exponent () {
0088 last_parameter_bin_.reset();
0089 root_cell_.subtree_hash (exponent_selector_,last_parameter_bin_);
0090 last_exponent_ = exponents_.find(last_parameter_bin_);
0091 if (last_exponent_ != exponents_.end())
0092 return;
0093 exponents_[last_parameter_bin_] = linear_interpolator();
0094 last_exponent_ = exponents_.find(last_parameter_bin_);
0095 double old_evo = last_point_[evolution_variable_];
0096 std::vector<double>::iterator exp_it = last_exponent_integrand_.begin();
0097 for (std::vector<double>::iterator esp = parameter_splits_[evolution_variable_].begin();
0098 esp < std::prev(parameter_splits_[evolution_variable_].end()); ++esp, ++exp_it) {
0099 last_point_[evolution_variable_] = (*esp + *std::next(esp))/2.;
0100 *exp_it = root_cell_.accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
0101 }
0102 exp_it = std::prev(last_exponent_integrand_.end());
0103 double total = 0.;
0104 for (std::vector<double>::iterator esp = std::prev(parameter_splits_[evolution_variable_].end());
0105 esp > parameter_splits_[evolution_variable_].begin(); --esp, --exp_it) {
0106 last_exponent_->second.set_interpolation(*esp,total);
0107 total += (*exp_it) * ((*esp) - (*std::prev(esp)));
0108 }
0109 last_exponent_->second.set_interpolation(parameter_splits_[evolution_variable_].front(),total);
0110 last_point_[evolution_variable_] = old_evo;
0111 }
0112
0113 template<class Function, class Random>
0114 std::set<std::vector<double> >
0115 exponential_generator<Function,Random>::parameter_points() {
0116 std::set<std::vector<double> > res;
0117 std::vector<double> pt(adaption_info_.dimension,0.);
0118 recursive_parameter_points(res,pt,0);
0119 return res;
0120 }
0121
0122 template<class Function, class Random>
0123 void exponential_generator<Function,Random>::
0124 recursive_parameter_points(std::set<std::vector<double> >& res,
0125 std::vector<double>& pt,
0126 size_t current) {
0127 if ( current == adaption_info_.dimension ) {
0128 res.insert(pt);
0129 return;
0130 }
0131 if ( sample_variables_[current] ) {
0132 recursive_parameter_points(res,pt,current+1);
0133 return;
0134 }
0135 for ( std::vector<double>::const_iterator sp =
0136 parameter_splits_[current].begin();
0137 sp != std::prev(parameter_splits_[current].end()); ++sp ) {
0138 pt[current] = (*sp + *std::next(sp))/2.;
0139 recursive_parameter_points(res,pt,current+1);
0140 }
0141 }
0142
0143 template<class Function, class Random>
0144 void exponential_generator<Function,Random>::compensate() {
0145 if (!did_split_ || !docompensate_) {
0146 assert(did_split_ || last_cell_ == root_cell_.begin());
0147 exponents_.clear();
0148 last_cell_->info().overestimate(last_value_,last_point_);
0149 last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
0150 last_point_ = function_->parameter_point();
0151 get_exponent();
0152 return;
0153 }
0154 std::vector<double> themaxpoint = last_point_;
0155 std::set<std::vector<double> > id_points
0156 = parameter_points();
0157 for ( std::set<std::vector<double> >::const_iterator id =
0158 id_points.begin(); id != id_points.end(); ++id ) {
0159 last_point_ = *id;
0160 get_exponent();
0161 }
0162 std::map<bit_container<parameter_hash_bits>,linear_interpolator >
0163 old_exponents = exponents_;
0164 double old_oe = last_cell_->info().overestimate();
0165 last_cell_->info().overestimate(last_value_,themaxpoint);
0166 last_cell_->integral(last_cell_->info().volume() * last_cell_->info().overestimate());
0167 exponents_.clear();
0168 for ( std::set<std::vector<double> >::const_iterator id =
0169 id_points.begin(); id != id_points.end(); ++id ) {
0170 last_point_ = *id;
0171 get_exponent();
0172 std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
0173 old_exp = old_exponents.find(last_parameter_bin_);
0174 std::map<bit_container<parameter_hash_bits>,linear_interpolator >::iterator
0175 new_exp = exponents_.find(last_parameter_bin_);
0176 assert(old_exp != old_exponents.end() && new_exp != exponents_.end());
0177 double old_norm = 1. - std::exp(-(old_exp->second)(adaption_info_.lower_left[evolution_variable_]));
0178 double new_norm = 1. - std::exp(-(new_exp->second)(adaption_info_.lower_left[evolution_variable_]));
0179 for (binary_tree<cell>::iterator it = root_cell_.begin();
0180 it != root_cell_.end(); ++it) {
0181 if ( !it->info().contains_parameter(last_point_,sample_variables_) )
0182 continue;
0183 double old_int = 0.;
0184 double new_int = 0.;
0185 for ( std::vector<double>::const_iterator sp = parameter_splits_[evolution_variable_].begin();
0186 sp != std::prev(parameter_splits_[evolution_variable_].end()); ++sp ) {
0187 if ( *sp >= it->info().lower_left()[evolution_variable_] &&
0188 *sp < it->info().upper_right()[evolution_variable_] ) {
0189 double xl = *sp;
0190 double xxl = *std::next(sp);
0191 double old_al =
0192 (old_exp->second.interpolation()[xxl] - old_exp->second.interpolation()[xl]) /
0193 (xxl-xl);
0194 double old_bl =
0195 (xxl * old_exp->second.interpolation()[xl] -
0196 xl * old_exp->second.interpolation()[xxl]) /
0197 (xxl-xl);
0198 double new_al =
0199 (new_exp->second.interpolation()[xxl] - new_exp->second.interpolation()[xl]) /
0200 (xxl-xl);
0201 double new_bl =
0202 (xxl * new_exp->second.interpolation()[xl] -
0203 xl * new_exp->second.interpolation()[xxl]) /
0204 (xxl-xl);
0205 if ( std::abs(old_al) > std::numeric_limits<double>::epsilon() ) {
0206 old_int += (exp(-(old_al*xl+old_bl)) - exp(-(old_al*xxl+old_bl)))/old_al;
0207 } else {
0208 old_int += (xxl-xl)*exp(-old_bl);
0209 }
0210 if ( std::abs(new_al) > std::numeric_limits<double>::epsilon() ) {
0211 new_int += (exp(-(new_al*xl+new_bl)) - exp(-(new_al*xxl+new_bl)))/new_al;
0212 } else {
0213 new_int += (xxl-xl)*exp(-new_bl);
0214 }
0215 }
0216 }
0217 double scaling;
0218 if (it != last_cell_) {
0219 if (old_int > std::numeric_limits<double>::epsilon() &&
0220 new_int > std::numeric_limits<double>::epsilon())
0221 scaling = ((old_norm * new_int) /
0222 (new_norm * old_int)) - 1.;
0223 else
0224 scaling = 0.;
0225 } else {
0226 if (old_int > std::numeric_limits<double>::epsilon() &&
0227 new_int > std::numeric_limits<double>::epsilon())
0228 scaling = ((last_value_ * old_norm * new_int) /
0229 (old_oe * new_norm * old_int)) - 1.;
0230 else
0231 scaling = 0.;
0232 }
0233 it->info().parametric_missing(last_parameter_bin_,
0234 it->info().parametric_missing(last_parameter_bin_) +
0235 static_cast<int>(round(scaling * it->info().attempted())));
0236 if (it->info().parametric_missing(last_parameter_bin_) != 0) {
0237 compensating_ = true;
0238 }
0239 }
0240 }
0241 last_point_ = function_->parameter_point();
0242 }
0243
0244 template<class Function, class Random>
0245 double exponential_generator<Function,Random>::generate(double enhance) {
0246 if ( enhance == 0.0 )
0247 return 0.;
0248 if (compensating_) {
0249 compensating_ = false;
0250 for (binary_tree<cell>::iterator it = root_cell_.begin();
0251 it != root_cell_.end(); ++it)
0252 if (it->info().parametric_compensating()) {
0253 compensating_ = true;
0254 break;
0255 }
0256 parametric_sampler_.compensate(compensating_);
0257 }
0258 last_point_ = function_->parameter_point();
0259 if (last_point_[evolution_variable_] < evolution_cutoff_) {
0260 return 0.;
0261 }
0262 unsigned long n_hit_miss = 0;
0263 unsigned long n_select = 0;
0264 double minus_log_r;
0265 root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
0266 get_exponent();
0267 while (true) {
0268 n_select = 0;
0269 minus_log_r = -std::log(rnd_gen_())/enhance/detuning_ +
0270 last_exponent_->second(last_point_[evolution_variable_]);
0271 if (!last_exponent_->second.invertible(minus_log_r)) {
0272 return 0.;
0273 }
0274 try {
0275 last_point_[evolution_variable_] = last_exponent_->second.unique_inverse(minus_log_r);
0276 } catch (constant_interpolation& c) {
0277 last_point_[evolution_variable_] = rnd_gen_(c.range.first,c.range.second);
0278 }
0279 assert(std::isfinite(last_point_[evolution_variable_]));
0280 if (last_point_[evolution_variable_] < evolution_cutoff_) {
0281 return 0.;
0282 }
0283 ++attempts_;
0284 if (compensating_) {
0285 root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
0286 }
0287 if (parameter_splits_[evolution_variable_].size() > 2)
0288 root_cell_.tree_accumulate(parametric_selector_,integral_accessor_,std::plus<double>());
0289 if (did_split_)
0290 while ((last_cell_ = root_cell_.select(parametric_sampler_)) == root_cell_.end()) {
0291 root_cell_.tree_accumulate(missing_accessor_,std::plus<int>());
0292 if(++n_select > adaption_info_.maxtry)
0293 throw selection_maxtry();
0294 }
0295 else
0296 last_cell_ = root_cell_.begin();
0297 last_cell_->info().select(rnd_gen_,last_point_,sample_other_variables_);
0298 last_value_ = function_->evaluate(last_point_);
0299 assert(last_value_ >= 0.);
0300 last_cell_->info().selected(last_point_,last_value_,adaption_info_);
0301 if (last_value_ > last_cell_->info().overestimate()) {
0302 if ( std::abs(last_value_)/last_cell_->info().overestimate() > 2. ) {
0303 last_value_ =
0304 last_cell_->info().overestimate()*
0305 (1.+exp(2.*(2.-std::abs(last_value_)/last_cell_->info().overestimate())));
0306 }
0307 compensate();
0308 throw exponential_regenerate();
0309 }
0310 if (last_cell_->info().attempted() % check_events_ == 0) {
0311 if (split()) {
0312 throw exponential_regenerate();
0313 }
0314 }
0315 if (last_value_/last_cell_->info().overestimate()/detuning_ > rnd_gen_()) {
0316 function_->accept(last_point_,enhance*last_value_,enhance*detuning_*last_cell_->info().overestimate());
0317 break;
0318 }
0319 if ( last_value_ != 0.0 ) {
0320 function_->veto(last_point_,enhance*last_value_,enhance*detuning_*last_cell_->info().overestimate());
0321 }
0322 if(++n_hit_miss > adaption_info_.maxtry)
0323 throw hit_and_miss_maxtry();
0324 }
0325 if (last_value_ == 0.)
0326 return 0.;
0327 ++accepts_;
0328 ++check_events_;
0329 last_cell_->info().accept();
0330 return 1.;
0331 }
0332
0333
0334 template<class Function, class Random>
0335 template<class OStream>
0336 void exponential_generator<Function,Random>::put (OStream& os) const {
0337 os << check_events_; ostream_traits<OStream>::separator(os);
0338 adaption_info_.put(os);
0339 root_cell_.put(os);
0340 os << did_split_; ostream_traits<OStream>::separator(os);
0341 os << initialized_; ostream_traits<OStream>::separator(os);
0342 os << evolution_variable_; ostream_traits<OStream>::separator(os);
0343 os << evolution_cutoff_; ostream_traits<OStream>::separator(os);
0344 os << sample_variables_; ostream_traits<OStream>::separator(os);
0345 os << sample_other_variables_; ostream_traits<OStream>::separator(os);
0346 os << parameter_splits_; ostream_traits<OStream>::separator(os);
0347
0348 os << last_point_; ostream_traits<OStream>::separator(os);
0349 os << last_value_; ostream_traits<OStream>::separator(os);
0350 last_parameter_bin_.put(os);
0351 os << exponents_.size(); ostream_traits<OStream>::separator(os);
0352 for ( std::map<bit_container<parameter_hash_bits>,linear_interpolator >::const_iterator
0353 ex = exponents_.begin(); ex != exponents_.end() ; ++ex ) {
0354 ex->first.put(os);
0355 ex->second.put(os);
0356 }
0357 os << last_exponent_integrand_; ostream_traits<OStream>::separator(os);
0358 os << compensating_; ostream_traits<OStream>::separator(os);
0359 os << attempts_; ostream_traits<OStream>::separator(os);
0360 os << accepts_; ostream_traits<OStream>::separator(os);
0361 os << splits_; ostream_traits<OStream>::separator(os);
0362 os << docompensate_; ostream_traits<OStream>::separator(os);
0363 }
0364
0365 template<class Function, class Random>
0366 template<class IStream>
0367 void exponential_generator<Function,Random>::get (IStream& is) {
0368 is >> check_events_;
0369 adaption_info_.get(is);
0370 root_cell_.get(is);
0371 is >> did_split_ >> initialized_ >> evolution_variable_
0372 >> evolution_cutoff_ >> sample_variables_ >> sample_other_variables_
0373 >> parameter_splits_;
0374
0375 is >> last_point_ >> last_value_;
0376 last_parameter_bin_.get(is);
0377 size_t dim; is >> dim;
0378 for ( size_t k = 0; k < dim ; ++k ) {
0379 bit_container<parameter_hash_bits> key;
0380 key.get(is);
0381 exponents_[key].get(is);
0382 }
0383 is >> last_exponent_integrand_;
0384 last_exponent_ = exponents_.find(last_parameter_bin_);
0385 is >> compensating_ >> attempts_ >> accepts_ >> splits_ >> docompensate_;
0386 }
0387
0388 }
