PR: Implement support for *Python Array API Standard*.

.github/workflows/continuous-integration-quality-unit-tests.yml

@@ -47,6 +47,11 @@ jobs:
           uv sync --all-extras
           uv run python -c "import imageio;imageio.plugins.freeimage.download()"
         shell: bash
+      - name: Install Array API Backends
+        if: matrix.os == 'macOS-latest' && matrix.python-version == '3.13'
+        run: |
+          uv pip install jax jaxlib torch
+        shell: bash
       - name: Prek (All Files)
         if: matrix.os == 'macOS-latest'
         run: |

colour/__init__.py

@@ -190,6 +190,7 @@
     TVS_ILLUMINANTS_HUNTERLAB,
     WHITENESS_METHODS,
     YELLOWNESS_METHODS,
+    CIE_illuminant_D_series,
     MultiSpectralDistributions,
     SpectralDistribution,
     SpectralShape,
@@ -203,8 +204,11 @@
     luminous_efficacy,
     luminous_efficiency,
     luminous_flux,
+    msds_blackbody,
+    msds_CIE_illuminant_D_series,
     msds_constant,
     msds_ones,
+    msds_rayleigh_jeans,
     msds_to_XYZ,
     msds_zeros,
     sd_blackbody,
@@ -615,6 +619,7 @@
     "TVS_ILLUMINANTS_HUNTERLAB",
     "WHITENESS_METHODS",
     "YELLOWNESS_METHODS",
+    "CIE_illuminant_D_series",
     "MultiSpectralDistributions",
     "SpectralDistribution",
     "SpectralShape",
@@ -628,8 +633,11 @@
     "luminous_efficacy",
     "luminous_efficiency",
     "luminous_flux",
+    "msds_blackbody",
+    "msds_CIE_illuminant_D_series",
     "msds_constant",
     "msds_ones",
+    "msds_rayleigh_jeans",
     "msds_to_XYZ",
     "msds_zeros",
     "sd_blackbody",

colour/adaptation/cie1994.py

@@ -34,12 +34,14 @@
     Range100,
 )
 from colour.utilities import (
+    array_namespace,
     as_float_array,
     from_range_100,
     to_domain_100,
     tsplit,
     tstack,
     usage_warning,
+    xp_as_float_array,
 )
 
 __author__ = "Colour Developers"
@@ -151,10 +153,14 @@ def chromatic_adaptation_CIE1994(
 
     XYZ_1 = to_domain_100(XYZ_1)
     Y_o = as_float_array(to_domain_100(Y_o))
-    E_o1 = as_float_array(E_o1)
-    E_o2 = as_float_array(E_o2)
 
-    if np.any(Y_o < 18) or np.any(Y_o > 100):
+    xp = array_namespace(XYZ_1, Y_o, E_o1, E_o2)
+
+    Y_o = xp_as_float_array(Y_o, xp=xp, like=XYZ_1)
+    E_o1 = xp_as_float_array(E_o1, xp=xp, like=XYZ_1)
+    E_o2 = xp_as_float_array(E_o2, xp=xp, like=XYZ_1)
+
+    if xp.any(Y_o < 18) or xp.any(Y_o > 100):
         usage_warning(
             '"Y_o" luminance factor must be in [18, 100] domain, '
             "unpredictable results may occur!"
@@ -295,10 +301,13 @@ def effective_adapting_responses(
     """
 
     xez = as_float_array(xez)
-    Y_o = as_float_array(Y_o)
-    E_o = as_float_array(E_o)
 
-    return ((Y_o[..., None] * E_o[..., None]) / (100 * np.pi)) * xez
+    xp = array_namespace(xez, Y_o, E_o)
+
+    Y_o = xp_as_float_array(Y_o, xp=xp, like=xez)
+    E_o = xp_as_float_array(E_o, xp=xp, like=xez)
+
+    return ((Y_o[..., None] * E_o[..., None]) / (100 * xp.pi)) * xez
 
 
 @typing.overload
@@ -447,12 +456,15 @@ def K_coefficient(
     np.float64(1.0)
     """
 
+    xp = array_namespace(xez_1, xez_2, bRGB_o1, bRGB_o2, Y_o, n)
+
+    Y_o = xp_as_float_array(Y_o, xp=xp, like=xez_1)
+    n = xp_as_float_array(n, xp=xp, like=xez_1)
+
     xi_1, eta_1, _zeta_1 = tsplit(xez_1)
     xi_2, eta_2, _zeta_2 = tsplit(xez_2)
     bR_o1, bG_o1, _bB_o1 = tsplit(bRGB_o1)
     bR_o2, bG_o2, _bB_o2 = tsplit(bRGB_o2)
-    Y_o = as_float_array(Y_o)
-    n = as_float_array(n)
 
     return (
         spow((Y_o * xi_1 + n) / (20 * xi_1 + n), (2 / 3) * bR_o1)
@@ -523,14 +535,17 @@ def corresponding_colour(
     array([23.1636901..., 20.0211948..., 16.2001664...])
     """
 
+    xp = array_namespace(RGB_1, xez_1, xez_2, bRGB_o1, bRGB_o2, Y_o, K, n)
+
+    Y_o = xp_as_float_array(Y_o, xp=xp, like=RGB_1)
+    K = xp_as_float_array(K, xp=xp, like=RGB_1)
+    n = xp_as_float_array(n, xp=xp, like=RGB_1)
+
     R_1, G_1, B_1 = tsplit(RGB_1)
     xi_1, eta_1, zeta_1 = tsplit(xez_1)
     xi_2, eta_2, zeta_2 = tsplit(xez_2)
     bR_o1, bG_o1, bB_o1 = tsplit(bRGB_o1)
     bR_o2, bG_o2, bB_o2 = tsplit(bRGB_o2)
-    Y_o = as_float_array(Y_o)
-    K = as_float_array(K)
-    n = as_float_array(n)
 
     def RGB_c(
         x_1: NDArrayFloat,

colour/adaptation/cmccat2000.py

@@ -43,10 +43,11 @@
 from colour.utilities import (
     CanonicalMapping,
     MixinDataclassIterable,
-    as_float_array,
+    array_namespace,
     from_range_100,
     to_domain_100,
     validate_method,
+    xp_as_float_array,
 )
 
 __author__ = "Colour Developers"
@@ -193,20 +194,23 @@ def chromatic_adaptation_forward_CMCCAT2000(
     XYZ = to_domain_100(XYZ)
     XYZ_w = to_domain_100(XYZ_w)
     XYZ_wr = to_domain_100(XYZ_wr)
-    L_A1 = as_float_array(L_A1)
-    L_A2 = as_float_array(L_A2)
+
+    xp = array_namespace(XYZ, XYZ_w, XYZ_wr, L_A1, L_A2)
+
+    L_A1 = xp_as_float_array(L_A1, xp=xp, like=XYZ)
+    L_A2 = xp_as_float_array(L_A2, xp=xp, like=XYZ)
 
     RGB = vecmul(CAT_CMCCAT2000, XYZ)
     RGB_w = vecmul(CAT_CMCCAT2000, XYZ_w)
     RGB_wr = vecmul(CAT_CMCCAT2000, XYZ_wr)
 
     D = surround.F * (
-        0.08 * np.log10(0.5 * (L_A1 + L_A2))
+        0.08 * xp.log10(0.5 * (L_A1 + L_A2))
         + 0.76
         - 0.45 * (L_A1 - L_A2) / (L_A1 + L_A2)
     )
 
-    D = np.clip(D, 0, 1)
+    D = xp.clip(D, 0, 1)
     a = D * XYZ_w[..., 1] / XYZ_wr[..., 1]
 
     RGB_c = RGB * (a[..., None] * (RGB_wr / RGB_w) + 1 - D[..., None])
@@ -288,20 +292,23 @@ def chromatic_adaptation_inverse_CMCCAT2000(
     XYZ_c = to_domain_100(XYZ_c)
     XYZ_w = to_domain_100(XYZ_w)
     XYZ_wr = to_domain_100(XYZ_wr)
-    L_A1 = as_float_array(L_A1)
-    L_A2 = as_float_array(L_A2)
+
+    xp = array_namespace(XYZ_c, XYZ_w, XYZ_wr, L_A1, L_A2)
+
+    L_A1 = xp_as_float_array(L_A1, xp=xp, like=XYZ_c)
+    L_A2 = xp_as_float_array(L_A2, xp=xp, like=XYZ_c)
 
     RGB_c = vecmul(CAT_CMCCAT2000, XYZ_c)
     RGB_w = vecmul(CAT_CMCCAT2000, XYZ_w)
     RGB_wr = vecmul(CAT_CMCCAT2000, XYZ_wr)
 
     D = surround.F * (
-        0.08 * np.log10(0.5 * (L_A1 + L_A2))
+        0.08 * xp.log10(0.5 * (L_A1 + L_A2))
         + 0.76
         - 0.45 * (L_A1 - L_A2) / (L_A1 + L_A2)
     )
 
-    D = np.clip(D, 0, 1)
+    D = xp.clip(D, 0, 1)
     a = D * XYZ_w[..., 1] / XYZ_wr[..., 1]
 
     RGB = RGB_c / (a[..., None] * (RGB_wr / RGB_w) + 1 - D[..., None])

colour/adaptation/datasets/cat.py

@@ -276,7 +276,7 @@
 :cite:`Nayatani1995a`
 """
 
-CAT_XYZ_SCALING: NDArrayFloat = np.reshape(np.array(np.identity(3)), (3, 3))
+CAT_XYZ_SCALING: NDArrayFloat = np.reshape(np.array(np.eye(3)), (3, 3))
 """
 *XYZ Scaling* chromatic adaptation transform.
 

colour/adaptation/fairchild1990.py

@@ -29,12 +29,12 @@
     Range100,
 )
 from colour.utilities import (
+    array_namespace,
     as_float_array,
     from_range_100,
     ones,
-    row_as_diagonal,
     to_domain_100,
-    tstack,
+    xp_as_float_array,
 )
 
 __author__ = "Colour Developers"
@@ -134,7 +134,10 @@ def chromatic_adaptation_Fairchild1990(
     XYZ_1 = to_domain_100(XYZ_1)
     XYZ_n = to_domain_100(XYZ_n)
     XYZ_r = to_domain_100(XYZ_r)
-    Y_n = as_float_array(Y_n)
+
+    xp = array_namespace(XYZ_1, XYZ_n, XYZ_r, Y_n)
+
+    Y_n = xp_as_float_array(Y_n, xp=xp, like=XYZ_1)
 
     LMS_1 = XYZ_to_RGB_Fairchild1990(XYZ_1)
     LMS_n = XYZ_to_RGB_Fairchild1990(XYZ_n)
@@ -145,18 +148,14 @@ def chromatic_adaptation_Fairchild1990(
     a_LMS_1 = p_LMS / LMS_n
     a_LMS_2 = p_LMS / LMS_r
 
-    A_1 = row_as_diagonal(a_LMS_1)
-    A_2 = row_as_diagonal(a_LMS_2)
-
-    LMSp_1 = vecmul(A_1, LMS_1)
-
-    c = 0.219 - 0.0784 * np.log10(Y_n)
-    C = row_as_diagonal(tstack([c, c, c]))
+    # Diagonal matrix products collapse to elementwise scaling; their inverses
+    # are reciprocals. Avoids materialising ``(N, 3, 3)`` tensors at scale.
+    # The ``c`` matrix in *Fairchild (1990)* Eq. 21-22 cancels exactly between
+    # the forward ``LMS_a = c * LMSp_1`` and inverse ``LMSp_2 = LMS_a / c``
+    # steps, so the intermediate is elided.
+    LMSp_1 = a_LMS_1 * LMS_1
 
-    LMS_a = vecmul(C, LMSp_1)
-    LMSp_2 = vecmul(np.linalg.inv(C), LMS_a)
-
-    LMS_c = vecmul(np.linalg.inv(A_2), LMSp_2)
+    LMS_c = LMSp_1 / a_LMS_2
     XYZ_c = RGB_to_XYZ_Fairchild1990(LMS_c)
 
     return from_range_100(XYZ_c)
@@ -249,21 +248,25 @@ def degrees_of_adaptation(
 
     LMS = as_float_array(LMS)
 
+    xp = array_namespace(LMS, Y_n, v)
+
     if discount_illuminant:
-        return ones(LMS.shape)
+        return xp_as_float_array(ones(LMS.shape), xp=xp, like=LMS)
 
-    Y_n = as_float_array(Y_n)
-    v = as_float_array(v)
+    Y_n = xp_as_float_array(Y_n, xp=xp, like=LMS)
+    v = xp_as_float_array(v, xp=xp, like=LMS)
 
     # E illuminant.
     LMS_E = vecmul(CAT_VON_KRIES, ones(LMS.shape))
 
+    LMS_E = xp_as_float_array(LMS_E, xp=xp, like=LMS)
+
     Ye_n = spow(Y_n, v)
 
     def m_E(x: NDArrayFloat, y: NDArrayFloat) -> NDArrayFloat:
         """Compute the :math:`m_E` term."""
 
-        return (3 * x / y) / np.sum(x / y, axis=-1)[..., None]
+        return (3 * x / y) / xp.sum(x / y, axis=-1)[..., None]
 
     def P_c(x: NDArrayFloat) -> NDArrayFloat:
         """Compute the :math:`P_L`, :math:`P_M` or :math:`P_S` terms."""

colour/adaptation/fairchild2020.py

@@ -38,13 +38,15 @@
 from colour.utilities import (
     CanonicalMapping,
     MixinDataclassIterable,
+    array_namespace,
     as_float_array,
     from_range_1,
     get_domain_range_scale,
     optional,
     row_as_diagonal,
     to_domain_1,
     validate_method,
+    xp_as_float_array,
 )
 
 __author__ = "Colour Developers"
@@ -214,13 +216,15 @@ def matrix_chromatic_adaptation_vk20(
     XYZ_r = as_float_array(XYZ_r)
     XYZ_p = as_float_array(XYZ_p)
 
+    xp = array_namespace(XYZ_n, XYZ_r, XYZ_p)
+
     transform = validate_method(
         transform,
         tuple(CHROMATIC_ADAPTATION_TRANSFORMS),
         '"{0}" chromatic adaptation transform is invalid, it must be one of {1}!',
     )
 
-    M = CHROMATIC_ADAPTATION_TRANSFORMS[transform]
+    M = xp_as_float_array(CHROMATIC_ADAPTATION_TRANSFORMS[transform], xp=xp, like=XYZ_n)
 
     D_n, D_r, D_p = coefficients.values
 
@@ -231,9 +235,9 @@ def matrix_chromatic_adaptation_vk20(
     with sdiv_mode():
         D = row_as_diagonal(sdiv(1, (D_n * LMS_n + D_r * LMS_r + D_p * LMS_p)))
 
-    M_CAT = np.matmul(np.linalg.inv(M), D)
+    M_CAT = xp.matmul(xp.linalg.inv(M), D)
 
-    return np.matmul(M_CAT, M)
+    return xp.matmul(M_CAT, M)
 
 
 def chromatic_adaptation_vK20(
@@ -344,6 +348,10 @@ def chromatic_adaptation_vK20(
         )
     )
 
+    xp = array_namespace(XYZ, XYZ_p, XYZ_n, XYZ_r)
+
+    XYZ_r = xp_as_float_array(XYZ_r, xp=xp, like=XYZ)
+
     M_CAT = matrix_chromatic_adaptation_vk20(
         XYZ_p, XYZ_n, XYZ_r, transform, coefficients
     )