
% Include the results from p4_6:
p4_6;

% For the initial probability vector we choose a uniform distribution:
n=0;
nu=ones(1,6)/6;
epsilon=0.00001;
c=1/2;
% As a norm we choose the sum of the absolute values of the probability vector,
% as well as for the residual error 


last_residualerror = 0
residualerror = norm( nu(end,:) - nu(end,:) * P, 1 );
while residualerror >= epsilon

  % Geometric convergence
  if n >= 1

    if residualerror > c * last_residualerror
      nu(end,:) = nu(end,:) * P;
    end
     
  end

  n = n + 1;

  % Now calculate the macro-state transition probability matrix according
  % to Eq. (4.21):
  Gamma = zeros(3,3);
  Gamma_Ij = zeros( 3, 6 );

  for I = (0:2)
    for J = (0:2)

      % Determine the partitions
      switch I
        case 0
          SI=S0;
        case 1
          SI=S1;
        case 2
          SI=S2;
      end
      switch J
        case 0
          SJ=S0;
        case 1
          SJ=S1;
        case 2
          SJ=S2;
      end

      PIJ = P(SI,SJ);
      [ k, l ] = size( PIJ );
      row_sums = PIJ * ones( l, 1 );
      nu_macro = nu(end,SI)/sum(nu(end,SI));
      Gamma(I+1,J+1) = nu_macro * row_sums;

       % Now the transition probabilities from macro-state I to micro-states j in S can easily
       % be derived in this iteration (Eq. (4.22))
       Gamma_Ij(I+1,:) = nu_macro * P(SI,:);
       Gamma_Ij(I+1,SI) = 0;

    end
  % Now the transition probabilities from macro-state I 
  end

  Gamma
  Gamma_Ij

  % Now calculate the equilibrium of Gamma by its left-sided eigenvector with the largest
  % eigenvalue (=1)
  [ evec, eval ] = eig(Gamma');
  gamma = evec(:, find( max(eval') == max(max(eval)) ) );
  gamma = gamma' / sum(gamma)

  % Add a new approximation vector
  nu = [ nu; zeros(1,6) ];

  % Now for every subset S_I (0 <= I < 3) the linear equations from Eq.~(4.23) need to be solved:
  for I = (0:2)

    % Determine the partition
    switch I
      case 0
        SI=S0;
      case 1
        SI=S1;
      case 2
        SI=S2;
    end

    % We have to solve the linear equations system x*A = b <=> x = b * inv(A), where
    A = P(SI,SI)
    A = A - eye(size(P(SI,SI)));
    K = find( (0:2) != I);
    b = -gamma(K) * Gamma_Ij( K, SI )

    nu(end,SI) = b * inv(A)

  end

  % Normalize nu
  nu(end,:) = nu(end,:) / sum(nu(end,:))


  % Calculate the residual error
  last_residualerror = residualerror;
  residualerror = norm( nu(end,:) - nu(end,:) * P, 1 );

end

nu
