Fix pfft gradients (#34)

* update jaxdecomp version and test gradients

* Prepare for DTO tests

* format

.github/workflows/python-publish.yml

@@ -0,0 +1,55 @@
+name: Upload Python Package
+
+on:
+  release:
+    types: [published]
+
+permissions:
+  contents: read
+
+jobs:
+  release-build:
+    runs-on: ubuntu-latest
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.x"
+
+      - name: Build release distributions
+        run: |
+          # NOTE: put your own distribution build steps here.
+          python -m pip install build
+          python -m build
+
+      - name: Upload distributions
+        uses: actions/upload-artifact@v4
+        with:
+          name: release-dists
+          path: dist/
+
+  pypi-publish:
+    runs-on: ubuntu-latest
+    needs:
+      - release-build
+    permissions:
+      # IMPORTANT: this permission is mandatory for trusted publishing
+      id-token: write
+
+    environment:
+      name: pypi
+      url: https://pypi.org/p/jaxpm
+
+    steps:
+      - name: Retrieve release distributions
+        uses: actions/download-artifact@v4
+        with:
+          name: release-dists
+          path: dist/
+
+      - name: Publish release distributions to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          packages-dir: dist/

.github/workflows/tests.yml

@@ -34,7 +34,8 @@ jobs:
         pip install git+https://github.com/MP-Gadget/pfft-python
         pip install git+https://github.com/MP-Gadget/pmesh
         pip install git+https://github.com/ASKabalan/fastpm-python --no-build-isolation
-        pip install .[test]
+        pip install -r requirements-test.txt
+        pip install .
 
     - name: Run Single Device Tests
       run: |

.gitignore

@@ -132,3 +132,6 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+# Hide version file
+_version.py

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2021 Differentiable Universe Initiative
+Copyright (c) 2021-2025 Differentiable Universe Initiative
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

MANIFEST.in

@@ -0,0 +1,2 @@
+prune notebooks
+prune tests

README.md

@@ -1,9 +1,26 @@
 # JaxPM
-[![Tests](https://github.com/DifferentiableUniverseInitiative/JaxPM/actions/workflows/tests.yml/badge.svg)](https://github.com/DifferentiableUniverseInitiative/JaxPM/actions/workflows/tests.yml) <!-- ALL-CONTRIBUTORS-BADGE:START - Do not remove or modify this section -->
+[![Notebook](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/DifferentiableUniverseInitiative/JaxPM/blob/main/notebooks/01-Introduction.ipynb)
+[![PyPI version](https://img.shields.io/pypi/v/jaxpm)](https://pypi.org/project/jaxpm/) [![Tests](https://github.com/DifferentiableUniverseInitiative/JaxPM/actions/workflows/tests.yml/badge.svg)](https://github.com/DifferentiableUniverseInitiative/JaxPM/actions/workflows/tests.yml) <!-- ALL-CONTRIBUTORS-BADGE:START - Do not remove or modify this section -->
 [![All Contributors](https://img.shields.io/badge/all_contributors-5-orange.svg?style=flat-square)](#contributors-)
 <!-- ALL-CONTRIBUTORS-BADGE:END -->
 JAX-powered Cosmological Particle-Mesh N-body Solver
 
+> ### Note
+> **The new JaxPM v0.1.xx** supports multi-GPU model distribution while remaining compatible with previous releases. These significant changes are still under development and testing, so please report any issues you encounter.
+> For the older but more stable version, install:
+> ```bash
+> pip install jaxpm==0.0.2
+> ```
+
+## Install
+
+Basic installation can be done using pip:
+```bash
+pip install jaxpm
+```
+For more advanced installation for optimized distribution on gpu clusters, please install jaxDecomp first. See instructions [here](https://github.com/DifferentiableUniverseInitiative/jaxDecomp).
+
+
 ## Goals
 
 Provide a modern infrastructure to support differentiable PM N-body simulations using JAX:

pyproject.toml

@@ -3,7 +3,7 @@ requires = ["setuptools", "wheel", "setuptools-scm"]
 build-backend = "setuptools.build_meta"
 
 [project]
-name = "JaxPM"
+name = "jaxpm"
 dynamic = ["version"]
 description = "A simple Particle-Mesh implementation in JAX"
 authors = [{ name = "JaxPM developers" }]
@@ -11,20 +11,7 @@ readme = "README.md"
 requires-python = ">=3.9"
 license = { file = "LICENSE" }
 urls = { "Homepage" = "https://github.com/DifferentiableUniverseInitiative/JaxPM" }
-dependencies = ["jax_cosmo", "jax>=0.4.30", "jaxdecomp>=0.2.2"]
+dependencies = ["jax_cosmo", "jax>=0.4.35", "jaxdecomp>=0.2.3"]
-
-[project.optional-dependencies]
-test = [
-    "jax>=0.4.30",
-    "numpy",
-    "jax_cosmo",
-    "jaxdecomp>=0.2.2",
-    "pytest>=8.0.0",
-    "pfft-python @ git+https://github.com/MP-Gadget/pfft-python",
-    "pmesh @ git+https://github.com/MP-Gadget/pmesh",
-    "fastpm @ git+https://github.com/ASKabalan/fastpm-python",
-    "diffrax"
-]
 
 [tool.setuptools]
 packages = ["jaxpm"]

requirements-test.txt

@@ -0,0 +1,5 @@
+pytest>=8.0.0
+diffrax
+pfft-python @ git+https://github.com/MP-Gadget/pfft-python
+pmesh @ git+https://github.com/MP-Gadget/pmesh
+fastpm @ git+https://github.com/ASKabalan/fastpm-python

tests/test_gradients.py

@@ -0,0 +1,87 @@
+import jax
+import pytest
+from diffrax import (BacksolveAdjoint, Dopri5, ODETerm, PIDController,
+                     RecursiveCheckpointAdjoint, SaveAt, diffeqsolve)
+from helpers import MSE
+from jax import numpy as jnp
+
+from jaxpm.distributed import uniform_particles
+from jaxpm.painting import cic_paint, cic_paint_dx
+from jaxpm.pm import lpt, make_diffrax_ode
+
+
+@pytest.mark.single_device
+@pytest.mark.parametrize("order", [1, 2])
+@pytest.mark.parametrize("absolute_painting", [True, False])
+@pytest.mark.parametrize("adjoint", ['DTO', 'OTD'])
+def test_nbody_grad(simulation_config, initial_conditions, lpt_scale_factor,
+                    nbody_from_lpt1, nbody_from_lpt2, cosmo, order,
+                    absolute_painting, adjoint):
+
+    mesh_shape, _ = simulation_config
+    cosmo._workspace = {}
+
+    if adjoint == 'OTD':
+        pytest.skip("OTD adjoint not implemented yet (needs PFFT3D JVP)")
+
+    adjoint = RecursiveCheckpointAdjoint(
+    ) if adjoint == 'DTO' else BacksolveAdjoint(solver=Dopri5())
+
+    @jax.jit
+    @jax.grad
+    def forward_model(initial_conditions, cosmo):
+
+        # Initial displacement
+        if absolute_painting:
+            particles = uniform_particles(mesh_shape)
+            dx, p, _ = lpt(cosmo,
+                           initial_conditions,
+                           particles,
+                           a=lpt_scale_factor,
+                           order=order)
+            ode_fn = ODETerm(make_diffrax_ode(cosmo, mesh_shape))
+            y0 = jnp.stack([particles + dx, p])
+
+        else:
+            dx, p, _ = lpt(cosmo,
+                           initial_conditions,
+                           a=lpt_scale_factor,
+                           order=order)
+            ode_fn = ODETerm(
+                make_diffrax_ode(cosmo, mesh_shape, paint_absolute_pos=False))
+            y0 = jnp.stack([dx, p])
+
+        solver = Dopri5()
+        controller = PIDController(rtol=1e-7,
+                                   atol=1e-7,
+                                   pcoeff=0.4,
+                                   icoeff=1,
+                                   dcoeff=0)
+
+        saveat = SaveAt(t1=True)
+
+        solutions = diffeqsolve(ode_fn,
+                                solver,
+                                t0=lpt_scale_factor,
+                                t1=1.0,
+                                dt0=None,
+                                y0=y0,
+                                adjoint=adjoint,
+                                stepsize_controller=controller,
+                                saveat=saveat)
+
+        if absolute_painting:
+            final_field = cic_paint(jnp.zeros(mesh_shape), solutions.ys[-1, 0])
+        else:
+            final_field = cic_paint_dx(solutions.ys[-1, 0])
+
+        return MSE(final_field,
+                   nbody_from_lpt1 if order == 1 else nbody_from_lpt2)
+
+    bad_initial_conditions = initial_conditions + jax.random.normal(
+        jax.random.PRNGKey(0), initial_conditions.shape) * 0.5
+    best_ic = forward_model(initial_conditions, cosmo)
+    bad_ic = forward_model(bad_initial_conditions, cosmo)
+
+    assert jnp.max(best_ic) < 1e-5
+    assert jnp.max(bad_ic) > 1e-5