%-------------------------------------------------------------------------%
%--------------------------------- hmc.m ---------------------------------%
%
%--------------- hmc settings ---------------%
numstep=1000000; % number of states generated
pvar=1; % step size in momentum space
qvar=0.1; % step size in position space
nint=10; % number of jumps in each leapfrog run
stepsize=qvar/nint; % step size in each jump

%--------------- spect decomp settings ---------------%
d=4; % dimension of Hilbert space
% nt=d*(d+1)/2-1; % number of theta parameters in HMC
% nf=d*(d-1)/2; % number of phi parameters in HMC
num=d^2-1; % total number of indpendent parameters

%-------------------------------------------------------------------------%
%--------------- POVM settings ---------------%
sigma=zeros(2, 2, 4); % Pauli matrices
sigma(:,:,1)=[1 0; 0 1]; % identity
sigma(:,:,2)=[1 0; 0 -1]; % sigma_z
sigma(:,:,3)=[0 1; 1 0]; % sigma_x
sigma(:,:,4)=[0 -1i; 1i 0]; % sigma_y

T=zeros(2,2,4); % Tetrahedron POVM
T(:,:,1)=sigma(:,:,1)+(sigma(:,:,2)-sigma(:,:,3)-sigma(:,:,4))/sqrt(3);
T(:,:,2)=sigma(:,:,1)+(-sigma(:,:,2)+sigma(:,:,3)-sigma(:,:,4))/sqrt(3);
T(:,:,3)=sigma(:,:,1)+(-sigma(:,:,2)-sigma(:,:,3)+sigma(:,:,4))/sqrt(3);
T(:,:,4)=sigma(:,:,1)+(sigma(:,:,2)+sigma(:,:,3)+sigma(:,:,4))/sqrt(3);

Q=zeros(d,d,d^2); % overall POVM: Two tetrahedrons
for j=1:d^2
    Q(:,:,j)=kron(T(:,:,ceil(j/d)),T(:,:,1+rem(j-1,d)));
end

pom=Q(:,:,1:num); % reduced POVM: independent POVM elements only,
% corresponding to the same number of independent parameters used in HMC

%-------------------------------------------------------------------------%
%--------------- hmc starts ---------------%
tic
t0=cputime;

q=zeros(num,numstep);
q(:,1)=2*pi*rand(num,1);
[A,JacDet,u]=spect(q(:,1),d,pom);

%
rho=zeros(d,d,numstep);
rho(:,:,1)=A'*A;
purity=zeros(1,numstep);
purity(1)=trace(rho(:,:,1)*rho(:,:,1)');

%
accept=0;
for i=2:numstep
    p=normrnd(0,pvar,num,1);
    
    pstar=p;
    qstar=q(:,i-1);
    
    u_ini=u;
    JacDet_ini=JacDet;
    for j=1:nint 
        phalf=pstar+(stepsize/2)*u';
        qstar=qstar+stepsize*phalf;
        [A, JacDetstar, u]=spect(qstar,d,pom);
        pstar=phalf+(stepsize/2)*u';
    end
    pstar=-pstar;
    
    aprob=exp((p'*p-pstar'*pstar)/2)*JacDetstar/JacDet;
        
    U=rand;
    if aprob>U
        q(:,i)=qstar;
        rho(:,:,i)=A'*A;
        JacDet=JacDetstar;
        accept=accept+1;
    else
        q(:,i)=q(:,i-1);
        rho(:,:,i)=rho(:,:,i-1);
        u=u_ini;
        JacDet=JacDet_ini;
    end
    
    purity(i)=trace(rho(:,:,i)*rho(:,:,i)');

end
acceptrate=accept/numstep;

%
t=cputime-t0;
fprintf('Elapsed CPU time in seconds: %f\n', t);
toc