tests+fix: starry comparison (flux)

tests+fix: starry comparison

.github/workflows/ci.yml

@@ -12,18 +12,66 @@ on:
   workflow_dispatch:
 
 jobs:
+  tests:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: ["ubuntu-latest"]
+        python-version: ["3.9", "3.11"]
+        session:
+          - "test"
+        include:
+          - os: "ubuntu-latest"
+            python-version: "3.11"
+            session: "test"
+          - os: "macos-latest"
+            python-version: "3.11"
+            session: "test"
+          - os: "ubuntu-latest"
+            python-version: "3.9"
+            session: "comparison"
+
+    steps:
+      - name: "Init: checkout"
+        uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+          submodules: true
+
+      - name: "Init: Python"
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ matrix.python-version }}
+
+      - name: "Install: dependencies"
+        run: |
+          python -m pip install -U pip
+          python -m pip install -U nox
+
+      - name: "Tests: run"
+        run: |
+          python -m nox --non-interactive --error-on-missing-interpreter \
+            --session "${{matrix.session}}-${{matrix.python-version}}"
+
   build:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
-      - name: Set up Python 3.9
-        uses: actions/setup-python@v2
+      - uses: actions/checkout@v4
         with:
-          python-version: 3.9
+          fetch-depth: 0
+      - uses: actions/setup-python@v5
+        name: Install Python
+        with:
+          python-version: "3.11"
       - name: Install dependencies
         run: |
-          python -m pip install -U pip poetry
-          python -m poetry install --with dev
-      - name: Test with pytest
-        run: |
-          python -m poetry run python -m pytest
+          python -m pip install -U pip
+          python -m pip install -U build twine
+      - name: Build the distribution
+        run: python -m build .
+      - name: Check the distribution
+        run: python -m twine check --strict dist/*
+      - uses: actions/upload-artifact@v4
+        with:
+          path: dist/*

noxfile.py

@@ -0,0 +1,29 @@
+# mypy: ignore-errors
+
+import nox
+
+ALL_PYTHON_VS = ["3.9", "3.11"]
+
+
+@nox.session(python=ALL_PYTHON_VS)
+def test(session):
+    session.install(".[test]")
+    session.run("pytest", *session.posargs)
+
+
+@nox.session(python=["3.9"])
+def comparison(session):
+    session.install(".[test,comparison]", "numpy<1.22")
+    session.run("python", "-c", "import starry")
+    session.run("python", "-c", "import theano")
+    if session.posargs:
+        args = session.posargs
+    else:
+        args = ("tests/starry_comparison",)
+    session.run(
+        "pytest",
+        "-n",
+        "auto",
+        *args,
+        env={"JAX_ENABLE_X64": "True"},
+    )

pyproject.toml

@@ -1,36 +1,30 @@
-[tool.poetry]
+[project]
 name = "spotter"
-version = "0.0.5-beta"
+version = "0.0.6-beta"
 description = "Stellar contamination estimates from rotational light curves"
-authors = ["Lionel Garcia", "Benjamin Rackham"]
+authors = [{name="Lionel Garcia"}, {name="Benjamin Rackham"}]
 license = "MIT"
 readme = "readme.md"
+requires-python = ">=3.9"
 packages = [{ include = "spotter" },]
+dependencies = ["numpy", "healpy", "matplotlib", "jax", "jaxlib"]
 
-[tool.poetry.dependencies]
-python = "^3.9"
-numpy = "^1.25.0"
-healpy = "*"
-matplotlib = "*"
-jax = "*"
-jaxlib = "*"
-
-[tool.poetry.group.dev.dependencies]
-pytest = "*"
-black = "*"
-
-[tool.poetry.group.docs.dependencies]
-sphinx = "*"
-docutils = "*"
-jupyterlab = "*"
-myst-parser = "*"
-sphinx-book-theme = "^1.0.0"
-myst-nb = "*"
-sphinx-copybutton = "*"
-toml = "*"
-ipywidgets = "*"
-tqdm = "*"
+[project.optional-dependencies]
+dev = ["black", "pytest", "nox"]
+test = ["pytest", "pytest-xdist"]
+comparison = ["starry", "tqdm", "xarray<2023.10.0", "numpy<1.22"]
+docs = [
+    "sphinx",
+    "docutils",
+    "jupyterlab",
+    "myst-parser", 
+    "sphinx-book-theme", 
+    "myst-nb", 
+    "sphinx-copybutton",
+    "toml", 
+    "ipywidgets"
+    ]
 
 [build-system]
-requires = ["poetry-core"]
-build-backend = "poetry.core.masonry.api"
+requires = ["hatchling", "hatch-vcs"]
+build-backend = "hatchling.build"

readme.md

@@ -11,7 +11,7 @@
   <br>
   <p align="center">
     <a href="https://github.com/lgrcia/spotter">
-      <img src="https://img.shields.io/badge/github-lgrcia/spotter-indianred.svg?style=flat" alt="github"/></a>
+      <img src="https://img.shields.io/badge/github-lgrcia/spotter-e3a8a1.svg?style=flat" alt="github"/></a>
     <a href="LICENCE">
       <img src="https://img.shields.io/badge/license-MIT-lightgray.svg?style=flat" alt="license"/>
     </a>

spotter/star.py

@@ -338,7 +338,9 @@ def jax_flux(self, phases):
         def flux(spot_map):
             _spot = (1 - spot_map) * limb_darkening
             _geometry = mask * projected_area
-            return (_spot * _geometry).sum(1) / _geometry.sum(1)
+            return (
+                np.pi * (_spot * _geometry).sum(1) / (_geometry * limb_darkening).sum(1)
+            )
 
         return flux
 
@@ -446,17 +448,25 @@ def flux(self, phases):
         projected_area = np.vectorize(
             core.projected_area(self._thetas, self._phis), signature="()->(n)"
         )(phases)
-        limb_darkening = np.vectorize(
-            core.polynomial_limb_darkening(self._thetas, self._phis),
-            signature="()->(n)",
-            excluded={0},
-        )(self.u, phases)
+        limb_darkening = (
+            np.vectorize(
+                core.polynomial_limb_darkening(self._thetas, self._phis),
+                signature="()->(n)",
+                excluded={0},
+            )(self.u, phases)
+            if len(self.u) > 0
+            else 1
+        )
         _spot = (1 - self.map_spot) * limb_darkening
         _geometry = mask * projected_area
         # faculae contribution, with same ld for now (TODO)
-        _faculae = self.map_faculae * limb_darkening
+        _faculae = 0  # self.map_faculae * limb_darkening
 
-        return ((_spot + _faculae) * _geometry).sum(1) / _geometry.sum(1)
+        return (
+            np.pi
+            * ((_spot + _faculae) * _geometry).sum(1)
+            / (_geometry * limb_darkening).sum(1)
+        )
 
     def map(self, phase=None, limb_darkening=False):
         """
@@ -487,7 +497,8 @@ def map(self, phase=None, limb_darkening=False):
         faculae_limb_brightening = 1
         m = (1 - self.map_spot) * mask * spot_limb_darkening
         spots = self.map_spot == 0.0
-        m[spots] = m[spots] + (self.map_faculae * faculae_limb_brightening)[spots]
+        if np.any(spots):
+            m[spots] = m[spots] + (self.map_faculae * faculae_limb_brightening)[spots]
         return m
 
     def show(

tests/starry_comparison/test_flux.py

@@ -0,0 +1,79 @@
+from collections import defaultdict
+
+import healpy as hp
+import numpy as np
+import pytest
+
+from spotter import Star
+
+
+@pytest.mark.parametrize("deg", (3, 10))
+@pytest.mark.parametrize("u", ([], [0.1, 0.4]))
+def test_starry(deg, u):
+    starry = pytest.importorskip("starry")
+    starry.config.lazy = False
+
+    # starry map with random coefficients
+    np.random.seed(deg + len(u))
+    y = np.random.randn((deg + 1) ** 2)
+    y[0] = 1.0
+    y[1:] *= 1e-2
+    inc = np.pi / 2
+    ms = starry.Map(ydeg=deg, udeg=len(u), inc=np.rad2deg(inc))
+    if len(u) > 0:
+        ms[1:] = u
+    ms[:, :] = y
+
+    # rotation to map healpix
+    _m = starry._core.core.OpsYlm(deg, 0, 0, 1)
+    ry = _m.dotR([y], 0.0, 1.0, 0.0, np.pi / 2)[0]
+    ry = _m.dotR([ry], 1.0, 0.0, 0.0, -np.pi / 2)[0]
+    ry = _m.dotR([ry], 0.0, 0.0, 1.0, np.pi)[0]
+
+    def starry2healpy(y):
+        # Edmonds to Condon-Shortley phase convention
+        lmax = int(np.floor(np.sqrt(len(y))))
+
+        _hy = defaultdict(lambda: 0 + 0j)
+
+        i = 0
+
+        for l in range(0, lmax):
+            for m in range(-l, l + 1):
+                j = hp.sphtfunc.Alm.getidx(lmax, l, np.abs(m))
+                if m < 0:
+                    _hy[j] += 1j * y[i] / (np.sqrt(2) * (-1) ** m)
+                elif m == 0:
+                    _hy[j] += y[i]
+                else:
+                    _hy[j] += y[i] / (np.sqrt(2) * (-1) ** m)
+                i += 1
+
+        hn = hp.sphtfunc.Alm.getsize(lmax)
+        hy = np.zeros(hn, dtype=np.complex128)
+        for i in range(hn):
+            hy[i] = _hy[i]
+
+        return hy
+
+    # starry to healpix to spotter
+    N = 2**7
+    y2 = starry2healpy(ry)
+    mh = hp.alm2map(y2, nside=N)
+
+    # mh with same ptp as ims
+    ims = ms.render(projection="moll")
+    mh = mh - np.nanmin(mh)
+    mh = mh / np.nanmax(mh)
+    mh = mh * (np.nanmax(ims) - np.nanmin(ims))
+    mh = mh + np.nanmin(ims)
+
+    star = Star(N=N, u=u)
+    star.map_spot = 1 - mh
+
+    # comparison
+    phases = np.linspace(0, 2 * np.pi, 100)
+    expected = np.array(ms.flux(theta=np.rad2deg(phases)))
+    calc = star.flux(phases)
+
+    np.testing.assert_allclose(calc, expected, atol=1e-4)