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
0014 template<class Random>
0015 void cell_info::select (Random& rnd_gen,
0016 std::vector<double>& p) {
0017 std::transform(lower_left_.begin(),lower_left_.end(),
0018 upper_right_.begin(),p.begin(),
0019 rnd_gen);
0020 ++attempted_;
0021 }
0022
0023
0024 template<class Random>
0025 void cell_info::select (Random& rnd_gen,
0026 std::vector<double>& p,
0027 const std::vector<bool>& sample) {
0028 conditional_transform(lower_left_.begin(),lower_left_.end(),
0029 upper_right_.begin(),sample.begin(),
0030 p.begin(),rnd_gen);
0031 ++attempted_;
0032 }
0033
0034
0035
0036 template<class Random, class Function, class SlaveStatistics>
0037 void cell_info::explore(Random& rnd_gen,
0038 const adaption_info& ainfo,
0039 Function* function, statistics* stats,
0040 SlaveStatistics& opt) {
0041 function->start_presampling();
0042 unsigned long n_sampled = 0;
0043 std::vector<double> ll = lower_left_;
0044 std::vector<double> ur = upper_right_;
0045 double val = 0.;
0046 std::vector<double> pos (ll.size());
0047 std::vector< std::pair<double,std::vector<double> > > vals;
0048 unsigned long ivalnonzero = 0;
0049 while (n_sampled < ainfo.presampling_points) {
0050 std::transform(ll.begin(),ll.end(),
0051 ur.begin(),pos.begin(),
0052 rnd_gen);
0053 val = function->evaluate(pos);
0054 vals.push_back( std::pair<double,std::vector<double> > (val,pos) );
0055 if ( val != 0 ) ivalnonzero++;
0056 ++n_sampled;
0057 }
0058 while ( ivalnonzero > 0 ) {
0059 double avg = 0;
0060 double err = 0;
0061 double maxval = 0;
0062 std::vector<double> maxpos (ll.size());
0063 unsigned long imax(0);
0064 for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
0065 val = std::abs(vals[ival].first);
0066 if ( val == 0 ) continue;
0067 avg += val;
0068 err += sqr(val);
0069 if ( val > maxval ) {
0070 maxval = val;
0071 maxpos = vals[ival].second;
0072 imax = ival;
0073 }
0074 }
0075 avg /= ivalnonzero;
0076 err /= ivalnonzero;
0077 err = sqrt(err-sqr(avg));
0078 if ( maxval <= avg+sqrt(ivalnonzero/2.)*err ) {
0079 overestimate_ = maxval;
0080 last_max_position_ = maxpos;
0081 break;
0082 }
0083 vals.erase(vals.begin()+imax);
0084 ivalnonzero--;
0085 }
0086
0087 for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
0088 val = vals[ival].first;
0089 stats->presampled(val);
0090 opt.select(val);
0091 selected(pos,std::abs(val),ainfo);
0092 }
0093
0094 function->stop_presampling();
0095 }
0096
0097
0098 template<class Random, class Function>
0099 void cell_info::explore (Random& rnd_gen,
0100 const adaption_info& ainfo, Function* function) {
0101 function->start_presampling();
0102 unsigned long n_sampled = 0;
0103 std::vector<double> ll = lower_left_;
0104 std::vector<double> ur = upper_right_;
0105 double val = 0.;
0106 std::vector<double> pos (ll.size());
0107 std::vector< std::pair<double,std::vector<double> > > vals;
0108 while (n_sampled < ainfo.presampling_points) {
0109 std::transform(ll.begin(),ll.end(),
0110 ur.begin(),pos.begin(),
0111 rnd_gen);
0112 val = function->evaluate(pos);
0113 if ( std::abs(val) > 0 )
0114 vals.push_back( std::pair<double,std::vector<double> > (std::abs(val),pos) );
0115 ++n_sampled;
0116 }
0117 while ( vals.size() > 0 ) {
0118 double avg = 0;
0119 double err = 0;
0120 double maxval = 0;
0121 std::vector<double> maxpos (ll.size());
0122 unsigned long imax(0);
0123 for ( unsigned long ival=0; ival < vals.size(); ival++ ) {
0124 double thisval = vals[ival].first;
0125 avg += thisval;
0126 err += sqr(thisval);
0127 if ( thisval > maxval ) {
0128 maxval = thisval;
0129 maxpos = vals[ival].second;
0130 imax = ival;
0131 }
0132 }
0133 avg /= vals.size();
0134 err /= vals.size();
0135 err = sqrt(err-sqr(avg));
0136 if ( maxval <= avg+sqrt(vals.size()/2.)*err ) {
0137 overestimate_ = maxval;
0138 last_max_position_ = maxpos;
0139 break;
0140 }
0141 vals.erase(vals.begin()+imax);
0142 }
0143
0144 function->stop_presampling();
0145 }
0146
0147
0148 template<class OStream>
0149 void cell_info::put (OStream& os) const {
0150 os << overestimate_;
0151 ostream_traits<OStream>::separator(os);
0152 os << volume_;
0153 ostream_traits<OStream>::separator(os);
0154 os << lower_left_.size();
0155 ostream_traits<OStream>::separator(os);
0156 for (std::size_t k = 0; k < lower_left_.size(); ++k) {
0157 os << lower_left_[k];
0158 ostream_traits<OStream>::separator(os);
0159 }
0160 for (std::size_t k = 0; k < upper_right_.size(); ++k) {
0161 os << upper_right_[k];
0162 ostream_traits<OStream>::separator(os);
0163 }
0164 for (std::size_t k = 0; k < mid_point_.size(); ++k) {
0165 os << mid_point_[k];
0166 ostream_traits<OStream>::separator(os);
0167 }
0168 for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
0169 os << last_max_position_[k];
0170 ostream_traits<OStream>::separator(os);
0171 }
0172 for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
0173 os << avg_weight_[k].first;
0174 ostream_traits<OStream>::separator(os);
0175 os << avg_weight_[k].second;
0176 ostream_traits<OStream>::separator(os);
0177 }
0178 os << attempted_;
0179 ostream_traits<OStream>::separator(os);
0180 os << accepted_;
0181 ostream_traits<OStream>::separator(os);
0182 os << parametric_missing_map_.size();
0183 ostream_traits<OStream>::separator(os);
0184 for ( std::map<bit_container<parameter_hash_bits>,int>::const_iterator p =
0185 parametric_missing_map_.begin(); p != parametric_missing_map_.end(); ++p ) {
0186 p->first.put(os);
0187 os << p->second;
0188 ostream_traits<OStream>::separator(os);
0189 }
0190 }
0191
0192
0193 template<class IStream>
0194 void cell_info::get (IStream& is) {
0195 std::size_t dim;
0196 is >> overestimate_ >> volume_ >> dim;
0197 lower_left_.resize(dim);
0198 for (std::size_t k = 0; k < lower_left_.size(); ++k) {
0199 is >> lower_left_[k];
0200 }
0201 upper_right_.resize(dim);
0202 for (std::size_t k = 0; k < upper_right_.size(); ++k) {
0203 is >> upper_right_[k];
0204 }
0205 mid_point_.resize(dim);
0206 for (std::size_t k = 0; k < mid_point_.size(); ++k) {
0207 is >> mid_point_[k];
0208 }
0209 last_max_position_.resize(dim);
0210 for (std::size_t k = 0; k < last_max_position_.size(); ++k) {
0211 is >> last_max_position_[k];
0212 }
0213 avg_weight_.resize(dim);
0214 for (std::size_t k = 0; k < avg_weight_.size(); ++k) {
0215 is >> avg_weight_[k].first >> avg_weight_[k].second;
0216 }
0217 is >> attempted_ >> accepted_ >> dim;
0218 for ( size_t k = 0; k < dim; ++k ) {
0219 bit_container<parameter_hash_bits> in;
0220 in.get(is);
0221 is >> parametric_missing_map_[in];
0222 }
0223 }
0224
0225 template<class Random, class Function>
0226 std::pair<cell,cell >
0227 cell::split (std::pair<std::size_t,double> split_d,
0228 Random& rnd_gen,
0229 Function* function,
0230 const adaption_info& ainfo,
0231 const std::vector<bool>& sampled) {
0232 assert(!missing_events() && !info().parametric_compensating());
0233 split_dimension_ = split_d.first;
0234 split_point_ = split_d.second;
0235 std::vector<double> lower_left1 = info().lower_left();
0236 std::vector<double> upper_right1 = info().upper_right();
0237 std::vector<double> lower_left2 = info().lower_left();
0238 std::vector<double> upper_right2 = info().upper_right();
0239 upper_right1[split_dimension_] = split_point_;
0240 lower_left2[split_dimension_] = split_point_;
0241 std::pair<cell,cell> children;
0242 if (sampled.empty())
0243 children = std::pair<cell,cell>(cell(lower_left1,upper_right1,ainfo),
0244 cell(lower_left2,upper_right2,ainfo));
0245 else
0246 children = std::pair<cell,cell> (cell(lower_left1,upper_right1,sampled,ainfo),
0247 cell(lower_left2,upper_right2,sampled,ainfo));
0248 if (info().last_max_position()[split_dimension_] <= split_point_) {
0249 children.first.info().overestimate(info().overestimate(),info().last_max_position());
0250 children.second.info().explore(rnd_gen,ainfo,function);
0251 } else {
0252 children.second.info().overestimate(info().overestimate(),info().last_max_position());
0253 children.first.info().explore(rnd_gen,ainfo,function);
0254 }
0255 cell_info_.reset(0);
0256 children.first.integral(children.first.info().volume() * children.first.info().overestimate());
0257 children.second.integral(children.second.info().volume() * children.second.info().overestimate());
0258 return children;
0259 }
0260
0261
0262 template<class OStream>
0263 void cell::put (OStream& os) const {
0264 os << split_dimension_;
0265 ostream_traits<OStream>::separator(os);
0266 os << split_point_;
0267 ostream_traits<OStream>::separator(os);
0268 os << integral_;
0269 ostream_traits<OStream>::separator(os);
0270 os << missing_events_;
0271 ostream_traits<OStream>::separator(os);
0272 if (cell_info_) {
0273 os << "has_cell_info";
0274 ostream_traits<OStream>::separator(os);
0275 cell_info_->put(os);
0276 } else {
0277 os << "has_no_cell_info";
0278 ostream_traits<OStream>::separator(os);
0279 }
0280 }
0281
0282
0283 template<class IStream>
0284 void cell::get (IStream& is) {
0285 std::string info_tag;
0286 is >> split_dimension_ >> split_point_
0287 >> integral_ >> missing_events_
0288 >> info_tag;
0289 if (info_tag == "has_cell_info") {
0290 cell_info_.reset(new cell_info());
0291 cell_info_->get(is);
0292 }
0293 }
0294
0295 inline cell_info::cell_info()
0296 : overestimate_(0.), volume_(0.),
0297 lower_left_(), upper_right_(), mid_point_(),
0298 last_max_position_(), avg_weight_(),
0299 attempted_(0), accepted_(0) {}
0300
0301 inline cell_info::cell_info(const std::vector<double>& ll,
0302 const std::vector<double>& ur,
0303 const adaption_info& ainfo)
0304 : overestimate_(0.), volume_(),
0305 lower_left_(ll), upper_right_(ur), mid_point_(),
0306 last_max_position_(),
0307 avg_weight_(std::vector<std::pair<double,double> >
0308 (ainfo.dimension,std::make_pair(0.,0.))),
0309 attempted_(0), accepted_(0) {
0310
0311 std::vector<double> delta;
0312 std::transform(ur.begin(),ur.end(),
0313 ll.begin(),std::back_inserter(delta),
0314 std::minus<double>());
0315
0316 volume_ =
0317 std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
0318
0319 std::transform(ur.begin(),ur.end(),
0320 ll.begin(),std::back_inserter(mid_point_),
0321 std::plus<double>());
0322
0323 for (std::size_t k = 0; k < ainfo.dimension; ++k)
0324 mid_point_[k] /= 2.;
0325
0326 }
0327
0328 inline cell_info::cell_info(const std::vector<double>& ll,
0329 const std::vector<double>& ur,
0330 const std::vector<bool>& sampled_variables,
0331 const adaption_info& ainfo)
0332 : overestimate_(0.), volume_(),
0333 lower_left_(ll), upper_right_(ur), mid_point_(),
0334 last_max_position_(),
0335 avg_weight_(std::vector<std::pair<double,double> >
0336 (ainfo.dimension,std::make_pair(0.,0.))),
0337 attempted_(0), accepted_(0) {
0338
0339 std::vector<double> delta;
0340 conditional_transform(ur.begin(),ur.end(),
0341 ll.begin(),sampled_variables.begin(),
0342 std::back_inserter(delta),
0343 std::minus<double>());
0344
0345 volume_ =
0346 std::accumulate(delta.begin(),delta.end(),1.,std::multiplies<double>());
0347
0348 std::transform(ur.begin(),ur.end(),
0349 ll.begin(),std::back_inserter(mid_point_),
0350 std::plus<double>());
0351
0352 for (std::size_t k = 0; k < ainfo.dimension; ++k)
0353 mid_point_[k] /= 2.;
0354
0355 }
0356
0357
0358 inline int cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id) const {
0359 std::map<bit_container<parameter_hash_bits>,int>::const_iterator mit
0360 = parametric_missing_map_.find(id);
0361 if (mit == parametric_missing_map_.end())
0362 return 0;
0363 return mit->second;
0364 }
0365
0366 inline void cell_info::parametric_missing(const bit_container<parameter_hash_bits>& id, int n) {
0367 if (n == 0) {
0368 std::map<bit_container<parameter_hash_bits>,int>::iterator mit
0369 = parametric_missing_map_.find(id);
0370 if (mit != parametric_missing_map_.end())
0371 parametric_missing_map_.erase(mit);
0372 return;
0373 }
0374 parametric_missing_map_[id] = n;
0375 }
0376
0377 inline void cell_info::increase_parametric_missing(const bit_container<parameter_hash_bits>& id) {
0378 std::map<bit_container<parameter_hash_bits>,int>::iterator mit
0379 = parametric_missing_map_.find(id);
0380 if (mit != parametric_missing_map_.end()) {
0381 mit->second += 1;
0382 if (mit->second == 0) parametric_missing_map_.erase(mit);
0383 } else parametric_missing_map_[id] = 1;
0384 }
0385
0386 inline void cell_info::decrease_parametric_missing(const bit_container<parameter_hash_bits>& id) {
0387 std::map<bit_container<parameter_hash_bits>,int>::iterator mit
0388 = parametric_missing_map_.find(id);
0389 if (mit != parametric_missing_map_.end()) {
0390 mit->second -= 1;
0391 if (mit->second == 0) parametric_missing_map_.erase(mit);
0392 } else assert(false);
0393 }
0394
0395 inline void cell_info::selected(const std::vector<double>& p,
0396 double weight,
0397 const adaption_info& ainfo) {
0398 for (std::size_t k = 0; k < p.size(); ++k) {
0399 if (ainfo.adapt[k]) {
0400 if (p[k] < mid_point_[k])
0401 avg_weight_[k].first += weight;
0402 else
0403 avg_weight_[k].second += weight;
0404 }
0405 }
0406 }
0407
0408 inline std::pair<std::size_t,double> cell_info::get_split (const adaption_info& ainfo,
0409 bool& worth) const {
0410 std::size_t split_d = 0;
0411 double gain = 0.;
0412 for (std::size_t k = 0; k < ainfo.dimension; ++k) {
0413 double xgain = 0.;
0414 double left = avg_weight_[k].first;
0415 double right = avg_weight_[k].second;
0416 if (left+right > 0.) {
0417 xgain = std::abs(left-right)/(left+right);
0418 }
0419 if (xgain > gain) {
0420 gain = xgain;
0421 split_d = k;
0422 }
0423 }
0424 worth = (gain >= ainfo.gain_threshold);
0425 return std::make_pair(split_d,mid_point_[split_d]);
0426 }
0427
0428 inline bool cell_info::contains_parameter (const std::vector<double>& point,
0429 const std::vector<bool>& sampled) const {
0430 std::vector<double>::const_iterator p = point.begin();
0431 std::vector<double>::const_iterator l = lower_left_.begin();
0432 std::vector<double>::const_iterator u = upper_right_.begin();
0433 std::vector<bool>::const_iterator f = sampled.begin();
0434 for (; p < point.end(); ++p, ++f, ++l, ++u)
0435 if (!(*f)) {
0436 if (((*l) > (*p)) ||
0437 ((*u) < (*p)))
0438 return false;
0439 }
0440 return true;
0441 }
0442
0443
0444 inline cell::cell()
0445 : split_dimension_(0), split_point_(0.),
0446 integral_(0.), missing_events_(0),
0447 cell_info_(nullptr) {}
0448
0449 inline cell::cell(const std::vector<double>& ll,
0450 const std::vector<double>& ur,
0451 const adaption_info& ainfo)
0452 : split_dimension_(0), split_point_(0.),
0453 integral_(0.), missing_events_(0),
0454 cell_info_(new cell_info(ll,ur,ainfo)) {}
0455
0456 inline cell::cell(const std::vector<double>& ll,
0457 const std::vector<double>& ur,
0458 const std::vector<bool>& sampled_variables,
0459 const adaption_info& ainfo)
0460 : split_dimension_(0), split_point_(0.),
0461 integral_(0.), missing_events_(0),
0462 cell_info_(new cell_info(ll,ur,sampled_variables,ainfo)) {}
0463
0464 inline cell::cell(const cell& x)
0465 : split_dimension_(x.split_dimension_),
0466 split_point_(x.split_point_),
0467 integral_(x.integral_),
0468 missing_events_(x.missing_events_),
0469 cell_info_(nullptr) {
0470 if (x.cell_info_)
0471 cell_info_.reset(new cell_info(*x.cell_info_));
0472 }
0473
0474 inline cell& cell::operator=(const cell& x) {
0475 if (this == &x)
0476 return *this;
0477 split_dimension_ = x.split_dimension_;
0478 split_point_ = x.split_point_;
0479 integral_ = x.integral_;
0480 missing_events_ = x.missing_events_;
0481 if (x.cell_info_)
0482 cell_info_.reset(new cell_info(*x.cell_info_));
0483 return *this;
0484 }
0485
0486
0487 }
