retrieval_fpi.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_RETRIEVAL_RETRIEVAL_FPI_H
#define LINE_API_RETRIEVAL_RETRIEVAL_FPI_H
 
/**
 * @file
 * @ingroup api_retrieval
 * Fixed-point heuristic for a delayed-hit (list-based) cache.
 *
 * Templated port of matlab/src/api/retrieval/retrieval_fpi.m, cross-checked
 * against jar/src/main/java/jline/api/retrieval/Retrieval_fpi.java.
 *
 * The exact recursion of retrieval_metrics costs O(2^n); the heuristic
 * truncates the perturbation expansion at zeroth order,
 * pi_{i,l}(m-1_j) ~ pi_{i,l}(m), and solves the resulting nonlinear system by
 * successive substitution (paper eqs. theta/xi/pi/pi_0/phi_FPI). One sweep is
 *
 *   F_{s,i}   = 1 + sum_{k!=i} phi_{s,k}
 *   D_i       = 1 + lambda_i eta_{0,i} + sum_s lambda_i eta_{s,i} F_{s,i}
 *   theta_ij  = gamma_ij / D_i
 *   xi        : solved JOINTLY from m_j = sum_i theta_ij xi_j / (1 + sum_l theta_il xi_l),
 *               i.e. cache_xi_fp(theta, m), the no-delay cache fixed point with gamma := theta
 *   A_i       = 1 / (1 + sum_l theta_il xi_l)
 *   pi_ij     = theta_ij xi_j A_i
 *   pi_i0     = A_i / D_i
 *   phi_{0,i} = A_i lambda_i eta_{0,i} / D_i,  phi_{s,i} = A_i lambda_i eta_{s,i} F_{s,i} / D_i
 *
 * starting from phi = 0 and pi from the cache WITHOUT delayed hits
 * (cache_xi_fp(gamma, m), pi_i0 = 1 - sum_j pi_ij), and stopping when the
 * largest relative change of the miss, hit and delayed-hit ratios falls below
 * tol. Every iterate is feasible, sum_i pi_ij = m_j. When r = 0 (no PS
 * station) D does not depend on phi, so one sweep is the fixed point and the
 * iteration stops after it.
 *
 * ARITHMETIC: the iteration stops on a tolerance, so its answer is the fixed
 * point only to within tol whatever the arithmetic; it is gated on
 * has_transcendental so nobody instantiates it at exact arithmetic expecting
 * an exact result. Use retrieval_metrics or retrieval_mva for that.
 */
 
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <iostream>
#include <string>
#include <vector>
 
#include "line/api/cache/cache_xi_fp.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/line_console.h"
#include "line/util/matrix.h"
 
namespace line {
namespace retrieval {
 
/** Options mirroring the trailing (max_iter, tol) arguments of the MATLAB function. */
struct FpiOptions {
    std::size_t max_iter = 1000;
    double tol = 1e-6;
};
 
/** Mirrors the [pmiss, phit, pdh] return list, plus the iteration diagnostics. */
template <class T>
struct RetrievalFpiResult {
    std::vector<T> pmiss;  ///< (n) miss ratios pi_{i,0}
    Matrix<T> phit;        ///< (h x n) hit ratios pi_{i,j}
    Matrix<T> pdh;         ///< ((r+1) x n) delayed-hit probabilities phi_{s,i}, s = 0..r
    std::size_t iterations = 0;
    bool converged = false;
    /// True when the iterate left the finite range; MATLAB warns and breaks here.
    bool diverged = false;
    /// The reference's line_warning text (divergence or iteration cap), empty when it did not warn.
    std::string warning;
};
 
namespace detail {
 
/** MATLAB's local reldiff: max|a-b| over max|b|, with an all-zero b read as 1. */
template <class T>
double fpi_reldiff(const std::vector<T>& a, const std::vector<T>& b) {
    double num = 0.0, den = 0.0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        const double d = num_traits<T>::to_double(num_abs(T(a[i] - b[i])));
        const double e = num_traits<T>::to_double(num_abs(T(b[i])));
        if (d > num) num = d;
        if (e > den) den = e;
    }
    if (den == 0.0) den = 1.0;
    return num / den;
}
 
}  // namespace detail
 
/**
 * @brief Fixed-point heuristic for a delayed-hit (list-based) cache.
 *
 * @param m      (h) cache list capacities
 * @param lambda (n) per-item arrival rates
 * @param eta    (n x (r+1)) fetching demands, column 0 = IS station, columns 1..r = PS
 * @param gamma  (n x h) access factors
 * @param options fixed-point options (tolerance, iteration cap, damping)
 */
template <class T>
RetrievalFpiResult<T> retrieval_fpi(const std::vector<int>& m, const std::vector<T>& lambda,
                                    const Matrix<T>& eta, const Matrix<T>& gamma,
                                    const FpiOptions& options = FpiOptions()) {
    static_assert(num_traits<T>::has_transcendental,
                  "retrieval_fpi requires transcendental arithmetic: it is a successive "
                  "substitution stopped on a relative tolerance, so its answer is the fixed "
                  "point only to within tol whatever the arithmetic");
    const std::size_t n = lambda.size();
    const std::size_t h = m.size();
    if (eta.rows() != n || gamma.rows() != n)
        throw InputError("retrieval_fpi: eta/gamma and lambda disagree on the item count");
    if (gamma.cols() != h)
        throw InputError("retrieval_fpi: gamma and m disagree on the number of lists");
    if (eta.cols() == 0) throw InputError("retrieval_fpi: eta has no columns");
    const std::size_t r = eta.cols() - 1;
 
    const T one = num_traits<T>::from_int(1);
    const T zero = num_traits<T>::from_int(0);
 
    // initial guess of the paper: phi = 0 and pi from the cache without delayed
    // hits, so that sum_i pi_ij = m_j holds from the start.
    // phi is stored (r+1) x n and pij is stored h x n, as MATLAB returns them.
    Matrix<T> phi(r + 1, n, zero);
    Matrix<T> pij(h, n, zero);
    std::vector<T> pi0(n, one);
    {
        const cache::CacheXiFpResult<T> c0 = cache::cache_xi_fp(gamma, m);
        for (std::size_t i = 0; i < n; ++i) {
            T sum_pij = zero;
            for (std::size_t j = 0; j < h; ++j) {
                pij(j, i) = c0.pij(i, j);
                sum_pij += pij(j, i);
            }
            pi0[i] = one - sum_pij;
        }
    }
 
    RetrievalFpiResult<T> out;
    double last_delta = 0.0;
    for (std::size_t t = 1; t <= options.max_iter; ++t) {
        out.iterations = t;
 
        // F(s,i) = 1 + sum_{k != i} phi_{s,k}
        Matrix<T> F(r, n, one);
        for (std::size_t s = 0; s < r; ++s) {
            T tot = zero;
            for (std::size_t i = 0; i < n; ++i) tot += phi(s + 1, i);
            for (std::size_t i = 0; i < n; ++i) F(s, i) = one + (tot - phi(s + 1, i));
        }
 
        std::vector<T> D(n);
        for (std::size_t i = 0; i < n; ++i) {
            D[i] = one + lambda[i] * eta(i, 0);
            for (std::size_t s = 0; s < r; ++s) D[i] += lambda[i] * eta(i, s + 1) * F(s, i);
        }
 
        Matrix<T> theta(n, h);
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < h; ++j) theta(i, j) = gamma(i, j) / D[i];
 
        // xi solves the capacity constraints jointly: the no-delay cache fixed point with gamma := theta
        const std::vector<T> xi = cache::cache_xi_fp(theta, m).xi;
 
        Matrix<T> pij_new(h, n);
        std::vector<T> pi0_new(n);
        Matrix<T> phi_new(r + 1, n);
        for (std::size_t i = 0; i < n; ++i) {
            T denom = one;
            for (std::size_t l = 0; l < h; ++l) denom += theta(i, l) * xi[l];
            const T A = one / denom;
            for (std::size_t j = 0; j < h; ++j) pij_new(j, i) = theta(i, j) * xi[j] * A;
            pi0_new[i] = A / D[i];
            phi_new(0, i) = A * lambda[i] * eta(i, 0) / D[i];
            for (std::size_t s = 0; s < r; ++s)
                phi_new(s + 1, i) = A * lambda[i] * eta(i, s + 1) * F(s, i) / D[i];
        }
 
        const std::vector<T> pij_flat(pij_new.data(), pij_new.data() + pij_new.size());
        const std::vector<T> pij_old(pij.data(), pij.data() + pij.size());
        const std::vector<T> phi_flat(phi_new.data(), phi_new.data() + phi_new.size());
        const std::vector<T> phi_old(phi.data(), phi.data() + phi.size());
        double delta = detail::fpi_reldiff(pi0_new, pi0);
        const double d2 = detail::fpi_reldiff(pij_flat, pij_old);
        const double d3 = detail::fpi_reldiff(phi_flat, phi_old);
        if (d2 > delta) delta = d2;
        if (d3 > delta) delta = d3;
 
        pij = pij_new;
        pi0 = pi0_new;
        phi = phi_new;
        last_delta = delta;
 
        if (!std::isfinite(delta)) {
            out.diverged = true;
            out.warning = "FPI diverged (non-finite iterate) at iteration " + std::to_string(t) + ".";
            break;
        }
        // r == 0: D does not depend on phi, so one sweep is the fixed point
        if (delta < options.tol || r == 0) {
            out.converged = true;
            break;
        }
    }
 
    if (!out.converged && !out.diverged && options.max_iter > 0) {
        char buf[64];
        std::snprintf(buf, sizeof(buf), "%.2e", last_delta);
        out.warning = "FPI did not converge within " + std::to_string(options.max_iter) +
                      " iterations (last rel. change " + buf + ").";
    }
    // line_warning(mfilename, ...) in the reference: printed unless the session is SILENT
    if (!out.warning.empty() && util::LineConsole::get_verbose() != util::VerboseLevel::SILENT)
        std::cerr << "[LINE] Warning: retrieval_fpi: " << out.warning << std::endl;
 
    out.pmiss = pi0;
    out.phit = pij;
    out.pdh = phi;
    return out;
}
 
}  // namespace retrieval
}  // namespace line
 
#endif  // LINE_API_RETRIEVAL_RETRIEVAL_FPI_H