getSensitivity.m

function [S, SS, dpi, pi] = getSensitivity(self, param, reward)
% [S, SS, DPI, PI] = GETSENSITIVITY(PARAM, REWARD)
%
% Parametric sensitivity of a steady-state reward to a scalar model
% parameter, following Trivedi and Bobbio (2017), Sec. 9.7.
%
% PARAM describes the parameter theta and how to set it on the model:
%   param.name    identifier used in reports
%   param.value   nominal value theta
%   param.set     handle (model, value) -> void, applying theta to the model
%   param.step    optional finite-difference step, default value*1e-6
%
% The generator derivative dQ/dtheta is obtained by central differences on
% the rate with the state space held fixed. This is exact to O(step^2) and
% requires no symbolic differentiation of the rate assembly; the state space
% is unaffected because it depends on the topology and the cutoff, not on
% rate values. The steady-state sensitivity then follows from one linear
% solve, see ctmc_sens.
%
% REWARD is either a function handle mapping the state space to a reward rate
% vector, or a numeric reward rate vector over the states. If omitted, dpi is
% returned and S is empty.
%
% Note: this returns d(E[r])/dtheta with dr/dtheta = 0, i.e. it assumes the
% reward rates do not themselves depend on theta. Rewards that depend on
% theta need the second term of Eq. (9.83) and are not handled here.
%
% @param param Struct describing the parameter, see above
% @param reward Reward rate vector or handle over the state space (optional)
% @return S Unscaled sensitivity d(E[r])/dtheta, Eq. (9.79)
% @return SS Scaled sensitivity (theta/E[r]) d(E[r])/dtheta, Eq. (9.80)
% @return dpi Sensitivity of the steady-state distribution (1 x n)
% @return pi Steady-state distribution (1 x n)
%
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
if ~isstruct(param) || ~isfield(param, 'set') || ~isfield(param, 'value')
    line_error(mfilename, 'param must be a struct with fields value and set');
end
if nargin < 3
    reward = [];
end
 
theta = param.value;
if isfield(param, 'step') && ~isempty(param.step)
    h = param.step;
else
    h = max(abs(theta), 1) * 1e-6;
end
 
% Nominal generator and state space
[Q, ~, ~] = self.getGenerator();
Q = full(Q);
space = self.getStateSpace();
n = size(Q, 1);
 
% Central differences on theta with the state space fixed
Qp = perturbedGenerator(self, param, theta + h);
Qm = perturbedGenerator(self, param, theta - h);
if size(Qp, 1) ~= n || size(Qm, 1) ~= n
    line_error(mfilename, ['Perturbing the parameter changed the state space size, so the ', ...
        'generators cannot be differenced. This happens when the parameter switches a ', ...
        'transition on or off (e.g. a zero rate or an immediate transition).']);
end
dQ = (Qp - Qm) / (2 * h);
 
pi = ctmc_solve(Q);
dpi = ctmc_sens(Q, dQ, pi);
 
S = [];
SS = [];
if isempty(reward)
    return;
end
 
if isa(reward, 'function_handle')
    r = reward(space);
else
    r = reward;
end
r = r(:);
if length(r) ~= n
    line_error(mfilename, 'reward must have one entry per state');
end
 
% Eq. (9.83) with dr/dtheta = 0
S = dpi * r;
Er = pi * r;
if abs(Er) > GlobalConstants.Zero
    SS = (theta / Er) * S;
else
    SS = NaN;
end
end
 
function Q = perturbedGenerator(self, param, value)
% Q = PERTURBEDGENERATOR(SELF, PARAM, VALUE)
% Rebuild the generator with theta set to VALUE, on a copy of the model so
% the caller's model is left untouched.
%
% The hard refresh is required, not defensive: setService and setArrival
% deliberately leave the cached struct in place, so a copy that inherited a
% built struct would report the old rate and the difference quotient would
% silently come out as zero.
 
modelCopy = self.model.copy();
param.set(modelCopy, value);
modelCopy.refreshStruct(true);
solverCopy = SolverCTMC(modelCopy, self.getOptions());
Q = full(solverCopy.getGenerator());
end