# autoparam3

Automatic choice of the regularization parameters alpha and beta. (This is the third version.)

Warning: This function clears variables, creates a directory and creates files.

## Contents

- Syntax
- Description
- Example
- Input "Arguments"
- In-place Input
- Output "Arguments"
- In-place Output
- More About
- References
- See Also
- Code

## Syntax

    init;
    autoparam3;

## Description

The routine autoparam3 conatins a heuristic algorithm for an automatic, data-driven choice of the regularization parameters (sparsity) and (total variation) in dependence of the relative noise level for the inverse scattering problem. Note that the ratio is fixed.

Related routines:

- For an automatic choice of `alpha` using `beta = 0`, see autoparam3alpha.html.
- For an automatic choice of `beta` using `alpha = 0`, see autoparam3beta.html.

## Example

    init; autoparam3;

    init; inseti = 'example'; autoparam3;

The parameter `inseti` is used to compute with parameters beyond the default values. See start.html in section "Quick Start" for details.

## Input "Arguments"

There are no classical input arguments because autoparam3 is not a function.

The input is analogous to varalpha.html, i.e. from default deviating parameters in struct `seti` are set as described in the section "Quick Start" of start.html. Some of these input parameters are described in example.html.

## In-place Input

Define in this file.

- alpha0 : Regularization parameter alpha, e.g. 1E5, to start the heuristic algorithm.
- N : Maximal number of steps in bisection method, e.g. 30.
- q : Quotient of alpha/beta to set beta in dependence of alpha. (e.g. 500/1E-5 using default seti.alpha = 500 and seti.beta = 1E-5).
- dummy : 0: default; 1: to test the heuristic instead of the scattering scheme the discrepancy is computed by a simple function.

## Output "Arguments"

There are no classical output arguments because autoparam3 is not a function.

Output results of the reconstructions with various regularization parameters are saved in a folder created after `start` like this

    output/20180907T171151_example

For details see start.html in the section "Output Arguments".

Furthermore, a plot containing the discrepancy of the reconstruction over the used regularization parameter alpha is saved in the folder "output", e.g.

    20180907T171151_autoparam3

The discrepancy principle border is also marked in this plot.

## In-place Output

    * disAlphaBeta       : Vector of size N x 1 storing the discrepancy with
                           the reconstruction result using |alpha(i)|,
                           |beta(i)| for each alpha-beta-combination |i|.
    * errAlphaBeta       : Vector of size N x 1 storing the error.
    * disPrincAlphaBeta  : Vector of size N x 1 storing if the discrepancy
                           principle was fulfilled or not at the i-th
                           reconstruction.
                           (1: if $|disAlphaBeta(i)| \leq \tau_\mathrm{dis}*\delta$; 0: otherwise)

## More About

As the ratio is fixed we only consider keeping in mind that changes correspondingly.

We want to choose as high as possible, such that the stabilizing penalty terms have a high effect, but as low as necessary, such that the relative discrepancy is under the discrepancy principle border . Therefore the idea of the heuristic algorithm is:

1. Find an called `a` that fulfills the discrepancy principle.
2. Find a higher called `b` that does not.
3. Halve the interval [a, b] as in bisection method.

A detailed description is given in [1].

Note that the routine `autoparam3` calls the routine `start` with various input parameters for the regularization parameters (sparsity) and (total variation) to tackle the inverse scattering problem with these parameters.

As this routine is similar to autoparam3alpha.html and autoparam3beta.html the internal scripts inAuto01init.html, inAuto02eval.html, inAuto03decAlpha.html (or inAuto03decBeta.html in the case of autoparam3beta.html) and inAuto04plot.html are used. The use of these scripts is improveable. It should be written in functions.

## References

- [1] Florian Bürgel. Effective and Efficient Reconstruction Schemes for the Inverse Medium Problem in Scattering. PhD thesis, Universität Bremen, 2019.

## See Also

- setInput.html
- varalpha.html
- varbeta.html
- varalphabeta.html
- vardelta.html
- vartol.html
- start.html
- autoparam3alpha.html
- autoparam3beta.html

## Code

    disp(' ')
    disp('----- Automatic choice of regularization parameters alpha and beta... -----')
    disp(' ')

    close all;
    if ~exist('inseti','var')
        inseti = '';
    end
    clearvars -except inseti;
    ticAutoParam = tic;

In-place settings:

    % ----------------------------------------------------------------------
    % Set initial regularization parameter alpha:

    switch inseti
        case 'example3Dalg'
            alpha0 = 1000;
        case 'fresnelPda'
            alpha0 = 1E6;
        case 'shipBorePda'
            alpha0 = 1000;
            % alpha0 = 10000; % 1000 is high enough. However, 10000 was used in the example in the phd thesis
        case 'fresnelShrink'
            alpha0 = 1E3;
        case 'fresnelSynthShrink'
            alpha0 = 0.01;
        otherwise
            alpha0 = 1E5;
    end

    % ----------------------------------------------------------------------

    % Further in-place settings (changes are not recommended):
    %
    N = 30; % maximal number of steps in bisection method
    q = 500/1E-5; % quotient of alpha/beta (default: seti.alpha = 500; seti.beta = 1E-5;)

    dummy = 0; % 0: no dummy function (use IPscatt); 1: dummy function to test this routine (result \approx 500)

Initialization

    inAuto01init;
    alpha(1) = alpha0;

    % -- diary on does not work because of evalc later on...:
    % disp(' ')
    % disp('-- start diary log ---------------------------------------')
    % disp(' ')
    % diary(sprintf('output/%s_autoparam3_diary.log',dirDatetime));
    % diary on; % later diary off

Process

    for i = 1:N
        % Evaluate the current alpha/beta

        beta = alpha./q; % compute beta (here, beta is a vector(!))
        inAuto02eval; % Reconstruction with alpha and beta; check, whether discrepancy principle border reached; output result

        inAuto03decAlpha; % decision (try alpha/beta and decide which alpha is next tried)
        if breakfor % breakfor is set in inAuto03decAlpha;
            break;
        end
    end

Output

    if startbis == 0
        disp('Attention: startbis is still 0, i.e. the choice was stopped although bisection method was not started yet.')
        disp('Suggestions:');
        disp('  - Increase the number of iterations N to search a parameter.')
        disp('  - Increase seti.tau to reach the discrepancy principle.')
        disp('  - Increase the numer of iterations of the reconstruction scheme to (maybe) get a better result.')
    end

    if startbis == 1
        disp('Choose last a as alpha, because it fulfills the discrepancy principle.')
        disp('Result of automatic parameter choice:')
        % fprintf('alpha = %g, beta = %g\n',alpha(i),beta(i));
        % fprintf('results in dis = %g\n',disAlphaBeta(i)');

        fprintf('alpha = %g, beta = %g\n',a,a./q); % is changed in case of autoparam3alpha and autoparam3beta
        fprintf('results in dis = %g\n',da);
    end
    tocAutoParam = toc(ticAutoParam);
    fprintf('Elapsed time of automatic parameter choice is %05.1f min.\n',tocAutoParam/60);

    % diary off;

    % -- Plot results --
    regplot = alpha(1:stopind); % regplot is alpha(...) or beta(...)
    inAuto04plot; % plot

    usevaralphabeta = 0; % take it out of service for further computations