pymor.operators

Extended functionality for pyMOR operators

Imports

import xarray as xr
from IPython.display import HTML
# HTML('''<style>
# html,body,address,blockquote,div,form,fieldset,caption,h1,h2,h3,h4,h5,h6,hr,ul,ol,ul,li,table,tr,td,th,p,img {margin:0; padding:0;}
# ul, ol, li {list-style-type: none;}
# img,fieldset{border: none;}
# ul {margin-left: -2em}
# </style>''')
# HTML('''foo<ul><li>foo</li></ul>foo''')

operator_html_tree


def operator_html_tree(
    op:Operator, opened:bool=False
)->str:

op = LincombOperator(
    [
        NumpyMatrixOperator(np.ones([1])), 
        LincombOperator([NumpyMatrixOperator(np.ones([1]))] * 2, [.1] * 2)
    ],
    [1.] * 2
)
HTML(operator_html_tree(op, True))
    • LincombOperator
    • NumpyVectorSpace(1) → NumpyVectorSpace(1)
    • []
      • NumpyMatrixOperator
      • NumpyVectorSpace(1) → NumpyVectorSpace(1)
      • []
        [[1.]]
      • +
      • LincombOperator
      • NumpyVectorSpace(1) → NumpyVectorSpace(1)
      • []
        • NumpyMatrixOperator
        • NumpyVectorSpace(1) → NumpyVectorSpace(1)
        • []
          [[1.]]
        • +
        • NumpyMatrixOperator
        • NumpyVectorSpace(1) → NumpyVectorSpace(1)
        • []
          [[1.]] 
NumpyMatrixOperator(np.ones([1]))
    • NumpyMatrixOperator
    • NumpyVectorSpace(1) → NumpyVectorSpace(1)
    • []
      [[1.]] 
op
    • LincombOperator
    • NumpyVectorSpace(1) → NumpyVectorSpace(1)
    • []
      • NumpyMatrixOperator
      • NumpyVectorSpace(1) → NumpyVectorSpace(1)
      • []
        [[1.]]
      • +
      • LincombOperator
      • NumpyVectorSpace(1) → NumpyVectorSpace(1)
      • []
        • NumpyMatrixOperator
        • NumpyVectorSpace(1) → NumpyVectorSpace(1)
        • []
          [[1.]]
        • +
        • NumpyMatrixOperator
        • NumpyVectorSpace(1) → NumpyVectorSpace(1)
        • []
          [[1.]]

XarrayMatrixOperator

# #| export
# def _get_space_dims(dims, matrix, default):
#     if dims is None:
#         if hasattr(matrix, 'dims'): return list(matrix.dims)
#         else: return [default]
#     elif isinstance(dims, str): return [dims]
#     elif hasattr(dims, 'dims'): return dims.dims
#     return list(dims)
space_a = XarrayVectorSpace({'A': np.linspace(0, 1, 2)})
_label_dims_with(['A'], Lbl.SRC)
['A (source)']
_label_dims_with(validate_coord_data(space_a.coords), Lbl.SRC)
{'A (source)': (array([0., 1.]), {})}
_label_dims_with(space_a, Lbl.SRC)
A (source)(2) 
<xarray.DataArray 'XarrayVectorSpace' (A (source): 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
matrix = DataArray(
        np.arange(6).reshape(2, 3), 
        {'A': np.linspace(0, 1, 2), 'B': np.linspace(0, 1, 3)}
    )
matrix = _label_dims(matrix, ['A'], ['B'])
matrix
<xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
crds = _label_dims(space_a.coords, ['A'], ['B'])
crds
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
spc = _label_dims(space_a, ['A'], ['B'])
spc
A (source)(2) 
<xarray.DataArray 'XarrayVectorSpace' (A (source): 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
matrix
<xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
unlabel_dims('A (source)')
'A'
unlabel_dims(matrix)
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A        (A) float64 16B 0.0 1.0
  * B        (B) float64 24B 0.0 0.5 1.0
unlabel_dims(crds)
Coordinates:
  * A        (A) float64 16B 0.0 1.0
spc
A (source)(2) 
<xarray.DataArray 'XarrayVectorSpace' (A (source): 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
unlabel_dims(spc)
A(2) 
<xarray.DataArray 'XarrayVectorSpace' (A: 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A        (A) float64 16B 0.0 1.0
_swap_dim_labels(matrix)
<xarray.DataArray (A (range): 2, B (source): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (range)   (A (range)) float64 16B 0.0 1.0
  * B (source)  (B (source)) float64 24B 0.0 0.5 1.0

coordinates_to_list


def coordinates_to_list(
    coords, dims
):

XarrayMatrixOperator


def XarrayMatrixOperator(
    matrix:DataArray | ndarray | SparseArray | sparray, # N-dimensional matrix
    source:XarrayVectorSpace | DataArray | Coordinates | dict | list | str | None=None, # Source vector space or coordinates
    range:XarrayVectorSpace | DataArray | Coordinates | dict | list | str | None=None, # Range vector space or coordinates
    solver_options:dict=None, # Options for matrix solver
    name:str=None, # Operator name
):

An Operator backed by an xarray DataArray.

A pyMOR Operator maps a VectorArray from a source VectorSpace to another VectorArray in a range VectorSpace. The XarrayMatrixOperator acts on an XarrayVectorArray from an XarrayVectorSpace by means of multiplication by an xarray DataArray.

The html representation of an XarrayMatrixOperator is operator name{source space} → {range space}, where X is used as the name if the operator has not been given an explict name.

The matrix can be supplied as a DataArray, or as an unlabeled numpy ndarray, scipy.sparse sparray, or sparse SparseArray.

Create from a numpy array and specified range and source spaces:

space_a = XarrayVectorSpace({'A': np.linspace(0, 1, 2)})
space_b = XarrayVectorSpace({'B': ['c', 'd', 'e']})
op = XarrayMatrixOperator(np.arange(6).reshape(2, 3), source=space_b, range=space_a, name='foo')
op
    • foo
    • XarrayMatrixOperator
    • B(3) → A(2)
    • []
     
      <xarray.DataArray 'foo' (A (range): 2, B (source): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (range)   (A (range)) float64 16B 0.0 1.0
  * B (source)  (B (source)) <U1 12B 'c' 'd' 'e'

The DataArray backing the Operator is stored in the matrix attribute:

op.matrix
<xarray.DataArray 'foo' (A (range): 2, B (source): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (range)   (A (range)) float64 16B 0.0 1.0
  * B (source)  (B (source)) <U1 12B 'c' 'd' 'e'

The source and range attributes store the corresponding VectorArrays:

op.source
B(3) 
<xarray.DataArray 'XarrayVectorSpace' (B: 3)> Size: 24B
array([0., 0., 0.])
Coordinates:
  * B        (B) <U1 12B 'c' 'd' 'e'
op.range
A(2) 
<xarray.DataArray 'XarrayVectorSpace' (A: 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A        (A) float64 16B 0.0 1.0

Create from a DataArray with defined coords:

da = DataArray(
        np.arange(6).reshape(2, 3), 
        {'A (source)': np.linspace(0, 1, 2), 'B (range)': np.linspace(0, 1, 3)}
    )
op = XarrayMatrixOperator(da)
op
    • XarrayMatrixOperator
    • A(2) → B(3)
    • []
     
      <xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
op.matrix
<xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
op.range
B(3) 
<xarray.DataArray 'XarrayVectorSpace' (B: 3)> Size: 24B
array([0., 0., 0.])
Coordinates:
  * B        (B) float64 24B 0.0 0.5 1.0
da = DataArray(
        np.arange(6).reshape(2, 3), 
        {'A': np.linspace(0, 1, 2), 'B': np.linspace(0, 1, 3)}
    )
_label_dims(da, src_dims=['A'], rng_dims=['B'])
<xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
op = XarrayMatrixOperator(da, source=['A'], range=['B'])
op
    • XarrayMatrixOperator
    • A(2) → B(3)
    • []
     
      <xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
op.matrix
<xarray.DataArray (A (source): 2, B (range): 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A (source)  (A (source)) float64 16B 0.0 1.0
  * B (range)   (B (range)) float64 24B 0.0 0.5 1.0
XarrayMatrixOperator(sparse2d_rand([4, 5], density=.5, random_state=42), source='A', range='B')
    • XarrayMatrixOperator
    • A(5) → B(4)
    • []
     
      <xarray.DataArray (B (range): 4, A (source): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (range)   (B (range)) int64 32B 0 1 2 3
  * A (source)  (A (source)) int64 40B 0 1 2 3 4
op = XarrayMatrixOperator(sparse2d_rand([4, 5], density=.5, random_state=42), source={'A': np.arange(5)}, range={'B': np.arange(4)})
op
    • XarrayMatrixOperator
    • A(5) → B(4)
    • []
     
      <xarray.DataArray (B (range): 4, A (source): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (range)   (B (range)) int64 32B 0 1 2 3
  * A (source)  (A (source)) int64 40B 0 1 2 3 4
op.matrix
<xarray.DataArray (B (range): 4, A (source): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (range)   (B (range)) int64 32B 0 1 2 3
  * A (source)  (A (source)) int64 40B 0 1 2 3 4

XarrayMatrixOperator.H


def H(
    
):

The adjoint operator.

op
    • XarrayMatrixOperator
    • A(5) → B(4)
    • []
     
      <xarray.DataArray (B (range): 4, A (source): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (range)   (B (range)) int64 32B 0 1 2 3
  * A (source)  (A (source)) int64 40B 0 1 2 3 4
op.H
    • XarrayMatrixOperator_adjoint
    • XarrayMatrixOperator
    • B(4) → A(5)
    • []
     
      <xarray.DataArray (B (source): 4, A (range): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (source)  (B (source)) int64 32B 0 1 2 3
  * A (range)   (A (range)) int64 40B 0 1 2 3 4

contract_matvec


def contract_matvec(
    a, b, axes
):

Contract one axis of the first (2D) array with one axis of the second (ND) array.

contract_matvec(np.ones([2,5]), np.ones([4,5,6]), [1, 1]).shape
(4, 2, 6)
contract_matvec(np.ones([6,2]), np.ones([4,5,6]), [0, 2]).shape
(4, 5, 2)

dims_matvec


def dims_matvec(
    a, b
):

Contract one axis of the first (2D) array with the axis of the second (ND) array of the same name.

dims_matvec([np.ones([2,5]), ["A", "B"]], [np.ones([4,5,6]), ["C", "B", "D"]])
(array([[[5., 5., 5., 5., 5., 5.],
         [5., 5., 5., 5., 5., 5.]],
 
        [[5., 5., 5., 5., 5., 5.],
         [5., 5., 5., 5., 5., 5.]],
 
        [[5., 5., 5., 5., 5., 5.],
         [5., 5., 5., 5., 5., 5.]],
 
        [[5., 5., 5., 5., 5., 5.],
         [5., 5., 5., 5., 5., 5.]]]),
 ['C', 'A', 'D'])
dims_matvec([np.ones([2]), ['A']], [np.ones([2, 3]), ['A', 'B']])
(array([2., 2., 2.]), ['B'])

XarrayVectorSpace.__truediv__


def __truediv__(
    other
):

replace_item


def replace_item(
    dct, old_key, new_pair
):

XarrayVectorSpace.range_under


def range_under(
    op:XarrayMatrixOperator
):

Find the space that this space is mapped to by an XarrayMatrixOperator. Assumes that the operator has one-dimensional source and range.

U = op.source.zeros({'len': [0]})
U.space.range_under(op)
B(4) 
<xarray.DataArray 'XarrayVectorSpace' (B: 4)> Size: 32B
array([0., 0., 0., 0.])
Coordinates:
  * B        (B) int64 32B 0 1 2 3
S = XarrayVectorSpace({'A': np.arange(3)})
op1 = XarrayMatrixOperator(np.ones(3), source=S)
S.range_under(op1)
<xarray.DataArray 'XarrayVectorSpace' ()> Size: 8B
array(0.)
S = XarrayVectorSpace({'A': np.arange(3), 'B': np.arange(2)})
S.range_under(op1)
B(2) 
<xarray.DataArray 'XarrayVectorSpace' (B: 2)> Size: 16B
array([0., 0.])
Coordinates:
  * B        (B) int64 16B 0 1

XarrayMatrixOperator.apply


def apply(
    U, # `XarrayVectorArray` of vectors to which the operator is applied
    mu:NoneType=None, # The parameter values for which to evaluate the operator
    method:str='numpy'
)->XarrayVectorArray: # `XarrayVectorArray` in the range `XarrayVectorSpace`

Apply the operator to an XarrayVectorArray in the source vector space.

U1 = S.ones()
op1.apply(U1)
∅ ⛒ B(2) 
<xarray.DataArray (B: 2)> Size: 16B
array([3., 3.])
Coordinates:
  * B        (B) int64 16B 0 1
op1 = XarrayMatrixOperator(np.ones((2, 2)), source={'A': np.arange(2)}, range={'B': np.arange(2)})
U1 = XarrayVectorSpace({'A': np.arange(2), 'C': np.arange(3)}).ones()
op1.apply(U1)
∅ ⛒ B(2) ⛒ C(3) 
<xarray.DataArray (B: 2, C: 3)> Size: 48B
array([[2., 2., 2.],
       [2., 2., 2.]])
Coordinates:
  * B        (B) int64 16B 0 1
  * C        (C) int64 24B 0 1 2
op1.apply(U1, method='xarray')
∅ ⛒ B(2) ⛒ C(3) 
<xarray.DataArray (B: 2, C: 3)> Size: 48B
array([[2., 2., 2.],
       [2., 2., 2.]])
Coordinates:
  * B        (B) int64 16B 0 1
  * C        (C) int64 24B 0 1 2
from pymor.tools.random import new_rng
with new_rng(42):
    U = op.source.random(extended_coords={'foo': np.arange(2)})
U
foo(2) ⛒ A(5) 
<xarray.DataArray (foo: 2, A: 5)> Size: 80B
array([[0.77395605, 0.85859792, 0.09417735, 0.7611397 , 0.12811363],
       [0.43887844, 0.69736803, 0.97562235, 0.78606431, 0.45038594]])
Coordinates:
  * foo      (foo) int64 16B 0 1
  * A        (A) int64 40B 0 1 2 3 4
op.apply(U)
foo(2) ⛒ B(4) 
<xarray.DataArray (foo: 2, B: 4)> Size: 64B
array([[0.12601957, 0.62030984, 0.00597851, 0.8288864 ],
       [0.44580162, 0.73140252, 0.02101756, 1.17285242]])
Coordinates:
  * foo      (foo) int64 16B 0 1
  * B        (B) int64 32B 0 1 2 3
xr.testing.assert_allclose(op.apply(U, method='numpy').array, op.apply(U, method='xarray').array)
from pymor.basic import Parameters
dof = Coordinates({'dof': ('dof', ['x', 'v'], {'long_name': "Degrees of freedom"})})
space = XarrayVectorSpace(dof)
A = XarrayMatrixOperator(np.array([[0., -1.], [1., 0]]), source=space, range=space)
B = XarrayMatrixOperator(np.array([[0., 0.], [0., 1.]]), source=space, range=space)

operator = LincombOperator(
    [A, B], 
    [
        ExpressionParameterFunctional("w0**2", Parameters({'w0': 1})), 
        ExpressionParameterFunctional("c", Parameters({'c': 1}))
    ]
)
mu = dict(c=.1, w0="sin(t)")
mu = Parameters(dict(c=1, w0=1)).parse(mu)

y = operator.source.from_numpy(np.array([1., 0.]))

operator.apply(y, mu=mu.at_time(t=.1))
∅ ⛒ dof(2) 
<xarray.DataArray (dof: 2)> Size: 16B
array([0.        , 0.00996671])
Coordinates:
  * dof      (dof) <U1 8B 'x' 'v'
CPU times: user 278 μs, sys: 10 μs, total: 288 μs
Wall time: 284 μs
∅ ⛒ dof(2) 
<xarray.DataArray (dof: 2)> Size: 16B
array([0.        , 0.00996671])
Coordinates:
  * dof      (dof) <U1 8B 'x' 'v'

XarrayMatrixOperator.apply_adjoint


def apply_adjoint(
    V, # `XarrayVectorArray` of vectors to which the operator is applied
    mu:NoneType=None, # The parameter values for which to evaluate the operator
)->XarrayVectorArray: # `XarrayVectorArray` in the source `XarrayVectorSpace`

Apply the adjoint of the operator to an XarrayVectorArray in the range vector space.

V = op.range.zeros({'len': [0]})
op.source
A(5) 
<xarray.DataArray 'XarrayVectorSpace' (A: 5)> Size: 40B
array([0., 0., 0., 0., 0.])
Coordinates:
  * A        (A) int64 40B 0 1 2 3 4
op.range
B(4) 
<xarray.DataArray 'XarrayVectorSpace' (B: 4)> Size: 32B
array([0., 0., 0., 0.])
Coordinates:
  * B        (B) int64 32B 0 1 2 3
op.H
    • XarrayMatrixOperator_adjoint
    • XarrayMatrixOperator
    • B(4) → A(5)
    • []
     
      <xarray.DataArray (B (source): 4, A (range): 5)> Size: 160B
<COO: shape=(4, 5), dtype=float64, nnz=10, fill_value=0.0>
Coordinates:
  * B (source)  (B (source)) int64 32B 0 1 2 3
  * A (range)   (A (range)) int64 40B 0 1 2 3 4
op.H.range
A(5) 
<xarray.DataArray 'XarrayVectorSpace' (A: 5)> Size: 40B
array([0., 0., 0., 0., 0.])
Coordinates:
  * A        (A) int64 40B 0 1 2 3 4
V.space.range_under(op.H)
A(5) 
<xarray.DataArray 'XarrayVectorSpace' (A: 5)> Size: 40B
array([0., 0., 0., 0., 0.])
Coordinates:
  * A        (A) int64 40B 0 1 2 3 4
op.apply_adjoint(V)
len(1) ⛒ A(5) 
<xarray.DataArray (len: 1, A: 5)> Size: 40B
array([[0., 0., 0., 0., 0.]])
Coordinates:
  * len      (len) int64 8B 0
  * A        (A) int64 40B 0 1 2 3 4
# #| export
# @patch
# def _assemble_lincomb(self:XarrayMatrixOperator, operators, coefficients, identity_shift=0., solver_options=None, name=None):
#     if not all(isinstance(op, XarrayMatrixOperator) for op in operators):
#         return None

#     common_mat_dtype = reduce(np.promote_types,
#                               (op.matrix.dtype for op in operators if hasattr(op, 'matrix')))
#     common_coef_dtype = reduce(np.promote_types, (type(c) for c in coefficients + [identity_shift]))
#     common_dtype = np.promote_types(common_mat_dtype, common_coef_dtype)

#     if coefficients[0] == 1:
#         matrix = operators[0].matrix.astype(common_dtype)
#     else:
#         matrix = operators[0].matrix * coefficients[0]
#         if matrix.dtype != common_dtype:
#             matrix = matrix.astype(common_dtype)

#     for op, c in zip(operators[1:], coefficients[1:]):
#         if c == 1:
#             try:
#                 matrix += op.matrix
#             except NotImplementedError:
#                 matrix = matrix + op.matrix
#         elif c == -1:
#             try:
#                 matrix -= op.matrix
#             except NotImplementedError:
#                 matrix = matrix - op.matrix
#         else:
#             try:
#                 matrix += (op.matrix * c)
#             except NotImplementedError:
#                 matrix = matrix + (op.matrix * c)

#     if identity_shift: raise NotImplementedError

#     return XarrayMatrixOperator(matrix, source=self.source, range=self.range, solver_options=solver_options)

XarrayMatrixOperator.to_numpy


def to_numpy(
    
):

Return the operator in dense ndarray format.

The returned numpy array will be of shape \(r_1\times \ldots \times r_n \times s_1 \times \ldots \times s_m\) where \(r_i\) are the dimensions of the \(n\) vector spaces making up the range product vector space, and \(s_j\) are the dimensions of the \(m\) vector spaces making up the source product vector space.

op.to_numpy()
array([[7.78765841e-04, 0.00000000e+00, 7.06630522e-03, 0.00000000e+00,
        9.73755519e-01],
       [6.11653160e-01, 2.30624250e-02, 0.00000000e+00, 0.00000000e+00,
        9.92211559e-01],
       [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
        4.66656632e-02],
       [3.99860972e-01, 0.00000000e+00, 5.24774660e-01, 6.17481510e-01,
        0.00000000e+00]])

densify


def densify(
    obj:Model | Operator, # Object to densify
)->Model | Operator: # Densified object

Convert sparse operators to dense.

An operator backed by a sparse DataArray:

op.sparse
True

Convert the operator to one backed by a dense DataArray:

densify(op).sparse
False

XarrayFunctionalOperator

XarrayFunctionalOperator


def XarrayFunctionalOperator(
    functional:ParameterFunctional, range, source
):

An Operator described by a ParameterFunctional that assembles to a XarrayMatrixOperator.

from pymor.basic import GenericParameterFunctional
func = GenericParameterFunctional(lambda mu: DataArray(np.arange(4).reshape(2, 2)**mu['c'], dims=['A', 'B']), {'c': 1})
func.evaluate(Mu({'c': 2}))
<xarray.DataArray (A: 2, B: 2)> Size: 32B
array([[0, 1],
       [4, 9]])
Dimensions without coordinates: A, B
op = XarrayFunctionalOperator(
        GenericParameterFunctional(lambda mu: DataArray(np.arange(4).reshape(2, 2)**mu['c'], dims=['A (source)', 'B (range)']), {'c': 1}),
        range=XarrayVectorSpace({'A': [0, 1]}),
        source=XarrayVectorSpace({'B': [0, 1]}),
    )
op
    • XarrayFunctionalOperator
    • B(2) → A(2)
    • ['c']
op.assemble(Mu({'c': 2}))
    • XarrayMatrixOperator
    • A(2) → B(2)
    • []
     
      <xarray.DataArray (A (source): 2, B (range): 2)> Size: 32B
array([[0, 1],
       [4, 9]])
Dimensions without coordinates: A (source), B (range)

IdentityOperator

# @patch
# def apply(self:IdentityOperator, U, mu=None):
#     if isinstance(U, XarrayVectorArray):

# space_a = XarrayVectorSpace({'A (source)': np.linspace(0, 1, 2)})
# space_b = XarrayVectorSpace({'B': ['c', 'd', 'e']})

# space_a * space_b

# space

# IdentityOperator(space).source.dims

SumOperator

SumOperator


def SumOperator(
    source_coords:Coordinates | dict, # Coordinates to sum over
    name:NoneType=None, # Operator name
):

An Operator that sums over one or more dimensions of a XarrayVectorSpace.

space = XarrayVectorSpace({'A': [1, 2], 'B': ['a', 'b', 'c']}, name='foo')
np.random.seed(42)
U = space.from_numpy(np.random.rand(space.dim))
U.array.A
<xarray.DataArray 'A' (A: 2)> Size: 16B
array([1, 2])
Coordinates:
  * A        (A) int64 16B 1 2
op = SumOperator(U.array.A)
op
    • SumOperator
    • A(2) → Sum(1)
    • []
op.source
A(2) 
<xarray.DataArray 'A' (A: 2)> Size: 16B
array([0., 0.])
Coordinates:
  * A        (A) int64 16B 1 2
op.range
Sum(1) 
<xarray.DataArray 'XarrayVectorSpace' (Sum: 1)> Size: 8B
array([0.])
Coordinates:
  * Sum      (Sum) <U3 12B 'Sum'
op.apply(U).array
<xarray.DataArray (Sum: 1, B: 3)> Size: 24B
array([[0.9731986 , 1.10673295, 0.88798846]])
Coordinates:
  * Sum      (Sum) <U3 12B 'Sum'
  * B        (B) <U1 12B 'a' 'b' 'c'

ScaleOperator

ScaleOperator


def ScaleOperator(
    array, space:NoneType=None, name:NoneType=None
):

A scaling operator for XarrayVectorSpaces.

space = XarrayVectorSpace({'A': [1, 2], 'B': ['a', 'b', 'c']}, name='foo')
U = space.from_numpy(np.arange(6))
U.array
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[0, 1, 2],
       [3, 4, 5]])
Coordinates:
  * A        (A) int64 16B 1 2
  * B        (B) <U1 12B 'a' 'b' 'c'
arr = U.array[0].drop_vars('A')
arr
<xarray.DataArray (B: 3)> Size: 24B
array([0, 1, 2])
Coordinates:
  * B        (B) <U1 12B 'a' 'b' 'c'
scale_op = ScaleOperator(arr, name='bar')
scale_op
    • bar
    • ScaleOperator
    • B(3) → B(3)
    • []
scale_op.apply(U)
∅ ⛒ A(2) ⛒ B(3) 
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[ 0,  1,  4],
       [ 0,  4, 10]])
Coordinates:
  * A        (A) int64 16B 1 2
  * B        (B) <U1 12B 'a' 'b' 'c'
space.ones().array.data
array([[1., 1., 1.],
       [1., 1., 1.]])
ScaleOperator(space.ones().array.data, space=space)
    • ScaleOperator
    • A(2) ⛒ B(3) → A(2) ⛒ B(3)
    • []

LincombOperator

LincombOperator.__init__


def __init__(
    operators, coefficients, solver:NoneType=None, name:NoneType=None
):

Initialize self. See help(type(self)) for accurate signature.

# #|export
# @match_class(LincombOperator)
# def action_LincombOperator(self, op):
#     op_coefficients = op.evaluate_coefficients(self.mu)
#     res = op_coefficients[0] * self.apply(op.operators[0])
#     for i in range(1, len(op.operators)):
#         res = res + op_coefficients[i] * self.apply(op.operators[i])
# a = XarrayMatrixOperator(np.ones([3, 2]), source=space_a, range=space_b)
# b = XarrayMatrixOperator(np.eye(2), source='source a', range='range a')
# LincombOperator([a, b], [1., 1.])
_read_sparse_dataarray('test')
<xarray.DataArray 'test' (dim_1: 2, another dim: 3, third dim: 4)> Size: 64B
<COO: shape=(2, 3, 4), dtype=float64, nnz=2, fill_value=0.0>
Coordinates:
  * dim_1        (dim_1) <U8 64B 'A**2' 'A*sin(B)'
  * another dim  (another dim) <U3 36B 'foo' 'bar' 'baz'
  * third dim    (third dim) <U3 48B '$a$' '$b$' '$c$' '$d$'
Attributes:
    parameters:  ['A', 'B']
with importlib.resources.path("pylgs.systems.NaD1_Toy", "Flux.mtxn") as path:
    out = _read_sparse_dataarray(path)
out
<xarray.DataArray 'Flux' (coefficients: 1, Transition: 1,
                          Density matrix (source): 4)> Size: 32B
<COO: shape=(1, 1, 4), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * coefficients             (coefficients) <U1 4B '1'
  * Transition               (Transition) <U33 132B '3P<sub>1/2</sub>→3S<sub>...
  * Density matrix (source)  (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...
Attributes:
    parameters:  ['']

LincombOperator.from_file


def from_file(
    file_name:str | Path, # File name
    solver:NoneType=None, name:NoneType=None
):

Read a LincombOperator (list of sparse arrays and corresponding ExpressionParameterFunctionals) from a .mtxn file.

with importlib.resources.path("pylgs.systems.NaD1_Toy", "b.mtxn") as path:
    op = LincombOperator.from_file(path)
op
    • LincombOperator
    • → Density matrix(4)
    • ['BeamTransitRatePerS']
      • b
      • XarrayMatrixOperator
      • → Density matrix(4)
      • []
       
        <xarray.DataArray 'b' (Density matrix (range): 4, none: 1)> Size: 24B
<COO: shape=(4, 1), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * Density matrix (range)  (Density matrix (range)) <U50 800B 'ρ<sub>Re, 3S<...
  * none                    (none) <U4 16B 'none'
op.as_range_array(Mu({'BeamTransitRatePerS': 1}))
∅ ⛒ Density matrix(4) 
<xarray.DataArray (Density matrix: 4)> Size: 32B
array([1., 0., 0., 0.])
Coordinates:
  * Density matrix  (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, 3...
with importlib.resources.path("pylgs.systems.NaD1_Toy", "Flux.mtxn") as path:
    op = LincombOperator.from_file(path)
op
    • LincombOperator
    • Density matrix(4) →
    • []
      • Flux
      • XarrayMatrixOperator
      • Density matrix(4) →
      • []
       
        <xarray.DataArray 'Flux' (Transition: 1, Density matrix (source): 4)> Size: 24B
<COO: shape=(1, 4), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * Transition               (Transition) <U33 132B '3P<sub>1/2</sub>→3S<sub>...
  * Density matrix (source)  (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...
op.source
Density matrix(4) 
<xarray.DataArray 'XarrayVectorSpace' (Density matrix: 4)> Size: 32B
array([0., 0., 0., 0.])
Coordinates:
  * Density matrix  (Density matrix) <U50 800B 'ρ<sub>Re, 3S<sub>1/2</sub>, 3...
op.range
<xarray.DataArray 'XarrayVectorSpace' ()> Size: 8B
array(0.)
op.assemble()
    • Flux
    • XarrayMatrixOperator
    • Density matrix(4) →
    • []
     
      <xarray.DataArray 'Flux' (Transition: 1, Density matrix (source): 4)> Size: 24B
<COO: shape=(1, 4), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * Transition               (Transition) <U33 132B '3P<sub>1/2</sub>→3S<sub>...
  * Density matrix (source)  (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...
# op = LincombOperator.from_file('test.mtxn')
# op

ExpressionParameterFunctional.partial_evaluate


def partial_evaluate(
    mu:dict | Mu, # Parameter values
)->float | ExpressionParameterFunctional:

Substitute parameter values, returning a new ExpressionParameterFunctional if expression does not evaluate to a number.

# op.coefficients[1].partial_evaluate({'A': [.5]})
# op.coefficients[1].partial_evaluate({'A': [.5], 'B': [2.]})

LincombOperator.partial_evaluate_coefficients


def partial_evaluate_coefficients(
    mu:Mu, # Parameter values to substitute
):

Substitute parameter values into linear coefficients, returning a new ExpressionParameterFunctional if expression does not evaluate to a number.

# op.partial_evaluate_coefficients(Mu({'A': .5}))

LincombOperator.partial_assemble


def partial_assemble(
    mu:NoneType=None, # Parameter values to substitute
)->Operator:

Substitute parameter values into the linear coefficients, returning a new operator with fewer (or no) parameters.

op
    • LincombOperator
    • Density matrix(4) →
    • []
      • Flux
      • XarrayMatrixOperator
      • Density matrix(4) →
      • []
       
        <xarray.DataArray 'Flux' (Transition: 1, Density matrix (source): 4)> Size: 24B
<COO: shape=(1, 4), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * Transition               (Transition) <U33 132B '3P<sub>1/2</sub>→3S<sub>...
  * Density matrix (source)  (Density matrix (source)) <U50 800B 'ρ<sub>Re, 3...
# op.partial_assemble(Mu({'A': .5}))
# op.partial_assemble(Mu({'A': .5, 'B': 2.}))
# op.partial_assemble(Mu({'A': 0}))

LincombOperator.terms


def terms(
    
):

# #| export
# @patch
# def scalar_part(self:LincombOperator, mu):
#     mu_vector = Mu({k: v for k, v in mu.items() if v.size > 1})
#     return np.sum([o for o in self.terms if not set(mu_vector).intersection(o.parameters)])
# #| export
# @patch
# def vector_part(self:LincombOperator, mu):
#     mu_vector = Mu({k: v for k, v in mu.items() if v.size > 1})
#     return np.sum([o for o in self.terms if set(mu_vector).intersection(o.parameters)])

ProductOperator

ProductOperator


def ProductOperator(
    operators, # Sequence of operators that each assemble to XarrayMatrixOperator
    name:NoneType=None, # Optional name
):

An Operator given by the direct product of operators operators.

A = XarrayMatrixOperator(array([[1, 0], [0, 1]]), source={'source a': [0, 1]}, range={'range a': [0, 1]})
B = XarrayMatrixOperator(array([[0, 1], [1, 0]]), source={'source b': [0, 1]}, range={'range b': [0, 1]})

Operator.__mul__


def __mul__(
    other
):

ProductOperator.__str__


def __str__(
    
):

Return str(self).

IdentityOperator.__str__


def __str__(
    
):

Return str(self).

A * IdentityOperator(B.source)
    • ProductOperator
    • source a(2) ⛒ source b(2) → range a(2) ⛒ source b(2)
    • []
      • XarrayMatrixOperator
      • source a(2) → range a(2)
      • []
       
        <xarray.DataArray (range a (range): 2, source a (source): 2)> Size: 32B
array([[1, 0],
       [0, 1]])
Coordinates:
  * range a (range)    (range a (range)) int64 16B 0 1
  * source a (source)  (source a (source)) int64 16B 0 1
      • ⛒
      • IdentityOperator
      • source b(2) → source b(2)
      • []
IdentityOperator(B.source) * A
    • ProductOperator
    • source b(2) ⛒ source a(2) → source b(2) ⛒ range a(2)
    • []
      • IdentityOperator
      • source b(2) → source b(2)
      • []
      • ⛒
      • XarrayMatrixOperator
      • source a(2) → range a(2)
      • []
       
        <xarray.DataArray (range a (range): 2, source a (source): 2)> Size: 32B
array([[1, 0],
       [0, 1]])
Coordinates:
  * range a (range)    (range a (range)) int64 16B 0 1
  * source a (source)  (source a (source)) int64 16B 0 1
def vec(a):
    return a.T.ravel()
s=1
A = XarrayMatrixOperator(
    sparse2d_rand([3*s, 2*s], density=.1, random_state=42), 
    source={'source a': np.arange(2*s)}, 
    range={'range a': np.arange(3*s)}
)
A
    • XarrayMatrixOperator
    • source a(2) → range a(3)
    • []
     
      <xarray.DataArray (range a (range): 3, source a (source): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (range)    (range a (range)) int64 24B 0 1 2
  * source a (source)  (source a (source)) int64 16B 0 1
mu = Mu(A=.5)
B = ExpressionParameterFunctional('A**2', {'A': 1}) * XarrayMatrixOperator(
    sparse2d_rand([4*s, 5*s], density=.1, random_state=42),
    source={'source b': np.arange(5*s)}, 
    range={'range b': np.arange(4*s)}
)
B = B.assemble(mu)
op = A * B
op
    • ProductOperator
    • source a(2) ⛒ source b(5) → range a(3) ⛒ range b(4)
    • []
      • XarrayMatrixOperator
      • source a(2) → range a(3)
      • []
       
        <xarray.DataArray (range a (range): 3, source a (source): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (range)    (range a (range)) int64 24B 0 1 2
  * source a (source)  (source a (source)) int64 16B 0 1
      • ⛒
      • XarrayMatrixOperator
      • source b(5) → range b(4)
      • []
       
        <xarray.DataArray (range b (range): 4, source b (source): 5)> Size: 32B
<COO: shape=(4, 5), dtype=float64, nnz=2, fill_value=0.0>
Coordinates:
  * range b (range)    (range b (range)) int64 32B 0 1 2 3
  * source b (source)  (source b (source)) int64 40B 0 1 2 3 4
op.source
source a(2) ⛒ source b(5) 
<xarray.DataArray 'XarrayVectorSpace' (source a: 2, source b: 5)> Size: 80B
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])
Coordinates:
  * source a  (source a) int64 16B 0 1
  * source b  (source b) int64 40B 0 1 2 3 4
op.range
range a(3) ⛒ range b(4) 
<xarray.DataArray 'XarrayVectorSpace' (range a: 3, range b: 4)> Size: 96B
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])
Coordinates:
  * range a  (range a) int64 24B 0 1 2
  * range b  (range b) int64 32B 0 1 2 3
B.source._array
<xarray.DataArray 'XarrayVectorSpace' (source b: 5)> Size: 40B
array([0., 0., 0., 0., 0.])
Coordinates:
  * source b  (source b) int64 40B 0 1 2 3 4
sum_op = SumOperator(B.source._array['source b'])
sum_op * A
    • ProductOperator
    • source b(5) ⛒ source a(2) → Sum(1) ⛒ range a(3)
    • []
      • SumOperator
      • source b(5) → Sum(1)
      • []
      • ⛒
      • XarrayMatrixOperator
      • source a(2) → range a(3)
      • []
       
        <xarray.DataArray (range a (range): 3, source a (source): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (range)    (range a (range)) int64 24B 0 1 2
  * source a (source)  (source a (source)) int64 16B 0 1
A * sum_op
    • ProductOperator
    • source a(2) ⛒ source b(5) → range a(3) ⛒ Sum(1)
    • []
      • XarrayMatrixOperator
      • source a(2) → range a(3)
      • []
       
        <xarray.DataArray (range a (range): 3, source a (source): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (range)    (range a (range)) int64 24B 0 1 2
  * source a (source)  (source a (source)) int64 16B 0 1
      • ⛒
      • SumOperator
      • source b(5) → Sum(1)
      • []
np.random.seed(42)
U = op.source.from_numpy(np.random.rand(op.source._array.size))
V = op.range.from_numpy(np.random.rand(op.range._array.size))
U.space
source a(2) ⛒ source b(5) 
<xarray.DataArray 'XarrayVectorSpace' (source a: 2, source b: 5)> Size: 80B
array([[0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0.]])
Coordinates:
  * source a  (source a) int64 16B 0 1
  * source b  (source b) int64 40B 0 1 2 3 4
An = A.matrix.data
Bn = B.assemble(mu).matrix.data
Un = U.to_numpy().T
Vn = V.to_numpy().T
An = A.matrix.data
Bn = B.assemble(mu).matrix.data
Un = U.to_numpy().T
out1 = np.kron(An, Bn).dot(vec(Un))
CPU times: user 96.3 ms, sys: 2.56 ms, total: 98.8 ms
Wall time: 98.2 ms
An = A.matrix.data
Bn = B.assemble(mu).matrix.data
Un = U.array.T.data
out3 = An.dot(Bn.dot(Un).T).ravel()
CPU times: user 682 μs, sys: 35 μs, total: 717 μs
Wall time: 692 μs
B.apply(U)
∅ ⛒ source a(2) ⛒ range b(4) 
<xarray.DataArray (source a: 2, range b: 4)> Size: 64B
array([[7.29197627e-05, 3.87008746e-02, 0.00000000e+00, 0.00000000e+00],
       [3.03708010e-05, 1.75639449e-01, 0.00000000e+00, 0.00000000e+00]])
Coordinates:
  * source a  (source a) int64 16B 0 1
  * range b   (range b) int64 32B 0 1 2 3
out4 = A.apply(B.apply(U)).to_numpy()
CPU times: user 694 μs, sys: 19 μs, total: 713 μs
Wall time: 702 μs
test_close(out1.ravel(), out3.ravel())
test_close(out1.ravel(), out4.ravel())
test_close(
    np.kron(An, Bn).dot(vec(Un)),
    An.dot(Bn.dot(Un).T).ravel()
)

N.B.: to_numpy on a NumpyVectorArray effectively takes the transpose.

ProductOperator.apply


def apply(
    U, # |VectorArray| of vectors to which the operator is applied.
    mu:NoneType=None, # The |parameter values| for which to evaluate the operator.
):

Apply the operator to a |VectorArray|.

W = space.from_numpy(np.arange(12), extended_coord_data='len')
C = ScaleOperator(W.array[0, 0].drop_vars(['A', 'len']), name='C')
C.apply(W).name
'C'
D = SumOperator(W.array.A, name='D')
D.apply(W).name
'D'
C = ScaleOperator(W.array[0, 0].drop_vars(['A', 'len']))
E = (C * D)
E.apply(W).name
'XarrayVectorArray'
E.name
'ProductOperator'
E = C * (C * D)
E.apply(W).array.name
K = A * B.assemble(mu)
out2 = K.apply(U).to_numpy()
CPU times: user 1.08 ms, sys: 58 μs, total: 1.14 ms
Wall time: 1.11 ms
# test_close(out1.ravel(), out2.ravel())
test_close(
    op.apply(U, mu).to_numpy().ravel(), 
    np.kron(An, Bn).dot(vec(Un))
)

\((A\bigotimes B)^T=A^T\bigotimes B^T\)

test_close(np.kron(An, Bn).T, np.kron(An.T, Bn.T))
C = XarrayMatrixOperator(np.ones([2, 2]), source='A', range='A')
D = XarrayMatrixOperator(np.ones([3, 3]), source='B', range='B')
U1 = (C * D).source.ones()
op1 = C + C * D
op1
    • LincombOperator
    • A(2) ⛒ B(3) → A(2) ⛒ B(3)
    • []
      • XarrayMatrixOperator
      • A(2) → A(2)
      • []
       
        <xarray.DataArray (A (range): 2, A (source): 2)> Size: 32B
array([[1., 1.],
       [1., 1.]])
Coordinates:
  * A (range)   (A (range)) int64 16B 0 1
  * A (source)  (A (source)) int64 16B 0 1
      • +
      • ProductOperator
      • A(2) ⛒ B(3) → A(2) ⛒ B(3)
      • []
        • XarrayMatrixOperator
        • A(2) → A(2)
        • []
         
          <xarray.DataArray (A (range): 2, A (source): 2)> Size: 32B
array([[1., 1.],
       [1., 1.]])
Coordinates:
  * A (range)   (A (range)) int64 16B 0 1
  * A (source)  (A (source)) int64 16B 0 1
        • ⛒
        • XarrayMatrixOperator
        • B(3) → B(3)
        • []
         
          <xarray.DataArray (B (range): 3, B (source): 3)> Size: 72B
array([[1., 1., 1.],
       [1., 1., 1.],
       [1., 1., 1.]])
Coordinates:
  * B (range)   (B (range)) int64 24B 0 1 2
  * B (source)  (B (source)) int64 24B 0 1 2
C.apply(U1)
∅ ⛒ A(2) ⛒ B(3) 
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[2., 2., 2.],
       [2., 2., 2.]])
Coordinates:
  * A        (A) int64 16B 0 1
  * B        (B) int64 24B 0 1 2
D.apply(U1)
∅ ⛒ A(2) ⛒ B(3) 
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[3., 3., 3.],
       [3., 3., 3.]])
Coordinates:
  * A        (A) int64 16B 0 1
  * B        (B) int64 24B 0 1 2
(C * D).apply(U1)
∅ ⛒ A(2) ⛒ B(3) 
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[6., 6., 6.],
       [6., 6., 6.]])
Coordinates:
  * A        (A) int64 16B 0 1
  * B        (B) int64 24B 0 1 2
(C + (C * D)).apply(U1)
∅ ⛒ A(2) ⛒ B(3) 
<xarray.DataArray (A: 2, B: 3)> Size: 48B
array([[8., 8., 8.],
       [8., 8., 8.]])
Coordinates:
  * A        (A) int64 16B 0 1
  * B        (B) int64 24B 0 1 2

ProductOperator.H


def H(
    
):

op.H
    • ProductOperator
    • range a(3) ⛒ range b(4) → source a(2) ⛒ source b(5)
    • []
      • XarrayMatrixOperator_adjoint
      • XarrayMatrixOperator
      • range a(3) → source a(2)
      • []
       
        <xarray.DataArray (range a (source): 3, source a (range): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (source)  (range a (source)) int64 24B 0 1 2
  * source a (range)  (source a (range)) int64 16B 0 1
      • ⛒
      • XarrayMatrixOperator_adjoint
      • XarrayMatrixOperator
      • range b(4) → source b(5)
      • []
       
        <xarray.DataArray (range b (source): 4, source b (range): 5)> Size: 32B
<COO: shape=(4, 5), dtype=float64, nnz=2, fill_value=0.0>
Coordinates:
  * range b (source)  (range b (source)) int64 32B 0 1 2 3
  * source b (range)  (source b (range)) int64 40B 0 1 2 3 4

\(\text{vec}(V)(A\bigotimes B)=(A\bigotimes B)^T\text{vec}(V)\).

op.H.source
range a(3) ⛒ range b(4) 
<xarray.DataArray 'XarrayVectorSpace' (range a: 3, range b: 4)> Size: 96B
array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])
Coordinates:
  * range a  (range a) int64 24B 0 1 2
  * range b  (range b) int64 32B 0 1 2 3
test_close(
    vec(Vn).dot(np.kron(An, Bn).todense()),
    (np.kron(An, Bn).T).dot(vec(Vn))
)

ProductOperator.apply_adjoint


def apply_adjoint(
    V, # |VectorArray| of vectors to which the adjoint operator is applied.
    mu:NoneType=None, # The |parameter values| for which to apply the adjoint operator.
):

Apply the adjoint operator.

For any given linear |Operator| op, |parameter values| mu and |VectorArrays| U, V in the :attr:~pymor.operators.interface.Operator.source resp. :attr:~pymor.operators.interface.Operator.range we have::

op.apply_adjoint(V, mu).dot(U) == V.inner(op.apply(U, mu))

Thus, when op is represented by a matrix M, apply_adjoint is given by left-multiplication of (the complex conjugate of) M with V.

test_close(
    op.apply_adjoint(V, mu).to_numpy().ravel(),
    vec(Vn).dot(np.kron(An, Bn).todense())
)

ProductOperator.assemble


def assemble(
    mu:NoneType=None, # The |parameter values| for which to assemble the operator.
):

Assemble the operator for given |parameter values|.

The result of the method strongly depends on the given operator. For instance, a matrix-based operator will assemble its matrix, a |LincombOperator| will try to form the linear combination of its operators, whereas an arbitrary operator might simply return a :class:~pymor.operators.constructions.FixedParameterOperator. The only assured property of the assembled operator is that it no longer depends on a |Parameter|.

If the operator has a |Solver|, the assembled |Operator| will be equipped with the same |Solver|.

op.assemble(mu)
    • ProductOperator
    • source a(2) ⛒ source b(5) → range a(3) ⛒ range b(4)
    • []
      • XarrayMatrixOperator
      • source a(2) → range a(3)
      • []
       
        <xarray.DataArray (range a (range): 3, source a (source): 2)> Size: 16B
<COO: shape=(3, 2), dtype=float64, nnz=1, fill_value=0.0>
Coordinates:
  * range a (range)    (range a (range)) int64 24B 0 1 2
  * source a (source)  (source a (source)) int64 16B 0 1
      • ⛒
      • XarrayMatrixOperator
      • source b(5) → range b(4)
      • []
       
        <xarray.DataArray (range b (range): 4, source b (source): 5)> Size: 32B
<COO: shape=(4, 5), dtype=float64, nnz=2, fill_value=0.0>
Coordinates:
  * range b (range)    (range b (range)) int64 32B 0 1 2 3
  * source b (source)  (source b (source)) int64 40B 0 1 2 3 4
zero = (0 * B).assemble()
(A * zero).assemble()
    • ZeroOperator
    • source a(2) ⛒ source b(5) → range a(3) ⛒ range b(4)
    • []
ProductOperator([IdentityOperator(XarrayVectorSpace({k:v})) for k,v in space.coord_data.items()])
    • ProductOperator
    • A(2) ⛒ B(3) → A(2) ⛒ B(3)
    • []
      • IdentityOperator
      • A(2) → A(2)
      • []
      • ⛒
      • IdentityOperator
      • B(3) → B(3)
      • []
expand(IdentityOperator(space))
    • ProductOperator
    • A(2) ⛒ B(3) → A(2) ⛒ B(3)
    • []
      • IdentityOperator
      • A(2) → A(2)
      • []
      • ⛒
      • IdentityOperator
      • B(3) → B(3)
      • []
# #| hide
# #| exporti
# @match_class(ProductOperator)
# def to_matrix_ProductOperator(self, op):
#     return self.apply(contract(op))
    
# ToMatrixRules.insert_rule(-2, to_matrix_ProductOperator)
# #| hide
# to_matrix(op, mu=mu)
to_matrix(op, mu=mu)
<Compressed Sparse Row sparse array of dtype 'float64'
    with 2 stored elements and shape (12, 10)>
test_array('operators', 'to_matrix_ProductOperator_1', _.todense())
op1 = densify(op.operators[0]) * op.operators[1]
to_matrix(op1, mu=mu)
<Compressed Sparse Row sparse array of dtype 'float64'
    with 2 stored elements and shape (12, 10)>
test_array('operators', 'to_matrix_ProductOperator_1', _.todense())

ProductOperator.as_source_array


def as_source_array(
    mu:NoneType=None, # The |parameter values| for which to return the |VectorArray|
representation.
): # The |VectorArray| defined above.

Return a |VectorArray| representation of the operator in its source space.

In the case of a linear operator with |NumpyVectorSpace| as :attr:~pymor.operators.interface.Operator.range, this method returns for given |parameter values| mu a |VectorArray| V in the operator’s :attr:~pymor.operators.interface.Operator.source, such that ::

self.range.make_array(V.inner(U)) == self.apply(U, mu)

for all |VectorArrays| U.

# op.as_source_array(mu)
# test_close(
#     np.kron(An, Bn).todense().ravel(),
#     op.as_source_array(mu).to_numpy().T.ravel()
# )

ProductOperator.as_range_array


def as_range_array(
    mu:NoneType=None, # The |parameter values| for which to return the |VectorArray|
representation.
): # The |VectorArray| defined above.

Return a |VectorArray| representation of the operator in its range space.

In the case of a linear operator with |NumpyVectorSpace| as :attr:~pymor.operators.interface.Operator.source, this method returns for given |parameter values| mu a |VectorArray| V in the operator’s :attr:~pymor.operators.interface.Operator.range, such that ::

V.lincomb(U.to_numpy()) == self.apply(U, mu)

for all |VectorArrays| U.

# test_close(
#     np.kron(An, Bn).todense().ravel(),
#     op.as_range_array(mu).to_numpy().ravel()
# )