From ab6c122dc5402dc21e289ab3da0f1e732e92bdef Mon Sep 17 00:00:00 2001 From: Joseph Abbott Date: Fri, 22 May 2026 15:12:44 +0200 Subject: [PATCH] Finish off `spherical_to_cartesian` --- .../featomic/clebsch_gordan/__init__.py | 5 +- .../clebsch_gordan/_cartesian_spherical.py | 384 ++++++++++++++++++ .../featomic/clebsch_gordan/_coefficients.py | 97 +++++ .../featomic/clebsch_gordan/_dispatch.py | 14 + .../clebsch_gordan/cartesian_spherical.py | 60 ++- 5 files changed, 555 insertions(+), 5 deletions(-) diff --git a/python/featomic/featomic/clebsch_gordan/__init__.py b/python/featomic/featomic/clebsch_gordan/__init__.py index 058329b83..5aa0a162f 100644 --- a/python/featomic/featomic/clebsch_gordan/__init__.py +++ b/python/featomic/featomic/clebsch_gordan/__init__.py @@ -1,4 +1,7 @@ -from ._cartesian_spherical import cartesian_to_spherical # noqa: F401 +from ._cartesian_spherical import ( # noqa: F401 + cartesian_to_spherical, + spherical_to_cartesian, +) from ._cg_product import ClebschGordanProduct # noqa: F401 from ._coefficients import calculate_cg_coefficients # noqa: F401 from ._density_correlations import DensityCorrelations # noqa: F401 diff --git a/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py b/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py index 7fd7c96c3..c730b5298 100644 --- a/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py +++ b/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py @@ -308,6 +308,226 @@ def cartesian_to_spherical( ) +def spherical_to_cartesian( + tensor: TensorMap, + new_component_names: Optional[List[str]] = None, + cg_backend: Optional[str] = None, + cg_coefficients: Optional[TensorMap] = None, +) -> TensorMap: + """ + Transform a ``tensor`` from spherical form back to cartesian form. + + This is the exact inverse of :py:func:`cartesian_to_spherical`. Starting from a + tensor expressed on a basis of spherical harmonics ``Y^M_L`` (as produced by + :py:func:`cartesian_to_spherical` with ``keep_l_in_keys=True`` and + ``remove_blocks_threshold=None``), this function reconstructs the corresponding + cartesian product-basis tensor. + + The input keys must contain ``o3_lambda`` and ``o3_sigma``, and the single ``o3_mu`` + component will be replaced by ``tensor_rank`` cartesian components (each with values + ``[0, 1, 2]``). For tensors of rank 3 and above the keys must also contain the + ``l_{i}`` and ``k_{i}`` dimensions describing the coupling history (these are added + automatically by :py:func:`cartesian_to_spherical`). For rank 1 and 2 tensors, the + missing ``l_{i}`` dimensions are re-created internally and assumed to be 1. + + :param tensor: input tensor, using a spherical harmonics basis + :param new_component_names: names to give to the reconstructed cartesian components. + There should be exactly ``tensor_rank`` of them. If ``None``, the components are + named ``xyz_1``, ``xyz_2``, ... + :param cg_backend: Backend to use for Clebsch-Gordan calculations. This can be + ``"python-dense"`` or ``"python-sparse"``. If ``None``, this is automatically + determined. + :param cg_coefficients: Cache containing Clebsch-Gordan coefficients. This is + optional except when using this function from TorchScript. The coefficients + should be computed with :py:func:`calculate_cg_coefficients`, using the same + ``cg_backend`` as this function. + + :return: :py:class:`TensorMap` containing cartesian components instead of spherical + components. + """ + if len(tensor) == 0: + # nothing to do + return tensor + + if torch_jit_is_scripting(): + use_torch = True + else: + if isinstance(tensor.block(0).values, TorchTensor): + use_torch = True + elif isinstance(tensor.block(0).values, np.ndarray): + use_torch = False + else: + raise TypeError( + f"unknown array type in tensor ({type(tensor.block(0).values)}), " + "only numpy and torch are supported" + ) + + if cg_backend is None: + # TODO: benchmark & change the default? + if use_torch: + cg_backend = "python-dense" + else: + cg_backend = "python-sparse" + + key_names = tensor.keys.names + if "o3_lambda" not in key_names or "o3_sigma" not in key_names: + raise ValueError( + "this tensor is missing one of `o3_lambda` or `o3_sigma` in the keys, " + "is it already in cartesian form?" + ) + + if "o3_mu" not in tensor.component_names: + raise ValueError( + "this tensor does not have an `o3_mu` component, " + "is it already in cartesian form?" + ) + + new_keys = tensor.keys + array_of_ones = _dispatch.int_array_like([1] * len(new_keys), new_keys.values) + + # add back l_1, l_2 keys if they are missing from a rank 1 or 2 tensor + lambda_min = int(_dispatch.min(tensor.keys.column("o3_lambda"))) + lambda_max = int(_dispatch.max(tensor.keys.column("o3_lambda"))) + if lambda_min < 0: + raise ValueError("expected all o3_lambda to be >= 0 in the keys") + elif lambda_max < 1: + raise ValueError("expected o3_lambda to got at least to 1 in the keys") + elif lambda_max == 1: + if "l_1" not in key_names: + new_keys = new_keys.insert(0, "l_1", array_of_ones) + elif lambda_max == 2: + if "l_1" not in key_names: + new_keys = new_keys.insert(0, "l_1", array_of_ones) + if "l_2" not in key_names: + new_keys = new_keys.insert(0, "l_2", array_of_ones) + + # check the l_x and k_x keys + l_keys = [] + for n in range(1, lambda_max + 1): + if f"l_{n}" not in key_names: + raise ValueError(f"expected a l_{n} dimension in the keys") + one = _dispatch.int_array_like([1], new_keys.values) + column = _dispatch.unique(new_keys.column(f"l_{n}")) + if not _dispatch.all(column == one): + raise ValueError(f"expected only 1 as value for the l_{n} dimension") + + if n > 2: + if f"k_{n - 2}" not in key_names: + raise ValueError(f"expected a k_{n - 2} dimension in the keys") + + column = _dispatch.unique(new_keys.column(f"k_{n - 2}")) + if not _dispatch.max(column) < lambda_max: + raise ValueError( + f"expected only all values for k_{n - 2} to be below {lambda_max}" + ) + + l_keys.insert(0, f"k_{n - 2}") + + l_keys.insert(0, f"l_{n}") + + # de-couple the basis back into a product of L=1 spherical harmonics + if cg_coefficients is None: + if torch_jit_is_scripting(): + raise ValueError( + "in TorchScript mode, `cg_coefficients` must be pre-computed " + "and given to this function explicitly" + ) + else: + cg_coefficients = _coefficients.calculate_cg_coefficients( + lambda_max=lambda_max, + cg_backend=cg_backend, + # use_torch=use_torch, + arrays_backend="torch" if use_torch else "numpy", + dtype=tensor.block(0).values.dtype, + device=tensor.block(0).values.device if use_torch else "cpu", + ) + + tensor_rank = lambda_max + + if new_component_names is None: + new_component_names = [f"xyz_{i}" for i in range(1, tensor_rank + 1)] + elif len(new_component_names) != tensor_rank: + raise ValueError( + f"`new_component_names` should have {tensor_rank} entries (the rank of the " + f"tensor), got {len(new_component_names)}" + ) + + # work on a tensor that has the (possibly re-created) `l_{i}` dimensions in its keys + tensor = TensorMap(new_keys, [block.copy() for block in tensor.blocks()]) + + # position of the (single) `o3_mu` component that will be expanded into the + # `tensor_rank` cartesian components + mu_position = tensor.component_names.index("o3_mu") + + # Step 1: de-couple the basis back into a product of L=1 spherical harmonics. + # + # This reverses the iterative coupling done by `cartesian_to_spherical`: at each + # step we split off the left-most L=1 spherical harmonic (`l_{m}`) from the coupled + # `o3_lambda`, leaving the remaining (still coupled) intermediate to the right. The + # intermediate is decoupled further on the next iteration, until only L=1 spherical + # harmonics remain. + component_position = mu_position + m = tensor_rank + while m >= 2: + l_left_name = f"l_{m}" + if m > 2: + lambda_right_name = f"k_{m - 2}" + else: + lambda_right_name = "l_1" + + tensor = _do_decoupling( + tensor=tensor, + component_position=component_position, + l_left_name=l_left_name, + lambda_right_name=lambda_right_name, + cg_coefficients=cg_coefficients, + cg_backend=cg_backend, + ) + component_position += 1 + m -= 1 + + # Step 2: transform the `tensor_rank` L=1 spherical harmonics back into cartesian + # components. Each L=1 component is stored as (m=-1, 0, 1), i.e. (y, z, x), so a + # `roll` by +1 brings it back to (x, y, z). The reconstructed components occupy the + # axes `mu_position, ..., mu_position + tensor_rank - 1`. + new_blocks: List[TensorBlock] = [] + for block in tensor.blocks(): + axes = [mu_position + 1 + i for i in range(tensor_rank)] + values = _dispatch.roll( + block.values, + shifts=[1] * tensor_rank, + axis=axes, + ) + + new_components: List[Labels] = [] + for idx, component in enumerate(block.components): + if mu_position <= idx < mu_position + tensor_rank: + new_components.append( + Labels( + names=[new_component_names[idx - mu_position]], + values=_dispatch.int_array_like( + [[0], [1], [2]], tensor.keys.values + ), + ) + ) + else: + new_components.append(component) + + new_blocks.append( + TensorBlock(values, block.samples, new_components, block.properties) + ) + + # remove the now-unused spherical keys + new_keys = tensor.keys + new_keys = new_keys.remove("o3_lambda") + new_keys = new_keys.remove("o3_sigma") + if "l_1" in new_keys.names: + # only left over for rank 1 tensors (no decoupling happened) + new_keys = new_keys.remove("l_1") + + return TensorMap(new_keys, new_blocks) + + def _do_coupling( tensor: TensorMap, component_1: int, @@ -446,3 +666,167 @@ def _do_coupling( _dispatch.int_array_like(new_keys_values, new_keys.values), ) return TensorMap(new_keys, new_blocks) + + +def _do_decoupling( + tensor: TensorMap, + component_position: int, + l_left_name: str, + lambda_right_name: str, + cg_coefficients: TensorMap, + cg_backend: str, +) -> TensorMap: + """ + Undo a single coupling step performed by :py:func:`_do_coupling`. + + The component at ``component_position`` (which transforms like a single spherical + harmonic ``Y^M_L`` with ``L = o3_lambda``) is split back into a product of two + spherical harmonics: an ``L=1`` term (given by ``l_left_name``, always equal to 1) + at ``component_position`` and an intermediate term (given by ``lambda_right_name``) + at ``component_position + 1``. This is the inverse of coupling those two terms into + ``o3_lambda`` via Clebsch-Gordan coefficients. + + All blocks that share the same keys except for ``o3_lambda`` (and the dependent + ``o3_sigma``) are grouped together and summed through :py:func:`cg_uncouple`, so the + multiple ``o3_lambda`` blocks created by the forward coupling collapse back into a + single block. + + :param tensor: input :py:class:`TensorMap`, in coupled form + :param component_position: index (in ``block.components``) of the component to split + :param l_left_name: name of the key holding the degree of the left-hand spherical + harmonic to recover (``l_{m}``); its value is always 1 + :param lambda_right_name: name of the key holding the degree of the intermediate + (right-hand) spherical harmonic to recover (``k_{m - 2}``, or ``l_1`` for the + last decoupling step) + :param cg_coefficients: pre-computed Clebsch-Gordan coefficients + :param cg_backend: which backend to use for the calculation + + :return: :py:class:`TensorMap` with ``o3_lambda`` set to the recovered intermediate + degree, ``l_left_name`` and ``lambda_right_name`` removed from the keys, and the + single component at ``component_position`` replaced by two components. + """ + keys = tensor.keys + names = keys.names + blocks = tensor.blocks() + + lambda_values = _dispatch.to_int_list(keys.column("o3_lambda")) + sigma_values = _dispatch.to_int_list(keys.column("o3_sigma")) + left_values = _dispatch.to_int_list(keys.column(l_left_name)) + right_values = _dispatch.to_int_list(keys.column(lambda_right_name)) + + # group blocks by all key dimensions except `o3_lambda` and `o3_sigma` + group_names: List[str] = [] + for name in names: + if name != "o3_lambda" and name != "o3_sigma": + group_names.append(name) + group_columns = [_dispatch.to_int_list(keys.column(name)) for name in group_names] + + # output keys: keep everything except the two dimensions we are recovering, and put + # back `o3_lambda` and `o3_sigma` at the front + kept_names: List[str] = [] + for name in group_names: + if name != l_left_name and name != lambda_right_name: + kept_names.append(name) + out_names = ["o3_lambda", "o3_sigma"] + kept_names + + group_signatures: List[List[int]] = [] + group_members: List[List[int]] = [] + for block_idx in range(len(blocks)): + signature = [group_columns[c][block_idx] for c in range(len(group_names))] + found = -1 + for g in range(len(group_signatures)): + if group_signatures[g] == signature: + found = g + break + if found == -1: + group_signatures.append(signature) + group_members.append([block_idx]) + else: + group_members[found].append(block_idx) + + new_keys_values: List[List[int]] = [] + new_blocks: List[TensorBlock] = [] + mu_axis = component_position + 1 + for g in range(len(group_signatures)): + members = group_members[g] + first = members[0] + + l1 = left_values[first] # always 1 + l2 = right_values[first] + + ref_block = blocks[first] + shape = ref_block.values.shape + n_before = 1 + for axis in range(mu_axis): + n_before *= shape[axis] + n_after = 1 + for axis in range(mu_axis + 1, len(shape)): + n_after *= shape[axis] + + o3_lambdas: List[int] = [] + arrays: List[Array] = [] + for block_idx in members: + o3_lambda = lambda_values[block_idx] + o3_lambdas.append(o3_lambda) + arrays.append( + blocks[block_idx].values.reshape(n_before, 2 * o3_lambda + 1, n_after) + ) + + uncoupled = _coefficients.cg_uncouple( + arrays, l1, l2, o3_lambdas, cg_coefficients, cg_backend + ) + + new_shape = ( + list(shape[:mu_axis]) + [2 * l1 + 1, 2 * l2 + 1] + list(shape[mu_axis + 1:]) + ) + values = uncoupled.reshape(new_shape) + + new_components: List[Labels] = [] + for idx in range(len(ref_block.components)): + if idx == component_position: + new_components.append( + Labels( + names=[f"__sph_{component_position}"], + values=_dispatch.int_array_like( + [[mu] for mu in range(-l1, l1 + 1)], keys.values + ), + ) + ) + new_components.append( + Labels( + names=[f"__sph_{component_position + 1}"], + values=_dispatch.int_array_like( + [[mu] for mu in range(-l2, l2 + 1)], keys.values + ), + ) + ) + else: + new_components.append(ref_block.components[idx]) + + new_blocks.append( + TensorBlock( + values=values, + samples=ref_block.samples, + components=new_components, + properties=ref_block.properties, + ) + ) + + # recover the `o3_sigma` of the intermediate. The forward coupling set + # `o3_sigma = old_sigma * (-1)^(l1 + l2 + o3_lambda)`, and `(-1)^x` is its own + # inverse, so we recover `old_sigma` the same way. + old_sigma = sigma_values[first] * ((-1) ** (l1 + l2 + o3_lambdas[0])) + + signature = group_signatures[g] + kept_values: List[int] = [] + for c in range(len(group_names)): + if group_names[c] != l_left_name and group_names[c] != lambda_right_name: + kept_values.append(signature[c]) + + new_keys_values.append([l2, old_sigma] + kept_values) + + new_keys = Labels( + out_names, + _dispatch.int_array_like(new_keys_values, keys.values), + ) + return TensorMap(new_keys, new_blocks) \ No newline at end of file diff --git a/python/featomic/featomic/clebsch_gordan/_coefficients.py b/python/featomic/featomic/clebsch_gordan/_coefficients.py index 85d0c3b85..1a6465f50 100644 --- a/python/featomic/featomic/clebsch_gordan/_coefficients.py +++ b/python/featomic/featomic/clebsch_gordan/_coefficients.py @@ -582,6 +582,103 @@ def _cg_couple_dense( return _dispatch.tensordot(array, cg_l1l2lam[0, ..., 0], axes=([2, 1], [1, 0])) +def cg_uncouple( + arrays: List[Array], + l1: int, + l2: int, + o3_lambdas: List[int], + cg_coefficients: TensorMap, + cg_backend: str, +) -> Array: + """ + Exact inverse of :py:func:`cg_couple`. + + Go from a coupled basis that behaves like a single spherical harmonic back to the + uncoupled product basis that behaves like a product of two spherical harmonics of + degrees ``l1`` and ``l2``. + + The coupling performed by :py:func:`cg_couple` is the orthogonal transformation + + ``T^{L}_{M} = \\sum_{m1 m2} (-1)^{l1 + l2 + L} T^{l1}_{m1} + T^{l2}_{m2}`` + + Because the (real) Clebsch-Gordan coefficients are orthogonal and complete over the + full set of ``L`` in ``range(|l1 - l2|, l1 + l2 + 1)``, this can be inverted exactly + by applying the transpose: + + ``T^{l1}_{m1} T^{l2}_{m2} = \\sum_{L M} (-1)^{l1 + l2 + L} + T^{L}_{M}`` + + The same sign convention and the same cached coefficients as :py:func:`cg_couple` + are used, so ``cg_uncouple(cg_couple(...))`` recovers the original array (summing + over whichever ``o3_lambdas`` are provided). + + :param arrays: list of arrays, one for each entry in ``o3_lambdas``, each with shape + ``[n_s, 2 * o3_lambda + 1, n_q]``. ``n_s`` and ``n_q`` must be the same for all + the arrays. + :param l1: degree of the first spherical harmonic to recover + :param l2: degree of the second spherical harmonic to recover + :param o3_lambdas: degrees of the coupled spherical harmonics given in ``arrays``, + in the same order as ``arrays``. For an exact inversion this should contain all + of ``range(|l1 - l2|, l1 + l2 + 1)``; any missing degree is treated as a block of + zeros. + :param cg_coefficients: CG coefficients as returned by + :py:func:`calculate_cg_coefficients` with the same ``cg_backend`` given to this + function + :param cg_backend: ``"python-dense"`` or ``"python-sparse"`` + + :return: array with shape ``[n_s, 2 * l1 + 1, 2 * l2 + 1, n_q]`` + """ + assert len(arrays) == len(o3_lambdas) + assert len(arrays) > 0 + + n_s = arrays[0].shape[0] + n_q = arrays[0].shape[2] + + output = _dispatch.zeros_like(arrays[0], (n_s, 2 * l1 + 1, 2 * l2 + 1, n_q)) + + if cg_backend == "python-sparse": + for o3_lambda, array in zip(o3_lambdas, arrays): + cg_l1l2lam = cg_coefficients.block( + {"l1": l1, "l2": l2, "lambda": o3_lambda} + ) + sign = (-1) ** (l1 + l2 + o3_lambda) + samples = cg_l1l2lam.samples + for i in range(len(samples)): + m1m2mu = samples.entry(i) + m1 = int(m1m2mu[0]) + m2 = int(m1m2mu[1]) + mu = int(m1m2mu[2]) + output[:, m1, m2, :] += ( + array[:, mu, :] * sign * cg_l1l2lam.values[i, 0] + ) + + return output + + elif cg_backend == "python-dense": + for o3_lambda, array in zip(o3_lambdas, arrays): + # cg block has shape (1, 2*l1 + 1, 2*l2 + 1, 2*o3_lambda + 1, 1) + cg_l1l2lam = (-1) ** (l1 + l2 + o3_lambda) * cg_coefficients.block( + {"l1": l1, "l2": l2, "lambda": o3_lambda} + ).values + # contract the `mu` axis of `array` with the `mu` axis of the coefficients + # array (n_s, 2*o3_lambda + 1, n_q) . (2*l1+1, 2*l2+1, 2*o3_lambda+1) + # -> (n_s, n_q, 2*l1 + 1, 2*l2 + 1) + contribution = _dispatch.tensordot( + array, cg_l1l2lam[0, ..., 0], axes=([1], [2]) + ) + # -> (n_s, 2*l1 + 1, 2*l2 + 1, n_q) + output = output + _dispatch.permute(contribution, [0, 2, 3, 1]) + + return output + + else: + raise ValueError( + f"invalid `cg_backend`, got '{cg_backend}', " + "only 'python-dense', or 'python-sparse' are supported" + ) + + # ======================================================================= # # =============== Functions for performing CG tensor products =========== # # ======================================================================= # diff --git a/python/featomic/featomic/clebsch_gordan/_dispatch.py b/python/featomic/featomic/clebsch_gordan/_dispatch.py index 2bc1e8d7d..6c4253b6c 100644 --- a/python/featomic/featomic/clebsch_gordan/_dispatch.py +++ b/python/featomic/featomic/clebsch_gordan/_dispatch.py @@ -383,6 +383,20 @@ def max(array): raise TypeError(UNKNOWN_ARRAY_TYPE) +def min(array): + """ + Takes the smallest value of the array. + This function has the same behavior as + ``np.min(array)`` or ``torch.min(array)``. + """ + if isinstance(array, TorchTensor): + return torch.min(array) + elif isinstance(array, np.ndarray): + return np.min(array) + else: + raise TypeError(UNKNOWN_ARRAY_TYPE) + + def any(array): """Test whether any array elements along a given axis evaluate to True. diff --git a/python/featomic/tests/clebsch_gordan/cartesian_spherical.py b/python/featomic/tests/clebsch_gordan/cartesian_spherical.py index 16d805d74..b0ba57f4a 100644 --- a/python/featomic/tests/clebsch_gordan/cartesian_spherical.py +++ b/python/featomic/tests/clebsch_gordan/cartesian_spherical.py @@ -3,7 +3,7 @@ import pytest from metatensor import Labels, TensorBlock, TensorMap -from featomic.clebsch_gordan import cartesian_to_spherical +from featomic.clebsch_gordan import cartesian_to_spherical, spherical_to_cartesian from .rotations import WignerDReal, cartesian_rotation @@ -146,18 +146,70 @@ def test_cartesian_to_spherical_rank_3(cartesian): @pytest.mark.parametrize("cg_backend", ["python-dense", "python-sparse"]) @pytest.mark.parametrize( - "components", [["xyz_1"], ["xyz_1", "xyz_12"], ["xyz_1", "xyz_2", "xyz_3"]] + "components", [["xyz_1"], ["xyz_1", "xyz_2"], ["xyz_1", "xyz_2", "xyz_3"]] ) def test_cartesian_to_spherical_and_back(cartesian, components, cg_backend): spherical = cartesian_to_spherical( cartesian, - components=["xyz_1", "xyz_2", "xyz_3"], + components=components, keep_l_in_keys=True, + # keep every block so the transformation stays exactly invertible + remove_blocks_threshold=None, cg_backend=cg_backend, ) assert "o3_lambda" in spherical.keys.names - # TODO: check for identity after spherical_to_cartesian + + # going back to cartesian should recover the original tensor exactly + back = spherical_to_cartesian(spherical, cg_backend=cg_backend) + + assert mts.equal_metadata(back, cartesian) + assert mts.allclose(back, cartesian) + + +@pytest.mark.parametrize("cg_backend", ["python-dense", "python-sparse"]) +@pytest.mark.parametrize("rank", [1, 2, 3, 4]) +def test_spherical_to_cartesian_roundtrip(rank, cg_backend): + """ + Self-contained round-trip test: a random cartesian tensor of the given ``rank`` is + transformed to spherical form and back, and must be recovered exactly. + """ + rng = np.random.default_rng(0x5C1E) + + n_samples, n_properties = 6, 3 + shape = (n_samples,) + (3,) * rank + (n_properties,) + components = [Labels.range(f"xyz_{i}", 3) for i in range(1, rank + 1)] + + cartesian = TensorMap( + keys=Labels.single(), + blocks=[ + TensorBlock( + values=rng.random(shape), + samples=Labels.range("system", n_samples), + components=components, + properties=Labels.range("p", n_properties), + ) + ], + ) + + component_names = [f"xyz_{i}" for i in range(1, rank + 1)] + spherical = cartesian_to_spherical( + cartesian, + components=component_names, + keep_l_in_keys=True, + remove_blocks_threshold=None, + cg_backend=cg_backend, + ) + + # the spherical representation should no longer contain cartesian components + assert "o3_mu" in spherical.component_names + for name in component_names: + assert name not in spherical.component_names + + back = spherical_to_cartesian(spherical, cg_backend=cg_backend) + + assert mts.equal_metadata(back, cartesian) + assert mts.allclose(back, cartesian) def test_cartesian_to_spherical_errors(cartesian):