% Prop8InvMultf: Computation of M^{-1}f according to Proposition 8 
function [TPsi,M1,M2,M1invf,M2invf,n] = Prop8InvMultf(c,r,TPhi,TG,m,f,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

% c frame vectors and each one - with r coordinates; r rows, c columns;
% rxc - the size of the synthesis matrices of the given frames
% m - positive semi-normalized sequence

%% 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('gen',TPhi); %generates the frame FramePhi
[APhi,BPhi]=framebounds(FramePhi); % gives frame bounds

a=min(m);
b=max(m);

%calculating mubound - the number which should be upper bound for mu
a1=a*APhi; %needed below as well
mubound=((a1)^2)/(BPhi*(b^2));

% creating a frame G based on which we will create FramePsi=FramePhi + FrameG
FrameG=frame('gen',TG);  
[AG,BG]=framebounds(FrameG); % gives frame bounds

% 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/4; % we use a number p from the interval (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;
%FrameG=frame('gen',TG);
end

% taking FramePsi=FramePhi + FrameG - it now satisfies the assumpitons of Prop. 8
TPsi=TPhi+TG;
%FramePsi=frame('gen',TPsi);

%create a row with the square roots of the elements of m
sqrtm=sqrt(m);
%create a matrix with r rows, each one being equal to sqrtm, 
%(this we need to get pointwise m_n phi_n below)
sqrtmmatrix=repmat(sqrtm,r,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=b*sqrt(BG*BPhi);
a1minusb1=a1-b1;
b1overa1=b1/a1;

% to determine n - the needed number of iterations
n=0;
normf=norm(f);
n1=log(a1minusb1*e/(normf))/(log(b1overa1));
if n1>=1 
    n=floor(n1);
end


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

M1invf=SsqrtmPhiInv*f;  % initialization (the zero iteration)
M2invf=M1invf;  % 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');

if n>0
 %if M1 is not needed as output 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=M1invf; % initialization (the zero iteration)
    q2=M1invf; % initialization (the zero iteration)
for i=1:n
    q1=P1*q1;% the i-th iteration for Q1
    M1invf=M1invf + q1; % the i-th iteration for M1inv
    q2=P2*q2;
    M2invf=M2invf + q2;
end
end



end