% Prop8InvMultOpGaborForFigure: 
% Inversion of Frame multiplier according to Prop. 8 for Gabor frames
% this function is to be used in Prop8InvMultOpGaborPlotFigure.m
% with AIM to plot the Convergence rate (on Log-scale) of the algorithm.
function [TPhi,TPsi,M1,M2,M1inv,M2inv,n,ListErrorbounds,ListAbsErrors] = Prop8InvMultOpGaborForFigure(L,a,M,gPhi,gG,m,e)
%
%
% Author:
% Diana Stoeva, dtstoeva@kfs.oeaw.at
% last modification 11.01.2018
%
% This work is licensed under the Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
% To view a copy of this license, visit http://creativecommons.org/licenses/by-nc-sa/3.0/

%% internal options

% L - the length of the transform 
% gPhi - the window of the Gabor frame Phi
% a - the time-shift, should be divisor of L
% M - the number of channels, should be divisor of L
% gG - the window of the Gabor frame G (G is with the meaning of Psi - Phi)
% m - positive semi-normalized sequence
% e - the desired error 

%denoting N=L/a, we have MN=ML/a number of frame elements 
%(i.e. ML/a columns in the synthesis matrix)
% and L corresponds to the number of coordinates 
%(the number or rows in the synthesis matrix)
%% intialization of outside code
%ltfatstart; % intialize LTFAT


%% frame Phi 
% using the frame concept in LTFAT 
% generating a general frame via the synthesis matrix T
%(the frame vectors are the column-vectors in T)
FramePhi=frame('dgt',gPhi,a,M); %generates the frame FramePhi being a Gabor frame with window g, time-shift a and M channels
[APhi,BPhi]=framebounds(FramePhi); % gives frame bounds of FramePhi
%TPhi - the synthesis matrix of FramePhi
TPhi=frsynmatrix(FramePhi,L);

am=min(m);
bm=max(m);

%calculating the number which should be upper bound for mu
a1=am*APhi; %used later as well
mubound=((a1)^2)/(BPhi*(bm^2));


FrameG=frame('dgt',gG,a,M);  
[AG,BG]=framebounds(FrameG); % gives frame bounds
%the synthesis matrix of FrameG
TG=frsynmatrix(FrameG,L);


% if the upper frame bound of G is bigger or equal then mubound, we can
% adjust G by deleting its element to appropriate numbers to get new frame
% G with the desired limited upper frame bound
if BG>=mubound
    warning('The entered frame G does not satisfy the assumptions of the proposition! The program adjusts G by an appropriate constant in order to satisfy the assumptions and continues the algorithm.')
    p=1/10; % we use a number p from (0,1) to adjust G
    k=(sqrt(mubound/BG))*p; %the constant for adjusting g_n, we will use g_n*k
    TG=TG*k;
    BG=BG*k^2;
%[AG,BG]=framebounds(frame('gen',TG));
end



% now the frame FramePsi=FramePhi + FrameG will satisfy the assumpitons of Prop. 8
%TPsi - the synthesis matrix of FramePsi
TPsi=TPhi+TG;

%create a row with the square roots of the elements of m
sqrtm=sqrt(m);
%create a matrix with rows being equal to sqrtm 
%(this we need to get pointwise m_n phi_n below)
sqrtmmatrix=repmat(sqrtm,L,1);



% to get (sqrt(m_n) phi_n), multiply element-wise sqrtmmatrix and TPhi
TsqrtmPhi=sqrtmmatrix.*TPhi;

%SsqrtmPhi - frame operator of (sqrt(m_n) phi_n)_n; and SsqrtmPhiInv - its inverse
SsqrtmPhi=TsqrtmPhi*(TsqrtmPhi');
SsqrtmPhiInv=inv(SsqrtmPhi);

%pieces of calculations
b1=bm*sqrt(BG*BPhi);
a1minusb1=a1-b1;
b1overa1=b1/a1;

n=0;
n1=log(a1minusb1*e)/(log(b1overa1));
if n1>=1 
    n=floor(n1);
end

%Multiplier M1 is M_{m,\Phi,\Psi}, M2 is M_{m,\Psi,\Phi}

M1inv=SsqrtmPhiInv;  % initialization (the zero iteration)
M2inv=SsqrtmPhiInv;  % initialization (the zero iteration)


%if M1,M2 are in the output - necessary to define them here if it happens to be the case n=0
mmatrixdiag=diag(m);
M1=TPhi*mmatrixdiag*(TPsi');
M2=TPsi*mmatrixdiag*(TPhi');

M1invMatlab = inv(M1); % for the plot of the absolute errors ||M1inv_i-(M1)^{-1}||

if n>0
  %if M1,M2 are not needed as outputs of the program, 
  %then the def of M1,M2 can be moved here as it would not be needed in the case n=0
    %mmatrixdiag=diag(m);
    %M1=TPhi*mmatrixdiag*(TPsi');
    %M2=TPsi*mmatrixdiag*(TPhi');
    P1=SsqrtmPhiInv*(SsqrtmPhi-M1);
    P2=SsqrtmPhiInv*(SsqrtmPhi-M2);
    Q1=SsqrtmPhiInv; % initialization (the zero iteration)
    Q2=SsqrtmPhiInv; % initialization (the zero iteration)
    errorbound=b1overa1/a1minusb1;
    ListErrorbounds=[errorbound];
    ListAbsErrors=[norm(M1inv-M1invMatlab)];
for i=1:n
    Q1=P1*Q1;% the i-th iteration for Q1
    M1inv=M1inv + Q1; % the i-th iteration for M1inv
    Q2=P2*Q2;
    M2inv=M2inv + Q2;
    errorbound=errorbound*b1overa1;
    ListErrorbounds=[ListErrorbounds errorbound];
    ListAbsErrors=[ListAbsErrors norm(M1inv-M1invMatlab)];
end
end



end