ctmc_takahashi.m

function [p,p_1,pcourt,Qperm,eps,epsMAX]=ctmc_takahashi(Q,MS,numSteps)
% CTMC_TAKAHASHI - Takahashi's aggregation-disaggregation method
% [p,p_1,pcourt,Qperm,eps,epsMAX] = CTMC_TAKAHASHI(Q,MS,numSteps)
% -- Input
% Q       : infinitesimal generator matrix
% MS      : cell array where MS{i} is the set of rows of Q in macrostate i
% numSteps: number of iterative steps
% -- Output
% p       : estimated steady-state probability vector
% p_1     :
% pcourt  : steady-state probability vector estimated by ctmc_courtois
% Qperm   : permuted Q matrix w.r.t MS as returned by ctmc_courtois
% eps     : NCD index as returned by ctmc_courtois
% epsMAX  : max acceptable value for eps (otherwise Q is not NCD)
% -- Remarks
% * The initial approximate solutions is obtained by calling CTMC_COURTOIS(Q,MS)
% * No convergence stop criterion is currently implented
 
%% INIT
nMacroStates = size(MS,1); % Number of macro-states
nStates=size(Q,1);
%% START FROM COURTOIS DECOMPOSITION SOLUTION
[pcourt,Qperm,Qdec,eps,epsMAX,P,B,C,q]=ctmc_courtois(Q,MS);
% Randomize the UNPERMUTED Q. The rate is passed explicitly rather than left to
% the ctmc_randomization default, which is max|Q|+rand and therefore unseeded:
% the iteration below is invariant to q (the aggregation and disaggregation
% equations are homogeneous in P-I = Q/q), so the random rate only injected
% irreproducibility. 1.05*max|Q| is the same rate ctmc_courtois derives.
P=ctmc_randomization(Q,1.05*max(max(abs(Q))));
%% STEP 0
pn=pcourt;
% P and pn are both in the ORIGINAL state ordering, so each macro-state is
% addressed by the state indices MS{I} carries. The former contiguous offsets
% (procRows+i) only coincide with MS{I} when every macro-state happens to be a
% contiguous range of Q in original order: otherwise they read the CONTENTS of
% MS{I} nowhere and silently solve for a different partition of the same block
% sizes.
%% MAIN LOOP
for n=1:numSteps
pn_1=pn;
p_1=pn_1;
%% AGGREGATION STEP
G=zeros(nMacroStates,nMacroStates);
for I = 1:nMacroStates  % for each source macro-state
    idxI=MS{I}(:)';
    S=sum(pn_1(idxI));
    for J = 1:nMacroStates  % for dest macro-state
        if I~=J
            idxJ=MS{J}(:)';
            for i=1:length(idxI)
                for j=1:length(idxJ)
                    if S>1e-14
                    G(I,J)=G(I,J)+P(idxI(i),idxJ(j))*pn_1(idxI(i))/S;
                    end
                end
            end
        end
    end
end
for i = 1:nMacroStates  % for each source macro-state
    G(i,i)=1-sum(G(i,:));
end
gamma=dtmc_solve(G); %compute macroprobabilities
 
%% DISAGGREGATION STEP
GI=zeros(nMacroStates,nStates);
for I = 1:nMacroStates  % for each source macro-state
    idxI=MS{I}(:)';
    S=sum(pn_1(idxI));
    for j=1:nStates
        GI(I,j)=GI(I,j)+sum(P(idxI,j).*pn_1(idxI)')/S;
    end
end
for I = 1:nMacroStates  % for each source macro-state
    idxI=MS{I}(:)';
    A=eye(length(idxI),length(idxI));
    b=zeros(length(idxI),1);
    for i=1:length(idxI)
        for j=1:length(idxI)
        A(i,j)=A(i,j)-P(idxI(j),idxI(i));
        end
        for K=1:nMacroStates
            if K~=I
                b(i)=b(i)+gamma(K)*GI(K,idxI(i));
            end
        end
    end
    xI=[];
    if size(A,1) > 6000
        [xg,gflag]=ctmc_gmres(sparse(A),b);
        if gflag==0
            xI=xg;
        end
    end
    if isempty(xI)
        xI=A\b;
    end
    pn(idxI)=xI;
end
%% END LOOP
pn=pn/sum(pn);
end
%% OUTPUT
p=pn(:)';
end