solver_fld_cacheqn_analyzer.m

function [QN, UN, RN, TN, CN, XN, t, QNt, UNt, TNt, xvec_iter, hitprob, missprob, runtime, it] = solver_fld_cacheqn_analyzer(sn, options)
% SOLVER_FLD_CACHEQN_ANALYZER Fluid solver for integrated caching-queueing networks
%
% Delegates the decomposition-aggregation alternation between the isolated
% caches and the fluid ODE solution of the surrounding queueing network to
% da_cacheqn, until the cache arrival rates converge. RANDOM(m) caches use
% the refined mean field miss algorithm (cache_miss_rmf, 1/N-accurate);
% other replacement strategies use the FPI approximation (cache_miss_fpi).
 
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
 
T0 = tic;
K = sn.nclasses;
 
[res, hitprob, missprob, it, sn] = da_cacheqn(sn, @miss_isolated, @netsolve, options);
QN = res.QN; UN = res.UN; RN = res.RN; TN = res.TN; XN = res.XN;
t = res.t; QNt = res.QNt; UNt = res.UNt; TNt = res.TNt; xvec_iter = res.xvec_iter;
 
% Compute CN
CN = zeros(1, K);
for k = 1:K
    if sn.refstat(k) > 0
        CN(k) = sn.njobs(k) ./ XN(k);
    end
end
 
% xvec_iter is already a cell array from solver_fluid_matrix
runtime = toc(T0);
 
    function missrate = miss_isolated(gamma, m, lambda_cache, ch)
        if ch.replacestrat == ReplacementStrategy.RR
            % RANDOM(m) replacement: refined mean field (1/N-accurate)
            [~, missrate] = cache_miss_rmf(gamma, m, lambda_cache);
        else
            [~, missrate] = cache_miss_fpi(gamma, m, lambda_cache);
        end
    end
 
    function res = netsolve(snit)
        % Solve the queueing network using the fluid matrix method
        res = struct();
        fluid_options = options;
        fluid_options.method = 'matrix';
        fluid_options.init_sol = solver_fluid_initsol(snit, fluid_options);
        [res.QN, res.UN, res.RN, res.TN, res.xvec_iter, res.QNt, res.UNt, res.TNt, ~, res.t] = solver_fluid_matrix(snit, fluid_options);
 
        % Compute system throughputs
        res.XN = zeros(1, K);
        for k = 1:K
            if snit.refstat(k) > 0
                res.XN(k) = res.TN(snit.refstat(k), k);
            end
        end
    end
end