Remove julia code and dependencies

.github/CHANGELOG.md

@@ -13,6 +13,7 @@
   [(#506)](https://github.com/XanaduAI/MrMustard/pull/506)
 
 ### Breaking changes
+* Remove julia code and dependencies. PRECISION_BITS_HERMITE_POLY setting removed.
 
 ### Improvements
 * Switch from the `julia` Python package to `juliacall` for easier installation and usage.

.github/workflows/builds.yml

@@ -45,10 +45,5 @@ jobs:
       - name: Install only test dependencies
         run: poetry install --no-root --extras "ray" --with dev
 
-      - name: Setup Julia
-        uses: julia-actions/setup-julia@v1
-        with:
-          version: 1.9.3
-
       - name: Run tests
         run: python -m pytest tests -p no:warnings --tb=native --backend=${{ matrix.backend }}

.github/workflows/tests_numpy.yml

@@ -29,21 +29,13 @@ jobs:
         with:
           python-version: ${{ matrix.python-version }}
 
-      - name: Setup Julia
-        uses: julia-actions/setup-julia@v1
-        with:
-          version: 1.9.3
-
       - name: Install dependencies
         run: |
           python -m pip install --no-cache-dir --upgrade pip
           pip install --no-cache-dir poetry==1.8.0
           poetry config virtualenvs.create false
           poetry install --extras "ray" --with dev
 
-      - name: Initialize Julia project  # to keep test times stable
-        run: python -c 'from juliacall import Main as jl'
-
       - name: Create durations file name
         id: record_file
         run: echo "filename=durations_$(date +%s)_$(git rev-parse --short HEAD).txt" >> "$GITHUB_OUTPUT"

.github/workflows/tests_tensorflow.yml

@@ -29,21 +29,13 @@ jobs:
         with:
           python-version: ${{ matrix.python-version }}
 
-      - name: Setup Julia
-        uses: julia-actions/setup-julia@v1
-        with:
-          version: 1.9.3
-
       - name: Install dependencies
         run: |
           python -m pip install --no-cache-dir --upgrade pip
           pip install --no-cache-dir poetry==1.8.0
           poetry config virtualenvs.create false
           poetry install --extras "ray" --with dev
 
-      - name: Initialize Julia project  # to keep test times stable
-        run: python -c 'from juliacall import Main as jl'
-
       - name: Create durations file name
         id: record_file
         run: echo "filename=durations_$(date +%s)_$(git rev-parse --short HEAD).txt" >> "$GITHUB_OUTPUT"

.github/workflows/upload.yml

@@ -33,11 +33,6 @@ jobs:
       - name: Install only test dependencies
         run: poetry install --no-root --extras "ray" --with dev
 
-      - name: Setup Julia
-        uses: julia-actions/setup-julia@v1
-        with:
-          version: 1.9.3
-
       - name: Run tests
         run: python -m pytest tests -p no:warnings --tb=native --backend=tensorflow
 

.readthedocs.yml

@@ -10,8 +10,6 @@ build:
     - graphviz
   jobs:
     post_install:
-      - "curl -fsSL https://install.julialang.org | sh -s -- --yes"
-      - "ln -s /home/docs/.juliaup/bin/julia /home/docs/.asdf/bin/julia"
       - pip install --no-cache-dir poetry==1.4.0
       - poetry config virtualenvs.create false
       - poetry install --with doc

README.md

@@ -31,13 +31,6 @@ Mr Mustard supports:
 - Plug-and-play backends (`numpy` as default)
 - An abstraction layer `XPTensor` for seamless symplectic algebra (experimental)
 
-# Increased numerical stability using Julia [optional]
-
-Converting phase space objects to Fock space can be numerically unstable due to accumulating floating point errors.
-To resolve this, the conversion can be performed with extended-precision arithmetic. To use this feature,
-an installation of [Julia](https://github.com/JuliaLang/juliaup#installation) is required (version 1.9.3 recommended).
-If no valid version of Julia is found, it will be installed automatically before trying to run any Julia code.
-
 # The lab module
 The lab module contains things you'd find in a lab: states, transformations, measurements, circuits. States can be used at the beginning of a circuit as well as at the end, in which case a state is interpreted as a measurement (a projection onto that state). Transformations are usually parametrized and map states to states. The action on states is differentiable with respect to the state and to the gate parameters.
 

doc/development/development_guide.rst

@@ -37,13 +37,6 @@ using development mode:
 The ``-e`` flag ensures that edits to the source code will be reflected when
 importing Mr Mustard in Python.
 
-Increased numerical stability using Julia [optional]
-------------------
-Converting phase space objects to Fock space can be numerically unstable due to accumulating floating point errors.
-To resolve this, the conversion can be performed with extended-precision arithmetic. To use this feature,
-an installation of `Julia <https://github.com/JuliaLang/juliaup#installation>`_ is required (version 1.9.3 recommended).
-If no valid version of Julia is found, it will be installed automatically before trying to run any Julia code.
-
 Development environment
 -----------------------
 

env.Dockerfile

@@ -24,10 +24,6 @@ RUN python -m pip install --no-cache-dir --upgrade pip
 # Install all dependencies
 RUN poetry install --no-root --all-extras --with dev,doc
 
-# Install Julia
-COPY juliapkg.json .
-RUN poetry run python -c 'from juliacall import Main'
-
 ENV DEBIAN_FRONTEND=dialog
 
 # Add source code, tests and configuration

mrmustard/math/backend_numpy.py

@@ -518,23 +518,11 @@ def hermite_renormalized(
             The renormalized Hermite polynomial of given shape.
         """
 
-        precision_bits = settings.PRECISION_BITS_HERMITE_POLY
         shape = tuple(int(i) for i in shape)  # ensure each item in the tuple is of the same dtype
-        if precision_bits == 128:  # numba
-            if settings.STABLE_FOCK_CONVERSION:
-                G = vanilla_stable(tuple(shape), A, b, c)
-            else:
-                G = vanilla(tuple(shape), A, b, c)
-        else:  # julia (with precision_bits = 512)
-            # The following import must come after running "jl = Julia(compiled_modules=False)" in settings.py
-            from juliacall import Main as jl  # pylint: disable=import-outside-toplevel
-
-            A, b, c = (
-                np.array(A, dtype=np.complex128),
-                np.array(b, dtype=np.complex128),
-                np.array(c, dtype=np.complex128),
-            )
-            G = jl.Vanilla.vanilla(A, b, c.item(), np.array(shape, dtype=np.int64), precision_bits)
+        if settings.STABLE_FOCK_CONVERSION:
+            G = vanilla_stable(shape, A, b, c)
+        else:
+            G = vanilla(shape, A, b, c)
 
         return G
 

mrmustard/math/backend_tensorflow.py

@@ -454,28 +454,12 @@ def hermite_renormalized(
             The renormalized Hermite polynomial of given shape.
         """
 
-        precision_bits = settings.PRECISION_BITS_HERMITE_POLY
-
         A, b, c = self.asnumpy(A), self.asnumpy(b), self.asnumpy(c)
 
-        if precision_bits == 128:  # numba
-            if settings.STABLE_FOCK_CONVERSION:
-                G = strategies.vanilla_stable(tuple(shape), A, b, c)
-            else:
-                G = strategies.vanilla(tuple(shape), A, b, c)
-        else:  # julia
-            # The following import must come after settings settings.PRECISION_BITS_HERMITE_POLY
-            from juliacall import Main as jl  # pylint: disable=import-outside-toplevel
-
-            A, b, c = (
-                A.astype(np.complex128),
-                b.astype(np.complex128),
-                c.astype(np.complex128),
-            )
-
-            G = self.astensor(
-                jl.Vanilla.vanilla(A, b, c.item(), np.array(shape, dtype=np.int64), precision_bits)
-            )
+        if settings.STABLE_FOCK_CONVERSION:
+            G = strategies.vanilla_stable(tuple(shape), A, b, c)
+        else:
+            G = strategies.vanilla(tuple(shape), A, b, c)
 
         def grad(dLdGconj):
             dLdA, dLdB, dLdC = strategies.vanilla_vjp(G, c, np.conj(dLdGconj))
@@ -587,35 +571,17 @@ def hermite_renormalized_diagonal_reorderedAB(
             The renormalized Hermite polynomial.
         """
         A, B, C = self.asnumpy(A), self.asnumpy(B), self.asnumpy(C)
-        precision_bits = settings.PRECISION_BITS_HERMITE_POLY
 
-        if precision_bits == 128:  # numba (complex128)
-            poly0, poly2, poly1010, poly1001, poly1 = tf.numpy_function(
-                hermite_multidimensional_diagonal, [A, B, C, cutoffs], [A.dtype] * 5
-            )
-        else:  # julia (higher precision than complex128)
-            # The following import must come after settings settings.PRECISION_BITS_HERMITE_POLY
-            from juliacall import Main as jl  # pylint: disable=import-outside-toplevel
-
-            (poly0, poly2, poly1010, poly1001, poly1) = (
-                self.asnumpy(val)
-                for val in jl.DiagonalAmps.fock_diagonal_amps(
-                    A, B, C.item(), tuple(cutoffs), precision_bits
-                )
-            )
+        poly0, poly2, poly1010, poly1001, poly1 = tf.numpy_function(
+            hermite_multidimensional_diagonal, [A, B, C, cutoffs], [A.dtype] * 5
+        )
 
         def grad(dLdpoly):
-            if precision_bits == 128:  # numba (complex128)
-                dpoly_dC, dpoly_dA, dpoly_dB = tf.numpy_function(
-                    grad_hermite_multidimensional_diagonal,
-                    [A, B, C.item(), poly0, poly2, poly1010, poly1001, poly1],
-                    [poly0.dtype] * 3,
-                )
-            else:  # julia (higher precision than complex128)
-                dpoly_dC = poly0 / C.item()
-                dpoly_dA, dpoly_dB = jl.DiagonalGrad.fock_diagonal_grad(
-                    A, B, poly0, poly2, poly1010, poly1001, poly1, precision_bits
-                )
+            dpoly_dC, dpoly_dA, dpoly_dB = tf.numpy_function(
+                grad_hermite_multidimensional_diagonal,
+                [A, B, C.item(), poly0, poly2, poly1010, poly1001, poly1],
+                [poly0.dtype] * 3,
+            )
 
             ax = tuple(range(dLdpoly.ndim))
             dLdA = self.sum(dLdpoly[..., None, None] * self.conj(dpoly_dA), axis=ax)
@@ -685,43 +651,19 @@ def hermite_renormalized_1leftoverMode_reorderedAB(
             The renormalized Hermite polynomial.
         """
         A, B, C = self.asnumpy(A), self.asnumpy(B), self.asnumpy(C)
-        precision_bits = settings.PRECISION_BITS_HERMITE_POLY
 
-        if precision_bits == 128:  # numba (complex128)
-            poly0, poly2, poly1010, poly1001, poly1 = tf.numpy_function(
-                hermite_multidimensional_1leftoverMode,
-                [A, B, C.item(), cutoffs],
-                [A.dtype] * 5,
-            )
-        else:  # julia (higher precision than complex128)
-            # The following import must come after settings settings.PRECISION_BITS_HERMITE_POLY
-            from juliacall import Main as jl  # pylint: disable=import-outside-toplevel
-
-            (
-                poly0,
-                poly2,
-                poly1010,
-                poly1001,
-                poly1,
-            ) = (
-                self.asnumpy(val)
-                for val in jl.LeftoverModeAmps.fock_1leftoverMode_amps(
-                    A, B, C.item(), tuple(cutoffs), precision_bits
-                )
-            )
+        poly0, poly2, poly1010, poly1001, poly1 = tf.numpy_function(
+            hermite_multidimensional_1leftoverMode,
+            [A, B, C.item(), cutoffs],
+            [A.dtype] * 5,
+        )
 
         def grad(dLdpoly):
-            if precision_bits == 128:  # numba (complex128)
-                dpoly_dC, dpoly_dA, dpoly_dB = tf.numpy_function(
-                    grad_hermite_multidimensional_1leftoverMode,
-                    [A, B, C, poly0, poly2, poly1010, poly1001, poly1],
-                    [poly0.dtype] * 3,
-                )
-            else:  # julia (higher precision than complex128)
-                dpoly_dC = poly0 / C.item()
-                (dpoly_dA, dpoly_dB) = jl.LeftoverModeGrad.fock_1leftoverMode_grad(
-                    A, B, poly0, poly2, poly1010, poly1001, poly1, precision_bits
-                )
+            dpoly_dC, dpoly_dA, dpoly_dB = tf.numpy_function(
+                grad_hermite_multidimensional_1leftoverMode,
+                [A, B, C, poly0, poly2, poly1010, poly1001, poly1],
+                [poly0.dtype] * 3,
+            )
 
             ax = tuple(range(dLdpoly.ndim))
             dLdA = self.sum(dLdpoly[..., None, None] * self.conj(dpoly_dA), axis=ax)