function [A lambda] = spec ( M )
% SPEC generates the spectral expansion matrices of the supplied matrix
% [A lambda] = SPEC(MATRIX)
% A is a cell array of the spectral expansion matrices
% lambda is a vector of the corresponding eigenvalues
[h v] = size(M);
if ( h ~= v ), error 'supplied matrix must be square!'; end

A = {};
[V, D] = eig(M);
rV = inv(V);
lambda = diag(D);

for k = 1:h
    cE = V(:, k);
    rE = rV(k, :);
    A{k} = cE * rE;
end;


function B = sylv ( M, f, args )
% SYLV calculates a function of a matrix using Sylvester's Theorem
% B = SYLV(MATRIX, FUNC, ARGS)
% FUNC must be of the form RESULT = F(MATRIX, ARG)
% B is a cell array containing the results for all ARG in ARGS
[h v] = size(M);
if ( h ~= v ), error 'supplied matrix must be square!'; end

[A lam] = spec(M);
B = {};

for r = 1:length(args)
    Br = f(lam(1), args(r)) * A{1};
    for k = 2:h
        Br = Br + f(lam(k), args(r)) * A{k};
    end
    B{r} = Br;
end


function [D_o D_s D_b] = qideal ( Q, kA, kB, kC, t )
% QIDEAL calculates the ideal open, shut and burst time distributions for a Q matrix
% [f_open f_shut f_burst] = QIDEAL( Q, kA, kB, kC, t )
%
% Q must be in canonical form, with A (open) states before B (fast shut) before C (slow shut)
% kA, kB and kC are the number of each kind of state
% t is a vector of time points for which to calculate the f values
%
% no account is taken of the distortion caused by missed events, and error checking is minimal
[h v] = size(Q);
if ( h ~= v ), error 'Q matrix must be square'; end

kF = kB + kC;
kE = kA + kB;
k = kA + kF;

if ( k ~= h ), error 'kA + kB + kC is not the same as k'; end

% get the partition matrices
Q_AA = Q(1:kA, 1:kA);
Q_AF = Q(1:kA, (kA+1):k);
Q_FA = Q((kA+1):k, 1:kA);
Q_FF = Q((kA+1):k, (kA+1):k);

Q_BB = Q((kA+1):(kA+kB), (kA+1):(kA+kB));
Q_CB = Q((kE+1):k, (kA+1):kE);
Q_CA = Q((kE+1):k, 1:kA);
Q_BA = Q((kA+1):kE, 1:kA);
Q_AB = Q(1:kA, (kA+1):kE);
Q_EE = Q(1:kE, 1:kE);

% define the (time independent) G matrices
G_AF = inv(-Q_AA) * Q_AF;
G_AB = inv(-Q_AA) * Q_AB;
G_BA = inv(-Q_BB) * Q_BA;

% calculate equilibrium occupancies
S = [Q, ones(k,1)];
p_inf = ones(1,k) * inv(S * S');

% calculate initial state vectors
pF_inf = p_inf((kA+1):k);
phi_o = pF_inf * Q_FA;
phi_o = phi_o / sum(phi_o);

phi_s = phi_o * G_AF;

pC_inf = p_inf((kE+1):k);
phi_b = pC_inf * (Q_CB * G_BA + Q_CA);
phi_b = phi_b / sum(phi_b);

% calculate burst-end vector
e_b = (eye(kA) - G_AB * G_BA) * ones(kA, 1);

% f_open(t) = sum(phi_o * exp(Q_AA * t) * Q_AF)
% f_shut(t) = sum(phi_s * exp(Q_FF * t) * Q_FA)
% f_b(t) = phi_b * exp(Q_EE * t)_AA * (-Q_AA) * e_b
f = @(x,t) exp(x*t);
expQA = sylv(Q_AA, f, t);
expQF = sylv(Q_FF, f, t);
expQE = sylv(Q_EE, f, t);

D_o = zeros(size(t));
D_s = zeros(size(t));
D_b = zeros(size(t));

for ti = 1:length(t)
    D_o(ti) = sum(phi_o * expQA{ti} * Q_AF);
    D_s(ti) = sum(phi_s * expQF{ti} * Q_FA);
    
    tmp = expQE{ti};
    D_b(ti) = phi_b * tmp(1:kA,1:kA) * (-Q_AA) * e_b;
end