spm_A_reduce.m 1.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960
  1. function [J,z,v,s] = spm_A_reduce(J,x,T,N)
  2. % FORMAT [J,z,v,s] = spm_A_reduce(J,x,T,N)
  3. % reduction of Markovian partition
  4. % J - Jacobian (x)
  5. % x - {3 x n} particular partition of states
  6. % T - eigenvalue threshold to retain eigenvectors [default: 8]
  7. % N - maxiumum number to retain [default: 8]
  8. %
  9. % J - Jacobian (z)
  10. % z - {1 x n} partition of states at the next level
  11. % v - {1 x n} eigenvector (adiabatic) operator
  12. % s - {1 x n} eigenvalues
  13. %
  14. % Adiabatic reduction operator (R)
  15. %__________________________________________________________________________
  16. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  17. % Karl Friston
  18. % $Id: spm_A_reduce.m 7655 2019-08-25 20:10:20Z karl $
  19. % preliminaries
  20. %--------------------------------------------------------------------------
  21. nx = size(x,2); % number of partitions
  22. if nargin < 3
  23. T = 8; % adiabatic threshold
  24. end
  25. if nargin < 4
  26. N = 8; % maxiumum number
  27. end
  28. % reduction
  29. %--------------------------------------------------------------------------
  30. for i = 1:nx
  31. % Lyapunov exponents (eigensolution) for this partition
  32. %----------------------------------------------------------------------
  33. y{i} = spm_vec(x(1:2,i));
  34. Jii = full(J(y{i},y{i}));
  35. [e,r] = eig(Jii);
  36. r = diag(r);
  37. [d,j] = sort(real(r),'descend');
  38. % Adiabatic threshold
  39. %----------------------------------------------------------------------
  40. n(i) = sum(d > -T);
  41. n(i) = min(n(i),N);
  42. s{i} = r( j(1:n(i)));
  43. v{i} = e(:,j(1:n(i)));
  44. u{i} = pinv(v{i});
  45. end
  46. for i = 1:nx
  47. for j = 1:nx
  48. Jij = full(J(spm_vec(x(1:2,i)),spm_vec(x(1:2,j))));
  49. A{i,j} = u{i}*Jij*v{j};
  50. end
  51. z{i} = sum(n(1:(i - 1))) + (1:n(i));
  52. end
  53. J = spm_cat(A);