function x=BinaryMatrixTransformation(x,action,B)
%INPUT:
% x         real symmetric matrix in IR^(1+n)x(1+n)
%
% action    a string(case sensitive) with one of the following values
%           'notranspose'           -- Evaluates CurvyT_B(x)
%           'invnotranspose'        -- Evaluates CurvyT_B^(-1)(x)
%           'transpose'             -- Evaluates CurvyT_B^(*)(x)
%           'invtranspose'          -- Evaluates CurvyT_B^(-*)(x)
% 
% B         logical column vector in {0,1}^n
%           
%
%OUTPUT:
% x         real symmetric matrix in IR^(1+n)x(1+n)
%
%
%
%Documentation:
% This function evaluates the Operator CurvyT_B^(action)(x) which is a gen-
% alized version of the operator CurvyT^(action) given in:
%
% Lieder, Felix, Fatemeh Bani Asadi Rad, and Florian Jarre. "Unifying semi-
% definite and set-copositive relaxations of binary problems and randomiza-
% tion techniques." Computational Optimization and Applications (2014):1-20
%
% which maps a matrix variable arising in SDP relaxations of {0,1} problems
% to  their counterpart arising from the SDP-Relaxations of {-1,1} problems
% The main difference here is that we work with a minors (specified by B)
% of the matrix x. Specifically for the partition 
% 
%       | y_0  y_B^T y_N^T|  
%   y = | y_B  Y_BB  Y_BN |                                      
%       | y_N  Y_NB  Y_NN |
%
% the application of CurvyT_B results in
%
%              | y_0               y_B^T-y_0 e^T              y_N^T      |   
%CurvyT_B(y) = | y_B-y_0e  4Y_BB-2x_Be^T-2ex_B+X_00ee^T    2Y_BN -ey_N^T |                                       |
%              | y_N               Y_NB-x_Ne^T                Y_NN       |
%
% or equivalently for
%
%              |  1      0          0    |   
%         Q  = | -e    2I_BB        0    |                                
%              |  0               I_NN   |
%
% we have CurvyT_B(y)= QyQ^T
%
% with the following properties for all symmetric matrices X Y
%          CurvyT_B^(-1)(CurvyT_B(X))=X= CurvyT_B(CurvyT_B^(-1)(X))
%          CurvyT_B^(-*)(CurvyT_B^(*)(X))=X= CurvyT_B^(*)(CurvyT_B^(-*)(X))
%          X \bullet Y = CurvyT_B^(-*)(X) \bullet  CurvyT_B(X)
%          X \bullet Y = CurvyT_B^(*)(X) \bullet  CurvyT_B^(-1)(X)
%
%          X is SPD if and only if CurvyT_B(X) is SPD
%
% Copyright (c) 2015 by
%   Felix Lieder
%
% Tested with MATLAB 2014a
%Disclaimer: no warranty no guarantee... If things go wrong: shit happens!

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Validating Input %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    error(['Specify action to be performed on input argument! Valid',...
        ' actions are "transpose" "notranspose" "invtranspose"',...
         '"invnotranspose"']);
else
   
    [mx,nx]=size(x);
    if mx~=nx 
    error('Expects Square Matrix')
    elseif issymmetric(x)==0
    error('Matrix is not symmetric!')
    elseif isreal(x)==0
        error('Matrix should have real entries')
    elseif nx==1
    error('Dimension of input argument to small');
    end
    
    if nargin==3
        [mB,nB]=size(B);
        if nB~=1 || mB~= mx-1 || islogical(B)==0
    error('Last Input is expected to be a logical column vector of size n')
        end
        N=[false(1);B==0];
        B=[false(1);B];
        
    else
        %If B is not specified we assume that all the variables were binary
        disp('Assuming B=true(n,1)');
        N=false(nx,1);
        B=[false(1);true(nx-1,1)];
    end
    
end
%%%%%%%%%%%%%%%%%%%%%% x is about to be transformed%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(action,'notranspose')
%% x=CurvyT_B^(x) [(0,1) is mapped to  (+-1)]
sB=sum(B);
x(B,B)=4*x(B,B)-2*x(B,1)*ones(1,sB)-2*ones(sB,1)*x(B,1).'+x(1,1)*ones(sB);
x(B,B)=0.5*(x(B,B)+x(B,B).');
x(B,1)=2*x(B,1)-ones(sB,1)*x(1,1);
x(1,B)=x(B,1).';
x(N,B)=2*x(N,B)-x(N,1)*ones(1,sB);
x(B,N)=x(N,B).';
elseif strcmp(action,'invnotranspose')
%% x=CurvyT_B^(-1)(x) [(+-1) is mapped to  (0,1)]
sB=sum(B);
x(B,B)=0.25*(x(B,B)+x(B,1)*ones(1,sB)+ones(sB,1)*x(B,1).'+x(1,1)*ones(sB));
x(B,B)=0.5*(x(B,B)+x(B,B).');
x(B,1)=0.5*(x(B,1)+ones(sB,1)*x(1,1));
x(1,B)=x(B,1).';
x(N,B)=0.5*(x(N,B)+x(N,1)*ones(1,sB));
x(B,N)=x(N,B).';
elseif strcmp(action,'transpose')
%% x=CurvyT_B^(*)(x) [Adjoint Operator of CurvyT_B]
x(1,1)=x(1,1)-2*sum(x(B,1))+sum(sum(x(B,B)),2);
x(B,1)=2*x(B,1)-2*sum(x(B,B),2);
x(1,B)=x(B,1);
x(B,B)=4*x(B,B);
x(N,B)=2*x(N,B);
x(B,N)=x(N,B).';
x(N,1)=x(N,1)-0.5*sum(x(N,B),2);
x(1,N)=x(N,1).';
elseif strcmp(action,'invtranspose')
%% x=CurvyT_B^(-*)(x) [Adjoint Operator of CurvyT_B^(-1)] 
x(1,1)=x(1,1)+sum(x(B,1))+0.25*sum(sum(x(B,B)),2);
x(B,1)=0.5*x(B,1)+0.25*sum(x(B,B),2);
x(1,B)=x(B,1);
x(B,B)=0.25*x(B,B);
x(N,B)=0.5*x(N,B);
x(B,N)=x(N,B).';
x(N,1)=x(N,1)+sum(x(N,B),2);
x(1,N)=x(N,1).';
else
error(['Error: Wrongly specified action to be performed on input'... 
    ' argument! Valid actions are "transpose"'...
    ' "notranspose" "invtranspose" "invnotranspose"']);
end