cache_miss_rmf.m

%{ @file cache_miss_rmf.m
 %  @brief Computes miss rates using the refined mean field (RMF) method
 %
 %  @author LINE Development Team
%}
 
%{
 % @brief Computes cache miss rates for RANDOM(m) replacement via refined mean field
 %
 % @details
 % This function computes global, per-user, and per-item miss rates for
 % multi-list caches with RANDOM(m) replacement using the DDPP mean field
 % approximation with 1/N correction (refined mean field). The cache
 % occupancy is computed from the aggregate request stream; per-user miss
 % rates weight the per-item miss probabilities by each user's rates.
 %
 % Reference:
 %   N. Gast, "Expected Values Estimated via Mean-Field Approximation are
 %   1/N-Accurate", Proc. ACM Meas. Anal. Comput. Syst., 2017.
 %
 % @par Syntax:
 % @code
 % [M, MU, MI, pi0] = cache_miss_rmf(gamma, m, lambda)
 % @endcode
 %
 % @par Parameters:
 % <table>
 % <tr><th>Name<th>Description
 % <tr><td>gamma<td>Item access factors (used for sizing only)
 % <tr><td>m<td>Cache capacity vector
 % <tr><td>lambda<td>Arrival rates per user per item per list
 % </table>
 %
 % @par Returns:
 % <table>
 % <tr><th>Name<th>Description
 % <tr><td>M<td>Global miss rate
 % <tr><td>MU<td>Per-user miss rate
 % <tr><td>MI<td>Per-item miss rate
 % <tr><td>pi0<td>Per-item miss probability
 % </table>
%}
function [M,MU,MI,pi0,tout,pi0_t,MU_t,xtraj] = cache_miss_rmf(gamma, m, lambda, tspan, x0init) %#ok<INUSL>
% Optional TSPAN = [t0, t1] and initial occupancy X0INIT request the
% transient mean-field trajectory: integrating the same drift F(x) that
% RMF_FIXED_POINT drives to steady state, over the finite window. Returns
% the time grid TOUT, the per-item list-0 occupancy PI0_T (n_items x nt)
% and the per-user miss-rate trajectory MU_T (u x nt). Omitting TSPAN
% preserves the original steady-state-only contract.
if nargin < 4, tspan = []; end
if nargin < 5, x0init = []; end
u = size(lambda,1);
n_items = size(lambda,2);
h = length(m);
m = m(:)';
 
% aggregate per-item request rates over users
lam_i = zeros(1, n_items);
for v = 1:u
    row = lambda(v,:,1);
    row(~isfinite(row)) = 0;
    lam_i = lam_i + row;
end
p = lam_i / sum(lam_i);
 
model_dim = n_items * (h + 1);
 
% Build initial state: first m(1) items in list 1, next m(2) in
% list 2, etc.; remaining items outside cache (list 0)
x0 = zeros(model_dim, 1);
obj_idx = 0;
for k = 1:h
    for jj = 1:m(k)
        obj_idx = obj_idx + 1;
        if obj_idx <= n_items
            x0(rmf_index(obj_idx, k, n_items)) = 1.0;
        end
    end
end
for i = (obj_idx + 1):n_items
    x0(rmf_index(i, 0, n_items)) = 1.0;
end
 
% Mean field fixed point, then 1/N refinement when computable. The
% refinement solves a reduced Jacobian/Lyapunov system that turns singular
% for small or non-hyperbolic fixed points; MATLAB signals this with a
% (non-error) warning and returns non-finite entries rather than throwing,
% so a bare try/catch does not catch it. Accept the 1/N correction only
% when it is finite, otherwise keep the plain mean-field fixed point.
xss = rmf_fixed_point(x0, p, m, n_items, h, model_dim);
% The refinement's reduced linear system is expected to be singular for
% small/non-hyperbolic fixed points; silence only that specific warning here
% since the non-finite result is detected and rejected below.
ws1 = warning('off', 'MATLAB:singularMatrix');
ws2 = warning('off', 'MATLAB:nearlySingularMatrix');
restoreWarn = onCleanup(@() warning([ws1 ws2]));
try
    [pi_mf, V] = rmf_expansion_steady_state(x0, p, m, n_items, h, model_dim);
    xref = pi_mf + V / n_items;
    if all(isfinite(xref))
        xss = xref;
    end
catch
    % Fall back to plain mean field
end
clear restoreWarn
 
% Per-item miss probability (occupancy of list 0), clipped to [0,1]
pi0 = zeros(n_items, 1);
for i = 1:n_items
    pi0(i) = max(0, min(1, xss(rmf_index(i, 0, n_items))));
end
 
MI = lam_i(:) .* pi0;
MU = zeros(1, u);
for v = 1:u
    row = lambda(v,:,1);
    row(~isfinite(row)) = 0;
    MU(v) = row * pi0;
end
M = sum(MI);
 
% Transient mean-field trajectory (optional): integrate the drift over the
% requested window from the supplied (or default) initial occupancy. XTRAJ is
% the full DDPP occupancy trajectory (model_dim x nt), used to carry the cache
% mean occupancy across environment switches.
tout = []; pi0_t = []; MU_t = []; xtraj = [];
if ~isempty(tspan)
    if isempty(x0init)
        x0init = x0;
    else
        x0init = x0init(:);
    end
    ode_func = @(t, x) rmf_drift(x, p, m, n_items, h, model_dim);
    odeopt = odeset('AbsTol', 1e-10, 'RelTol', 1e-8);
    [tout, xtraj] = ode15s(ode_func, tspan, x0init, odeopt);
    xtraj = xtraj'; % model_dim x nt
    nt = numel(tout);
    pi0_t = zeros(n_items, nt);
    for i = 1:n_items
        pi0_t(i,:) = max(0, min(1, xtraj(rmf_index(i, 0, n_items), :)));
    end
    MU_t = zeros(u, nt);
    for v = 1:u
        row = lambda(v,:,1);
        row(~isfinite(row)) = 0;
        MU_t(v,:) = row * pi0_t;
    end
end
end
 
function idx = rmf_index(i, k, n_items)
% RMF_INDEX Map (item i, list k) to flat state index
%   i: item index (1-based)
%   k: list index (0 = outside cache, 1..h = cache lists)
%   n_items: total number of items
idx = i + k * n_items;
end
 
function hr = rmf_hit_rate(x, p, list_number, n_items)
% RMF_HIT_RATE Compute hit rate contribution from a specific list
%   hr = sum_i p(i) * x(index(i, list_number))
hr = 0.0;
for i = 1:n_items
    hr = hr + p(i) * x(rmf_index(i, list_number, n_items));
end
end
 
function dX = rmf_drift(x, p, m, n_items, h, model_dim)
% RMF_DRIFT Compute mean field drift F(x) for RANDOM(m) replacement
%
% dx[i,k]/dt = -p(i)*x[i,k] + hitRate[k]*x[i,k+1]/m(k)  (promotion from k to k+1)
% dx[i,k+1]/dt = p(i)*x[i,k] - hitRate[k]*x[i,k+1]/m(k)
 
hit_rates = zeros(1, h + 1);
for k = 0:h
    hit_rates(k + 1) = rmf_hit_rate(x, p, k, n_items);
end
 
dX = zeros(model_dim, 1);
for i = 1:n_items
    for k = 0:(h - 1)
        flow = p(i) * x(rmf_index(i, k, n_items)) ...
             - hit_rates(k + 1) * x(rmf_index(i, k + 1, n_items)) / m(k + 1);
        dX(rmf_index(i, k, n_items))     = dX(rmf_index(i, k, n_items))     - flow;
        dX(rmf_index(i, k + 1, n_items)) = dX(rmf_index(i, k + 1, n_items)) + flow;
    end
end
end
 
function Fp = rmf_jacobian(x, p, m, n_items, h, model_dim)
% RMF_JACOBIAN Compute Jacobian dF/dx at state x
 
hit_rates = zeros(1, h + 1);
for k = 0:h
    hit_rates(k + 1) = rmf_hit_rate(x, p, k, n_items);
end
 
Fp = zeros(model_dim, model_dim);
for i = 1:n_items
    for k = 0:(h - 1)
        ik  = rmf_index(i, k, n_items);
        ik1 = rmf_index(i, k + 1, n_items);
 
        % Direct rate terms
        Fp(ik,  ik)  = Fp(ik,  ik)  - p(i);
        Fp(ik1, ik)  = Fp(ik1, ik)  + p(i);
        Fp(ik,  ik1) = Fp(ik,  ik1) + hit_rates(k + 1) / m(k + 1);
        Fp(ik1, ik1) = Fp(ik1, ik1) - hit_rates(k + 1) / m(k + 1);
 
        % Indirect terms via hit rate dependence on x(j,k)
        for j = 1:n_items
            jk  = rmf_index(j, k, n_items);
            jk1 = rmf_index(j, k + 1, n_items);
            Fp(ik,  jk1) = Fp(ik,  jk1) - p(i) * x(ik) / m(k + 1);
            Fp(ik1, jk1) = Fp(ik1, jk1) + p(i) * x(ik) / m(k + 1);
            Fp(ik,  jk)  = Fp(ik,  jk)  + p(j) * x(ik1) / m(k + 1);
            Fp(ik1, jk)  = Fp(ik1, jk)  - p(j) * x(ik1) / m(k + 1);
        end
    end
end
end
 
function Fpp = rmf_hessian(~, p, m, n_items, h, model_dim)
% RMF_HESSIAN Compute Hessian d^2F/dx^2 (constant for this quadratic drift)
 
Fpp = zeros(model_dim, model_dim, model_dim);
for i = 1:n_items
    for k = 0:(h - 1)
        ik  = rmf_index(i, k, n_items);
        ik1 = rmf_index(i, k + 1, n_items);
        for j = 1:n_items
            if j ~= i
                jk  = rmf_index(j, k, n_items);
                jk1 = rmf_index(j, k + 1, n_items);
                % d^2 F[ik] / (d x[jk] d x[ik1])
                Fpp(ik,  jk,  ik1) = Fpp(ik,  jk,  ik1) + p(j) / m(k + 1);
                Fpp(ik,  ik1, jk)  = Fpp(ik,  ik1, jk)  + p(j) / m(k + 1);
                % d^2 F[ik] / (d x[jk1] d x[ik])
                Fpp(ik,  jk1, ik)  = Fpp(ik,  jk1, ik)  - p(i) / m(k + 1);
                Fpp(ik,  ik,  jk1) = Fpp(ik,  ik,  jk1) - p(i) / m(k + 1);
                % Symmetric for ik1
                Fpp(ik1, jk,  ik1) = Fpp(ik1, jk,  ik1) - p(j) / m(k + 1);
                Fpp(ik1, ik1, jk)  = Fpp(ik1, ik1, jk)  - p(j) / m(k + 1);
                Fpp(ik1, jk1, ik)  = Fpp(ik1, jk1, ik)  + p(i) / m(k + 1);
                Fpp(ik1, ik,  jk1) = Fpp(ik1, ik,  jk1) + p(i) / m(k + 1);
            end
        end
    end
end
end
 
function Q = rmf_noise_matrix(x, p, m, n_items, h, model_dim)
% RMF_NOISE_MATRIX Compute noise intensity matrix Q(x) for the DDPP
%
% Q[a,b] = sum_ell ell[a]*ell[b]*beta_ell(x)
% Each transition swaps items i and j across lists k and k+1.
 
Q = zeros(model_dim, model_dim);
signs = [-1, 1, 1, -1];
for i = 1:n_items
    for k = 0:(h - 1)
        for j = 1:n_items
            rate = p(i) * x(rmf_index(i, k, n_items)) ...
                 * x(rmf_index(j, k + 1, n_items)) / m(k + 1);
            indices = [rmf_index(i, k, n_items), ...
                       rmf_index(j, k, n_items), ...
                       rmf_index(i, k + 1, n_items), ...
                       rmf_index(j, k + 1, n_items)];
            for ia = 1:4
                for ib = 1:4
                    Q(indices(ia), indices(ib)) = Q(indices(ia), indices(ib)) ...
                        + rate * signs(ia) * signs(ib);
                end
            end
        end
    end
end
end
 
function pi = rmf_fixed_point(x0, p, m, n_items, h, model_dim)
% RMF_FIXED_POINT Compute mean field fixed point by ODE integration
%
% Integrates dx/dt = F(x) until steady state using ode15s.
 
tmax = 10000;
ode_func = @(t, x) rmf_drift(x, p, m, n_items, h, model_dim);
odeopt = odeset('AbsTol', 1e-10, 'RelTol', 1e-8);
[~, xvec] = ode15s(ode_func, [0, tmax], x0, odeopt);
%[~, xvec] = lsoda_solve(ode_func, [0, tmax], x0, odeopt);
pi = xvec(end, :)';
end
 
function [C, Cinv, rk] = rmf_dimension_reduction(Fp, n_items, h, model_dim)
% RMF_DIMENSION_REDUCTION Compute change-of-basis for singular Jacobian
%
% The Jacobian is singular because item populations are conserved
% (sum over lists for each item = 1). Returns matrices to project
% onto the non-singular subspace.
 
rk = rank(Fp);
 
C = zeros(model_dim, model_dim);
d = 0;
for l_idx = 0:h
    for i = 1:(n_items - 1)
        d = d + 1;
        C(d, rmf_index(i, l_idx, n_items)) = 1.0;
    end
end
 
[U, ~, ~] = svd(Fp);
C((rk + 1):model_dim, :) = U(:, (rk + 1):model_dim)';
Cinv = inv(C);
end
 
function [pi, V] = rmf_expansion_steady_state(x0, p, m, n_items, h, model_dim)
% RMF_EXPANSION_STEADY_STATE Compute refined mean field steady-state expansion
%
% Computes the mean field fixed point pi and the 1/N correction V
% using the Lyapunov equation approach with dimension reduction.
%
% The refined approximation for a system of N items is:
%   E[X] ~ pi + V/N + O(1/N^2)
 
pi = rmf_fixed_point(x0, p, m, n_items, h, model_dim);
 
Fp  = rmf_jacobian(pi, p, m, n_items, h, model_dim);
Fpp = rmf_hessian(pi, p, m, n_items, h, model_dim);
Q   = rmf_noise_matrix(pi, p, m, n_items, h, model_dim);
 
% Dimension reduction: project onto non-singular subspace
[C, Cinv, rk] = rmf_dimension_reduction(Fp, n_items, h, model_dim);
 
Fp_r = (C * Fp * Cinv);
Fp_r = Fp_r(1:rk, 1:rk);
 
% Reduce Hessian: Fpp_r(a,b,c) = sum_{i,j,k} C(a,i)*Fpp(i,j,k)*Cinv(j,b)*Cinv(k,c)
% First contraction: tmp1(a,j,k) = sum_i C(a,i)*Fpp(i,j,k)
tmp1 = zeros(model_dim, model_dim, model_dim);
for a = 1:rk
    for j = 1:model_dim
        for k = 1:model_dim
            tmp1(a, j, k) = C(a, :) * Fpp(:, j, k);
        end
    end
end
% Second contraction: tmp2(a,b,k) = sum_j tmp1(a,j,k)*Cinv(j,b)
tmp2 = zeros(rk, rk, model_dim);
for a = 1:rk
    for b = 1:rk
        for k = 1:model_dim
            tmp2(a, b, k) = tmp1(a, :, k) * Cinv(:, b);
        end
    end
end
% Third contraction: Fpp_r(a,b,c) = sum_k tmp2(a,b,k)*Cinv(k,c)
Fpp_r = zeros(rk, rk, rk);
for a = 1:rk
    for b = 1:rk
        Fpp_r(a, b, :) = reshape(tmp2(a, b, :), 1, []) * Cinv(:, 1:rk);
    end
end
 
Q_r = C * Q * C';
Q_r = Q_r(1:rk, 1:rk);
 
% Solve Lyapunov equation: Fp_r * W_r + W_r * Fp_r' + Q_r = 0
W_r = lyap(Fp_r, Q_r);
 
% First-order correction: V_r = -Fp_r \ (C_r / 2)
% where C_r = sum_{b,c} Fpp_r(:,b,c) * W_r(b,c)
C_r = zeros(rk, 1);
for a = 1:rk
    for b = 1:rk
        for c = 1:rk
            C_r(a) = C_r(a) + Fpp_r(a, b, c) * W_r(b, c);
        end
    end
end
V_r = -Fp_r \ (C_r / 2.0);
 
% Expand back to full dimension
V = Cinv(:, 1:rk) * V_r;
end