Merge pull request #327 from adrn/fix-solout

Fix compiler issue with solout in dop853 integrator

.github/workflows/tests.yml

@@ -85,8 +85,8 @@ jobs:
           token: ${{ secrets.GITHUB_TOKEN }}
         id: setup-ffmpeg
 
-      - name: Set up Python ${{ matrix.python-version }} on ${{ matrix.os }}
-        if: !startsWith(matrix.os, 'windows')
+      - name: Set up Python ${{ matrix.python }} on ${{ matrix.os }}
+        if: ${{ !startsWith(matrix.os, 'windows') }}
         uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python }}
@@ -128,12 +128,12 @@ jobs:
 
       # Any *nix:
       - name: Install Python dependencies - nix
-        if: "!startsWith(matrix.os, 'windows')"
+        if: ${{ !startsWith(matrix.os, 'windows') }}
         run: python -m pip install --upgrade tox codecov
 
       - name: Run tests - nix
-        if: "!startsWith(matrix.os, 'windows')"
-        run: tox ${{ matrix.toxargs }} -e ${{ matrix.toxenv }} -- ${{ matrix.toxposargs }}
+        if: ${{ !startsWith(matrix.os, 'windows') }}
+        run: tox -e ${{ matrix.toxenv }} -- ${{ matrix.toxposargs }}
 
       # Coverage:
       - name: Upload coverage report to codecov

CHANGES.rst

@@ -22,6 +22,7 @@ Bug fixes
 - Fixed a bug in ``Orbit.estimate_period()`` that would cause the method to fail with a
   ``UnitsError`` if one orbit returned a nan value for the period.
 
+- Fixed a bug when compiling the ``dop853`` integrator.
 
 API changes
 -----------

gala/integrate/cyintegrators/dop853.pyx

@@ -55,16 +55,15 @@ cdef extern from "stdio.h":
     ctypedef struct FILE
     FILE *stdout
 
-cdef void solout(long nr, double xold, double x, double* y, unsigned n, int* irtrn):
-    # TODO: see here for example in FORTRAN: http://www.unige.ch/~hairer/prog/nonstiff/dr_dop853.f
-    pass
 
 cdef void dop853_step(CPotential *cp, CFrameType *cf, FcnEqDiff F,
                       double *w, double t1, double t2, double dt0,
                       int ndim, int norbits, int nbody, void *args,
                       double atol, double rtol, int nmax) except *:
 
-    cdef int res
+    cdef:
+        int res
+        SolTrait solout = NULL
 
     res = dop853(ndim*norbits, F,
                  cp, cf, norbits, nbody, args, t1, w, t2,

gala/potential/potential/tests/test_interop_galpy.py

@@ -17,15 +17,12 @@
 
 # Set these globally!
 ro = 8.122 * u.kpc
-vo = 245 * u.km/u.s
+vo = 245 * u.km / u.s
 
 if HAS_GALPY:
     import galpy.potential as galpy_gp
 
-    from gala.potential.potential.interop import (
-        _gala_to_galpy,
-        _galpy_to_gala
-    )
+    from gala.potential.potential.interop import _gala_to_galpy, _galpy_to_gala
 
 
 def pytest_generate_tests(metafunc):
@@ -39,17 +36,17 @@ def pytest_generate_tests(metafunc):
     # Test the Gala -> Galpy direction
     for Potential in _gala_to_galpy.keys():
         init = {}
-        len_scale = 1.
+        len_scale = 1.0
         for k, par in Potential._parameters.items():
-            if k == 'm':
+            if k == "m":
                 val = 1.43e10 * u.Msun
-            elif par.physical_type == 'length':
+            elif par.physical_type == "length":
                 val = 5.12 * u.kpc * len_scale
                 len_scale *= 0.5
-            elif par.physical_type == 'dimensionless':
-                val = 1.
-            elif par.physical_type == 'speed':
-                val = 201.41 * u.km/u.s
+            elif par.physical_type == "dimensionless":
+                val = 1.0
+            elif par.physical_type == "speed":
+                val = 201.41 * u.km / u.s
             else:
                 continue
 
@@ -85,25 +82,27 @@ def pytest_generate_tests(metafunc):
         gala_pots.append(pot)
         galpy_pots.append(galpy_pot)
 
-    test_names = [f'{g1.__class__.__name__}:{g2.__class__.__name__}'
-                  for g1, g2 in zip(gala_pots, galpy_pots)]
+    test_names = [
+        f"{g1.__class__.__name__}:{g2.__class__.__name__}"
+        for g1, g2 in zip(gala_pots, galpy_pots)
+    ]
 
-    metafunc.parametrize(['gala_pot', 'galpy_pot'],
-                         list(zip(gala_pots, galpy_pots)),
-                         ids=test_names)
+    metafunc.parametrize(
+        ["gala_pot", "galpy_pot"], list(zip(gala_pots, galpy_pots)), ids=test_names
+    )
 
 
-@pytest.mark.skipif(not HAS_GALPY,
-                    reason="must have galpy installed to run these tests")
+@pytest.mark.skipif(
+    not HAS_GALPY, reason="must have galpy installed to run these tests"
+)
 class TestGalpy:
-
     def setup_method(self):
         # Test points:
         rng = np.random.default_rng(42)
         ntest = 4
 
         Rs = rng.uniform(1, 15, size=ntest) * u.kpc
-        phis = rng.uniform(0, 2*np.pi, size=ntest) * u.radian
+        phis = rng.uniform(0, 2 * np.pi, size=ntest) * u.radian
         zs = rng.uniform(1, 15, size=ntest) * u.kpc
 
         cyl = CylindricalRepresentation(Rs, phis, zs)
@@ -121,68 +120,83 @@ def setup_method(self):
         Jac[:, 0, 1] = xyz[1] / cyl.rho
         Jac[:, 1, 0] = (-xyz[1] / cyl.rho**2).to_value(1 / ro)
         Jac[:, 1, 1] = (xyz[0] / cyl.rho**2).to_value(1 / ro)
-        Jac[:, 2, 2] = 1.
+        Jac[:, 2, 2] = 1.0
         self.Jac = Jac
 
     def test_density(self, gala_pot, galpy_pot):
         if isinstance(gala_pot, gp.LogarithmicPotential):
             pytest.skip()
 
         gala_val = gala_pot.density(self.xyz).to_value(u.Msun / u.pc**3)
-        galpy_val = np.array([galpy_gp.evaluateDensities(galpy_pot,
-                                                         R=RR, z=zz, phi=pp)
-                              for RR, pp, zz in self.Rpz_iter])
+        galpy_val = np.array(
+            [
+                galpy_gp.evaluateDensities(galpy_pot, R=RR, z=zz, phi=pp)
+                for RR, pp, zz in self.Rpz_iter
+            ]
+        )
         assert np.allclose(gala_val, galpy_val)
 
     def test_energy(self, gala_pot, galpy_pot):
         gala_val = gala_pot.energy(self.xyz).to_value(u.km**2 / u.s**2)
-        galpy_val = np.array([galpy_gp.evaluatePotentials(galpy_pot,
-                                                          R=RR, z=zz, phi=pp)
-                              for RR, pp, zz in self.Rpz_iter])
+        galpy_val = np.array(
+            [
+                galpy_gp.evaluatePotentials(galpy_pot, R=RR, z=zz, phi=pp)
+                for RR, pp, zz in self.Rpz_iter
+            ]
+        )
 
         if isinstance(gala_pot, gp.LogarithmicPotential):
             # Logarithms are weird
-            gala_val -= (0.5 * gala_pot.parameters['v_c']**2 *
-                         np.log(ro.value**2)).to_value((u.km / u.s)**2)
+            gala_val -= (
+                0.5 * gala_pot.parameters["v_c"] ** 2 * np.log(ro.value**2)
+            ).to_value((u.km / u.s) ** 2)
 
         assert np.allclose(gala_val, galpy_val)
 
     def test_gradient(self, gala_pot, galpy_pot):
         gala_grad = gala_pot.gradient(self.xyz)
-        gala_grad = gala_grad.to_value(u.km/u.s/u.Myr)
+        gala_grad = gala_grad.to_value(u.km / u.s / u.Myr)
 
         # TODO: Starting with galpy 1.7, this has been failing because of a
         # units issue with dPhi/dphi
         if isinstance(gala_pot, gp.LongMuraliBarPotential):
             pytest.skip()
 
-        galpy_dR = np.array([-galpy_gp.evaluateRforces(galpy_pot,
-                                                       R=RR, z=zz, phi=pp)
-                            for RR, pp, zz in self.Rpz_iter])
-        galpy_dp = np.array([-galpy_gp.evaluatephiforces(galpy_pot,
-                                                         R=RR, z=zz, phi=pp)
-                            for RR, pp, zz in self.Rpz_iter])
-        galpy_dp = (galpy_dp*(u.km/u.s)**2).to_value(vo**2)
-
-        galpy_dz = np.array([-galpy_gp.evaluatezforces(galpy_pot,
-                                                       R=RR, z=zz, phi=pp)
-                            for RR, pp, zz in self.Rpz_iter])
-        galpy_dRpz = np.stack((galpy_dR, galpy_dp, galpy_dz),
-                              axis=1)
-
-        galpy_grad = np.einsum('nij,ni->nj', self.Jac, galpy_dRpz).T
+        galpy_dR = np.array(
+            [
+                -galpy_gp.evaluateRforces(galpy_pot, R=RR, z=zz, phi=pp)
+                for RR, pp, zz in self.Rpz_iter
+            ]
+        )
+        galpy_dp = np.array(
+            [
+                -galpy_gp.evaluatephitorques(galpy_pot, R=RR, z=zz, phi=pp)
+                for RR, pp, zz in self.Rpz_iter
+            ]
+        )
+        galpy_dp = (galpy_dp * (u.km / u.s) ** 2).to_value(vo**2)
+
+        galpy_dz = np.array(
+            [
+                -galpy_gp.evaluatezforces(galpy_pot, R=RR, z=zz, phi=pp)
+                for RR, pp, zz in self.Rpz_iter
+            ]
+        )
+        galpy_dRpz = np.stack((galpy_dR, galpy_dp, galpy_dz), axis=1)
+
+        galpy_grad = np.einsum("nij,ni->nj", self.Jac, galpy_dRpz).T
 
         assert np.allclose(gala_grad, galpy_grad)
 
     def test_vcirc(self, gala_pot, galpy_pot):
         tmp = self.xyz.copy()
-        tmp[2] = 0.
+        tmp[2] = 0.0
 
-        if (not hasattr(galpy_pot, 'vcirc')
-                or isinstance(gala_pot, gp.LongMuraliBarPotential)):
+        if not hasattr(galpy_pot, "vcirc") or isinstance(
+            gala_pot, gp.LongMuraliBarPotential
+        ):
             pytest.skip()
 
-        gala_vcirc = gala_pot.circular_velocity(tmp).to_value(u.km/u.s)
-        galpy_vcirc = np.array([galpy_pot.vcirc(R=RR)
-                                for RR, *_ in self.Rpz_iter])
+        gala_vcirc = gala_pot.circular_velocity(tmp).to_value(u.km / u.s)
+        galpy_vcirc = np.array([galpy_pot.vcirc(R=RR) for RR, *_ in self.Rpz_iter])
         assert np.allclose(gala_vcirc, galpy_vcirc)