# setFuncsPda

Set functions required in primal-dual algorithm (pda.m).

## Contents

- Syntax
- Description
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code: setFuncsPda
- Code: subfunction: potPotNum
- Code: subfunction: setminfParts
- Code: subfunction: setminfPartsWavelet

## Syntax

    seti = setFuncsPda(seti)

## Description

`seti = setFuncsPda(seti)` sets functions in structure array `seti`, that are required in primal-dual algorithm, see pda.html. Therefore this function is called in setRecon.m, see setRecon.html.

## Input Arguments

- seti : structure array

## Output Arguments

- seti : structure array

For the fields of the structrual array `seti` see the section "More About".

## More About

We only give a overview of the functions. See also pda.html. Detailed information are in Section 4.5 in [1].

Note that grid scaling is respected in the code, but for simplicitiy is not presented in the formulas.

Auxiliary quantities

For and , see "More About" in pda.html.

- ,
- ,
- .

Note that the names `Kg` and `vg` were choosen because the appearance of the gradient.

Parts of functional to minimize

- is the discrepancy (linearized problem),
- is the sparsity penalty,
- is the TV-penalty,
- is the penalty for physical bounds.

Note that is a indicator function, i.e. is 0, if all entries of the vector are between and and all entries of the vector are between and ; and is , otherwise.

Definition of parts of Tikhonov functional

The Tikhonov functional is `MT = M1 + M2`. Function is , function is (with and from PDA) This is done in function `minPda` (minPda.html) when function `pda` (pda.html) is called.

Further comments (in particular, interesting in not public version)

- In case of inversion by pda, it is always (automatically) set `seti.pNorm = 2` and `seti.qNorm = 1` (except for invNo = 7), see setInvType.html.
- Exponent p of the discrepancy term, see setInvType.html, is usually `seti.p = 2`. In case of shrinkage or pda with wavelets it may differ, see setInvType.html.
- `vs` and `vp` in term : If you choose `vs` and `vp` manually and `fs` and `fp` belong to term , then `vs = vp` is required.

## References

- [1] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. A sparsity regularization and total variation based computational framework for the inverse medium problem in scattering. Journal of Computational Physics, 339:1-30, 2017. URL: https://doi.org/10.1016/j.jcp.2017.03.011.

## See Also

- setRecon.html
- pda.html
- setInvType.html

## Code: setFuncsPda

    function seti = setFuncsPda(seti)

    % Definition of parts of function to minimize (f = minf)
    % parts for reconstruction process
    seti = setminfParts(seti);
    if seti.useWavelet == 1 % not available in public version
        seti = setminfPartsWavelet(seti);
    end

    end

## Code: subfunction: potPotNum

Numerical potential pot function.

    function [d1,d2] = potPotNum(seti)
    infty = 1E5;
    dp = @(a,b,x) 0.* and((a<=x),(x<=b))+infty.*or((x<a),(b<x)); % potential pot function
    d =  @(a,b,x) max(dp(a,b,x)); % d = 0 if ALL entries are between a and b, otherwise d = \infty
    % max because d = \infty if ONE entry of x \notin [a,b]
    clear infty

    l = seti.physBounds;
    reMin = l(1); reMax = l(2); imMin = l(3); imMax = l(4);
    tol = seti.tolDelta;
    d1 = @(x) d(reMin-tol,reMax+tol,x); % limits for real part
    d2 = @(x) d(imMin-tol,imMax+tol,x); % limits for imag part
    % tol important because e.g. (min(real(qROI))+1) = -2.2204e-16
    clear l a b c
    end

## Code: subfunction: setminfParts

- Y: depends on function `Ki` (`Ki` i.e. `Kd`, `Kg`, ...)
- X: xnRVD (this means real vector downscaled): size of elements: 2*seti.nInv^seti.dim
- `hs` was `xnRVU`
- `q` is standard in qCVU (is fixed in pda), but you can use qCVD etc. too, if you have them as argument

    function seti = setminfParts(seti)

f_d

- Remember that `DFFq = @(xnRVD) JA*diag(seti.GU(seti.T(xnRVD)))*JB;`
- The difference to DFFq is that finally `seti.S` is used.

    if strcmp(seti.model,'helmholtz2D') || strcmp(seti.model,'helmholtz3D')
        seti.Kd = @(xnRVD,JA,JB) seti.S(JA*diagsparse(seti.GU(seti.T(xnRVD)))*JB); % K_d(h) = S(\F'(q)[GU(T(h_S))]) %real
    elseif strcmp(seti.model,'helmholtzHMode2D')
        seti.Kd = @(xnRVD,JA,JB) seti.S( JA(:,:,1)*diagsparse(seti.GU(seti.T(xnRVD)))*JB(:,:,1) ...
            + JA(:,:,2)*diagsparse(seti.GU(seti.T(xnRVD)))*JB(:,:,2) ); % K_d(h) = S(\F'(q)[GU(T(h_S))]) %real
    else
        fprintf('Error in setFuncs - pda not implemented for model %s.\n',  seti.model)
    end
    % Finally: Kd(x): R^{2 nInv^dim} -> R^{2 measNb x incNb}
    % Kd = \F'(q) (from mimo \F'(q): JA*diag(h)*JB = \F'(q)[h])
    seti.KdAdj = @(yd,JA,JB) seti.S(seti.GD(ADFFqFast(seti.T(yd),JA,JB,seti)));
    seti.vd = @(FFq) seti.S(FFq-seti.FmeasDelta); % 2* measNb x incNb
    seti.fd = @(DFFqh,FFq) (1/2)*normws(seti.S(DFFqh)+seti.vd(FFq),seti)^2;

f_s

    seti.vs = @(qCVU) seti.S(seti.GD(qCVU));
    seti.fs = @(qCVU,xnRVD) seti.alpha*normroi(xnRVD+seti.vs(qCVU),seti);

f_g

    seti.Kg = @(xnRVD) seti.beta*gradientNeumann(xnRVD,seti.hInv,seti.nInv,seti.GInv,seti); % K_g
    %seti.KgAdj = @(yg) -seti.beta*div(yg,seti); % Note that this is wrong(!) (although analytically correct)
    seti.KgAdj = @(yg) seti.beta*gradientNeumannAdj(yg,seti); %yg is 2*dim x nInv x nInv (x nInv)
    seti.vg = @(qCVU) seti.Kg(seti.S(seti.GD(qCVU)));
    seti.fg = @(qCVU,xnRVD) normTVinv1(seti.Kg(xnRVD)+seti.vg(qCVU),seti); %f_g(x) = ||beta*\nabla x + v_g)||_\tv (gradient sparsity)

f_p

- Note that seti.vp = seti.vs is important because of prox Gsp (prox Gsp does mean that the term G consists of fs and fp. To compute the proximal mapping of G the equality of vp and vs is required!)

    [d1,d2] = potPotNum(seti);
    %seti.fp = @(q,hz) d1(real(q)+real(hz)) + d2(imag(q)+imag(hz));
    % use seti.vs in fp:
    %seti.fp = @(q,hz) d1(seti.R(seti.vs(q))+seti.R(seti.vs(hz))) + d2(seti.I(seti.vs(q))+seti.I(seti.vs(hz)));
    seti.fp = @(qCVU,xnRVD) d1(real(seti.GD(qCVU))+seti.R(xnRVD)) + d2(imag(seti.GD(qCVU))+seti.I(xnRVD));

    seti.fp2 = seti.fp; % In fact it is not the same function, but because the pixels in the background are zero the result is the same.

    end

## Code: subfunction: setminfPartsWavelet

Definitions of functions in case of wavelets. (This is not supported in the public version.)

    function seti = setminfPartsWavelet(seti)

    % f_dw
    % seti.Kdw = seti.Kd;
    %seti.KdwAdj = seti.KdAdj;
    seti.vdw = @(FFq) seti.S(FFq-seti.FmeasDelta);
    %seti.fdw = @(DFFqh,FFq) 1/seti.p*normws(seti.S(DFFqh)+seti.vdw(FFq),seti)^seti.p;

    seti.fdw1 = @(DFFqh,FFq) normws(seti.S(DFFqh)+seti.vdw(FFq),seti); % seti.p = 1
    seti.fdw2 = @(DFFqh,FFq) 1/2*normws(seti.S(DFFqh)+seti.vdw(FFq),seti)^2; % seti.p = 2
    seti.fdw3 = @(DFFqh,FFq) 1/seti.pNorm*normLp(DFFqh+FFq,seti)^seti.pNorm; % seti.qNorm ~= 2
    %fswPart = @(q,hz) seti.wW(q+hz);
    %seti.fsw = @(q,hz) normroi(seti.alpha*[real(fswPart(q,hz)); imag(fswPart(q,hz))],seti);

    % f_sw
    if isfield(seti,'wavIsom') && strcmp(seti.wavIsom,'W1p') && isfield(seti,'omegaW1p')
        seti.Ksw = @(xnRVD) seti.alpha*seti.S(seti.wWi(seti.omegaW1p.*seti.T(xnRVD)));
        seti.KswAdj = @(ys) seti.alpha*seti.S(seti.omegaW1p.*seti.wWstari(seti.T(ys)));
        %seti.KswAdj = @(ys) seti.alpha*[seti.wWstar(seti.R(ys)); seti.wWstar(seti.I(ys))];
        %seti.vsw = @(q) seti.alpha*[real(wWq(q)); imag(wWq(q))];
        seti.vsw = @(qCVU) seti.alpha*seti.S(seti.wWi(seti.omegaW1p.*seti.GD(qCVU))); % down scaled
        seti.fsw = @(qCVU,xnRVD) norminv(seti.vsw(seti.GU(seti.T(xnRVD)))+seti.vsw(qCVU),seti);
        % T(xnRVD) is upscaled and in seti.vsw down scaled...
        % using qCVD (down scaled) instead of qCVU would be useful...
        % In theory normroi(seti.S(...),...), but it does not matter
    else
        seti.Ksw = @(xnRVD) seti.alpha*seti.S(seti.wWi(seti.T(xnRVD)));
        seti.KswAdj = @(ys) seti.alpha*seti.S(seti.wWstari(seti.T(ys)));
        seti.vsw = @(qCVU) seti.alpha*seti.S(seti.wWi(seti.GD(qCVU))); % down scaled
        seti.fsw = @(qCVU,xnRVD) norminv(seti.vsw(seti.GU(seti.T(xnRVD)))+seti.vsw(qCVU),seti);
    end

    end