% Prop9InvMultOp: Inversion of Frame multiplier according to Proposition 9  
%for part II. Computation of M_{m,\Phi,\Psi} and M_{m,\Psi,\Phi}
% using the n0-th approximation M0inv of M_{m,\Phi,\Phi}^{-1} obtained in part I. 
function [m,TPsi,M0,M1,M2,M0inv,n0,M1inv,M2inv,n] = Prop9InvMultOp(c,r,TPhi,TG,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

% 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



% lambda - the supremum of |m_n-1|
lambda=abs(m(1)-1);
if c>1
for i=2:c
    if lambda<abs(m(i)-1)
        lambda=abs(m(i)-1);
    end
end
end



% if lambda<APhi/BPhi, one can continue with this m,
% but if  lambda>= APhi/BPhi, we determine new m
% such that its lambda=max|m_n-1| would satisfy lambda<APhi/BPhi
  k=APhi/BPhi;
if lambda>= k  
    warning('The entered m does not satisfy the assumptions of the proposition! The program adjusts m in order to satisfy the assumptions and continues the algorithm.')
    epsil=1/10; %we take a number from (0,1) to adjust m
    q=epsil*k/sqrt(2); % we use this number to adjust the elements of m
    % so that the new lambda=sup_i |m_i-1| will satisfy lambda < APhi/BPhi
    % this is needed in order to have positivity of the expression 
    % APhi-lambda*BPhi 
    mreal=1+rand(c,1)*q; %in this case |mreal_n-1|<epsil*k/sqrt(2)
    mimag=rand(c,1)*q; %in this case |mimag|<epsil*k/sqrt(2)
    %so |m_n-1|<epsil*k and hence the new lambda=max|m_n-1| <epsil*k
    m=complex(mreal,mimag);
    m=transpose(m);
    %for the max element of |m_i-1| we can do the algorithm again:
    lambda=abs(m(1)-1);
    if c>1
     for i=2:c
       if lambda<abs(m(i)-1)
        lambda=abs(m(i)-1);
       end
     end
    end
   % or for simplicity we can take lambda=epsil*k;
end



% I. Calculations for M_{m,\Phi,\Phi} and its inverse

%needed pieces of calculations, mu is BG
lambdaBPhi=lambda*BPhi;
u=lambdaBPhi/APhi;
a1=APhi-lambdaBPhi;

n0=0;
n01=log(a1*e)/(log(u));
if n01>=1 
    n0=floor(n01);
end

%SPhi - frame operator of Phi and SPhiInv - its inverse
SPhi=TPhi*(TPhi');
SPhiInv=inv(SPhi);


%Multiplier M0 is M_{m,\Phi,\Phi}, M0inv - its inverse
M0inv=SPhiInv;  % initialization (the zero iteration)

%needed here for the output check1-check4 if it happens to be the case n=0
mmatrixdiag=diag(m);
M0=TPhi*mmatrixdiag*(TPhi');

if n0>0
     %if M1 is not needed as output of the program, 
     %then the above 3 lines can be moved here 
     % as the calculations will not be needed in the case n=0
    %mmatrixdiag=diag(m);
    %M1=TPhi*mmatrixdiag*(TPsi');
    %M2=TPsi*mmatrixdiag*(TPhi');
    P0=SPhiInv*(SPhi-M0);
    Q0=SPhiInv; % initialization (the zero iteration)
for i=1:n0
    Q0=P0*Q0;% the i-th iteration for Q1
    M0inv=M0inv + Q0; % the i-th iteration for M0inv
end
end



% II. Calculations for M_{m,\Phi,\Psi} and M_{m,\Psi,\Phi}
% using the n0-th approximation M0inv of M_{m,\Phi,\Phi}^{-1} obtained above

%calculating the number which should be upper bound for mu 
%(mu - Bessel bound for G)
mubound=((a1)^2)/((lambda+1)^2*BPhi);

% creating a frame G based on which we will create FramePsi=FramePhi + FrameG
FrameG=frame('gen',TG);  % or create random frame FrameG = rand(m,n)
[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/2; % 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);
%FrameG=frame('gen',TG);
end

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


%needed pieces of calculations, mu is BG
b1=(lambda+1)*sqrt(BG*BPhi);
b1overa1=b1/a1;

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

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

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

%needed here for the output check1-check4 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, 
     %the following 3 lines can be here not needed to be
     %calculated if it is the case n=0
    %mmatrixdiag=diag(m);
    %M1=TPhi*mmatrixdiag*(TPsi');
    %M2=TPsi*mmatrixdiag*(TPhi');
    P1=M0inv*(M0-M1);
    P2=M0inv*(M0-M2);
    Q1=M0inv; % initialization (the zero iteration)
    Q2=M0inv; % initialization (the zero iteration)
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;
end
end




end