# mimo

multiple-input-multiple-output

## Contents

- Syntax
- Description
- Example
- Input Arguments
- Output Arguments
- Best Practice
- More About
- See Also
- Code: mimo
- Code: mimocalc

## Syntax

    [FFqmF,ADFFq] = mimo(seti,qROI,FmeasDelta)
    [FFqmF,ADFFq] = mimo(seti,qROI,FmeasDelta,'adjOfDer')
    [FFqMeas,~,FFqROI] = mimo(seti,qROI,'simo')
    [FFqMeas,ADFFq,FFqROI] = mimo(seti,qROI,'simo')
    [JA,JB] = mimo(seti,qROI,'jacobian')

## Description

Adjoint of derivative ('adjOfDer')

    [FFqmF, ADFFq] = mimo(seti,qROI,FmeasDelta)

and

    [FFqmF, ADFFq] = mimo(seti,qROI,FmeasDelta,'adjOfDer')

- computes the difference of the evalutated forward operator and the measured data, `FFqmF`, and
- computes the adjoint of derivative, `ADFFq`.
- You can leave `'adjOfDer'` because it is used automatically if the last argument `op` is empty.

Single-input-multiple-output ('simo')

    [FFqMeas,~,FFqROI] = mimo(seti,qROI,'simo')

- computes the scattered field evaluated on receivers' positions for each transmitter, `FFqMeas`, and
- the scattered field on ROI for each transmitter, `FFqROI`.
- ~ is the adjoint of derivative, `ADFFq`, but this is uninteresting because it is 0. (This is expected, because `simo` is used with input of FmeasDelta = 0 in `mimocalc`.)

Auxilary matrices JA and JB to compute the Jacobian matrix ('jacobian')

    [JA,JB] = mimo(seti,qROI,'jacobian')

computes auxiliary matrices JA and JB, such that the Frechet derivative of at evaluated on , i.e. , that is in finite-dimensional spaces is represented by the Jacobian matrix of at evaluted on can easily computed by

## Example

We give one example to compute the difference of the evalutated forward operator and the measured data, `FFqmF`, and the adjoint of derivative, `ADFFq`.

    init;
    seti.contrast = 'corner2D';
    seti = setGeomSim(seti);
    qROI = seti.qROIexact;
    FmeasDelta = zeros(seti.measNb,seti.incNb);
    [FFqmF,ADFFq] = mimo(seti,qROI,FmeasDelta,'adjOfDer');
    figure(1); imagesc(real(FFqmF)); axis xy; colorbar;
    figure(2); imagesc(real(seti.G(ADFFq))); axis xy; colorbar;

## Input Arguments

- `seti` : struct, details see below.
- `qROI` : contrast in ROI (region of interest) (matrix as vector of size seti.nROI^seti.dim x 1)
- `FmeasDelta` : measurement data (on receivers) with noise level (matrix of size seti.measNb x seti.incNb)

struct `seti`

As in convenience functions adjOfDer.html, forward.html and derivative.html all necessary fields (and more) are set by

    seti = setGeomSim(seti);

A list of all necessary (and not more) fields is below.

seti in setGeomSim.html:

- `seti.mCD` : internal parameter, set seti.mCD = 0 (no coarse grid is used)
- `seti.tol` : internal parameter, set 1E-6, is used in simo and solveLippmannSchwinger

seti in setGrid.html:

- `seti.dim` : dimension of the problem (2 or 3)
- `seti.nCD` : number of discretization points for each dimension of computational domain (CD)
- `seti.nROI` : discretization points for each dimension of region of interest (ROI)
- `seti.gridROI` : grid of region of interest (ROI) (matrix of size seti.dim x seti.nROI^seti.dim)
- `seti.ROImask` : biggest square in 2D (cube in 3D) in ball with radius rCD/2 (logical array of size seti.nCD x seti.nCD)
- `seti.dV` : pixel/voxel area/volume (number)

seti in setKernel.html:

- `seti.model` : model (helmholtz2D or helmholtz3D)
- `seti.k` : wave number
- `seti.kHat` : Fourier coefficients of discretized integral operator (seti.nCD x seti.nCD)

seti in expSetup.html (subfiles pntsGeometry.html, setIncField.html, setMeasKer.html):

- `seti.measNb` : number of receivers (number of measurements)
- `seti.incNb` : number of transmitters (number of incident fields)
- `seti.dSInc` : approximation of the infinitesimal element of closed contour with control points in case of transmitters
- `seti.dSMeas` : approximation of the infinitesimal element of closed contour with control points in case of receivers
- `seti.incField` : incident field, setIncField.html (matrix of size seti.nROI^seti.dim x seti.incNb)
- `seti.measKer` : matrix to compute simulated measurement data from solution, setMeasKer.html (size seti.measNb x seti.nROI^seti.dim)

## Output Arguments

More details of output arguments are in section More About. type: vector, matrix, etc...

- `FFqMeas` : scattered field evaluated on the receivers' positions for each transmitter (complex matrix of size seti.measNb x seti.incNb)
- `FFqmF` : `FFqMeas` minus noisy data (complex matrix of size seti.measNb x seti.incNb)
- `ADFFq` : adjoint of derivative (complex matrix written as vector of size seti.nROI^seti.dim x 1)
- `FFqROI` : scattered field evaluated on ROI (region of interest) stored as vector for each transmitter (complex matrix of size seti.nROI^seti.dim x seti.incNb)
- `JA`, `JB` : Auxilliary matrices (JA is a complex matrix of size seti.measNb x seti.nROI^seti.dim, JB is a complex matrix of size seti.nROI^seti.dim x seti.incNb)

## Best Practice

If you are only interested in the scattered field evaluated on receivers' positions (i.e. measurements), `FFqMeas`, use

    [FFqMeas,~,~] = mimo(seti,qROI,'simo')

instead of

    [FFqMeas,~] = mimo(seti,qROI,0,'adjOfDer')

or

    [FFqMeas,~] = mimo(seti,qROI,zeros(seti.measNb,seti.incNb),'adjOfDer')

The latter with `'adjOfDer'` computes the adjoint of derivative additionally, which is more expensive than `'simo'`.

## More About

- : contrast-to-measurement operator
- is the matrix of TX-RX scattered fields (TX: transmitters, RX: receivers), i.e. multiple-input-multiple-output setting, i.e. the k-the column is the RX-field for TX(k).
- , i.e. with noise level noisy simulated data (measurement) .
- : discrepancy
- is the scattered field evaluated on receivers' positions (measurements) (also called ). Note that .
- .

- is the scattered field evaluated on ROI stored for each transmitter for each transmitter (for one transmitter it is also called ).
- Scattered field with volume potential .
- Remind: FFqROI is on ROI, but FFqmF is on RX (receivers positions).

- is the Frechet derivative with auxiliary matrices and .
- is the adjoint of derivative.

- `[JA,JB] = mimo(seti, qROI,'jacobian')` is used in case of primal-dual algorithm.
- `[FFqmF,ADFFq] = mimo(seti,qROI,seti.FmeasDelta)` is used in case of soft-shrinkage (not available in public version).

## See Also

For easy usage we offer some convenience functions for mimo:

- adjOfDer: adjOfDer.html
- simo: forward.html
- jacobian: derivative.html

## Code: mimo

    function varargout = mimo(seti, qROI, varargin)

    if ~ischar(varargin{end})
        op = 'adjOfDer';
    else
        op = varargin{end};
    end

    if strcmpi(op, 'adjOfDer') % old name: gradient
        [varargout{1}, varargout{2}] = mimocalc(seti, qROI, varargin{1}, op);
        % Note: varargin{1} = FmeasDelta.
        % Output documentation:
        % varargout{1} = FFqmF
        % varargout{2} = ADFFq

    elseif strcmpi(op, 'simo')
        [varargout{1}, varargout{2}, varargout{3}] = mimocalc(seti, qROI, 0, op);
        % Output documentation:
        % varargout{1} = uScattRX = FFqMeas
        %     (and not FFqMeas - FmeasDelta, because FmeasDelta was set to 0).
        % varargout{2} = ADFFq = 0 (in this case)
        % varargout{3} = uScattROI = FFqROI

    elseif strcmpi(op, 'jacobian')
        if strcmp(seti.model, 'helmholtz3D') ...
                || strcmp(seti.model, 'helmholtz2D') ...
                || strcmp(seti.model, 'helmholtzHMode2D')
            [~, ~, ~,varargout{1},varargout{2}] = mimocalc(seti, qROI, 0, op);
            % Output documentation:
            % varargout{1} = JA
            % varargout{2} = JB
        else
            disp('mimo.m: Error - Jacobian is not (yet ... ?) implemented');
            varargout{1} = [];
            varargout{2} = [];
        end
    end

    end

## Code: mimocalc

    function [FFqmF, ADFFq, FFqROI, JA, JB] = mimocalc(seti, qROI, FmeasDelta, op)
    %Note: HMode is not available in public version.

    nnROI = length(seti.gridROI(1,:));
    ADFFq = 0;
    FFqROI = zeros(nnROI, seti.incNb);
    FFqmF = zeros(seti.measNb, seti.incNb);

    if strcmpi(op, 'adjOfDer')
        ADFFq = zeros(size(qROI));
    end
    if strcmpi(op, 'jacobian')
        if strcmp(seti.model, 'helmholtz3D') || strcmp(seti.model, 'helmholtz2D')
            JA = zeros(seti.measNb, nnROI); JB = zeros(nnROI, seti.incNb);
        elseif strcmp(seti.model, 'helmholtzHMode2D')
            JA = zeros(seti.measNb, nnROI,2); JB = zeros(nnROI, seti.incNb,2);
        else
            disp('mimo.m: Error - Jacobian is not (yet ... ?) implemented');
        end
    end

    dSInc = seti.dSInc.*ones(seti.incNb,1);

    [~, S, VGStar, VStar, TStar, ~, QStarU, ~, ~] = intOpsFuncs(qROI,seti);

    for jj = 1:seti.incNb
        SL = dSInc(jj)*seti.incField(:,jj);
        % SL = uIncROI (i.e. single layer potential/Herglotz wave function)
        % SL \phi = uIncROI
        [uMeas, uScattROI] = S(SL); % uMeas is uScattRX (i.e. u^s on measPnts)

        h = uMeas;

        if strcmpi(op, 'adjOfDer')
            h = h - FmeasDelta(:,jj); % result: h = FFqmF = FF(q) - FmeasDelta
        end;
        FFqmF(:,jj) = h; % In case of adjOfDer (default): FFqmF = FF(q) - FmeasDelta
                         % otherwise so just FFq:         FFqmF = FF(q)
        h( isnan(h) ) = 0; % workaround
        uTotROI = uScattROI + SL; % uTotROI = uScattROI + uIncROI
        if strcmpi(op, 'adjOfDer')
            if strcmp(seti.model, 'helmholtz3D') || strcmp(seti.model, 'helmholtz2D')
                f = VGStar(h);
                v = conj(uTotROI).*(f + VStar(TStar(QStarU(f)))); %VStar NOT VqStar
                ADFFq = ADFFq + v;
                %figure(101); x = ADFFq; imagesc(real(extendROItoCD(x,seti.ROImask))); axis xy; colorbar; caxis([-3,3]); error('stop');
            elseif strcmp(seti.model, 'helmholtzHMode2D')
                ADFFq = ADFFhelmHMode2D(h,seti,uTotROI);
            else
                disp('mimo.m: Error - model is not (yet ... ?) implemented');
            end
        end
        if strcmpi(op, 'simo')
            FFqROI(:,jj) = reshape(uScattROI,[nnROI 1]);
            ADFFq = 0;
        end
        if strcmpi(op, 'jacobian')
            if strcmp(seti.model, 'helmholtz3D') || strcmp(seti.model, 'helmholtz2D')
                JB(:,jj) = uTotROI;
            elseif strcmp(seti.model, 'helmholtzHMode2D')
                [uTotROI1, uTotROI2] = gradientHMode(uTotROI,seti);
                JB(:,jj,1) = uTotROI1;
                JB(:,jj,2) = uTotROI2;
            else
                disp('mimo.m: Error - Jacobian is not (yet ... ?) implemented');
            end
        end
    end
    if strcmpi(op, 'jacobian')
        dSMeas = seti.dSMeas.*ones(seti.measNb,1);
        if strcmp(seti.model, 'helmholtz3D') || strcmp(seti.model, 'helmholtz2D')
            for jj = 1:seti.measNb
                f = VGStar( ((1:seti.measNb)==jj)' * seti.dV/dSMeas(jj) );
                JA(jj,:) = conj(f + VStar(TStar(QStarU(f)))); %VStar NOT VqStar
            end;
        elseif strcmp(seti.model, 'helmholtzHMode2D')
            for jj = 1:seti.measNb
                [f1,f2] = helmholtzHMode2DZStar( ((1:seti.measNb)==jj)' * seti.dV/dSMeas(jj), seti); % Z*(g_j)
                qf1 = conj(seti.qROIexact).*f1; qf2 = conj(seti.qROIexact).*f2;

    %             qf1Ext = 1.0i*pi/seti.rCD * fft2(extendROItoCD(qf1,seti.ROImask)) * seti.gradMat;
    %             qf2Ext = 1.0i*pi/seti.rCD * seti.gradMat * fft2(extendROItoCD(qf2,seti.ROImask));
    %             qf1 = restrictCDtoROI(ifft2(qf1Ext+qf2Ext),seti.ROImask);

                qf1 = divHMode(seti,qf1,qf2);

                qf1 = TStar(qf1);

                [qf1, qf2] = helmholtzHMode2DNablaVStar(qf1,seti);

                JA(jj,:,1) = conj(f1 + qf1); JA(jj,:,2) = conj(f2 + qf2);

            end;
        else
            disp('mimo.m: Error - Jacobian is not (yet ... ?) implemented');
        end
    end

    end