function sense = periodicSense(P,t,h,n0,cycle,Opt)
% This function calculates the transient sensitivity of the total
% population size at time t after an initial perturbation in the elements
% of the phase h (i.e. P{cycle(h)})
%
% P     = an array of cells, where each cell contains a transition matrix
%           for one phase of a cycle (the matrices in P can be different
%           sizes provided they tally correctly for the specified cycle)
% t     = the time after the initial perturbation
% h     = the sensitivities for a perturbation in phase h
% n0    = the initial population vector
% cycle = is a vector giving the cyclical order of each phase.
%         The default is [1,2,3, ...,H], where H=number of phases
%          e.g. cycle = [1 3 2 4] means that the order of the
%            transition matrices in a cycle is P{1} followed by P{3}
%            and then P{2} and P{4}, so that the transition matrix for the
%            whole cycle is P{4}*P{2}*P{3}*P{1}
% Opt    = a structure that gives different output options
%          Opt.elasticity = true  Calculate elasticities
%          Opt.immediate = true   Calculate immediate sensitivities
%          Opt.timescaled = true   Calculate time-scaled sensitivities
%            [default for all is false]
%
% OUTPUT
% sense(i,j) = the transient sensitivity of the total population size a
%                time t after a perturbation of element P{cycle(h)}(i,j)
%                Depending on output options this could be within-cycle,
%                elasticity, etc


% Written by Jon Yearsley and Shana Mertens, Jan 2006
%       jonathan.yearsley@unil.ch
%       shana.mertens@bbsrc.ac.uk

if nargin<4,
    % Check that there are the right number of inputs
    error('periodicSense requires at least 4 input arguments: P, t, h and n0')
elseif nargin<5,
    nPhases = length(P);     % Number of phases
    cycle = 1:nPhases;       % Set default cycle
    Opt.elasticity = false;  % Set default options
    Opt.immediate = false;
    Opt.timescaled = false;
elseif nargin<6,
    nPhases = length(cycle);
    Opt.elasticity = false;  % Set default options
    Opt.immediate = false;
    Opt.timescaled = false;
else
    nPhases = length(cycle);
    % Check the options are set properly
    if ~isfield(Opt,'elasticity')
        warning('Setting elasticity option to default')
        Opt.elasticity = false;
    end
    if ~isfield(Opt,'immediate')
        warning('Setting immediate option to default')
        Opt.immediate = false;
    end
    if ~isfield(Opt,'timescaled')
        warning('Setting timescaled option to default')
        Opt.timescaled = false;
    end
end


if h>length(cycle),
    % Check that the required phaseis consistent with the cycle information
    error(['The argument h is too large. There are only ' ...
        num2str(nPhases) ' in the current cycle'])
end

nStages = zeros(length(P),2);
for p=1:length(P),
    nStages(p,:) = size(P{p});  % Number of stages in each transition matrix
end

% Finally check that the cycle is possible given the dimension of the
% matrices in P
if any(nStages(cycle,1) ~= nStages(cycle([2:nPhases,1]),2)),
    error('The matrix dimensions of P are wrong for the given cyclic permutation')
end

% Initialise the matrix to contain the sensitivities
sense = zeros(nStages(h,:));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Calculate decomposition of sensitivity
%
% Now calculate D^{(+h)} and D^{(-h)} (equations 9 and 10 in the paper)
% which when multiplied together give Dh
% Dh multiplies the "cycle" sensitivities and gives us
% the sensitivity of a particular element in phase h

Dminus_h = eye(nStages(cycle(end),1));
for i=1:h-1,
    Dminus_h = P{cycle(i)} * Dminus_h;
end

Dplus_h = eye(nStages(cycle(h),1));
for i=h+1:nPhases,
    Dplus_h = P{cycle(i)} * Dplus_h;
end

% Finally create the transition matrix for an entire cycle
A = Dplus_h * P{cycle(h)} * Dminus_h;


% Calculate the deviation of n0 from the stable stage structure
% Calculate right eigenvectors and eigenvalues
optEig.disp = 0;
[w1,~] = eigs(A,1,'lm',optEig);
w1= abs(w1);  % Make sure sign of w1 is positive

% The equilibrium stage distribution
if length(n0)~=length(A),
    error('The number of stages in n0 does not agree with the model')
end
n0 = reshape(n0,length(A),1); % This ensures n0 is a column vector

n_infty = sum(n0)/sum(w1); % Use unstandardised right eigenvector because
% S_h is calculated with unstandardised
% eigenvectors
% Deviation of initial population structure from asymptotic stable stage
% structure
delta_n = n0 - sum(n0)/sum(w1)*w1;

if Opt.elasticity,
    %Calculating population size at time t
    n = A^t * n0;
    nTot=sum(n);  % Total population size
    
    if Opt.immediate,
        n_within = Dplus_h \ n;
        nTot = sum(n_within);
    end
end

for k=1:nStages(h,1),
    for l=1:nStages(h,2),
        % This is equation 8 in the Ecology paper
        Dh = Dplus_h(:,k) * Dminus_h(l,:);

        % Find the elements of Dh which are non-zero
        [iInd,jInd] = find(Dh);
        % Now do the same thing only break up the sensitivity into two parts, a la
        % Ecological modelling paper (Yearsley, 2004)
        S_h = zeros(length(n0),1);
        deltaS_h = zeros(length(n0));
        if ~isempty(iInd),
            for i=1:length(iInd);
                S_h = S_h + Dh(iInd(i),jInd(i)) * real(S(A,iInd(i),jInd(i),t));
                deltaS_h = deltaS_h + Dh(iInd(i),jInd(i)) * real(delta_S(A,iInd(i),jInd(i),t));
            end

            % And finally calculate the sensitivity of the population size
            % based upon the decomposition in Ecological Modelling paper
            % (there are other ways to do this without using a
            % decomposition)
            % Note that we're claculating the sensitivity of total
            % population size, hence the sum( * ,1)
            if Opt.immediate,
%                sense(k,l) = sum(Dplus_inv * (S_h * n_infty + deltaS_h * delta_n),1);
                sense(k,l) = sum(Dplus_h \ (S_h * n_infty + deltaS_h * delta_n),1);
            else
                sense(k,l) = sum(S_h * n_infty + deltaS_h * delta_n,1);
            end
            
        end
    end
end

if Opt.elasticity,
    % And elasticities, scaled for time
    sense = P{cycle(h)} .* sense / nTot;
end

if Opt.timescaled
    sense = sense / t;
end