%% Fragment 1 %%

p = zeros(1,nMax);
for n = 1:nMax
    p(n)=exp(-mu*T+n*log(mu*T)-sum(log(1:n)))*xf(n,b)/P_ab;
    T = T + p(n) * xb(:,n-1) * R;
    for k in 2:n
        T =  T + p(n) * ( xb(:,n-k) * xf(k-1,:) ) .* R ;
    end
end

% However, the number of stepping stones is too large to do the above, so
% we focus on the analysis of fold transitions

%% Fragment 2 %%

Qfolds = matrix of rates summed over all stepping stones in each fold

% e.g. if we have 9 stepping stones of three lengths
% divided into three folds in the form

   fold = [  1,2,3,   1,2,3,   1,2,3   ];
 %          length 1 length 2 length 3

folds = cell(nFolds,1);
for n=1:nFolds
    folds(n) = fold==n;
end

% Then Qfolds would be

  % Qfolds = [ sum(Q([1,3,7],[1,3,7])) sum(Q([1,3,7],[2,4,8])) sum(Q([1,3,7],[3,6,9])) ;
  %	       sum(Q([2,4,8],[1,3,7])) sum(Q([2,4,8],[2,4,8])) sum(Q([2,4,8],[3,6,9])) ;
  %            sum(Q([3,6,9],[1,3,7])) sum(Q([3,6,9],[2,4,8])) sum(Q([3,6,9],[3,6,9])) ]

% But that won't give a proper transition matrix -- need to do the marginalisation
% on Q I think, and then form Rfolds from this. Does this give the right probs though?'

% Where Rfolds(i,j) indicates the rate for transitioning between fold i and fold j
% The above would be slow in Matlab, but could be computed at the same time as the
% original R matrix.

% Similarly we need a fold level version of xf and xb

for k=1:nFolds
	xffolds(1,k) = sum(xf(1,folds(k)));
end
for k=1:nFolds
	xbfolds(1,k) = sum(xb(1,folds(k)));
end

for i = 2:nMax
xffolds(i,:) = xffolds(i-1,:) * Rfolds;
xbfolds(:,i) = Rfolds * xbfolds(:,i-1);
% Then to compute the expected number of such transitions, we use Corollary 3 of H&J

p = zeros(1,nMax);
Tfolds = zeros(nFolds);
for n = 1:nMax
    p(n)=exp(-mu*T+n*log(mu*T)-sum(log(1:n)))*xffolds(n,b)/P_ab;
    Tfolds = Tfolds + p(n) * xbfolds(:,n-1) * R;
    for k in 2:n
        Tfolds =  Tfolds + p(n) * ( xbfolds(:,n-k) * xffolds(k-1,:) ) .* Rfolds ;
    end
end