% Resampled Monte-Carlo to compute the generalized Lyapunov exponents l(q)
% See J. Vanneste, Phys. Rev. E 81, 036701 (2010).
K = 100; N = 100; gamma = 1; omega = 0; q = -20:20;
zeta = sqrt(gamma^2-omega^2);
Gamma = gamma/zeta; Omega = omega/zeta;
cz = cosh(zeta); sz = sinh(zeta);
rng('default');
ev = zeros(K,2); ev(:,1) = 1;  % initial vectors [1 0]
evh = zeros(size(ev)); normevh = zeros(1,K); beta = zeros(1,N);
ell = []; ell0 = []; ellth = []; A = zeros(K,2,2);
for iq = 1:length(q)
  for n = 1:N
    % random matrix
    th2 = 2*pi*rand(1,K); ct = cos(th2); st = sin(th2);
    A(:,1,1) = cz + Gamma*ct*sz;      A(:,1,2) = (Gamma*st - Omega)*sz;
    A(:,2,1) = (Gamma*st + Omega)*sz; A(:,2,2) = cz - Gamma*ct*sz;
    % multiply matrix ev by matrix A (vectorized over realisations K)
    evh(:,1) = A(:,1,1).*ev(:,1) + A(:,1,2).*ev(:,2);
    evh(:,2) = A(:,2,1).*ev(:,1) + A(:,2,2).*ev(:,2);
    % make unit vector evh
    normevh = sqrt(evh(:,1).^2 + evh(:,2).^2);
    evh(:,1) = evh(:,1)./normevh; evh(:,2) = evh(:,2)./normevh;
    % save q^th power of the norm
    alpha = normevh.^q(iq);
    % resampling
    gamma = cumsum(alpha); beta(n) = gamma(K);
    eps = beta(n)*rand(1,K);
    for k = 1:K, ev(k,:) = evh(find(gamma-eps(k) >= 0,1),:); end
  end
  % l(q) without resampling
  ell0 = [ell0 log(mean(alpha))];
  % l(q) with resampling
  ell = [ell mean(log(beta))-log(K)];
  % the analytic expression for l(q)
  aa = sqrt(1 - (Gamma^2*cosh(2*zeta) - Omega^2)^-2);
  ellth = [ellth log(((1-aa)/(1+aa))^(q(iq)/4)*...
		     hypergeom([.5 -q(iq)/2],1,-2*aa/(1-aa)))];
end

plot(q,ell0,'go','LineWidth',2), hold on
plot(q,ell,'-','LineWidth',2)
plot(q,ellth,'r--','LineWidth',2)
legend('MC','resampled MC','exact','Location','North')
xlabel('$q$','Interpreter','LaTeX','FontSize',22)
ylabel('$\ell(q)$','Interpreter','LaTeX','FontSize',22)
set(gca,'FontSize',18,'FontName','Times')
hold off
print -dpdf resampled_mc.pdf
