lqn_focus.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_LQN_LQN_FOCUS_H
#define LINE_API_LQN_LQN_FOCUS_H
 
/**
 * @file
 * @ingroup api_lqn
 * Focused layered performance model by aggregation (FSPT).
 *
 * Port of `matlab/src/api/lqn/lqn_focus.m`, on the raw model `LqnModel<double>`
 * (the form LqnBuilder and the .lqnx reader produce). Reduces a layered network
 * to a focused model (FM) that keeps a set of preserved tasks and hosts as they
 * are and replaces every other task by one aggregate task per dependency group,
 * following Islam, Petriu and Woodside, "Focused Layered Performance Modelling
 * by Aggregation", ACM TOMPECS 7(2-4), 2022.
 *
 * The original model (OM) is solved once with SolverLN to obtain the
 * saturations S = U/m of every task and host. The preserved set is chosen by
 * the `sat` (forced elements plus every resource with S > Smax/r), `acc`
 * (forced elements plus the i most saturated resources, i grown until the
 * reference-task throughput error is below x) or `manual` (forced elements,
 * the most saturated resource and `preserve`) strategy. A preserved task
 * preserves its host and a preserved host all of its tasks. Tasks using
 * features FSPT does not aggregate (priorities, AND forks/joins, second
 * phases, asynchronous calls, forwarding, open arrivals, replication, fan-in
 * or fan-out, setup/delay-off and activity think times) are always preserved.
 *
 * Non-preserved tasks with the same set Pi(T) of preserved tasks depending on
 * them form a group, replaced by one aggregate task with one entry and one
 * activity (eqs. 7-11). A host shared by several groups is split into shadow
 * hosts of speed s_H (1 - saturation by the other groups) (eq. 15), or
 * preserved when that capacity is not positive. Aggregate processors have unit
 * speed, the aggregate speed of eq. 9 being folded into the aggregate demand.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <map>
#include <set>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>
 
#include "line/lang/lang_types.h"
#include "line/lang/lqn/lqn_reader.h"
#include "line/lang/lqn/lqn_struct.h"
#include "line/solvers/ln/solver_ln.h"
 
namespace line {
namespace lqn {
 
/** Options of lqn_focus, as the MATLAB name/value pairs. */
struct LqnFocusOptions {
    std::string strategy = "sat";        ///< "sat", "acc" or "manual"
    double r = 3.3;                      ///< SAT saturation ratio threshold
    double x = 0.02;                     ///< ACC relative throughput error target
    std::vector<std::string> preserve;   ///< task or host names forced into the focus
    /** Solver of a finalized struct; empty uses ln::SolverLN with default options. */
    std::function<ln::LnSolution<double>(const LqnStruct<double>&)> solver;
    bool validate = true;                ///< solve the FM and report errors (always on for "acc")
};
 
/** One dependency group and its aggregate task. */
struct LqnFocusGroup {
    std::vector<std::string> tasks, pi;
    std::string aggTask, aggEntry, aggActivity, aggHost;
    double mult = 0, hostMult = 0, hostSpeed = 1, demand = 0, rate = 0;
    std::vector<std::string> callDests;
    std::vector<double> callMeans;
};
 
/** OM and FM metrics of the retained tasks, entries and hosts. */
struct LqnFocusRetained {
    std::vector<std::string> names, types;
    std::vector<double> X_om, X_fm, R_om, R_fm, Q_om, Q_fm, U_om, U_fm;
};
 
/** Focused model and diagnostics, fields named as the MATLAB info struct. */
struct LqnFocusResult {
    LqnModel<double> fm;
    std::vector<std::string> preservedTasks, preservedHosts, forcedTasks, droppedTasks;
    std::vector<bool> preservedMask;  ///< (ntasks + nhosts) tasks then hosts
    std::vector<std::string> taskNames, hostNames;
    std::vector<double> taskSat, hostSat;
    std::vector<LqnFocusGroup> groups;
    /** {original, focused}; focused is "" for a dropped element. */
    std::vector<std::pair<std::string, std::string>> taskMap, entryMap, activityMap, hostMap;
    double SR = 0;
    double trustedRange[2] = {std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN()};
    double FT = 0, F = 0;
    double err = std::numeric_limits<double>::quiet_NaN();
    double errResp = std::numeric_limits<double>::quiet_NaN();
    std::vector<std::string> refTasks;
    std::vector<double> errRef, errRespRef;
    LqnFocusRetained retained;
    std::string strategy;
    double r = 0, x = 0;
    std::vector<std::vector<double>> accHistory;  ///< rows [i, nPreserved, F, err]
    double Smax = 0;
};
 
namespace focus_detail {
 
using namespace ::line::lang;
 
inline double nan_() { return std::numeric_limits<double>::quiet_NaN(); }
 
inline bool is_prio(SchedStrategy s) {
    std::string n = sched_to_text(s);
    for (auto& c : n) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    return n.find("prio") != std::string::npos || n == "hol";
}
 
inline double dmean(const Distrib<double>& d) {
    if (d.disabled || d.type == ProcessType::IMMEDIATE || d.type == ProcessType::DISABLED) return 0.0;
    return std::isnan(d.mean) ? 0.0 : d.mean;
}
 
inline bool has_time(const Distrib<double>& d) { return dmean(d) > GlobalConstants::FineTol; }
 
struct Call {
    std::size_t src, dst;
    double mean;
};
 
struct Inv {
    std::size_t nH = 0, nT = 0, nE = 0, nA = 0;
    std::vector<std::string> hName, tName, eName, aName;
    std::vector<double> hMult, hSpeed, tMult, aDem;
    std::vector<std::size_t> tHost, eTask, aTask;
    std::vector<bool> tIsRef, forcedT, forcedH;
    std::vector<Call> sync, async, fwd;
    std::vector<std::vector<bool>> reach;
};
 
struct Om {
    std::vector<double> hQ, hU, tQ, tU, tX, eQ, eU, eR, eX, aX;
};
 
struct Sat {
    std::vector<std::vector<double>> SHT;
    std::vector<double> SH, ST;
};
 
[[noreturn]] inline void fail(const std::string& msg) { throw std::invalid_argument("lqn_focus: " + msg); }
 
inline std::size_t index_of(const std::vector<std::string>& v, const std::string& s) {
    for (std::size_t i = 0; i < v.size(); ++i)
        if (v[i] == s) return i;
    return v.size();
}
 
inline Inv inventory(const LqnModel<double>& m) {
    Inv inv;
    inv.nH = m.procs.size();
    inv.nT = m.tasks.size();
    inv.nE = m.entries.size();
    inv.nA = m.acts.size();
    for (const auto& p : m.procs) inv.hName.push_back(p.name);
    for (const auto& t : m.tasks) inv.tName.push_back(t.name);
    for (const auto& e : m.entries) inv.eName.push_back(e.name);
    for (const auto& a : m.acts) inv.aName.push_back(a.name);
    for (std::size_t h = 0; h < inv.nH; ++h) {
        if (m.proc_linconrows.count(h) || m.proc_lincon.count(h) || m.proc_lldscaling.count(h) ||
            m.proc_cdscaling.count(h) || m.proc_jdscaling.count(h) || m.proc_pools.count(h))
            fail("Element \"" + inv.hName[h] + "\" uses admission constraints, load-dependent scaling or "
                 "server pools, which lqn_focus does not support.");
    }
    inv.hMult.resize(inv.nH);
    inv.hSpeed.resize(inv.nH);
    inv.forcedH.assign(inv.nH, false);
    for (std::size_t h = 0; h < inv.nH; ++h) {
        const auto& p = m.procs[h];
        inv.hMult[h] = p.mult;
        inv.hSpeed[h] = p.speed_factor;
        inv.forcedH[h] = p.repl > 1 || is_prio(p.sched);
    }
    inv.tMult.resize(inv.nT);
    inv.tHost.resize(inv.nT);
    inv.tIsRef.assign(inv.nT, false);
    inv.forcedT.assign(inv.nT, false);
    for (std::size_t t = 0; t < inv.nT; ++t) {
        const auto& tk = m.tasks[t];
        if (tk.nitems > 0) fail("Task \"" + tk.name + "\" is a cache task, which lqn_focus does not support.");
        if (!tk.linconrows.empty() || tk.lincon_A.rows() > 0 || !tk.lldscaling.empty() ||
            static_cast<bool>(tk.cdscaling) || static_cast<bool>(tk.jdscaling) || !tk.pools.empty())
            fail("Element \"" + tk.name + "\" uses admission constraints, load-dependent scaling or "
                 "server pools, which lqn_focus does not support.");
        if (tk.proc_slot >= inv.nH) fail("Task \"" + tk.name + "\" has no host.");
        inv.tMult[t] = tk.mult;
        inv.tHost[t] = tk.proc_slot;
        inv.tIsRef[t] = tk.sched == SchedStrategy::REF;
        bool f = inv.tIsRef[t] || tk.repl > 1 || tk.priority != 0 || is_prio(tk.sched) ||
                 has_time(tk.setuptime) || has_time(tk.delayofftime) || !tk.fanin.empty() || !tk.fanout.empty();
        for (const auto& pr : tk.precedences)
            if (pr.pretype == PrecedenceType::PRE_AND || pr.posttype == PrecedenceType::POST_AND) f = true;
        inv.forcedT[t] = f;
    }
    inv.eTask.resize(inv.nE);
    for (std::size_t e = 0; e < inv.nE; ++e) {
        const auto& en = m.entries[e];
        if (en.cardinality > 0) fail("Entry \"" + en.name + "\" is not a plain entry, which lqn_focus does not support.");
        inv.eTask[e] = en.task_slot;
    }
    auto dest = [&](const std::string& n) {
        std::size_t e2 = index_of(inv.eName, n);
        if (e2 == inv.nE) fail("Call destination \"" + n + "\" is not an entry of the model.");
        return e2;
    };
    inv.aTask.resize(inv.nA);
    inv.aDem.resize(inv.nA);
    for (std::size_t a = 0; a < inv.nA; ++a) {
        const auto& ac = m.acts[a];
        const std::size_t t = ac.task_slot;
        inv.aTask[a] = t;
        inv.aDem[a] = dmean(ac.hostdem);
        if (!ac.call_groups.empty())
            fail("Activity \"" + ac.name + "\" uses call groups, which lqn_focus does not support.");
        if (ac.phase > 1 || has_time(ac.thinktime)) inv.forcedT[t] = true;
        for (const auto& c : ac.sync_calls) inv.sync.push_back({a, dest(c.dest), c.mean});
        for (const auto& c : ac.async_calls) {
            const std::size_t e2 = dest(c.dest);
            inv.async.push_back({a, e2, c.mean});
            inv.forcedT[t] = true;
            inv.forcedT[inv.eTask[e2]] = true;
        }
    }
    for (std::size_t e = 0; e < inv.nE; ++e) {
        const auto& en = m.entries[e];
        if (en.has_arrival) inv.forcedT[inv.eTask[e]] = true;
        for (std::size_t k = 0; k < en.fwd_dest.size(); ++k) {
            const std::size_t e2 = dest(en.fwd_dest[k]);
            inv.fwd.push_back({e, e2, en.fwd_prob[k]});
            inv.forcedT[inv.eTask[e]] = true;
            inv.forcedT[inv.eTask[e2]] = true;
        }
    }
    auto& R = inv.reach;
    R.assign(inv.nT, std::vector<bool>(inv.nT, false));
    for (const auto& c : inv.sync) R[inv.aTask[c.src]][inv.eTask[c.dst]] = true;
    for (const auto& c : inv.async) R[inv.aTask[c.src]][inv.eTask[c.dst]] = true;
    for (const auto& c : inv.fwd) R[inv.eTask[c.src]][inv.eTask[c.dst]] = true;
    for (std::size_t k = 0; k < inv.nT; ++k)
        for (std::size_t i = 0; i < inv.nT; ++i)
            if (R[i][k])
                for (std::size_t j = 0; j < inv.nT; ++j)
                    if (R[k][j]) R[i][j] = true;
    for (std::size_t t = 0; t < inv.nT; ++t)
        if (R[t][t]) fail("The task call graph has a cycle, which FSPT does not support.");
    return inv;
}
 
inline double z(double v) { return std::isnan(v) ? 0.0 : v; }
 
inline Om solve(const LqnModel<double>& m, const LqnFocusOptions& opt, const Inv& inv) {
    const LqnStruct<double> l = lqn_finalize(m);
    const ln::LnSolution<double> s = opt.solver ? opt.solver(l) : ln::SolverLN<double>(l, ln::LnOptions()).get_ensemble_avg();
    Om om;
    for (std::size_t h = 0; h < inv.nH; ++h) {
        om.hQ.push_back(s.QN[1 + h]);
        om.hU.push_back(z(s.UN[1 + h]));
    }
    for (std::size_t t = 0; t < inv.nT; ++t) {
        const std::size_t i = l.tshift + 1 + t;
        om.tQ.push_back(z(s.QN[i]));
        om.tU.push_back(s.UN[i]);
        om.tX.push_back(z(s.TN[i]));
    }
    for (std::size_t e = 0; e < inv.nE; ++e) {
        const std::size_t i = l.eshift + 1 + e;
        om.eQ.push_back(s.QN[i]);
        om.eU.push_back(s.UN[i]);
        om.eR.push_back(s.RN[i]);
        om.eX.push_back(z(s.TN[i]));
    }
    for (std::size_t a = 0; a < inv.nA; ++a) om.aX.push_back(z(s.TN[l.ashift + 1 + a]));
    return om;
}
 
inline Sat saturation(const Inv& inv, const Om& om) {
    Sat s;
    s.SHT.assign(inv.nH, std::vector<double>(inv.nT, 0.0));
    s.SH.assign(inv.nH, 0.0);
    s.ST.assign(inv.nT, 0.0);
    for (std::size_t a = 0; a < inv.nA; ++a) {
        const std::size_t t = inv.aTask[a], h = inv.tHost[t];
        if (std::isfinite(inv.hMult[h])) s.SHT[h][t] += om.aX[a] * inv.aDem[a] / inv.hSpeed[h] / inv.hMult[h];
    }
    for (std::size_t h = 0; h < inv.nH; ++h)
        for (std::size_t t = 0; t < inv.nT; ++t) s.SH[h] += s.SHT[h][t];
    for (std::size_t t = 0; t < inv.nT; ++t)
        if (std::isfinite(inv.tMult[t]) && !inv.tIsRef[t]) s.ST[t] = om.tQ[t] / inv.tMult[t];
    return s;
}
 
inline void closure(const Inv& inv, std::vector<bool>& PT, std::vector<bool>& PH) {
    for (std::size_t t = 0; t < inv.nT; ++t)
        if (PT[t]) PH[inv.tHost[t]] = true;
    for (std::size_t t = 0; t < inv.nT; ++t)
        if (PH[inv.tHost[t]]) PT[t] = true;
}
 
inline std::string unique_name(std::set<std::string>& all, std::string nm) {
    while (all.count(nm)) nm += "_";
    all.insert(nm);
    return nm;
}
 
inline std::string map_name(int g, const std::string& name, const std::vector<LqnFocusGroup>& grp, bool task) {
    if (g == 0) return name;
    if (g < 0) return "";
    return task ? grp[g - 1].aggTask : grp[g - 1].aggActivity;
}
 
inline LqnFocusResult build(const LqnModel<double>& model, const Inv& inv, const Om& om, const Sat& sat,
                            std::vector<bool> PT, std::vector<bool> PH, const LqnFocusOptions& opt, bool validate) {
    const double tol = GlobalConstants::CoarseTol;
    const std::size_t nH = inv.nH, nT = inv.nT, nE = inv.nE, nA = inv.nA;
    std::vector<int> gid;
    std::vector<std::size_t> pidx;
    std::vector<std::vector<double>> shadow;
    std::size_t ng = 0;
    while (true) {
        closure(inv, PT, PH);
        pidx.clear();
        for (std::size_t t = 0; t < nT; ++t)
            if (PT[t]) pidx.push_back(t);
        gid.assign(nT, 0);
        std::vector<std::vector<bool>> keys;
        for (std::size_t t = 0; t < nT; ++t) {
            if (PT[t]) continue;
            std::vector<bool> key;
            bool any = false;
            for (std::size_t p : pidx) {
                key.push_back(inv.reach[p][t]);
                any = any || inv.reach[p][t];
            }
            if (!any) {
                gid[t] = -1;
                continue;
            }
            auto it = std::find(keys.begin(), keys.end(), key);
            if (it == keys.end()) {
                keys.push_back(key);
                gid[t] = static_cast<int>(keys.size());
            } else {
                gid[t] = static_cast<int>(it - keys.begin()) + 1;
            }
        }
        ng = keys.size();
        shadow.assign(ng + 1, std::vector<double>(nH, 1.0));
        bool changed = false;
        for (std::size_t h = 0; h < nH && !changed; ++h) {
            if (PH[h]) continue;
            std::set<int> gs;
            for (std::size_t t = 0; t < nT; ++t)
                if (inv.tHost[t] == h && gid[t] > 0) gs.insert(gid[t]);
            if (gs.size() < 2) continue;
            for (int g : gs) {
                double frac = 1.0;
                for (std::size_t t = 0; t < nT; ++t)
                    if (inv.tHost[t] == h && gid[t] > 0 && gid[t] != g) frac -= sat.SHT[h][t];
                if (frac <= tol) {
                    PH[h] = true;
                    changed = true;
                    break;
                }
                shadow[g][h] = frac;
            }
        }
        if (!changed) break;
    }
 
    std::set<std::string> all;
    for (const auto* v : {&inv.hName, &inv.tName, &inv.eName, &inv.aName}) all.insert(v->begin(), v->end());
    std::vector<LqnFocusGroup> grp(ng);
    for (std::size_t g = 1; g <= ng; ++g) {
        LqnFocusGroup& d = grp[g - 1];
        std::set<std::size_t> hg;
        std::size_t first = nT;
        for (std::size_t t = 0; t < nT; ++t) {
            if (gid[t] != static_cast<int>(g)) continue;
            if (first == nT) first = t;
            d.tasks.push_back(inv.tName[t]);
            d.mult += inv.tMult[t];
            hg.insert(inv.tHost[t]);
        }
        for (std::size_t p : pidx)
            if (inv.reach[p][first]) d.pi.push_back(inv.tName[p]);
        const std::string k = std::to_string(g);
        d.aggTask = unique_name(all, "AggT" + k);
        d.aggEntry = unique_name(all, "AggE" + k);
        d.aggActivity = unique_name(all, "AggA" + k);
        d.aggHost = unique_name(all, "AggP" + k);
        for (std::size_t h : hg) d.hostMult += inv.hMult[h];
        if (std::isfinite(d.hostMult)) {
            double s = 0;
            for (std::size_t h : hg) s += inv.hSpeed[h] * shadow[g][h] * inv.hMult[h];
            d.hostSpeed = s / d.hostMult;
        }
    }
    std::vector<std::string> fmEntry(nE);
    for (std::size_t e = 0; e < nE; ++e) {
        const int g = gid[inv.eTask[e]];
        fmEntry[e] = g == 0 ? inv.eName[e] : (g > 0 ? grp[g - 1].aggEntry : std::string());
    }
    for (const auto& c : inv.sync) {
        const int gs = gid[inv.aTask[c.src]], gd = gid[inv.eTask[c.dst]];
        if (gd > 0 && gs != gd && gs >= 0) grp[gd - 1].rate += om.aX[c.src] * c.mean;
    }
    for (std::size_t g = 1; g <= ng; ++g) {
        LqnFocusGroup& d = grp[g - 1];
        const double rateA = std::max(d.rate, GlobalConstants::Zero);
        double work = 0, time = 0;
        for (std::size_t a = 0; a < nA; ++a) {
            const std::size_t t = inv.aTask[a];
            if (gid[t] != static_cast<int>(g)) continue;
            const std::size_t h = inv.tHost[t];
            const double w = om.aX[a] * inv.aDem[a];
            work += w;
            time += w / (inv.hSpeed[h] * shadow[g][h]);
        }
        d.demand = std::isfinite(d.hostMult) ? work / rateA / d.hostSpeed : time / rateA;
        for (const auto& c : inv.sync) {
            if (gid[inv.aTask[c.src]] == static_cast<int>(g) && gid[inv.eTask[c.dst]] != static_cast<int>(g)) {
                const std::string& dn = fmEntry[c.dst];
                std::size_t k = index_of(d.callDests, dn);
                if (k == d.callDests.size()) {
                    d.callDests.push_back(dn);
                    d.callMeans.push_back(0.0);
                }
                d.callMeans[k] += om.aX[c.src] * c.mean / rateA;
            }
        }
    }
 
    // focused raw model: preserved elements copied verbatim, slots remapped
    LqnModel<double> fm;
    fm.name = model.name + "_focus";
    std::vector<std::size_t> hslot(nH, nH), tslot(nT, nT);
    for (std::size_t h = 0; h < nH; ++h) {
        if (!PH[h]) continue;
        hslot[h] = fm.procs.size();
        fm.procs.push_back(model.procs[h]);
    }
    for (std::size_t t = 0; t < nT; ++t) {
        if (!PT[t]) continue;
        tslot[t] = fm.tasks.size();
        auto tk = model.tasks[t];
        tk.proc_slot = hslot[inv.tHost[t]];
        fm.tasks.push_back(tk);
    }
    std::vector<std::size_t> aggTslot(ng);
    for (std::size_t g = 0; g < ng; ++g) {
        const LqnFocusGroup& d = grp[g];
        detail::RawProc p;
        p.name = d.aggHost;
        if (std::isfinite(d.hostMult)) {
            p.sched = SchedStrategy::PS;
            p.mult = d.hostMult;
        } else {
            p.sched = SchedStrategy::INF;
            p.mult = std::numeric_limits<double>::infinity();
        }
        fm.procs.push_back(p);
        detail::RawTask<double> tk;
        tk.name = d.aggTask;
        if (std::isfinite(d.mult)) {
            tk.sched = SchedStrategy::FCFS;
            tk.mult = d.mult;
        } else {
            tk.sched = SchedStrategy::INF;
            tk.mult = std::numeric_limits<double>::infinity();
        }
        tk.proc_slot = fm.procs.size() - 1;
        aggTslot[g] = fm.tasks.size();
        fm.tasks.push_back(tk);
    }
    for (std::size_t e = 0; e < nE; ++e) {
        if (!PT[inv.eTask[e]]) continue;
        auto en = model.entries[e];
        en.task_slot = tslot[inv.eTask[e]];
        fm.entries.push_back(en);
    }
    for (std::size_t g = 0; g < ng; ++g) {
        detail::RawEntry<double> en;
        en.name = grp[g].aggEntry;
        en.task_slot = aggTslot[g];
        en.reply_activities.push_back(grp[g].aggActivity);
        fm.entries.push_back(en);
    }
    for (std::size_t a = 0; a < nA; ++a) {
        if (!PT[inv.aTask[a]]) continue;
        auto ac = model.acts[a];
        ac.task_slot = tslot[inv.aTask[a]];
        ac.sync_calls.clear();
        for (const auto& c : model.acts[a].sync_calls) {
            const std::string& dn = fmEntry[index_of(inv.eName, c.dest)];
            auto it = std::find_if(ac.sync_calls.begin(), ac.sync_calls.end(),
                                   [&](const detail::RawCall<double>& x) { return x.dest == dn; });
            if (it == ac.sync_calls.end()) ac.sync_calls.push_back({dn, c.mean});
            else it->mean += c.mean;
        }
        fm.acts.push_back(ac);
    }
    for (std::size_t g = 0; g < ng; ++g) {
        const LqnFocusGroup& d = grp[g];
        detail::RawActivity<double> ac;
        ac.name = d.aggActivity;
        ac.hostdem = d.demand > GlobalConstants::FineTol ? Distrib<double>::exp_mean(d.demand)
                                                        : Distrib<double>::immediate();
        ac.bound_to_entry = d.aggEntry;
        ac.task_slot = aggTslot[g];
        for (std::size_t k = 0; k < d.callDests.size(); ++k) ac.sync_calls.push_back({d.callDests[k], d.callMeans[k]});
        fm.acts.push_back(ac);
    }
 
    LqnFocusResult res;
    std::size_t nPT = 0, nPH = 0;
    for (std::size_t t = 0; t < nT; ++t) {
        if (PT[t]) {
            res.preservedTasks.push_back(inv.tName[t]);
            ++nPT;
        }
        if (inv.forcedT[t]) res.forcedTasks.push_back(inv.tName[t]);
        if (gid[t] == -1) res.droppedTasks.push_back(inv.tName[t]);
    }
    for (std::size_t h = 0; h < nH; ++h)
        if (PH[h]) {
            res.preservedHosts.push_back(inv.hName[h]);
            ++nPH;
        }
    res.preservedMask = PT;
    res.preservedMask.insert(res.preservedMask.end(), PH.begin(), PH.end());
    res.taskNames = inv.tName;
    res.taskSat = sat.ST;
    res.hostNames = inv.hName;
    res.hostSat = sat.SH;
    for (std::size_t t = 0; t < nT; ++t) res.taskMap.push_back({inv.tName[t], map_name(gid[t], inv.tName[t], grp, true)});
    for (std::size_t e = 0; e < nE; ++e) res.entryMap.push_back({inv.eName[e], fmEntry[e]});
    for (std::size_t a = 0; a < nA; ++a)
        res.activityMap.push_back({inv.aName[a], map_name(gid[inv.aTask[a]], inv.aName[a], grp, false)});
    for (std::size_t h = 0; h < nH; ++h) {
        if (PH[h]) {
            res.hostMap.push_back({inv.hName[h], inv.hName[h]});
            continue;
        }
        std::set<int> gs;
        for (std::size_t t = 0; t < nT; ++t)
            if (inv.tHost[t] == h && gid[t] > 0) gs.insert(gid[t]);
        if (gs.empty()) res.hostMap.push_back({inv.hName[h], ""});
        for (int g : gs) res.hostMap.push_back({inv.hName[h], grp[g - 1].aggHost});
    }
    double soutMax = 0, sallMax = 0;
    bool anyOut = false;
    for (std::size_t t = 0; t < nT; ++t) {
        sallMax = std::max(sallMax, sat.ST[t]);
        if (!PT[t]) {
            soutMax = std::max(soutMax, sat.ST[t]);
            anyOut = true;
        }
    }
    for (std::size_t h = 0; h < nH; ++h) {
        sallMax = std::max(sallMax, sat.SH[h]);
        if (!PH[h]) {
            soutMax = std::max(soutMax, sat.SH[h]);
            anyOut = true;
        }
    }
    res.SR = (!anyOut || soutMax <= 0) ? std::numeric_limits<double>::infinity() : sallMax / soutMax;
    if (res.SR > 2) {
        res.trustedRange[0] = 1.0 / (res.SR - 2);
        res.trustedRange[1] = res.SR - 2;
    }
    res.FT = static_cast<double>(nPT) / nT;
    res.F = static_cast<double>(nPT + nPH + 2 * ng) / (nT + nH);
    std::vector<std::size_t> refs;
    for (std::size_t t = 0; t < nT; ++t)
        if (inv.tIsRef[t]) {
            refs.push_back(t);
            res.refTasks.push_back(inv.tName[t]);
        }
    res.errRef.assign(refs.size(), nan_());
    res.errRespRef.assign(refs.size(), nan_());
    if (validate) {
        const LqnStruct<double> l = lqn_finalize(fm);
        const ln::LnSolution<double> s = opt.solver ? opt.solver(l) : ln::SolverLN<double>(l, ln::LnOptions()).get_ensemble_avg();
        auto find = [&](const std::string& n, std::size_t lo, std::size_t hi) {
            for (std::size_t i = lo + 1; i <= hi; ++i)
                if (l.names[i] == n) return i;
            fail("element " + n + " not found in the focused model.");
        };
        LqnFocusRetained& ret = res.retained;
        auto add = [&](const std::string& n, const char* ty, double xo, double xf, double ro, double rf, double qo,
                       double qf, double uo, double uf) {
            ret.names.push_back(n);
            ret.types.push_back(ty);
            ret.X_om.push_back(xo);
            ret.X_fm.push_back(xf);
            ret.R_om.push_back(ro);
            ret.R_fm.push_back(rf);
            ret.Q_om.push_back(qo);
            ret.Q_fm.push_back(qf);
            ret.U_om.push_back(uo);
            ret.U_fm.push_back(uf);
        };
        std::map<std::size_t, std::size_t> taskRow;
        for (std::size_t t = 0; t < nT; ++t) {
            if (!PT[t]) continue;
            const std::size_t i = find(inv.tName[t], l.tshift, l.eshift);
            taskRow[t] = ret.names.size();
            add(inv.tName[t], "task", om.tX[t], s.TN[i], nan_(), nan_(), om.tQ[t], s.QN[i], om.tU[t], s.UN[i]);
        }
        for (std::size_t e = 0; e < nE; ++e) {
            if (!PT[inv.eTask[e]]) continue;
            const std::size_t i = find(inv.eName[e], l.eshift, l.ashift);
            add(inv.eName[e], "entry", om.eX[e], s.TN[i], om.eR[e], s.RN[i], om.eQ[e], s.QN[i], om.eU[e], s.UN[i]);
        }
        for (std::size_t h = 0; h < nH; ++h) {
            if (!PH[h]) continue;
            const std::size_t i = find(inv.hName[h], 0, l.nhosts);
            add(inv.hName[h], "host", nan_(), nan_(), nan_(), nan_(), om.hQ[h], s.QN[i], om.hU[h], s.UN[i]);
        }
        std::size_t kmax = refs.size(), kmaxR = refs.size();
        for (std::size_t k = 0; k < refs.size(); ++k) {
            const std::size_t t = refs[k], row = taskRow[t];
            const double xo = ret.X_om[row], xf = ret.X_fm[row];
            double zt = dmean(model.tasks[t].thinktime);
            if (zt <= GlobalConstants::FineTol) zt = 0;
            res.errRef[k] = (xf - xo) / xo;
            const double ro = inv.tMult[t] / xo - zt, rf = inv.tMult[t] / xf - zt;
            res.errRespRef[k] = (rf - ro) / ro;
            if (kmax == refs.size() || std::fabs(res.errRef[k]) > std::fabs(res.errRef[kmax])) kmax = k;
            if (kmaxR == refs.size() || std::fabs(res.errRespRef[k]) > std::fabs(res.errRespRef[kmaxR])) kmaxR = k;
        }
        if (kmax < refs.size()) {
            res.err = res.errRef[kmax];
            res.errResp = res.errRespRef[kmaxR];
        }
    }
    res.fm = std::move(fm);
    res.groups = std::move(grp);
    return res;
}
 
}  // namespace focus_detail
 
/**
 * Focused model of a layered network by FSPT aggregation.
 *
 * @param model raw original model (LqnBuilder::model() or read_lqnx_model).
 * @param opt   strategy, thresholds, user focus and solver.
 * @return the focused raw model (finalize it with lqn_finalize) and the
 *         diagnostics of the MATLAB `info` struct.
 */
inline LqnFocusResult lqn_focus(const LqnModel<double>& model, const LqnFocusOptions& opt = LqnFocusOptions()) {
    using namespace focus_detail;
    std::string strategy = opt.strategy;
    for (auto& c : strategy) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
    if (strategy != "sat" && strategy != "acc" && strategy != "manual")
        fail("unknown strategy \"" + opt.strategy + "\": use sat, acc or manual.");
    if (!(opt.r > 1)) fail("r must be greater than 1.");
    const Inv inv = inventory(model);
    const Om om = solve(model, opt, inv);
    const Sat sat = saturation(inv, om);
    std::vector<bool> baseT = inv.forcedT, baseH = inv.forcedH;
    for (const auto& nm : opt.preserve) {
        const std::size_t t = index_of(inv.tName, nm), h = index_of(inv.hName, nm);
        if (t == inv.nT && h == inv.nH)
            fail("Element \"" + nm + "\" in preserve is neither a task nor a host of the model.");
        if (t < inv.nT) baseT[t] = true;
        if (h < inv.nH) baseH[h] = true;
    }
    struct Cand {
        int kind;
        std::size_t i;
        double s;
    };
    std::vector<Cand> cand;
    for (std::size_t t = 0; t < inv.nT; ++t)
        if (sat.ST[t] > 0) cand.push_back({1, t, sat.ST[t]});
    for (std::size_t h = 0; h < inv.nH; ++h)
        if (sat.SH[h] > 0) cand.push_back({2, h, sat.SH[h]});
    std::stable_sort(cand.begin(), cand.end(), [](const Cand& a, const Cand& b) { return a.s > b.s; });
    const double Smax = cand.empty() ? 0.0 : cand[0].s;
    auto select = [&](std::size_t k, std::vector<bool>& T, std::vector<bool>& H) {
        T = baseT;
        H = baseH;
        for (std::size_t i = 0; i < k; ++i) (cand[i].kind == 1 ? T[cand[i].i] : H[cand[i].i]) = true;
    };
    LqnFocusResult res;
    std::vector<std::vector<double>> history;
    std::vector<bool> selT, selH;
    if (strategy == "sat") {
        std::size_t k = 0;
        while (k < cand.size() && Smax > 0 && cand[k].s > Smax / opt.r) ++k;
        select(k, selT, selH);
        res = build(model, inv, om, sat, selT, selH, opt, opt.validate);
    } else if (strategy == "manual") {
        select(std::min<std::size_t>(1, cand.size()), selT, selH);
        res = build(model, inv, om, sat, selT, selH, opt, opt.validate);
    } else {
        const std::size_t nc = cand.size();
        std::vector<bool> last;
        for (std::size_t i = 1; i <= std::max<std::size_t>(1, nc); ++i) {
            select(std::min(i, nc), selT, selH);
            closure(inv, selT, selH);
            std::vector<bool> key = selT;
            key.insert(key.end(), selH.begin(), selH.end());
            if (key == last && i < nc) continue;
            last = key;
            res = build(model, inv, om, sat, selT, selH, opt, true);
            std::size_t np = 0;
            for (bool b : res.preservedMask) np += b ? 1 : 0;
            history.push_back({static_cast<double>(i), static_cast<double>(np), res.F, res.err});
            if (std::fabs(res.err) < opt.x) break;
        }
    }
    res.strategy = strategy;
    res.r = opt.r;
    res.x = opt.x;
    res.accHistory = history;
    res.Smax = Smax;
    return res;
}
 
}  // namespace lqn
}  // namespace line
 
#endif  // LINE_API_LQN_LQN_FOCUS_H