% fisher   Quantum Fisher information. Usage: fisher(rho,A), where
%          rho is a density matrix and A is an operator.
%          Alternatively, it can also be used as fisher(rho,A,B).
%          Note that fisher(rho,A,A)=fisher(rho,A).
%          The form fisher(rho,A,B,threshold) makes it possible
%          to define the threshold below which an eigenvalue is
%          consdiered zero. The defaults value is 1e-20.

function f=fisher(rho,J,varargin)

% Make rho surely Hermitian
% A little bit of non-Hermitianicity can lead to
% non-orthogonal eigenvalues for eig
rho=(rho+rho')/2;

%Threshold=1e-10;
Threshold=1e-20;

if length(varargin)==0

   [v,d]=eig(rho);
   [sy,sx]=size(rho);
   f=0;
   for n=1:sx
       for m=1:n-1
           lambdan=d(n,n);
           lambdam=d(m,m);
           if abs(lambdam+lambdan)>Threshold;
               f=f+(lambdam-lambdan)^2/(lambdam+lambdan)*abs(braket(v(:,n),J,v(:,m)))^2;
           end %if      
       end %for
   end %for

   f=f*4; % Because of missing factor of 2, and because we count only half of the terms  
   
else
    
   if length(varargin)==2
      Threshold=varargin{2};
   end %if
      
   J2=varargin{1};
    
   [v,d]=eig(rho);
   [sy,sx]=size(rho);
   f=0;
   for n=1:sx
       for m=1:n-1
           lambdan=d(n,n);
           lambdam=d(m,m);
           if abs(lambdam+lambdan)>Threshold;
               f=f+(lambdam-lambdan)^2/(lambdam+lambdan)*...
               (braket(v(:,n),J,v(:,m)))*(braket(v(:,m),J2,v(:,n)));
           end %if      
       end %for
   end %for

   f=f*4; % Because of missing factor of 2, and because we count only half of the terms
   
end %if