function [zmean, zstd, z, L, R] = Delta_incremental(X,C,u,random_projections,num_trials)
global B;

% ==============
% X - (n,k)
% X = zB' + e;
% 
% X*inv(S_h) = z*(inv(S_h)*B)' + e*inv(S_h);
% 
% X_new = X*inv(S_h)
% u_new = inv(S_h)*u
% C_new = inv(S_h)*C*inv(S_h)';
% 
% S   = cov(X);
% S_h = sqrtm(S);
% X = X/S_h; 
% B = S_h\B;
% e = e/S_h;
% 
% ============

% %% Original 
% 
% A  = u;
% u1 = pinv(C)*A;
% 
% B1  = u1/(u1'*A);
% 
% X2 = X*(A*B1')';
% X1 = X - X2;

%% Normalisation

S       = cov(X);
S_h     = real(sqrtm(S));
S_h_inv = pinv(S_h);
X_new   = X*S_h_inv; 
u_new   = S_h_inv*u;
C_new   = S_h_inv*C*(S_h_inv');
A       = u_new;
u1      = pinv(C_new)*A;
B1      = u1/(u1'*A);
X2      = X_new*(A*B1')';
X1      = X_new - X2;

% X_new   = X/S_h; 
% u_new   = S_h\u;
% C_new   = (S_h\C)/(S_h');
% A       = u_new;
% u1      = pinv(C_new)*A;
% B1      = u1/(u1'*A);
% X2      = X_new*(A*B1')';
% X1      = X_new - X2;

%     X1=X(:,1); X2=X(:,2);
%keyboard
z = zeros(num_trials, 1);
k = size(X1,2);
u2_random_projections = random_projections(1:num_trials,:);
u3_random_projections = random_projections(num_trials+1:2*num_trials,:);

for i=1:num_trials
    %         u2 = randn(k,1); %u2=u2/norm(u2);
    %         u3 = randn(k,1); %u3=u3/norm(u3);
    u2 = u2_random_projections(i,:)';u2=u2/norm(u2);
    u3 = u3_random_projections(i,:)';u3=u3/norm(u3);
    %         t = randn(k,1);
    %         u2 = t;
    %         u3 = t;
    %keyboard
    x  = X1*u2;
    y  = X2*u3;

    Xn(:,1) = x;
    Xn(:,2) = y;
    S = cov(Xn);
    %         S = cov(X);

    % L = log(mean(exp(x+y)))-sum(sum(S))/2;
    % R = log(mean(exp(x)))+log(mean(exp(y)))-sum(diag(S))/2;

     % L(i) = log(mean(exp(sqrt(-1)*(x+y))))+sum(sum(S))/2;
     % R(i) = log(mean(exp(sqrt(-1)*x)))+log(mean(exp(sqrt(-1)*y)))+sum(diag(S))/2;
     
     L(i) = mean(exp(sqrt(-1)*(x+y)))*exp(-sum(diag(S))/2);
     R(i) = mean(exp(sqrt(-1)*x))*mean(exp(sqrt(-1)*y))*exp(-sum(sum(S))/2);
    %L(i) = mean(exp(sqrt(-1)*(x+y)))*exp(sum(sum(S))/2);
    %R(i) = mean(exp(sqrt(-1)*x))*mean(exp(sqrt(-1)*y))*exp(sum(diag(S))/2);

    %          L = log(mean(exp(x+y)))- t'*S*t/2;
    %          R = log(mean(exp(x)))+log(mean(exp(y)))-t'*diag(diag(S))*t/2;
    %         z(i) = min(abs(L-R)/(abs(L)+abs(R)), 1000);
    %         z(i) = min(abs(L-R)/k, 1000);

    z(i) = min(abs(L(i)-R(i)), 10000);
    %z(i) = min(abs(L/R-1),10000);
%     if abs(R(i))>10
%         keyboard
%     end
    % z(i) = min(abs(L/R - 1), 10000);

    %     if(isnan(L(i)))
    %         L(i) = 10000;
    %     end
    %     if isnan(R(i))
    %         R(i) = 10000;
    %     end
    % z(i) = abs(real(L-R));
    % if(z(i) >= 1000)
    % end
end
z = abs(L-R);
 % z=abs(abs(L./R)-1);
%z = abs(L-R)./max(abs(L),abs(R));
% if((sum(z>10^4)>0)||(sum(isnan(z))>0)||(sum(isinf(z))>0))
%     keyboard;
% end
z(isnan(z)|isinf(z)|z>10^10)=10^10;

%z=abs(L-R);
%zmean=mean(z);zstd=std(z)/sqrt(num_trials/log(num_trials));
zmean = mean(z);
% if(zmean > 100000000)
%     keyboard;
% end
% zstd=2*std(z)/sqrt(num_trials);
zstd  = 0;
% t=inv(B)*u;
% if max(abs(t/norm(t)))<.8 && mean(z)<.5
%     keyboard
% end
%[zmean zstd]

% if isnan(z) || isinf(z)||z>10^10
%     z=10^10;
%     %keyboard
% end
%z = abs(L./R-1);
%keyboard
%hist(abs(L-R))
end