code: def sequential_gaussian_filter_sample(
      init: Gaussian,
      trans: Gaussian,
      sample_shape: Tuple[int, ...] = (),
      noise: Optional[torch.Tensor] = None,
  ) -> torch.Tensor:
      """
      Draws a reparameterized sample from a Markov product of Gaussians via
      parallel-scan forward-filter backward-sample.
  
      :param Gaussian init: A Gaussian representing an initial state.
      :param Gaussian trans: A Gaussian representing as series of state transitions,
          with time as the rightmost batch dimension. This must have twice the event
          dim as ``init``: ``trans.dim() == 2 * init.dim()``.
      :param tuple sample_shape: An optional extra shape of samples to draw.
      :param torch.Tensor noise: An optional standard white noise tensor of shape
          ``sample_shape + batch_shape + (duration, state_dim)``, where
          ``duration = 1 + trans.batch_shape[-1]`` is the number of time points
          to be sampled, and ``state_dim = init.dim()`` is the state dimension.
          This is useful for computing the mean (pass zeros), varying temperature
          (pass scaled noise), and antithetic sampling (pass ``cat([z,-z])``).
      :returns: A reparametrized sample of shape
          ``sample_shape + batch_shape + (duration, state_dim)``.
      :rtype: torch.Tensor
      """
      assert isinstance(init, Gaussian)
      assert isinstance(trans, Gaussian)
      assert trans.dim() == 2 * init.dim()
      state_dim = trans.dim() // 2
      batch_shape = broadcast_shape(
          trans.batch_shape[:-1], init.batch_shape
      )
      if init.batch_shape != batch_shape:
          init = init.expand(batch_shape)
      dtype = trans.precision.dtype
      device = trans.precision.device
      perm = torch.cat(
          [
              torch.arange(1 * state_dim, 2 * state_dim, device=device),
              torch.arange(0 * state_dim, 1 * state_dim, device=device),
              torch.arange(2 * state_dim, 3 * state_dim, device=device),
          ]
      )
  
      # Forward filter, similar to sequential_gaussian_tensordot().
      tape = []
      shape = trans.batch_shape[:-1]  # Note trans may be unbroadcasted.
      gaussian = trans
      while gaussian.batch_shape[-1] > 1:
          time = gaussian.batch_shape[-1]
          even_time = time // 2 * 2
          even_part = gaussian[..., :even_time]
          x_y = even_part.reshape(shape + (even_time // 2, 2))
          x, y = x_y[..., 0], x_y[..., 1]
          x = x.event_pad(right=state_dim)
          y = y.event_pad(left=state_dim)
          joint = (x + y).event_permute(perm)
          tape.append(joint)
          contracted = joint.marginalize(left=state_dim)
          if time > even_time:
              contracted = Gaussian.cat(
                  (contracted, gaussian[..., -1:]), dim=-1
              )
          gaussian = contracted
      gaussian = gaussian[..., 0] + init.event_pad(right=state_dim)
  
      # Generate noise in batch, then allow blocks to be consumed incrementally.
      duration = 1 + trans.batch_shape[-1]
      shape = torch.Size(sample_shape) + init.batch_shape
      result_shape = shape + (duration, state_dim)
      noise_stride = (
          shape.numel() * state_dim
      )  # noise is consumed in time blocks
      noise_position: int = 0
      if noise is None:
          noise = torch.randn(result_shape, dtype=dtype, device=device)
      assert noise.shape == result_shape
  
      def rsample(
          g: Gaussian, sample_shape: Tuple[int, ...] = ()
      ) -> torch.Tensor:
          """Samples, extracting a time-block of noise."""
          nonlocal noise_position
          assert noise is not None
          numel = torch.Size(
              sample_shape + g.batch_shape + (g.dim(),)
          ).numel()
          assert numel % noise_stride == 0
          beg: int = noise_position
          end: int = noise_position + numel // noise_stride
          assert end <= duration, "too little noise provided"
          noise_position = end
          return g.rsample(sample_shape, noise=noise[..., beg:end, :])
  
      # Backward sample.
      result = rsample(gaussian, sample_shape).reshape(
          shape + (2, state_dim)
      )
      for joint in reversed(tape):
          # The following comments demonstrate two example computations, one
          # EVEN, one ODD.  Ignoring sample_shape and batch_shape, let each zn be
          # a single sampled event of shape (state_dim,).
          if joint.batch_shape[-1] == result.size(-2) - 1:  # EVEN case.
              # Suppose e.g. result = [z0, z2, z4]
              cond = result.repeat_interleave(
                  2, dim=-2
              )  # [z0, z0, z2, z2, z4, z4]
              cond = cond[..., 1:-1, :]  # [z0, z2, z2, z4]
              cond = cond.reshape(
                  shape + (-1, 2 * state_dim)
              )  # [z0z2, z2z4]
              sample = rsample(joint.condition(cond))  # [z1, z3]
              zipper = result.new_empty(
                  shape + (2 * result.size(-2) - 1, state_dim)
              )
              zipper[..., ::2, :] = result  # [z0, _, z2, _, z4]
              zipper[..., 1::2, :] = sample  # [_, z1, _, z3, _]
              result = zipper  # [z0, z1, z2, z3, z4]
          else:  # ODD case.
              assert joint.batch_shape[-1] == result.size(-2) - 2
              # Suppose e.g. result = [z0, z2, z3]
              cond = result[..., :-1, :].repeat_interleave(
                  2, dim=-2
              )  # [z0, z0, z2, z2]
              cond = cond[..., 1:-1, :]  # [z0, z2]
              cond = cond.reshape(shape + (-1, 2 * state_dim))  # [z0z2]
              sample = rsample(joint.condition(cond))  # [z1]
              zipper = result.new_empty(
                  shape + (2 * result.size(-2) - 2, state_dim)
              )
              zipper[..., ::2, :] = result[
                  ..., :-1, :
              ]  # [z0, _, z2, _]
              zipper[..., -1, :] = result[..., -1, :]  # [_, _, _, z3]
              zipper[..., 1:-1:2, :] = sample  # [_, z1, _, _]
              result = zipper  # [z0, z1, z2, z3]
  
      assert noise_position == duration, "too much noise provided"
      assert result.shape == result_shape
      return result
  
imports: import math
  from typing import Optional, Tuple, Union
  import torch
  from torch.distributions.utils import lazy_property
  from torch.nn.functional import pad
  from pyro.distributions.util import broadcast_shape
  from pyro.ops.tensor_utils import matmul, matvecmul, safe_cholesky, triangular_solve
label_desc: Use image node_id 6 for calling a global function with example usage: sequential_gaussian_filter_sample(init, trans, sample_shape, noise) and returns: result, g