Add traits for matrix interfaces

A lot of external packages are happy if operators and states are
abstract arrays. We can make that part of the interface.

Author: goerz

CLAUDE.md

@@ -5,21 +5,26 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 ## Development Commands
 
 ### Testing
+
 - `make test` - Run the full test suite
 - `julia --project=test -e 'include("test/runtests.jl")'` - Run tests directly
 
 ### Development Environment
+
 - `make devrepl` - Start interactive REPL with development environment (recommended for development)
 - `julia -i --banner=no devrepl.jl` - Alternative way to start development REPL
 
 ### Documentation
+
 - `make docs` - Build documentation
 
 ### Code Quality
+
 - `make codestyle` - Apply JuliaFormatter to entire project
 - `julia --project=test -e 'using JuliaFormatter; format(".", verbose=true)'` - Format code directly
 
 ### Cleanup
+
 - `make clean` - Clean build/doc/testing artifacts
 - `make distclean` - Restore to clean checkout state
 
@@ -28,6 +33,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 This is a Julia package for quantum system time propagation within the JuliaQuantumControl ecosystem.
 
 ### Core Structure
+
 - **src/QuantumPropagators.jl** - Main module definition with submodule includes
 - **src/interfaces/** - Interface definitions for operators, states, generators, propagators, etc.
 - **Submodules** organized by functionality:
@@ -43,27 +49,36 @@ This is a Julia package for quantum system time propagation within the JuliaQuan
   - `Generators` (generators.jl) - Hamiltonian/Liouvillian generators
 
 ### Key Propagation Methods
+
 - **Chebyshev propagation** (cheby_propagator.jl) - Polynomial expansion method
 - **Newton propagation** (newton_propagator.jl) - Newton interpolation method
 - **Matrix exponential** (exp_propagator.jl) - Direct exponentiation
 - **ODE integration** (ode_function.jl) - Interface to OrdinaryDiffEq.jl via extensions
 
 ### High-Level Interface
+
 - `propagate()` function for single time evolution
 - `propagate_sequence()` for sequential propagations with different parameters
 
 ### Extensions
+
 The package uses Julia's extension system for optional dependencies:
 - `QuantumPropagatorsODEExt` - OrdinaryDiffEq.jl integration
 - `QuantumPropagatorsRecursiveArrayToolsExt` - RecursiveArrayTools.jl support
 - `QuantumPropagatorsStaticArraysExt` - StaticArrays.jl optimization
 
 ### Testing Framework
+
 Uses SafeTestsets for isolated test execution. Tests are comprehensive and include:
 - Interface validation tests
 - Individual propagator method tests
 - Integration tests for complete propagation workflows
 
+Running `make test` prints out coverage data in a table.
+
+If necessary, detailed line-by-line coverage information can be obtained by running julia --project=test -e 'include("devrepl.jl"); generate_coverage_html()' after `make test`.
+This will produce html files inside the `coverage` subfolder, with `coverage/src` mirroring the structure of the `src` folder of `.jl` files. Lines with `<span class="tlaUNC">` are not covered. Ignore the raw tracefiles in the `.coverage` subfolder.
+
 ### Development Environment
 The project uses a sophisticated development setup:
 - Development REPL (devrepl.jl) with Revise.jl for hot reloading
@@ -72,6 +87,7 @@ The project uses a sophisticated development setup:
 - Code formatting with JuliaFormatter
 
 ### Documentation
+
 Uses Documenter.jl with comprehensive API documentation and examples. Documentation includes detailed method explanations.
 
 ## Docstrings

src/QuantumPropagators.jl

@@ -38,11 +38,13 @@ export init_prop, reinit_prop!, prop_step!
 
 #! format: off
 module Interfaces
-    export supports_inplace
+    export supports_inplace, supports_vector_interface, supports_matrix_interface
     export check_operator, check_state, check_tlist, check_amplitude
     export check_control, check_generator, check_propagator
     export check_parameterized_function, check_parameterized
     include("interfaces/supports_inplace.jl")
+    include("interfaces/supports_vector_interface.jl")
+    include("interfaces/supports_matrix_interface.jl")
     include("interfaces/utils.jl")
     include("interfaces/state.jl")
     include("interfaces/tlist.jl")

src/interfaces/operator.jl

@@ -17,7 +17,6 @@ verifies the given `op` relative to `state`. The `state` must pass
 An "operator" is any object that [`evaluate`](@ref) returns when evaluating a
 time-dependent dynamic generator. The specific requirements for `op` are:
 
-* `eltype(op)` must be defined and return a numeric type
 * `size(op)` must be defined and return a tuple of integers
 * `size(op, dim)` must be defined for each dimension and be consistent with
   `size(op)`
@@ -26,9 +25,11 @@ time-dependent dynamic generator. The specific requirements for `op` are:
 * `op * state` must be defined
 * The [`QuantumPropagators.Interfaces.supports_inplace`](@ref) method must be
   defined for `op`. If it returns `true`, it must be possible to evaluate a
-  generator in-place into the existing `op`. See [`check_generator`](@ref).
+  generator in-place into the existing `op`. See
+  [`QuantumPropagators.Interfaces.check_generator`](@ref).
 
-If [`QuantumPropagators.Interfaces.supports_inplace(state)`](@ref):
+If [`QuantumPropagators.Interfaces.supports_inplace(state)`](@ref
+QuantumPropagators.Interfaces.supports_inplace):
 
 * The 3-argument `LinearAlgebra.mul!` must apply `op` to the given `state`
 * The 5-argument `LinearAlgebra.mul!` must apply `op` to the given `state`
@@ -40,6 +41,30 @@ If `for_expval` (typically required for optimal control):
 * `LinearAlgebra.dot(state, op, state)` must return return a number
 * `dot(state, op, state)` must match `dot(state, op * state)`, if applicable
 
+If [`QuantumPropagators.Interfaces.supports_matrix_interface(op)`](@ref
+QuantumPropagators.Interfaces.supports_matrix_interface) is `true`, the
+operator must implement the
+[Abstract Array interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
+for two-dimensional arrays:
+
+* `eltype(op)` must be defined and return a numeric type
+* `getindex(op, i, j)` must be defined and return elements matching `eltype`
+* `length(op)` must equal `prod(size(op))`
+* `iterate(op)` must be defined
+* `similar(op)` must be defined and return a mutable array with the same shape
+  and element type
+* `similar(op, ::Type{S})` must return a mutable array with the same shape and
+  element type `S`
+* `similar(op, dims::Dims)` must return a mutable array with the same element
+  type and the given dimensions
+* `similar(op, ::Type{S}, dims::Dims)` must return a mutable array with the
+  given element type and dimensions
+
+If additionally [`QuantumPropagators.Interfaces.supports_inplace(op)`](@ref
+QuantumPropagators.Interfaces.supports_inplace) is `true` (read-write matrix):
+
+* `setindex!(op, v, i, j)` must be defined
+
 The function returns `true` for a valid operator and `false` for an invalid
 operator. Unless `quiet=true`, it will log an error to indicate which of the
 conditions failed.
@@ -74,20 +99,6 @@ function check_operator(
         success = false
     end
 
-    try
-        T = eltype(op)
-        if !(T isa Type && T <: Number)
-            quiet || @error "$(px)`eltype(op)` must return a numeric type, not $T"
-            success = false
-        end
-    catch exc
-        quiet || @error(
-            "$(px)`eltype(op)` must be defined.",
-            exception = (exc, catch_abbreviated_backtrace())
-        )
-        success = false
-    end
-
     try
         s = size(op)
         if !(s isa Tuple)
@@ -171,7 +182,6 @@ function check_operator(
         success = false
     end
 
-
     if supports_inplace(state)
 
         try
@@ -254,6 +264,197 @@ function check_operator(
 
     end
 
+    if supports_matrix_interface(op)
+
+        try
+            T = eltype(op)
+            if !(T isa Type && T <: Number)
+                quiet || @error "$(px)`eltype(op)` must return a numeric type, not $T"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`eltype(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            s = size(op)
+            if length(s) == 2 && all(d -> d > 0, s)
+                val = op[1, 1]
+                T = eltype(op)
+                if T isa Type && T <: Number && !(val isa T)
+                    quiet ||
+                        @error "$(px)`op[1, 1]` must return a value of type `eltype(op)=$T`, not $(typeof(val))"
+                    success = false
+                end
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`getindex(op, i, j)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            l = length(op)
+            s = size(op)
+            if l != prod(s)
+                quiet ||
+                    @error "$(px)`length(op)` must equal `prod(size(op))`: $l ≠ $(prod(s))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`length(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            itr = iterate(op)
+            s = size(op)
+            if isnothing(itr) && prod(s) > 0
+                quiet ||
+                    @error "$(px)`iterate(op)` must not return `nothing` for a non-empty operator"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`iterate(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            op2 = similar(op)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != size(op)
+                quiet ||
+                    @error "$(px)`similar(op)` must return an array with the same shape: size $(size(op2)) ≠ $(size(op))"
+                success = false
+            end
+            if eltype(op2) != eltype(op)
+                quiet ||
+                    @error "$(px)`similar(op)` must return an array with the same element type: $(eltype(op2)) ≠ $(eltype(op))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(op) == ComplexF64) ? ComplexF32 : ComplexF64
+            op2 = similar(op, S)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != size(op)
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return an array with the same shape: size $(size(op2)) ≠ $(size(op))"
+                success = false
+            end
+            if eltype(op2) != S
+                quiet ||
+                    @error "$(px)`similar(op, $S)` must return an array with element type $S, got $(eltype(op2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, ::Type{S})` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            dims = size(op)
+            op2 = similar(op, dims)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != dims
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return an array with size $dims, got $(size(op2))"
+                success = false
+            end
+            if eltype(op2) != eltype(op)
+                quiet ||
+                    @error "$(px)`similar(op, dims)` must return an array with the same element type: $(eltype(op2)) ≠ $(eltype(op))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        try
+            S = (eltype(op) == ComplexF64) ? ComplexF32 : ComplexF64
+            dims = size(op)
+            op2 = similar(op, S, dims)
+            if !supports_inplace(op2)
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return a mutable array (`supports_inplace` must be `true`), got $(typeof(op2))"
+                success = false
+            end
+            if size(op2) != dims
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return an array with size $dims, got $(size(op2))"
+                success = false
+            end
+            if eltype(op2) != S
+                quiet ||
+                    @error "$(px)`similar(op, $S, dims)` must return an array with element type $S, got $(eltype(op2))"
+                success = false
+            end
+        catch exc
+            quiet || @error(
+                "$(px)`similar(op, ::Type{S}, dims::Dims)` must be defined.",
+                exception = (exc, catch_abbreviated_backtrace())
+            )
+            success = false
+        end
+
+        if supports_inplace(op)
+
+            try
+                s = size(op)
+                if length(s) == 2 && all(d -> d > 0, s)
+                    v = op[1, 1]
+                    op[1, 1] = v  # write back the same value
+                end
+            catch exc
+                quiet || @error(
+                    "$(px)`setindex!(op, v, i, j)` must be defined for an in-place operator.",
+                    exception = (exc, catch_abbreviated_backtrace())
+                )
+                success = false
+            end
+
+        end
+
+    end
+
     return success
 
 end