Clang & Python 3.8 fix (#3)

* Replace std::sqrt in constexpr by a Babylon algorithm implementation to achieve compatibility with clang++.

* Relax typing requirements on Python version. Also relax type requirements where not needed.

* Change log.

README.md

@@ -92,6 +92,13 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to
 [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+### [1.0.1] - 2022-11-28
+#### Changed
+- Fix a C++ standard incompatibility that is compatible with g++ but
+  not with clang++.
+- Relax some typing requirements and make typing information compatible with
+  Python 3.8.
+
 ### [1.0.0] - 2022-11-27
 #### Added
 - First release version

docs/conf.py

@@ -9,7 +9,7 @@
 project = 'REHEATFUNQ'
 copyright = '2022, Deutsches GeoForschungsZentrum Potsdam & Malte J. Ziebarth'
 author = 'Malte J. Ziebarth'
-release = '1.0.0'
+release = '1.0.1'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

external/pdtoolbox/include/constexpr.hpp

@@ -0,0 +1,54 @@
+/*
+ * Constexpr implementations of mathematical functions.
+ *
+ * Author: Malte J. Ziebarth (ziebarth@gfz-potsdam.de)
+ *
+ * Copyright (C) 2022 Malte J. Ziebarth
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <cmath>
+
+namespace pdtoolbox {
+
+/*
+ * Square root using the Babylon algorithm.
+ * Adds some 'static asserts' by raising recursion depth errors.
+ */
+constexpr double const_sqrt_iter(const double S, const double d)
+{
+	/* A very hacky static assert, rasing a recursion depth error: */
+	if (S < 0)
+		return const_sqrt_iter(S,d);
+
+	double Snext = (S + d/S) / 2;
+	if (Snext == S)
+		return S;
+	return const_sqrt_iter(Snext, d);
+}
+
+constexpr double cnst_sqrt(const double d)
+{
+	/* Babylon: */
+	double S = const_sqrt_iter(d,d);
+
+	/* "static_assert" */
+	if (std::abs(S*S - d) > 1e-14*d)
+		return cnst_sqrt(d);
+
+	return S;
+}
+
+}

external/pdtoolbox/src/gamma_conjugate_prior.cpp

@@ -22,6 +22,7 @@
  * Miller (1980):
  */
 
+#include <constexpr.hpp>
 #include <gamma_conjugate_prior.hpp>
 #include <cmath>
 #include <boost/math/special_functions/digamma.hpp>
@@ -31,7 +32,8 @@
 #include <boost/math/quadrature/tanh_sinh.hpp>
 #include <boost/math/quadrature/exp_sinh.hpp>
 
-using pdtoolbox::GammaConjugatePriorBase;
+using pdtoolbox::GammaConjugatePriorBase,
+      pdtoolbox::cnst_sqrt;
 
 
 using std::abs, std::max, std::min, std::log, std::exp, std::sqrt, std::log1p,
@@ -289,7 +291,6 @@ static bool condition_warn(double S, double err, double L1, int where,
 
 
 
-
 /*
  * Compute the natural logarithm of the integration constant:
  */
@@ -367,7 +368,7 @@ static double ln_Phi_backend(double lp, double ls, double n, double v,
 	 * integration techniques: */
 	double err, L1;
 	double res = 0.0 ;
-	constexpr double term = std::sqrt(std::numeric_limits<double>::epsilon());
+	constexpr double term = cnst_sqrt(std::numeric_limits<double>::epsilon());
 
 	auto integrand1 = [&](double a, double distance_to_next_bound) -> double {
 		double lF = ln_F(a, P.lp, P.ls, P.n, P.v);
@@ -521,7 +522,7 @@ double GammaConjugatePriorBase::kullback_leibler(
 	};
 
 	exp_sinh<double> es;
-	constexpr double term = std::sqrt(std::numeric_limits<double>::epsilon());
+	constexpr double term = cnst_sqrt(std::numeric_limits<double>::epsilon());
 	double err, L1;
 	const double I = es.integrate(integrand, term, &err, &L1);
 

external/pdtoolbox/src/ziebarth2022a.cpp

@@ -41,6 +41,7 @@
 #include <algorithm>
 #include <iostream>
 
+#include <constexpr.hpp>
 #include <ziebarth2022a.hpp>
 #include <quantileinverter.hpp>
 
@@ -57,6 +58,7 @@ using boost::math::tools::eps_tolerance;
 using pdtoolbox::QuantileInverter;
 using pdtoolbox::OR;
 using pdtoolbox::AND;
+using pdtoolbox::cnst_sqrt;
 
 typedef gauss_kronrod<double, 15> GK;
 
@@ -244,7 +246,7 @@ static double outer_integrand(double z, const double lp,
                               const double w, const double log_scale)
 {
 	/* exp_sinh tolerance: */
-	constexpr double tol = std::sqrt(std::numeric_limits<double>::epsilon());
+	constexpr double tol = cnst_sqrt(std::numeric_limits<double>::epsilon());
 
 	// Set the inner parameters:
 	double l1p_kiz_sum = 0.0;
@@ -812,7 +814,7 @@ double a_integral_large_z(const double ym, const double h0, const double h1,
                  const double l1p_w, const double w)
 {
 	/* exp_sinh tolerance: */
-	constexpr double tol = std::sqrt(std::numeric_limits<double>::epsilon());
+	constexpr double tol = cnst_sqrt(std::numeric_limits<double>::epsilon());
 
 	/* Set the integrand's non-varying parameters: */
 	const double ly = std::log(ym);

reheatfunq/coverings/rgrdc.py

@@ -21,7 +21,7 @@
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 import numpy as np
-from typing import Optional
+from typing import Optional, List, Sequence, Iterable
 from numpy.typing import ArrayLike
 from pyproj import Proj, Geod
 from scipy.spatial import KDTree
@@ -30,11 +30,11 @@
 
 
 def random_global_R_disk_coverings(R: float, min_points: int, hf: ArrayLike,
-                                   buffered_poly_xy: list[ArrayLike],
+                                   buffered_poly_xy: Sequence[ArrayLike],
                                    proj_str: str,  N: int = 10000,
                                    MAX_DRAW: int = 100000, dmin: float = 0.0,
                                    seed: int = 982981,
-                                   used_points: Optional[list[int]] = None,
+                                   used_points: Optional[Iterable[int]] = None,
                                    a: float = 6378137.0):
     """
     Uses rejection sampling to draw a number of exclusive regional

reheatfunq/regional/prior.py

@@ -327,7 +327,7 @@ def maximum_likelihood_estimate(a: ArrayLike, b: ArrayLike,
 
 
     @staticmethod
-    def minimum_surprise_estimate(hf_samples: list[ArrayLike],
+    def minimum_surprise_estimate(hf_samples: Iterable[ArrayLike],
                                   pmin: float = 1.0, pmax: float = 1e5,
                                   smin: float = 0.0, smax: float = 1e3,
                                   vmin: float = 0.02, vmax: float = 1.0,
@@ -344,7 +344,7 @@ def minimum_surprise_estimate(hf_samples: list[ArrayLike],
 
         Parameters
         ----------
-        hf_samples : list[array_like]
+        hf_samples : Iterable[array_like]
             A set of heat flow data sets.
         pmin : float, optional
             Minimum value for the GCP :math:`p` parameter.
@@ -435,20 +435,21 @@ def cost(x):
         return GammaConjugatePrior(None, s, n, v, lp, amin)
 
 
-    def visualize(self, ax, distributions: Optional[list[ArrayLike]] = None,
+    def visualize(self, ax, distributions: Optional[Iterable[ArrayLike]] = None,
                   cax: Optional[object] = None, log_axes: bool = True,
                   cmap='inferno', color_scale: Literal['log','lin'] = 'log',
                   plot_mean: bool = True, q_mean: float = 68.3,
-                  q_plot: list[tuple[float,float,float,str] | float] = [],
-                  qstd_plot: list[tuple[float,float,float,str] | float] = []):
+                  q_plot: Iterable[Tuple[float,float,float,str] | float] = [],
+                  qstd_plot: Iterable[Tuple[float,float,float,str] | float] = []
+        ):
         """
         Visualize this GammaConjugatePrior instance on an axis.
 
         Parameters
         ----------
         ax : matplotlib.axes.Axes
            The :class:`matplotlib.axes.Axes` to visualize on.
-        distributions : list[array_like], optional
+        distributions : Iterable[array_like], optional
            A set of aggregate heat flow distributions, each given
            as a one-dimensional :class:`numpy.ndarray` of heat flow
            values in :math:`\\mathrm{mW}/\\mathrm{m}^2`. Each
@@ -472,16 +473,16 @@ def visualize(self, ax, distributions: Optional[list[ArrayLike]] = None,
         q_mean : float, optional
            The global mean heat flow in :math:`\\mathrm{mW}/\\mathrm{m}^2`.
            The default value is 68.3 from Lucazeau (2019).
-        q_plot : list[tuple[float,float,float,str] | float], optional
-           A list of additional average heat flow values to display.
+        q_plot : Iterable[Tuple[float,float,float,str] | float], optional
+           A set of additional average heat flow values to display.
            For each *q* a line through the :math:`(\\alpha,\\beta)` parameter
            space, enumerating parameter combinations whose distributions
            average to the given *q*. Each entry in `q_plot` needs to be
            either a float *q* or a tuple (*q*,*amin*,*amax*,*c*), where *amin*
            and *amax* denote the :math:`\\alpha`-interval within which the
            line should be plotted, and *c* is the color.
-        qstd_plot : list[tuple[float,float,float,str] | float], optional
-           A list of additional heat flow standard deviations to display.
+        qstd_plot : Iterable[Tuple[float,float,float,str] | float], optional
+           A set of additional heat flow standard deviations to display.
            For each *qstd* a line through the :math:`(\\alpha,\\beta)` parameter
            space, enumerating parameter combinations whose distributions
            are quantified by a standard deviation *qstd*. Each entry in

setup.py

@@ -30,7 +30,7 @@
 
 
 setup(name='REHEATFUNQ',
-      version='1.0.0',
+      version='1.0.1',
       author='Malte J. Ziebarth',
       description='',
       packages = ['reheatfunq','reheatfunq.regional','reheatfunq.anomaly',