% Prop11InvMultOp: Inversion of Frame multiplier according to Proposition 11  
% For part II. The computations for M_{m,\Phi,\Psi} and M_{m,\Psi,\Phi}
% use the inverse M0inv of M_{m,\Phi,\Phi}^{-1} obtained iteratively in part I. 
function [m,TPsi,M1,M2,M1inv,M2inv,n] = Prop11InvMultOp(c,r,TPhi,TPsi,m,e)
%
%
% Authors:
% Diana Stoeva, dtstoeva@kfs.oeaw.at
% last modification 12.10.2017
%
% 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

%%Check whether the entered frame Psi is an approximate dual of the entered frame Phi
Smixed=TPsi*(TPhi');
epsad=norm(Smixed-eye(r));
if epsad>=1
  % in this case Psi is not an approximate dual frame of Phi 
  % and we take Psi to be the canonical dual of Phi, i.e., framedual(FramePhi)
  %TPsi - the synthesis operator of the canonical dual of FramePhi; r is the number of the coordinates in a vector
  warning('The entered frame Psi is not an approximate dual of Phi. The program replaces the entered Psi with the canonical dual of Phi and continues the algorithm.')
  TPsi=frsynmatrix(framedual(FramePhi),r); 
  BPsi = 1/APhi;
  epsad=0; %because now Psi is a dual frame of Phi
else
  FramePsi=frame('gen',TPsi); %generates the frame FramePsi
  [APsi,BPsi]=framebounds(FramePsi); % gives frame bounds  
end



% 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<(1-epsad)/sqrt(BPhi*BPsi), one can continue with this m,
% but if  lambda>= (1-epsad)/sqrt(BPhi*BPsi), we determine new m
% such that its lambda=max|m_n-1| would satisfy lambda<(1-epsad)/sqrt(BPhi*BPsi)
sqB=sqrt(BPhi*BPsi);  % used also later
k=(1-epsad)/sqB;
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=sum |m_n-1| can satisfy lambda<k
    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



% Calculations for M_{m,\Phi,\Psi} and M_{m,\Psi,\Phi}


n=0;
u=lambda*sqB+epsad;
n1=log(e*(1-u))/(log(u));
if n1>=1 
    n=floor(n1);
end

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

M1inv=eye(r);  % initialization (the zero iteration)
M2inv=M1inv;  % initialization (the zero iteration)

mmatrixdiag=diag(m);
M1=TPhi*mmatrixdiag*(TPsi');
M2=TPsi*mmatrixdiag*(TPhi');

if n>0
    P1=eye(r)-M1;
    P2=eye(r)-M2;
    Q1=eye(r); % initialization (the zero iteration)
    Q2=eye(r); % 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