.. _concatenate:

==============================
Concatenating linear operators
==============================

It is possible to concatenate linear operators in a way that imitates the
``numpy.concatenate`` function. The only concatenations that are allowed are
in the vertical and horizontal directions.

The ``talon.concatenate`` function requires an iterable containing the linear
operators to concatenate and the axis along which they have to be concatenated.

The following code shows the correct syntax to concatenate two linear operators
:math:`A` and :math:`B` vertically and horizontally:

.. code:: python

    V = talon.concatenate((A, B), axis=0) # vertical (default)
    H = talon.concatenate((A, B), axis=1) # horizontal

which correspond to the following

.. math::
    V = \begin{bmatrix} A \\ B \end{bmatrix} \qquad
    H = \begin{bmatrix} A & B \end{bmatrix}.


Examples
------------------------
Build a tractogram with two crossing bundles of fibers and the corresponding
linear operator.

.. code:: python

    import numpy as np
    import talon

    from scipy.interpolate import interp1d

    # Set seed for reproducibility
    np.random.seed(1992)

    # The number of voxels in each dimension of the output image.
    image_size = 25
    center = image_size // 2

    n_points = int(image_size / 0.01)
    t = np.linspace(0, 1, n_points)
    
    streamlines_per_bundle = 50

    def generate_crossing_tractogram():
        tractogram = []

        horizontal_points = np.array([[0, center, center],
                                     [image_size - 1, center, center]])
        horizontal_streamline = interp1d([0, 1], horizontal_points, axis=0)(t)

        for k in range(streamlines_per_bundle):
            new_streamline = horizontal_streamline.copy()
            new_streamline[:,1] += (np.random.rand(1) - 0.5)
            tractogram.append(new_streamline)

        vertical_points = np.array([[center, 0, center],
                                   [center, image_size - 1, center]])
        vertical_streamline = interp1d([0, 1], vertical_points, axis=0)(t)

        for k in range(streamlines_per_bundle):
            new_streamline = vertical_streamline.copy()
            new_streamline[:,0] += (np.random.rand(1) - 0.5)
            tractogram.append(new_streamline)
        return tractogram

    cross_tractogram = generate_crossing_tractogram()
    directions = talon.utils.directions(1000)
    generators = np.ones((len(directions), 1))
    image_shape = (image_size,) * 3
    indices, lengths = talon.voxelize(cross_tractogram, directions, image_shape)

    A = talon.operator(generators, indices, lengths)


Vertical concatenation
^^^^^^^^^^^^^^^^^^^^^^
If multiple features for each streamline are encoded in different linear
operators we can concatenate different linear operators vertically. If :math:`A`
encodes the linear operator for the set of streamlines :math:`\alpha` and
generators :math:`G_1` and :math:`B` encodes the linear operator for the same
streamlines but with generators :math:`G_2`, instead of rebuilding the linear
operator from scratch we can concatenate :math:`A` and :math:`B` vertically
to obtain the same result.

.. code:: python

    G2 = np.random.rand(len(directions), 5) # New generators
    B = talon.operator(G2, indices, lengths)

    V = talon.concatenate((A,B), axis=0)

    print('Shape of A: {}'.format(A.shape))
    print('Shape of B: {}'.format(B.shape))
    print('Shape of V: {}'.format(V.shape))
    print('Check: {} + {} = {}'.format(A.shape[0], B.shape[0], A.shape[0] + B.shape[0]))

Notice that the :code:`axis=0` argument is redundant since it is the default.

The output should be the following:

.. code::

    Shape of A: (15625, 100)
    Shape of B: (78125, 100)
    Shape of V: (93750, 100)
    Check: 15625 + 78125 = 93750



Horizontal concatenation
^^^^^^^^^^^^^^^^^^^^^^^^
One (but not the only) reason to concatenate two linear operators horizontally
is to add a set of streamlines to the system. If :math:`A` encodes the linear
operator for the set of streamlines :math:`\alpha` and :math:`C` for set
:math:`\beta`, instead of rebuilding the linear operator from scratch we can
concatenate :math:`A` and :math:`C` horizontally to obtain the same result.

.. code:: python

    def generate_diagonal_tractogram():
        tractogram = []
        diagonal_points = np.array([[0, center, center],
                                   [center, image_size - 1, center]])
        diagonal_streamline = interp1d([0, 1], diagonal_points, axis=0)(t)

        for k in range(streamlines_per_bundle):
            new_streamline = diagonal_streamline.copy()
            new_streamline[:,0] += (np.random.rand(1) - 0.5)
            new_streamline[:,1] += (np.random.rand(1) - 0.5)
            tractogram.append(new_streamline)
        return tractogram

    diag_tractogram = generate_diagonal_tractogram()
    indices, lengths = talon.voxelize(diag_tractogram, directions, image_shape)

    C = talon.operator(generators, indices, lengths) # diagonal

The concatenation of the two linear operators is performed as follows:

.. code:: python

    H = talon.concatenate([A, C], axis=1)
    print('Shape of A: {}'.format(A.shape))
    print('Shape of C: {}'.format(C.shape))
    print('Shape of H: {}'.format(H.shape))

The output should be the following:

.. code::

    Shape of A: (15625, 100)
    Shape of C: (15625, 50)
    Shape of H: (15625, 150)


The matrix multiplication and transposition operations work as usual:

.. code:: python

    x = H @ np.random.rand(H.shape[1])
    y = H.T @ np.random.rand(H.shape[0])

    print('Shape of x: {}'.format(x.shape))
    print('Shape of y: {}'.format(y.shape))

The output should be the following:

.. code::

    Shape of x: (15625,)
    Shape of y: (150,)