ModelOperator

ModelOperator#

class torchsim.ModelOperator(acquisition, *unknown, bounds=None, scale=None, amplitude=True, subspace=None)[source]#

Bases: Module

A signal model over parameter maps, with its derivatives.

Maps are stacked on the last axis, as everywhere else in TorchSim: (..., channels) in, (..., contrasts) out, every leading axis a voxel axis. physics() is where that flips to the channel-first convention deepinv and mri-nufft use, and it is the only place it flips.

Parameters:

acquisition (Simulator) – The sequence being inverted, with every property that is not being solved for already fixed on it. A property bound as a map – a measured B1, a known T1 – is one value per voxel and rides along.
unknown (str, optional) – The property names being solved for, in the order their channels appear. At least one.
bounds (mapping, optional) – {name: (low, high)}, either end None for unbounded. A bound is kept by solving for a transformed variable, so no iterate is ever outside it – which matters more here than in a fit, because the model is evaluated at every voxel to predict every k-space sample and one voxel out of range corrupts the whole residual. A two-sided bound also puts the parameter on a scale of order one whatever its units, which is the preconditioning a mixed parameter set needs.
scale (mapping, optional) – {name: value}, the size of a step in a parameter left unbounded. Bounded parameters are already scaled by their interval.
amplitude (bool, optional) – Whether to carry a complex amplitude multiplying the model output. On, which is what a reconstruction wants; off for a fit whose model already exposes its own proton density.
subspace (Subspace, optional) – A Subspace the prediction is projected through, for a reconstruction that solves in the temporal basis rather than in the contrasts.

Examples

operator = ModelOperator(
    MultiEchoSimulator(TE=echo_times, T1=1000.0),
    "T2",
    bounds={"T2": (10.0, 300.0)},
)
maps = operator.initial((128, 128), T2=80.0)
images = operator.A(maps)          # (128, 128, echoes), complex

Notes

The operator holds nothing on a device of its own. It follows the maps it is called with, taking the sequence and any bound tissue along, so there is nothing here to move and to() moves nothing.

Equality constraints are written into the model, not declared here. A constraint that fixes one parameter in terms of the others removes a degree of freedom, so the way to impose it is to not have that freedom: write the model on the parameters that remain. For a fat-water separation whose two fractions must sum to one, make the fat fraction f the only unknown and write water as 1 - f inside the model. The constraint then holds identically at every Gauss-Newton iterate, to the bit, rather than being restored after each one.

Raises:: ValueError – If nothing is unknown, if a bound or a scale names something that is not unknown, or if a bound is not increasing.

Methods

`A`	The name a reconstruction knows the forward operator by.
`A_jvp`	The directional derivative `J d`, in one forward-mode pass.
`A_vjp`	The adjoint product `J^H v`, in one reverse-mode pass.
`initial`	Maps to start from, as the variables actually solved for.
`jacobian`	The derivative block of every voxel, built.
`physics`	This operator as a `deepinv.physics.Physics`.
`select`	This operator over a subset of the voxels.
`split`	The maps `x` stands for, in their own units.

property channels#: How many map channels x carries.

property names#: What each channel of x is, in order.

split(x)[source]#

The maps x stands for, in their own units.

Parameters:: x (torch.Tensor) – (..., channels), the variables being solved for.
Returns:: One entry per unknown, inside its bound, plus "amplitude" complex where one is carried.
Return type:: dict

initial(shape=(), **values)[source]#

Maps to start from, as the variables actually solved for.

Parameters:

shape (sequence of int, optional) – The voxel shape. Empty gives one voxel.
values (float or array-like, optional) – {name: value} in the property’s own units, one value or one per voxel – a map broadcasts against shape, which is what turns known parameters into a state and makes this the inverse of split(). An unknown left out starts at the middle of its bound, and a one-sided or absent bound has no middle, so it must be given. The amplitude starts at one unless amplitude= says otherwise.

Returns:

(*shape, channels).

Return type:

torch.Tensor

Raises:

ValueError – If a property with no two-sided bound is left out, or if a value sits on a bound – a value on a bound has no unconstrained image.

select(keep)[source]#

This operator over a subset of the voxels.

A property bound as a map – a measured B1, a known efficiency – has one value per voxel, so a loop that retires the voxels it has finished has to take the maps with it. Anything that is not one value per voxel is left alone.

Parameters:: keep (torch.Tensor) – An index or a boolean mask over the voxel axis.
Returns:: A copy over those voxels.
Return type:: ModelOperator

forward(x)[source]#

The contrasts these maps record.

Parameters:: x (torch.Tensor) – (..., channels).
Returns:: (..., contrasts), complex.
Return type:: torch.Tensor

A(x)#: The name a reconstruction knows the forward operator by.

A_jvp(x, d)[source]#

The directional derivative J d, in one forward-mode pass.

One pass whatever the channel count – the step itself is the tangent, so nothing is built column by column.

Parameters:

x (torch.Tensor) – (..., channels), the point and the direction.
d (torch.Tensor) – (..., channels), the point and the direction.

Returns:

(..., contrasts), complex.

Return type:

torch.Tensor

A_vjp(x, v)[source]#

The adjoint product J^H v, in one reverse-mode pass.

Real, because the maps are: what comes back is the real part of the conjugate product, which is the gradient direction of the squared modulus and the adjoint of A_jvp() under the real inner product.

Parameters:

x (torch.Tensor) – (..., channels), the point.
v (torch.Tensor) – (..., contrasts), the cotangent.

Returns:

(..., channels), real.

Return type:

torch.Tensor

jacobian(x)[source]#

The derivative block of every voxel, built.

Wanted where the blocks are solved outright, which is the voxel-wise fit; a Gauss-Newton step under an encoding operator wants A_jvp() and A_vjp() instead and never sees a block.

Parameters:: x (torch.Tensor) – (..., channels).
Returns:: (..., channels, contrasts), complex.
Return type:: torch.Tensor

physics(**kwargs)[source]#

This operator as a deepinv.physics.Physics.

The wrapper is where the map axis moves to the front, because that is the convention deepinv and mri-nufft share: (batch, channels, *xyz) in, (batch, contrasts, *xyz) out. Composing it with an encoding operator – encoding * operator.physics() – gives a ComposedPhysics every deepinv optimizer and prior applies to.

Parameters:: kwargs (dict, optional) – Passed to deepinv.physics.Physics, so a noise model can be given.
Return type:: deepinv.physics.Physics
Raises:: ImportError – If deepinv is not installed. TorchSim does not depend on it.

Notes

A ComposedPhysics takes its Jacobian products by automatic differentiation through the whole chain, so composing this way gives up the analytic derivative below. GaussNewton chains the two operators’ own products instead, and keeps it.

ModelOperator

Contents

ModelOperator#