function  CheckStateConstraints( CO_canonical,EV,N,M )
%CHECKSTATECONSTRAINTS Summary of this function goes here
%   Detailed explanation goes here


% CO_canonical: Constraint Operators in canonical basis 
% EV: Expectation Values associated with the linear constraints (either as
% a single number or as a range)
% N: dimension of system A
% M: dimension of system B
%
% The normalization constraint is imposed by the code
%
% This program REQUIRES the package SeDuMi, by J. Sturm, that can be downloaded
% from http://sedumi.ie.lehigh.edu/
%
% Code written F. Spedalieri.


% Determine the number of constraints


Number_of_constraints = size(CO_canonical,3);




% Check dimensions

if size(CO_canonical,1) ~= size(CO_canonical,2)
    error('Constraint operators are not square');
elseif size(CO_canonical,1) ~= N*M
    error('Constraints do not match systems dimensions');
elseif Number_of_constraints ~= size(EV,1)
    error('Expectation values do not match number of constraints');
end




% Compute the basis for the two subsystems
    
[basisA,basisB] = ComputeBasis(N,M);

alphaA=trace(basisA{1}*basisA{1});
alphaB=trace(basisB{1}*basisB{1});




% Add the normalization constraint matrix (the identity) to CO_canonical
% and recompute the number of constraints

CO_canonical(:,:,Number_of_constraints+1) = eye(size(CO_canonical,1));
EV(Number_of_constraints+1,1) = 1;

Number_of_constraints = size(CO_canonical,3);


% Compute the coefficients of the linear constraints in the tensor product
% basis

M_coeff = zeros(N^2,M^2,Number_of_constraints);

for k=1:Number_of_constraints
    for i=1:N^2
        for j=1:M^2
            M_coeff(i,j,k) = real(trace(kron(basisA{i},basisB{j})*CO_canonical(:,:,k)));
        end
    end
end


% Convert the coefficients of each constraint into the row of a matrix that
% will be used to solve the system of equations (note: the command vec(X)
% returns a column vector formed by the stacked columns of the matrix X)


LCM = zeros(Number_of_constraints,(N*M)^2);

for k=1:Number_of_constraints
    LCM(k,:) = vec(M_coeff(:,:,k)).';
end






if size(EV,2) == 1
    
    

    % Find a particular solution of the system represented by LCM and EV, and
    % determine a basis of its kernel


    sol_part = mldivide(LCM,EV); % This gives a column vector

    rho_part = zeros(N^2,M^2);



    

    for j=1:M^2
        rho_part(:,j) = (sol_part((j-1)*N^2+1:j*N^2,1));   % Transform the particular solution into a matrix of coefficients
    end                                                    % in the tensor product basis



    kernel_basis = null(LCM); % The columns are the vectorized form of the elements of the kernel basis


    dim_kernel = size(kernel_basis,2); % Dimension of LCM kernel, which is the number of extra F_i
                                       % matrices needed
    
                                       
                                       
                                       
                                       
    % If the kernel is non-empty, we do some manipulation with the kernel basis
    % vectors

    if dim_kernel > 0
        
        % Convert the kernel basis vectors into a matrix of coefficients in the
        % tensor product basis 
        % Compute the basis vectors of the kernel in the canonical basis
    
       
        
        V = zeros(N^2,M^2,dim_kernel);
        V_canonical = zeros(N*M,N*M,dim_kernel);
         
        for k=1:dim_kernel
            k
            for j=1:M^2
                V(:,j,k) = (kernel_basis((j-1)*N^2+1:j*N^2,k)); 
                for i=1:N^2    
                    V_canonical(:,:,k) = V_canonical(:,:,k) + V(i,j,k)*kron(basisA{i},basisB{j}); 
                end
                
            end
        end
    
      
    
        
        
    end



    





    % Compute the particular solution in the canonical basis

    rho_part_canonical = zeros(N*M);
    for i=1:N^2
        for j=1:M^2
            rho_part_canonical = rho_part_canonical + rho_part(i,j)*kron(basisA{i},basisB{j});
        end
    end
        

    % If the particular solution is unique we point out that fact and the
    % issues that it may bring

    if rank(LCM) == (N*M)^2
    
    
        if eig(rho_part_canonical) >= 10^(-10)
           disp('The state defined by the linear constraints seems to be unique.');
           disp('You may want to solve for it and apply PPTSE to check for entanglement');
        else
            error('The matrix defined by the linear constraints seems to be unique but not positive semidefinite. The linear constraints may be incompatible with a physical state.');
    
        end
    
        dim_kernel = 0;
    
    elseif (dim_kernel >0)
    
            % Sanity check: check whether there is a state (normalized and PSD)
            % satisfying these constraints.
    
            block_check=[2*M*N];
    
            f0 = [real(rho_part_canonical) imag(rho_part_canonical);-imag(rho_part_canonical) real(rho_part_canonical)];
            c_check = vec(f0);
    
            ct_check = 0;
    
            for k=1:dim_kernel
    
                ct_check = ct_check + 1;
    
                f_i = [real(V_canonical(:,:,k)) imag(V_canonical(:,:,k)); -imag(V_canonical(:,:,k)) real(V_canonical(:,:,k))];
    
                At_check(:,ct_check)=-vec(f_i);
            end
    
            ct_check = ct_check + 1;
    
            At_check(:,ct_check)  = -vec(eye(block_check));
    
            b_check = zeros(ct_check,1);
            b_check(ct_check) = -1;
    
            K_check.s=block_check;
            pars.fid=1;
            [x_check,y_check,info]=sedumi(At_check,b_check,c_check,K_check,pars);
    
            if (abs(y_check'*b_check - c_check'*x_check)< -c_check'*x_check)
                error('There is no state that satisfies these constraints.');
            else
                
                rho_extended = f0; % Add f0

                for i=1:ct_check-1 % Don't add the identity
                    rho_extended = rho_extended - y_check(i)*mat(At_check(1:2*(M*N)*2*(M*N),i));
                end
               
               
                rho(1:M*N,1:M*N) = rho_extended(1:M*N,1:M*N) + 1i*rho_extended(1:M*N,M*N +1:2*M*N);
            
                disp('The file RhoConstraint.mat has a state rho in the canonical basis that satisfies the constraints');
                save RhoConstraint rho CO_canonical EV
            end
    
    end



    
    
elseif size(find((EV(1:end-1,2)-EV(1:end-1,1))<0),1) ~= 0
    
    error('Check the ranges of the expectation values, some may have negative length')
    
    
    
    
    
elseif size(EV,2) == 2
    
 


    % Find particular solutions of the system represented by LCM and a homogeneous
    % term with only one nonzero component, that will be used to express any solution
    % of LCM with independent term within the ranges specified by EV.


    sol_part = zeros(N^2 * M^2,Number_of_constraints);


    for k=1:Number_of_constraints

        indep_term = zeros(Number_of_constraints,1);
        indep_term(k,1) = 1;                        % Define an independent term with only one nonzero element

        sol_part(:,k) = mldivide(LCM,indep_term); % Find a particular solution                                                         
    end




    % Compute the basis of the kernel of the linear system

    kernel_basis = null(LCM); % The columns are the vectorized form of the elements of the kernel basis


    dim_kernel = size(kernel_basis,2); % Dimension of LCM kernel, which is the number of extra F_i
                                       % matrices needed








    




    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %
    % Now we start setting up the matrices for the SDP
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    % Block structure. This is the block structure of the problem.
    % because the Re-Im trick screws up the block-diagonal structure, so we use only
    % the block structure induced by the partial transpose. We also include
    % single blocks associated with the ranges of the variables that control the expectaion
    % values

    block=[2*M*N  ones(1,2*(Number_of_constraints-1))]; 


    % Setup of the matrix F_0: this matrix is now associated with the ranges of
    % the variables that control the expectation values, and the particular
    % solution that is normalized (the last matrix in rho_part).

    G0=zeros(N*M);
    


    % The matrix G0 is constructed using the particular solution of the linear
    % constraint equations that is normalized


    for k=1:N^2
        for m=1:M^2    
            G0=G0+sol_part((m-1)*N^2+k,Number_of_constraints)*kron(basisA{k},basisB{m});
        end
    end


    F0 = zeros(2*(M*N)+2*(Number_of_constraints-1));



    F0(1:2*(M*N),1:2*(M*N))=[real(G0) imag(G0);-imag(G0) real(G0)];
    

    % Now add the part that will impose constraints on the expectation values

    for k=1:Number_of_constraints-1    
        F0(2*(M*N)+2*(k-1)+1,2*(M*N)+2*(k-1)+1) = -EV(k,1); % Constraint variable z_k to be greater than EV(k,1)
        F0(2*(M*N)+2*(k-1)+2,2*(M*N)+2*(k-1)+2) = EV(k,2);  % Constraint variable z_k to be smaller than EV(k,2)
    end




    % Vectorize F0 to pass later to Sedumi

    sr=1;
    fr=0;
    for k=1:length(block)
        fr=fr+block(k)^2;
        c(sr:fr)=vec(F0(sum(block(1:k-1))+1:sum(block(1:k)),sum(block(1:k-1))+1:sum(block(1:k))));
        sr=fr+1;
    end

    




    % Construct the F_i matrices associated with the kernel of the system and 
    % the solutions of the inhomogeneous system associated with the ranges of the 
    % expectation values.

    

    ct = 0; % Start the counter 
    
    if dim_kernel > 0


        for k=1:dim_kernel

            ct=ct+1;
            
            k

            GJ0=zeros(M*N,M*N);
            

            for i=2:N^2     
                for m=1:M^2
                    if (i==2 && m>=2)  
                        % This adds the terms corresponding to i=1, m>=2. We do
                        % this to avoid adding the term i=1, m=1 that
                        % is proportional to the identity and that will be
                        % included in At. But to reconstruct a solution of the
                        % system we need to include i=1. The (i==2) condition
                        % is to do this part only once.

                        GJ0=GJ0+kernel_basis((m-1)*N^2+1,k)*kron(basisA{1},basisB{m});
                    end

                    GJ0=GJ0+kernel_basis((m-1)*N^2+i,k)*kron(basisA{i},basisB{m});
                    
                end
            end




            F_i = zeros(2*(M*N)+2*(Number_of_constraints-1));



            F_i(1:2*(M*N),1:2*(M*N))=[real(GJ0) imag(GJ0);-imag(GJ0) real(GJ0)];
            
            % Add zeros for the new variables

            for kp=1:Number_of_constraints-1    
                F_i(2*(M*N)+2*(kp-1)+1,2*(M*N)+2*(kp-1)+1) = 0; 
                F_i(2*(M*N)+2*(kp-1)+2,2*(M*N)+2*(kp-1)+2) = 0;  
            end




            sr=1;
            fr=0;
            for i=1:length(block)
                fr=fr+block(i)^2;
                At(sr:fr,ct)=-vec(F_i(sum(block(1:i-1))+1:sum(block(1:i)),sum(block(1:i-1))+1:sum(block(1:i))));
                sr=fr+1;
            end

            

        end

    end


    % Construct the new F_i matrices associated with the varying expectation
    % values

    for k=1:Number_of_constraints-1

        ct=ct+1;




        GJ0=zeros(M*N,M*N);
        

        for i=2:N^2
            for m=1:M^2
                
                if (i==2 && m>=2)  
                    % This adds the terms corresponding to i=1, m>=2. We do
                    % this to avoid adding the term i=1, m=1 that
                    % is proportional to the identity and that will be
                    % included in At. But to reconstruct a solution of the
                    % system we need to include i=1. The (i==2) condition
                    % is to do this part only once.
                    
                    GJ0=GJ0+sol_part((m-1)*N^2+1,k)*kron(basisA{1},basisB{m});
                end


                GJ0=GJ0+sol_part((m-1)*N^2+i,k)*kron(basisA{i},basisB{m});
                
            end
        end

        F_i = zeros(2*(M*N)+2*(Number_of_constraints-1));


        F_i(1:2*(M*N),1:2*(M*N))=[real(GJ0) imag(GJ0);-imag(GJ0) real(GJ0)];
        
        % Add zeros for the new variables

        for kp=1:Number_of_constraints-1    
            F_i(2*(M*N)+2*(kp-1)+1,2*(M*N)+2*(kp-1)+1) = 0;  % Enlarge the F_i matrices
            F_i(2*(M*N)+2*(kp-1)+2,2*(M*N)+2*(kp-1)+2) = 0;  
        end

        F_i(2*(M*N)+2*(k-1)+1,2*(M*N)+2*(k-1)+1) = 1;  % To keep variable z_k in the required range
        F_i(2*(M*N)+2*(k-1)+2,2*(M*N)+2*(k-1)+2) = -1;


        sr=1;
        fr=0;
        for i=1:length(block)
            fr=fr+block(i)^2;
            At(sr:fr,ct)=-vec(F_i(sum(block(1:i-1))+1:sum(block(1:i)),sum(block(1:i-1))+1:sum(block(1:i))));
            sr=fr+1;
        end

    end


    % Add the identity to the block of rho
    
    F_i = zeros(2*M*N+2*(Number_of_constraints-1)); 
    
    F_i(1:2*M*N,1:2*M*N) = eye(2*M*N);
    
    
    ct = ct+1;
    sr=1;
    fr=0;
    for i=1:length(block)
        fr=fr+block(i)^2;
        At(sr:fr,ct)=-vec(F_i(sum(block(1:i-1))+1:sum(block(1:i)),sum(block(1:i-1))+1:sum(block(1:i))));
        sr=fr+1;
    end






    % Construct the vector b

    b = sparse(ct,1);   % this a larger vector to match the larger At
    b(ct) = -1;  % This the element corresponding to the identity matrix in the F_i
                            % Since we are using the Re-Im trick, the trace of
                            % mat(x)
                            % must be set to 2*dA if we want trace(Z)=1








    % Set parameters and call SeDuMi

   

    K.s=block;





    pars.fid=0;
    [x,y,info]=sedumi(At,b,c,K,pars);
    
    
    if (abs(y'*b - c*x)< -c*x)
                error('There is no state that satisfies these constraints.');
    else

        rho_extended = F0(1:2*N*M,1:2*N*M); % Add F0
        

        for i=1:ct-1 % Don't add the identity
            rho_extended = rho_extended - y(i)*mat(At(1:2*(M*N)*2*(M*N),i));
        end


        rho(1:M*N,1:M*N) = rho_extended(1:M*N,1:M*N) + 1i*rho_extended(1:M*N,M*N +1:2*M*N);

        disp('The file RhoConstraint.mat has a state rho in the canonical basis that satisfies the constraints');
        save RhoConstraint rho CO_canonical EV
    end
    

end
    
    
    
end


function [basisA,basisB] = ComputeBasis(N,M)

% Computes the basis of the two subsystems

d=N;
ct=0;

basisA=cell(d*(d+1)/2);

for m=1:d
    for n=1:m
        ct=ct+1;
        if (m==n) && (m==1)
            basisA{ct}=speye(d)/d;
        elseif m==n
            basisA{ct}=sparse(1:(m-1),1:(m-1),-ones(1,m-1)/sqrt(m*(m-1)),d,d);
            basisA{ct}(m,m)=(m-1)/sqrt(m*(m-1));
            norm=trace(basisA{ct}*basisA{ct});
            basisA{ct}=basisA{ct}/sqrt(d*norm);
        else
            basisA{ct}=sparse([m n],[n m],[1/sqrt(2) 1/sqrt(2)],d,d);
            norm=trace(basisA{ct}*basisA{ct});
            basisA{ct}=basisA{ct}/sqrt(d*norm);
            ct=ct+1;
            basisA{ct}=sparse([m n],[n m],[1i/sqrt(2) -1i/sqrt(2)],d,d);
            norm=trace(basisA{ct}*basisA{ct});
            basisA{ct}=basisA{ct}/sqrt(d*norm);
        end
    end
end





% Now compute a basis for party B, but only if M~=N.


if M==N
   basisB=basisA;
else
  d=M;
  ct=0;
  basisB=cell(d*(d+1)/2);
  
  for m=1:d
      for n=1:m
          ct=ct+1;
          if (m==n) && (m==1)
              basisB{ct}=speye(d)/d;
          elseif m==n
              basisB{ct}=sparse(1:(m-1),1:(m-1),-ones(1,m-1)/sqrt(m*(m-1)),d,d);
              basisB{ct}(m,m)=(m-1)/sqrt(m*(m-1));
              norm=trace(basisB{ct}*basisB{ct});
              basisB{ct}=basisB{ct}/sqrt(d*norm);
          else
              basisB{ct}=sparse([m n],[n m],[1/sqrt(2) 1/sqrt(2)],d,d);
              norm=trace(basisB{ct}*basisB{ct});
              basisB{ct}=basisB{ct}/sqrt(d*norm);
              ct=ct+1;
              basisB{ct}=sparse([m n],[n m],[1i/sqrt(2) -1i/sqrt(2)],d,d);
              norm=trace(basisB{ct}*basisB{ct});
              basisB{ct}=basisB{ct}/sqrt(d*norm);
          end
      end
  end
end




end