function w = whiteNoiseCalc(seti,W,G,iG,nCD,nROI)
% white noise calculation
% two cases: temp. definitions in this file (via calling)
% 1) nCD = seti.nCD, nROI = seti.nROI
% 2) nCD = seti.nCDinv, nROI = seti.nInv

% figure 4: weights from white noise
% figure 5: weights from white noise manually adapted
% figure 6: convergence of white noise

plotWhiteNoise = 0; % plot figures 2 and 3
plotConv = 0; % to plot figure 4

samples = seti.samples;

%numROI = seti.nROI.^seti.dim; % number of entries in ROI
%numInv = seti.nInv.^seti.dim; % number of entries in nInv
%n = zeros(numROI,1);
%n = zeros(numInv,1);
%numROI = nCD.^seti.dim;
numROI = nROI.^seti.dim;
n = zeros(numROI,1);
%wn = n; % weight from white noise
%w = n; % weight from white noise maybe manually adapted

state = rng();
rng(0)

%fileWhiteNoise = sprintf('mat/whiteNoise/wav_%s_ext_%s_smin_%d_nCD_%d_nROI_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,seti.nCD,seti.nROI,seti.dim,seti.samples);
%fileWhiteNoise = sprintf('mat/whiteNoise/wav_%s_ext_%s_smin_%d_nCDinv_%d_nInv_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,seti.nCDinv,seti.nInv,seti.dim,seti.samples);
fileWhiteNoise = sprintf('inmat/whiteNoise/wav_%s_ext_%s_smin_%d_nCD_%d_nROI_%d_dim_%d_samples_%d.mat',seti.wavelet,seti.extension,seti.smin,nCD,nROI,seti.dim,seti.samples);

if (exist(fileWhiteNoise, 'file') == 2 && seti.genWhiteNoiseNew ~= 1)
    %file with name exists
    disp('Loading white noise matrix.');
    load(fileWhiteNoise); % contains white noise matrix n
else
    disp('Generating white noise...');

    conv1vec = zeros(samples,1);
    conv2vec = zeros(samples,1);
    for i = 1:samples
        %noise = randn(size(n))+1i*randn(size(n));
        noise = randn(size(n));

        nBefore = n;
        %noise  = noise/normws(noise,g)*normws(n,g)*delta; % optional...
        %n = n + abs(ldwt(noise,'forward',wavelet,extension,smin)).^2;
        %W     = @(x) iG(ldwt(G(x),'forward',wavelet,extension,smin));
        n = n + abs(W(noise)).^2;

        % W     = @(x) iG(ldwt(G(x),'forward',wavelet,extension,smin));
        %n = n + abs(W(n)).^2;

        %n = n + abs(ldwt(randn(sizeq))).^2; %original

        % -- convergence... experimental
        nDiff = n-nBefore;
        conv1 = norm(nDiff)/norm(n);
        conv2 = max(nDiff(:))/norm(n); % maximal entry in matrix / norm
        conv1vec(i) = conv1;
        conv2vec(i) = conv2;
        breakWithConv = 0; % experimental
        if breakWithConv
            fprintf('%d | %g | %g\n',i,conv1,conv2);
            if (conv1 < 1E-4 && conv2 < 1E-5)
                disp('convergence criterium is true. break.')
                break;
            end
        end
        % --

    end

    if plotConv
        figure(6);
        xi = 100:samples;
        plot(xi,conv1vec(xi),xi,conv2vec(xi));
    end

    n = 1/samples.*n; %white noise matrix % now: vector!

    fprintf('finished.\n');
    rng(state);

    disp('Save white noise matrix.'); %now: vector!
    try
        mkdir inmat/whiteNoise/
    catch
    end
    save(fileWhiteNoise,'n','-mat')
end

disp('calculating weight matrix...') % now: vector!
wn = 1./n; %the resulting weights are
% wn: weights from noise

% now: use matrix...
wn = G(wn); % matrix!
w = zeros(size(wn));
%previously: nA = nCD... but here restricted to ROI

% -- adapt weights manually (only in 2D)

%nA = seti.nROI;
%nA = seti.nInv;
nA = nROI;
%nA: n of Area in 1 dimension (previously: in first step nA = nCD; now: nA = nROI)
if (nA < 2*seti.smin && seti.adaptWeightsMan)
    disp('ROI is too small to use manually adapted weights. So they are not used.')
end
if (nA >= 2*seti.smin && seti.adaptWeightsMan && ((strcmp(seti.wavelet,'legall53')) || (strcmp(seti.wavelet,'cdf97'))))
    % adapt weights manually AND use (wavelet legall53 OR wavelet cdf97)

    % ----- adapt weights manually

    mm = mean(mean(wn));
    iwdepth = 0;

    %while nA >= smin
    while nA >= 2*seti.smin % to avoid problems with boundary
    %for iwdepth = 1:wdepth
        iwdepth = iwdepth+1;
        %i1 = 1:floor(nA/2);
        %i2 = (floor(nA/2)+1:nA);
        i1 = 1:floor(nA/2);
        i2 = (floor(nA/2)+1:nA);

        wn1 = wn(i1,i1); %top left
        wn2 = wn(i1,i2); %top right
        wn3 = wn(i2,i1); %bottom left
        wn4 = wn(i2,i2); %bottom right

        %automatic weight
        m = zeros(4,1);
        m(1) = mean(mean(wn1));
        m(2) = mean(mean(wn2));
        m(3) = mean(mean(wn3)); % m(3) = m(2)
        m(4) = mean(mean(wn4));

        %wm1 = 2^iwdepth*m1/m2
        %wm2 = 2^iwdepth*m2/m2
        %wm3 = 2^iwdepth*m3/m2

        wm = m/mm; % compared with complete w-Matrix (mm = mean(mean(w)))
        %wm = wm.^20; % adapt weight manually
        % wm = wm;

        %wm: weight manual
        %wm1 = 2^1;
        %wm2 = 2^2;
        %wm3 = 2^3;

        %multiplication: manual weight and noise weight
        w1 = wm(1).*wn1;
        w2 = wm(2).*wn2;
        w3 = wm(3).*wn3;
        w4 = wm(4).*wn4;

        w(1:nA,1:nA) = [w1, w2; w3, w4]; % resulting weight matrix
        % next iteration will overwrite w1 (this is intended)
        % (except there is no next iteration step)

        %plot in every iteration (progress can be seen)
        %figure(2+iwdepth);
        %imagesc(w); colorbar;
        %axis xy;

        nA = floor(nA/2); % new nA for next iteration
    end

    % -----

else %case: wavelet = 'haar'; OR no manually adapted weights
    %do not use manually adapted weights (use just the weights from whiteNoise)
    disp('else...');
    w = wn;

end

% -----

if plotWhiteNoise
    figure(4);
    imagesc(wn); colorbar;
    axis xy;
    surf(wn);colorbar;shading interp;
    title('weights from white noise')

    figure(5);
    imagesc(w); colorbar;
    axis xy;
    surf(w);colorbar;shading interp;
    title('weights: manually adapted');

    error('stop')
end

% make a vector of the matrix
w = iG(w);
%size(w)

end