commit 0b974ba7cccf44d7fd62e6ef7e694a751c4a6bd2 Author: Guilhem Lavaux Date: Thu Nov 7 13:38:25 2024 +0200 initial import diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..e69de29 diff --git a/README.md b/README.md new file mode 100644 index 0000000..bd657a6 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# GLMath + diff --git a/ext_src/_cosmomath.pyx b/ext_src/_cosmomath.pyx new file mode 100644 index 0000000..fc4b3b6 --- /dev/null +++ b/ext_src/_cosmomath.pyx @@ -0,0 +1,45 @@ +#cython: language_level=3 +import numpy as np +cimport numpy as np + +np.import_array() +np.import_ufunc() + +cdef extern from "sys/types.h": + ctypedef np.int64_t int64_t + +cdef extern from "numpy/npy_common.h": + ctypedef npy_intp + +cdef extern from "special_math.hpp" namespace "CosmoTool": + T log_modified_bessel_first_kind[T](T v, T z) except + nogil + +cdef extern from "numpy_adaptors.hpp" namespace "CosmoTool": + void parallel_ufunc_dd_d[T,IT](char **args, IT* dimensions, IT* steps, void *func) + + +cdef np.PyUFuncGenericFunction loop_func[1] +cdef char input_output_types[3] +cdef void *elementwise_funcs[1] + +loop_func[0] = parallel_ufunc_dd_d[double,npy_intp] + +input_output_types[0] = np.NPY_DOUBLE +input_output_types[1] = np.NPY_DOUBLE +input_output_types[2] = np.NPY_DOUBLE + +elementwise_funcs[0] = log_modified_bessel_first_kind[double] + +log_modified_I = np.PyUFunc_FromFuncAndData( + loop_func, + elementwise_funcs, + input_output_types, + 1, # number of supported input types + 2, # number of input args + 1, # number of output args + 0, # `identity` element, never mind this + "log_modified_bessel_first_kind", # function name + "log_modified_bessel_first_kind(v: Float, z: Float) -> Float", # docstring + 0 # unused + ) + diff --git a/ext_src/numpy_adaptors.hpp b/ext_src/numpy_adaptors.hpp new file mode 100644 index 0000000..196206b --- /dev/null +++ b/ext_src/numpy_adaptors.hpp @@ -0,0 +1,30 @@ +#ifndef __COSMOTOOL_NUMPY_ADAPTOR_HPP +#define __COSMOTOOL_NUMPY_ADAPTOR_HPP + +namespace CosmoTool { + + template + void parallel_ufunc_dd_d(char **args, IT* dimensions, IT* steps, void *func) { + IT i; + IT n = dimensions[0]; + char *in = args[0], *in2 = args[1], *out = args[2]; + IT in_step = steps[0], in2_step = steps[1], out_step = steps[2]; + + double tmp; + typedef double (*F_t)(double,double); + + F_t f = (F_t)func; + +#pragma omp parallel for schedule(static) + for (i = 0; i < n; i++) { + T *out_t = (T *)(out + i * out_step); + T *in_t = (T *)(in + i * in_step); + T *in2_t = (T *)(in2 + i * in2_step); + *out_t = f(*in_t, *in2_t); + } + } + + +} + +#endif diff --git a/ext_src/special_math.hpp b/ext_src/special_math.hpp new file mode 100644 index 0000000..19e1414 --- /dev/null +++ b/ext_src/special_math.hpp @@ -0,0 +1,231 @@ +// Original code derived from Boost and is distributed here +// under the Boost license (licenses/boost-license.txt) +// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2007, 2017 John Maddock +// Secondary code copyright by its author and is distributed here +// under the BSD-3 license (LICENSE.md). Derived from +// stan/math/prim/fun/log_modified_bessel_first_kind.hpp + +#ifndef __COSMOTOOL_SPECIAL_MATH_HPP +#define __COSMOTOOL_SPECIAL_MATH_HPP + +#include "algo.hpp" +#include +#include +#include +#include + +// Taken and adapted from +// https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/log_modified_bessel_first_kind.hpp + +namespace CosmoTool { + +template T log1p_exp(T x) { + if (x > T(0)) { + return x + std::log1p(std::exp(-x)); + } + return std::log1p(std::exp(x)); +} + +template T multiply_log(T a, T b) { + if (a == 0 && b == 0) + return 0; + return a * std::log(b); +} + +template T inf() { return std::numeric_limits::infinity(); } + +template T log_sum_exp(T const a, T const b) { + if (a == -inf()) { + return b; + } + if (a == inf() && b == inf()) { + return inf(); + } + if (a > b) { + return a + log1p_exp(b - a); + } + return b + log1p_exp(a - b); +} + +/* Log of the modified Bessel function of the first kind, + * which is better known as the Bessel I function. See + * modified_bessel_first_kind.hpp for the function definition. + * The derivatives are known to be incorrect for v = 0 because a + * simple constant 0 is returned. + * + * @tparam T common type for calculation + * @param v Order, can be a non-integer but must be at least -1 + * @param z Real non-negative number + * @throws std::domain_error if either v or z is NaN, z is + * negative, or v is less than -1 + * @return log of Bessel I function + */ +template T log_modified_bessel_first_kind(T const v, T const z) { + using boost::math::tools::evaluate_polynomial; + using std::log; + using std::pow; + using std::sqrt; + static const double LOG_TWO = std::log(2.0); + static const double EPSILON = std::numeric_limits::epsilon(); + static const double TWO_PI = 2.0 * boost::math::constants::pi(); + + if (z == 0) { + if (v == 0) { + return 0.0; + } + if (v > 0) { + return -std::numeric_limits::infinity(); + } + return std::numeric_limits::infinity(); + } + if (std::isinf(z)) { + return z; + } + if (v == 0) { + // Modified Bessel function of the first kind of order zero + // we use the approximating forms derived in: + // "Rational Approximations for the Modified Bessel Function of the + // First Kind -- I0(x) for Computations with Double Precision" + // by Pavel Holoborodko, see + // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision + // The actual coefficients used are [Boost's] own, and extend + // Pavel's work to precisions other than double. + + if (z < 7.75) { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form : 3.042e-18 + // Max Error found at double precision = Poly : 5.106609e-16 + // Cheb : 5.239199e-16 + static const double P[] = { + 1.00000000000000000e+00, 2.49999999999999909e-01, + 2.77777777777782257e-02, 1.73611111111023792e-03, + 6.94444444453352521e-05, 1.92901234513219920e-06, + 3.93675991102510739e-08, 6.15118672704439289e-10, + 7.59407002058973446e-12, 7.59389793369836367e-14, + 6.27767773636292611e-16, 4.34709704153272287e-18, + 2.63417742690109154e-20, 1.13943037744822825e-22, + 9.07926920085624812e-25}; + return log1p_exp(multiply_log(2.0, z) - log(4.0) + + log(evaluate_polynomial(P, 0.25 * square(z)))); + } + if (z < 500) { + // Max error in interpolated form : 1.685e-16 + // Max Error found at double precision = Poly : 2.575063e-16 + // Cheb : 2.247615e+00 + static const double P[] = { + 3.98942280401425088e-01, 4.98677850604961985e-02, + 2.80506233928312623e-02, 2.92211225166047873e-02, + 4.44207299493659561e-02, 1.30970574605856719e-01, + -3.35052280231727022e+00, 2.33025711583514727e+02, + -1.13366350697172355e+04, 4.24057674317867331e+05, + -1.23157028595698731e+07, 2.80231938155267516e+08, + -5.01883999713777929e+09, 7.08029243015109113e+10, + -7.84261082124811106e+11, 6.76825737854096565e+12, + -4.49034849696138065e+13, 2.24155239966958995e+14, + -8.13426467865659318e+14, 2.02391097391687777e+15, + -3.08675715295370878e+15, 2.17587543863819074e+15}; + return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z); + } + // Max error in interpolated form : 2.437e-18 + // Max Error found at double precision = Poly : 1.216719e-16 + static const double P[] = {3.98942280401432905e-01, 4.98677850491434560e-02, + 2.80506308916506102e-02, 2.92179096853915176e-02, + 4.53371208762579442e-02}; + return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z); + } + if (v == 1) { // WARNING: will not autodiff for v = 1 correctly + // modified from Boost's bessel_i1_imp in the double precision case + // see credits above in the v == 0 case + if (z < 7.75) { + // Bessel I0 over[10 ^ -16, 7.75] + // Max error in interpolated form: 5.639e-17 + // Max Error found at double precision = Poly: 1.795559e-16 + + static const double P[] = { + 8.333333333333333803e-02, 6.944444444444341983e-03, + 3.472222222225921045e-04, 1.157407407354987232e-05, + 2.755731926254790268e-07, 4.920949692800671435e-09, + 6.834657311305621830e-11, 7.593969849687574339e-13, + 6.904822652741917551e-15, 5.220157095351373194e-17, + 3.410720494727771276e-19, 1.625212890947171108e-21, + 1.332898928162290861e-23}; + T a = square(z) * 0.25; + T Q[3] = {1, 0.5, evaluate_polynomial(P, a)}; + return log(z) + log(evaluate_polynomial(Q, a)) - LOG_TWO; + } + if (z < 500) { + // Max error in interpolated form: 1.796e-16 + // Max Error found at double precision = Poly: 2.898731e-16 + + static const double P[] = { + 3.989422804014406054e-01, -1.496033551613111533e-01, + -4.675104253598537322e-02, -4.090895951581637791e-02, + -5.719036414430205390e-02, -1.528189554374492735e-01, + 3.458284470977172076e+00, -2.426181371595021021e+02, + 1.178785865993440669e+04, -4.404655582443487334e+05, + 1.277677779341446497e+07, -2.903390398236656519e+08, + 5.192386898222206474e+09, -7.313784438967834057e+10, + 8.087824484994859552e+11, -6.967602516005787001e+12, + 4.614040809616582764e+13, -2.298849639457172489e+14, + 8.325554073334618015e+14, -2.067285045778906105e+15, + 3.146401654361325073e+15, -2.213318202179221945e+15}; + return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z); + } + // Max error in interpolated form: 1.320e-19 + // Max Error found at double precision = Poly: 7.065357e-17 + static const double P[] = { + 3.989422804014314820e-01, -1.496033551467584157e-01, + -4.675105322571775911e-02, -4.090421597376992892e-02, + -5.843630344778927582e-02}; + return z + log(evaluate_polynomial(P, 1 / z)) - multiply_log(0.5, z); + } + if (z > 100) { + // Boost does something like this in asymptotic_bessel_i_large_x + T lim = pow((square(v) + 2.5) / (2 * z), 3) / 24; + if (lim < (EPSILON * 10)) { + T s = 1; + T mu = 4 * square(v); + T ex = 8 * z; + T num = mu - 1; + T denom = ex; + s -= num / denom; + num *= mu - 9; + denom *= ex * 2; + s += num / denom; + num *= mu - 25; + denom *= ex * 3; + s -= num / denom; + s = z - log(sqrt(z * TWO_PI)) + log(s); + return s; + } + } + + T log_half_z = log(0.5 * z); + T lgam = (v > -1) ? lgamma(v + 1.0) : 0; + T lcons = (2.0 + v) * log_half_z; + T out; + if (v > -1) { + out = log_sum_exp(v * log_half_z - lgam, lcons - lgamma(v + 2)); + lgam += log1p(v); + } else { + out = lcons; + } + + double m = 2; + double lfac = 0; + T old_out; + do { + lfac += log(m); + lgam += log(v + m); + lcons += 2 * log_half_z; + old_out = out; + out = log_sum_exp(out, lcons - lfac - lgam); // underflows eventually + m++; + } while (out > old_out || out < old_out); + return out; +} + +} // namespace CosmoTool + +#endif diff --git a/new_release.sh b/new_release.sh new file mode 100644 index 0000000..28636c5 --- /dev/null +++ b/new_release.sh @@ -0,0 +1 @@ +semantic-release version --tag --changelog --vcs-release diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..e7766b7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,87 @@ +[project] +version = "1.0.0" +requires-python = ">= 3.8" +name = "pysphereproj" +dependencies = [ + "numpy" +] +authors = [ + { name = "Guilhem Lavaux", email = "guilhem.lavaux@iap.fr" } +] +readme = "README.md" +classifiers = [ + "Intended Audience :: Developers", + "License :: OSI Approved :: MIT License", + "Programming Language :: Python :: 3" +] + +[project.urls] +Homepage = "https://git.aquila-consortium.org/guilhem_lavaux/pysphereproj" +Issues = "https://git.aquila-consortium.org/guilhem_lavaux/pysphereproj/issues" + +[build-system] +requires = ["setuptools >= 61.0", "numpy", "wheel", "Cython", "python-semantic-release"] +build-backend = "setuptools.build_meta" + +[tool.semantic_release] +assets = [] +version_toml = [ + "pyproject.toml:project.version", +] +version_variables = [ + "src/sphereproj/__init__.py:__version__" +] +build_command_env = [] +commit_message = "{version}\n\nAutomatically generated by python-semantic-release" +commit_parser = "angular" +logging_use_named_masks = false +major_on_zero = true +allow_zero_version = true +tag_format = "v{version}" + +[tool.semantic_release.branches.main] +match = "(main|master)" +prerelease_token = "rc" +prerelease = false + +[tool.semantic_release.default_templates] +template_dir = "templates" +changelog_file = "CHANGELOG.md" +exclude_commit_patterns = [] + +[tool.semantic_release.changelog.environment] +block_start_string = "{%" +block_end_string = "%}" +variable_start_string = "{{" +variable_end_string = "}}" +comment_start_string = "{#" +comment_end_string = "#}" +trim_blocks = false +lstrip_blocks = false +newline_sequence = "\n" +keep_trailing_newline = false +extensions = [] +autoescape = true + +[tool.semantic_release.commit_author] +env = "GIT_COMMIT_AUTHOR" +default = "semantic-release " + +[tool.semantic_release.commit_parser_options] +allowed_tags = ["build", "chore", "ci", "docs", "feat", "fix", "perf", "style", "refactor", "test"] +minor_tags = ["feat"] +patch_tags = ["fix", "perf"] +default_bump_level = 0 + +[tool.semantic_release.remote] +name = "origin" +type = "gitea" +ignore_token_for_push = false +insecure = false +domain = "git.aquila-consortium.org" + +[tool.semantic_release.publish] +dist_glob_patterns = ["dist/*"] +upload_to_vcs_release = true + + diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..0c2faed --- /dev/null +++ b/setup.py @@ -0,0 +1,17 @@ +from setuptools import setup, Extension +from Cython.Build import cythonize +import numpy + +extensions = [ + Extension( + "glmath._cosmomath", + sources=[ + "ext_src/_cosmomath.pyx" + ], + include_dirs=["ext_src", numpy.get_include()], + language="c++" + ) +] + + +setup(ext_modules=cythonize(extensions)) diff --git a/src/glmath/__init__.py b/src/glmath/__init__.py new file mode 100644 index 0000000..e69de29