Merge branch 'tickets/DM-45473'

lsst-dm · Jul 31, 2024 · eb408df · eb408df
2 parents 0e535e6 + ebf3225
commit eb408df
Show file tree

Hide file tree

Showing 12 changed files with 145 additions and 73 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,7 +3,17 @@
 All notable changes to this project will be documented in this file.
 This project adheres to [Semantic Versioning](https://semver.org/).
 
-### [0.1.3] 2024-07-23
+### [0.1.4] 2024-07-30
+
+* Fixed: Put correct version in pkg-config file for standalone C++ builds
+* Added: -ffp-contract=off to default C++ compiler options
+* Fixed: Allow tiny (1e-25) deviations from equality in some Python tests
+* Fixed: Python formatting
+* Added: Black and ruff configurations to pyproject.toml
+* Fixed: Removed incorrect subdir for standalone Python build
+* See [DM-45473](https://rubinobs.atlassian.net/browse/DM-45473) for details.
+
+### [0.1.3] 2024-07-29
 
 * Fixed: test_channel is no longer sensitive to test order.
 * Fixed: test_model no longer calls compute_hessian with overly verbose print=True.
@@ -33,6 +43,8 @@ This project adheres to [Semantic Versioning](https://semver.org/).
 * Changed: Initial release, forked to https://github.com/lsst/gauss2d_fit.
 * See [DM-43906](https://rubinobs.atlassian.net/browse/DM-43906) for details.
 
+[0.1.4]: https://github.com/lsst-dm/gauss2d_fit/compare/0.1.3...0.1.4
+
 [0.1.3]: https://github.com/lsst-dm/gauss2d_fit/compare/0.1.2...0.1.3
 
 [0.1.2]: https://github.com/lsst-dm/gauss2d_fit/compare/0.1.1...0.1.2

diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-0.1.3
+0.1.4
diff --git a/include/lsst/gauss2d/fit/model.h b/include/lsst/gauss2d/fit/model.h
@@ -515,16 +515,14 @@ class Model : public ParametricModel {
         }
         auto map_extra_mut = expired ? std::make_shared<Indices>(n_gaussians_conv, 2, nullptr, coordsys)
                                      : map_extra_weak.lock();
-        auto map_grad_mut
-                = expired ? std::make_shared<Indices>(n_gaussians_conv, lsst::gauss2d::N_PARAMS_GAUSS2D,
-                                                      nullptr, coordsys)
-                          : map_grad_weak.lock();
+        auto map_grad_mut = expired ? std::make_shared<Indices>(
+                                    n_gaussians_conv, lsst::gauss2d::N_PARAMS_GAUSS2D, nullptr, coordsys)
+                                    : map_grad_weak.lock();
         auto factors_extra_mut = expired ? std::make_shared<Image>(n_gaussians_conv, 3, nullptr, coordsys)
                                          : factors_extra_weak.lock();
-        auto factors_grad_mut
-                = expired ? std::make_shared<Image>(n_gaussians_conv, lsst::gauss2d::N_PARAMS_GAUSS2D,
-                                                    nullptr, coordsys)
-                          : factors_grad_weak.lock();
+        auto factors_grad_mut = expired ? std::make_shared<Image>(
+                                        n_gaussians_conv, lsst::gauss2d::N_PARAMS_GAUSS2D, nullptr, coordsys)
+                                        : factors_grad_weak.lock();
         ;
 
         ExtraParamMap map_extra = {};
@@ -1405,12 +1403,17 @@ class Model : public ParametricModel {
      * @param findiff_add The minimum value of the finite difference increment.
      * @param rtol The allowed relative tolerance in the Jacobian as compared to the finite difference.
      * @param atol The allowed absolute tolerance in the Jacobian as compared to the finite difference.
+     * @param max_ll_diff The maximum allowed difference between equivalent log-likelihood evaluations.
+     *      Must be >= 0 and should only be a few orders of magnitude larger than machine epsilon.
      * @return A vector of error messages, one for each Parameter that failed verification.
      *
      * @note Verification is done using isclose(), which is modelled after Python's numpy.isclose.
      */
     std::vector<std::string> verify_jacobian(double findiff_frac = 1e-5, double findiff_add = 1e-5,
-                                             double rtol = 1e-3, double atol = 1e-3) {
+                                             double rtol = 1e-3, double atol = 1e-3, double max_ll_diff = 0) {
+        if (!(max_ll_diff >= 0)) {
+            throw std::invalid_argument("max_ll_diff=" + to_string_float(max_ll_diff) + " !> 0");
+        }
         if (_mode != EvaluatorMode::jacobian) {
             this->setup_evaluators(EvaluatorMode::jacobian);
         }
@@ -1466,9 +1469,11 @@ class Model : public ParametricModel {
 
             auto result = _make_evaluator(idx, EvaluatorMode::loglike_image);
             double loglike_img = result.first->loglike_pixel();
-            if (loglike_jac != loglike_img) {
-                errors.push_back("loglike_jac=" + to_string_float(loglike_jac) + "loglike_img="
-                                 + to_string_float(loglike_img) + "for obs[" + std::to_string(idx) + "]:");
+            auto is_close_ll = isclose(loglike_img, loglike_jac, 0., max_ll_diff);
+            if (!is_close_ll.isclose) {
+                errors.push_back("loglike_jac=" + to_string_float(loglike_jac) + ", loglike_img="
+                                 + to_string_float(loglike_img) + " !close (isclose=" + is_close_ll.str()
+                                 + ") for obs[" + std::to_string(idx) + "]:");
             }
 
             const auto& image_current = *(result.second.first);

diff --git a/meson.build b/meson.build
@@ -5,6 +5,7 @@ project(
   license : 'GPL-3.0-or-later',
   default_options : ['cpp_std=c++17',],
 )
+add_project_arguments('-ffp-contract=off', language : 'cpp')
 
 eups = find_program('eups', required: false)
 use_eups = eups.found()
@@ -50,7 +51,7 @@ if not use_eups
     pkg_mod = import('pkgconfig')
     if pkg_mod.found()
         pkg_mod.generate(libraries : gauss2d_fitlib,
-                         version : '0.1',
+                         version : meson.project_version(),
                          name : 'libgauss2d_fit',
                          filebase : 'lsst_gauss2d_fit',
                          description : 'Create, manipulate and evaluate 2D Gaussian mixtures and images thereof.')

diff --git a/python/meson.build b/python/meson.build
@@ -5,7 +5,7 @@ project(
     license : 'GPL-3.0-or-later',
     version : run_command('cat', meson.source_root() / 'VERSION', check: true).stdout().strip(),
 )
-
+add_project_arguments('-ffp-contract=off', language : 'cpp')
 add_project_arguments('-DVERSION=@0@'.format(meson.project_version()), language : 'cpp')
 
 eups = find_program('eups', required: false)

diff --git a/python/pyproject.toml b/python/pyproject.toml
@@ -41,3 +41,44 @@ where = ["src"]
 
 [tool.setuptools.dynamic]
 version = { attr = "lsst.gauss2d.__version__" }
+
+[tool.black]
+line-length = 110
+target-version = ["py311"]
+
+[tool.ruff]
+exclude = [
+    "src/lsst/gauss2d/fit/__init__.py",
+    "tests/*.py",
+]
+line-length = 110
+target-version = "py311"
+
+[tool.ruff.lint]
+ignore = [
+    "N802",
+    "N803",
+    "N806",
+    "N812",
+    "N815",
+    "N816",
+    "N999",
+    "D100",
+    "D102",
+    "D104",
+    "D105",
+    "D107",
+    "D200",
+    "D203",
+    "D205",
+    "D213",
+    "D400",
+    "UP007", # Allow UNION in type annotation
+]
+select = [
+    "E", # pycodestyle
+    "F", # pycodestyle
+    "N", # pep8-naming
+    "W", # pycodestyle
+    "D", # pydocstyle
+]
diff --git a/python/src/lsst/gauss2d/fit/__init__.py b/python/src/lsst/gauss2d/fit/__init__.py
@@ -1,4 +1,5 @@
 import pkgutil
+
 __path__ = pkgutil.extend_path(__path__, __name__)
 
 from ._gauss2d_fit import *

diff --git a/python/src/lsst/gauss2d/fit/model.cc b/python/src/lsst/gauss2d/fit/model.cc
@@ -82,7 +82,7 @@ void declare_model(py::module &m, std::string str_type) {
                  "outputs_prior"_a = std::vector<std::shared_ptr<Image>>{},
                  "residuals_prior"_a = std::shared_ptr<Image>{}, "force"_a = false, "print"_a = false)
             .def("verify_jacobian", &Model::verify_jacobian, "findiff_frac"_a = 1e-4, "findiff_add"_a = 1e-4,
-                 "rtol"_a = 1e-3, "atol"_a = 1e-3)
+                 "rtol"_a = 1e-3, "atol"_a = 1e-3, "max_diff_ll"_a = 0)
             .def("__repr__", [](const Model &self) { return self.repr(true); })
             .def("__str__", &Model::str);
 }

diff --git a/python/src/meson.build b/python/src/meson.build
@@ -80,7 +80,6 @@ else
     python.install_sources(
         python_sources,
         pure: false, # Will be installed next to binaries
-        subdir: 'lsst',  # Folder relative to site-packages to install to
         preserve_path: true,
     )
 endif
diff --git a/python/tests/test_channel.py b/python/tests/test_channel.py
@@ -17,7 +17,7 @@ def test_channel():
 
     name_new = "A long name that is very unlikely to have already been taken"
     channel_new = g2f.Channel(name_new)
-    assert (len(g2f.Channel.all) == (n_channels + 1))
+    assert len(g2f.Channel.all) == (n_channels + 1)
     with pytest.raises(ValueError):
         g2f.Channel.erase("None")
     with pytest.raises(ValueError):
@@ -27,7 +27,7 @@ def test_channel():
     del channel_new
     # ...only calling erase removes it from the registry, though
     g2f.Channel.erase(name_new)
-    assert (g2f.Channel.find(name_new) is None)
-    assert (len(g2f.Channel.all) == n_channels)
-    assert (len(g2f.Channel(name_new).all) == (n_channels + 1))
-    assert (len(g2f.Channel("Another unusually long name").all) == (n_channels + 2))
+    assert g2f.Channel.find(name_new) is None
+    assert len(g2f.Channel.all) == n_channels
+    assert len(g2f.Channel(name_new).all) == (n_channels + 1)
+    assert len(g2f.Channel("Another unusually long name").all) == (n_channels + 2)
diff --git a/python/tests/test_model.py b/python/tests/test_model.py
@@ -3,38 +3,44 @@
 import numpy as np
 import pytest
 
+# TODO: Investigate why clang/osx builds have non-zero values in DM-45308
+# GCC/linux are exactly 0. Hopefully it's not something like ffast-math
+max_diff_ll = 1e-25
 
-@pytest.fixture(scope='module')
+
+@pytest.fixture(scope="module")
 def channels():
-    return tuple(g2f.Channel(x) for x in 'rgb')
+    return tuple(g2f.Channel(x) for x in "rgb")
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture(scope="module")
 def data(channels):
     n_x, n_y = 5, 5
-    return g2f.DataD([
-        g2f.ObservationD(
-            g2d.ImageD(np.zeros((n_y, n_x))),
-            g2d.ImageD(np.full((n_y, n_x), 1e3)),
-            g2d.ImageB(np.ones((n_y, n_x))),
-            channel,
-        )
-        for channel in channels
-    ])
+    return g2f.DataD(
+        [
+            g2f.ObservationD(
+                g2d.ImageD(np.zeros((n_y, n_x))),
+                g2d.ImageD(np.full((n_y, n_x), 1e3)),
+                g2d.ImageB(np.ones((n_y, n_x))),
+                channel,
+            )
+            for channel in channels
+        ]
+    )
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture(scope="module")
 def psfmodels(data):
     psfmodels = tuple(
-        g2f.PsfModel([
-            g2f.GaussianComponent(
-                g2f.GaussianParametricEllipse(1., 1., 0.),
-                None,
-                g2f.LinearIntegralModel(
-                    [(g2f.Channel.NONE, g2f.IntegralParameterD(1.))]
-                ),
-            )
-        ])
+        g2f.PsfModel(
+            [
+                g2f.GaussianComponent(
+                    g2f.GaussianParametricEllipse(1.0, 1.0, 0.0),
+                    None,
+                    g2f.LinearIntegralModel([(g2f.Channel.NONE, g2f.IntegralParameterD(1.0))]),
+                )
+            ]
+        )
         for _ in range(len(data))
     )
     for psfmodel in psfmodels:
@@ -43,12 +49,12 @@ def psfmodels(data):
     return psfmodels
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture(scope="module")
 def interptype_sersic_default():
     return g2f.SersicMixComponentIndexParameterD().interptype
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture(scope="module")
 def sources(channels, interptype_sersic_default):
     last = None
     n_sources, n_components = 2, 2
@@ -61,21 +67,24 @@ def sources(channels, interptype_sersic_default):
             is_last = c == n_comp_max
             last = g2f.FractionalIntegralModel(
                 [
-                    (channel, g2f.ProperFractionParameterD(
-                        (is_last == 1) or 0.25 * (1 + idx_channel),
-                        fixed=is_last
-                    ))
+                    (
+                        channel,
+                        g2f.ProperFractionParameterD(
+                            (is_last == 1) or 0.25 * (1 + idx_channel), fixed=is_last
+                        ),
+                    )
                     for idx_channel, channel in enumerate(channels)
                 ],
-                g2f.LinearIntegralModel([
-                    (channel, g2f.IntegralParameterD(0.5 + 0.5 * (i + 1)))
-                    for channel in channels
-                ]) if (c == 0) else last,
+                g2f.LinearIntegralModel(
+                    [(channel, g2f.IntegralParameterD(0.5 + 0.5 * (i + 1))) for channel in channels]
+                )
+                if (c == 0)
+                else last,
                 is_last,
             )
             components[c] = g2f.GaussianComponent(
                 centroid=None,
-                ellipse=g2f.GaussianParametricEllipse(1., 1., 0.),
+                ellipse=g2f.GaussianParametricEllipse(1.0, 1.0, 0.0),
                 integral=last,
             )
         components.append(
@@ -90,15 +99,17 @@ def sources(channels, interptype_sersic_default):
                         transform=g2f.LogitLimitedTransformD(
                             limits=g2f.LimitsD(min=-0.999, max=0.999, name="ref_logit_sersic[0.5, 6.0]"),
                         ),
-                    )
+                    ),
+                ),
+                integral=g2f.LinearIntegralModel(
+                    [
+                        (
+                            channel,
+                            g2f.IntegralParameterD(1.0, transform=log10, label=f"Sersic {channel}-band"),
+                        )
+                        for channel in channels
+                    ]
                 ),
-                integral=g2f.LinearIntegralModel([
-                    (
-                        channel,
-                        g2f.IntegralParameterD(1.0, transform=log10, label=f"Sersic {channel}-band")
-                    )
-                    for channel in channels
-                ]),
                 # Linear interpolation fails at exact knot values
                 # Adding a small offset solves the problem
                 sersicindex=g2f.SersicMixComponentIndexParameterD(
@@ -107,15 +118,15 @@ def sources(channels, interptype_sersic_default):
                     transform=g2f.LogitLimitedTransformD(
                         limits=g2f.LimitsD(min=0.5, max=6.0, name="ref_logit_sersic[0.5, 6.0]"),
                     ),
-                )
+                ),
             ),
         )
         sources[i] = g2f.Source(components)
 
     return sources
 
 
-@pytest.fixture(scope='module')
+@pytest.fixture(scope="module")
 def model(data, psfmodels, sources):
     model = g2f.ModelD(data, list(psfmodels), list(sources))
     with pytest.raises(RuntimeError):
@@ -148,9 +159,9 @@ def test_model_eval_jacobian(model):
     result = np.array(model.evaluate())
     # TODO: Investigate why clang/osx builds have non-zero values in DM-45308
     # GCC/linux are exactly 0. Hopefully it's not something like ffast-math
-    assert (np.abs(result) < 1e-25).all()
+    assert (np.abs(result) <= max_diff_ll).all()
     # TODO: Investigate why this rtol needs to be so high
-    errors = model.verify_jacobian(atol=1e-3, rtol=5e-3)
+    errors = model.verify_jacobian(atol=1e-3, rtol=5e-3, max_diff_ll=max_diff_ll)
     # All of the IntegralParameters got double-counted - see DM-40674
     # ... but also, linear interpolators will fail the Sersic index params
     assert len(errors) == 6
@@ -160,16 +171,18 @@ def test_model_eval_loglike_grad(model):
     model.setup_evaluators(g2f.EvaluatorMode.loglike_grad, print=True)
     result = np.array(model.evaluate())
     # TODO: Investigate why clang/osx builds have non-zero values in DM-45308
-    assert (np.abs(result) < 1e-25).all()
+    assert (np.abs(result) <= max_diff_ll).all()
     dloglike = np.array(model.compute_loglike_grad(True))
-    assert (dloglike == 0).all()
+    assert (np.abs(dloglike) < max_diff_ll).all()
 
 
 def test_model_hessian(model):
     options = g2f.HessianOptions(return_negative=False)
     for idx in range(2):
-        hessians = {is_transformed: model.compute_hessian(transformed=is_transformed, options=options)
-                    for is_transformed in (False, True)}
+        hessians = {
+            is_transformed: model.compute_hessian(transformed=is_transformed, options=options)
+            for is_transformed in (False, True)
+        }
         # Note that since there is no noise in the image and none of the
         # param values have changed, the Hessian terms are not all negative -
         # the sign corrections are based on the sign of loglike_grad,