Solved bugs.

adolfomunoz · May 13, 2024 · d63372c · d63372c
1 parent f3cc24f
commit d63372c
Show file tree

Hide file tree

Showing 4 changed files with 33 additions and 9 deletions.
diff --git a/src/nested/error-heuristic.h b/src/nested/error-heuristic.h
@@ -21,21 +21,24 @@ class error_heuristic_default {
 template<typename ErrorMetric>
 class error_heuristic_size {
  ErrorMetric error_metric;
- double size_weight;
+ double size_weight; double min_size;
 public:
- error_heuristic_size(const ErrorMetric& em, double sw = 1.e-5) : 
- error_metric(em), size_weight(sw) {}
+ error_heuristic_size(const ErrorMetric& em, double sw = 1.e-5, double ms = 1.e-37) : 
+ error_metric(em), size_weight(sw), min_size(ms) {}
 
  template<typename R>
  auto operator()(const R& region) const {
  auto max_err = region.error(0,error_metric);
+ if ((region.range().max(0)-region.range().min(0))<min_size) max_err=0;
  max_err += size_weight*std::abs(region.range().max(0) - region.range().min(0)); 
  std::size_t max_dim = 0; 
  auto err = max_err;
  for (std::size_t d = 1; d<R::dimensions; ++d) {
  err = region.error(d,error_metric) + 
  size_weight*std::abs(region.range().max(d) - region.range().min(d));
- if (err>max_err) {
+
+ if ((region.range().max(d)-region.range().min(d))<min_size) err=0;
+ if (err>=max_err) {
  max_err = err; max_dim = d;
  }
  }

diff --git a/src/nested/error-metric.h b/src/nested/error-metric.h
@@ -27,9 +27,16 @@ class error_metric_absolute {
 };
 
 class error_metric_relative {
+ double min_val;
 public:
- float operator()(float a, float b) const { return std::abs(b-a)/(std::max(std::abs(a),std::abs(b))); }
- double operator()(double a, double b) const { return std::abs(b-a)/(std::max(std::abs(a),std::abs(b))); }
+ error_metric_relative(double mv = 1.e-10) : 
+ min_val(mv) {}
+ float operator()(float a, float b) const { 
+ return std::abs(b-a)/(std::max(std::max(std::abs(a),std::abs(b)),float(min_val))); 
+ }
+ double operator()(double a, double b) const { 
+ return std::abs(b-a)/(std::max(std::max(std::abs(a),std::abs(b)),min_val)); 
+ }
  template<typename V>
  auto operator()(const V& v1, const V& v2, typename std::enable_if_t<std::is_arithmetic_v<typename V::value_type>, int> = 0) const {
  auto i1 = v1.begin(); auto i2 = v2.begin();

diff --git a/src/newton-cotes/rules.h b/src/newton-cotes/rules.h
@@ -85,9 +85,14 @@ struct Simpson {
  template<typename Float, typename T> 
  constexpr T at(Float t, const std::array<T,samples>& p) const {
  auto c = coefficients(p);
- return (c[2]*t + c[1])*t + c[0]; 
+ return (c[2]*t + c[1])*t + c[0];
  }
 
+ static constexpr bool isnan(double d) { return std::isnan(d); }
+ static constexpr bool isnan(float d) { return std::isnan(d); }
+ template<typename V>
+ static constexpr bool isnan(const V& v) { return isnan(v[0]); }
+
  template<typename Float, typename T> 
  constexpr T subrange(Float a, Float b, const std::array<T,samples>& p) const {
  auto c = coefficients(p);

diff --git a/src/range.h b/src/range.h
@@ -43,8 +43,10 @@ class Range : public std::array<std::array<T,DIM>,2> {
  std::array<T,DIMSUB> pos_in_range(const std::array<T,DIMSUB>& pos) const {
  static_assert(DIMSUB<=DIM,"pos_in_range with too big dimensionsal index");
  std::array<T,DIMSUB> prange;
- for (std::size_t i = 0; i<DIMSUB; ++i) 
- prange[i] = (pos[i]-min(i))/(max(i) - min(i));
+ for (std::size_t i = 0; i<DIMSUB; ++i) {
+ if (min(i)>=max(i)) prange[i]=min(i);
+ else prange[i] = (pos[i]-min(i))/(max(i) - min(i));
+ }
  return prange;
  }
 
@@ -63,6 +65,13 @@ class Range : public std::array<std::array<T,DIM>,2> {
  return Range<T,DIM>(new_a, new_b);
  }
 
+ bool empty() const {
+ bool e = false;
+ for (std::size_t i = 0; (i<DIM) && (!e); ++i)
+ e = (min(i)>=max(i));
+ return e; 
+ }
+
 /*
  Range<T,DIM> intersection(const Range<T,DIM>& that) const {
  std::array<T,DIM> a, b;