483 return {
"srvn",
"srvn.ph",
"srvn.cs",
"flat",
"flat.cs",
"flat.ph",
"moment3",
"default"};
495 opt.layer_ssa.verbose =
false;
505 if (opt.method ==
"mwba.upper" || opt.method ==
"mwba.lower")
return box_bounds();
509 return is_ph_encoding() ? aggregate_ph() : aggregate();
525 if (lqn.nentries == 0)
return entrycdfrespt;
526 if (entrycdfrespt.size() <= lqn.nentries || entrycdfrespt[1].empty()) {
533 if (is_ph_encoding())
535 "getCdfRespT needs the routing encoding of the activity graph, which "
536 "method='srvn.ph' does not build. Rebuild the solver with "
537 "method='srvn.cs' or method='moment3'.");
544 const std::string saved = opt.method;
545 const std::string saved_resolved = lnmethod;
546 opt.method =
"moment3";
547 lnmethod =
"moment3";
551 lnmethod = saved_resolved;
553 return entrycdfrespt;
564 if (opt.ln_transient ==
"decoupled")
return tran_avg_decoupled();
565 if (opt.ln_transient ==
"coupled")
return tran_avg_coupled();
566 throw InputError(
"SolverLN: unknown ln_transient mode '" + opt.ln_transient +
567 "' (use 'coupled' or 'decoupled')");
583 if (results.empty()) iterate();
587 for (std::size_t e = 0; e < ensemble.size(); ++e) {
590 const bool exact_available =
591 opt.layer_solver ==
"mva" || opt.layer_solver ==
"nc";
593 ensemble[e], so, exact_available && !fj_tr[e].active(),
594 [
this, e]() {
return this->solve_layer(e); });
598 row.layer = ensemble[e].name;
605 out.
rows.push_back(row);
612 if (m.empty())
continue;
628 const bool upper = opt.method !=
"mwba.lower";
633 const std::size_t N = lqn.nidx;
634 auto blank = [&](std::vector<T>& v, std::vector<bool>& d) {
635 v.assign(N + 1, Tzero());
636 d.assign(N + 1,
false);
644 for (std::size_t i = 1; i <= N; ++i) {
645 if (b.defined_T[i]) {
646 s.
TN[i] = upper ? b.TN_up[i] : b.TN_lo[i];
649 if (b.defined_U[i]) {
650 s.
UN[i] = upper ? b.UN_up[i] : b.UN_lo[i];
660 std::size_t
nlayers()
const {
return ensemble.size(); }
661 const std::vector<qn::Layer<T>>&
layers()
const {
return ensemble; }
673 const std::size_t E = ensemble.size();
678 for (std::size_t e = 0; e < E; ++e) {
681 b.
msz[e] = ensemble[e].nstations;
682 b.
ksz[e] = ensemble[e].nclasses;
712 const std::size_t E = ensemble.size();
713 layer_tran_init.assign(E, std::vector<double>());
714 if (n.
empty())
return;
718 "SolverLN::init_from_marginal: the marginal is " + std::to_string(n.
rows()) +
"x" +
719 std::to_string(n.
cols()) +
" where the block-diagonal union of the layers is " +
720 std::to_string(b.
M) +
"x" + std::to_string(b.
K) +
721 "; the aggregate view is layerBlocks', not the LQN element count");
722 for (std::size_t e = 0; e < E; ++e) {
725 std::vector<double> y(lay.
nstates, 0.0);
727 for (std::size_t i = 0; i < L.
nstations && i < b.
msz[e]; ++i)
728 for (std::size_t r = 0; r < L.
nclasses && r < b.
ksz[e]; ++r) {
729 if (!lay.
enabled[i][r])
continue;
730 const double q = std::max(0.0, n(b.
roff[e] + i, b.
coff[e] + r));
731 y[lay.
qidx[i][r]] = q;
732 if (q > 0.0) any =
true;
738 if (any) layer_tran_init[e] = y;
765 std::vector<qn::Layer<T>> ensemble;
766 std::vector<long> idxhash;
767 std::vector<bool> ignore;
771 std::size_t idx, aidx, node, cls;
773 std::vector<UpdRow> servt_map, thinkt_map, call_map, actthinkt_map;
783 std::vector<UpdRow> arv_call_map;
786 std::size_t idx, tidx_caller, eidx, nodefrom, nodeto, cfrom, cto;
788 std::vector<RouteRow> route_map;
789 std::vector<std::size_t> unique_route_idx;
790 std::vector<std::size_t> route_reset, svc_reset;
792 Matrix<double> njobs;
795 std::vector<T> servt, residt, tput, util, thinkt;
796 std::vector<T> callservt, callresidt;
803 bool has_phase2 =
false;
804 std::vector<T> servt_ph1, servt_ph2, prOvertake;
806 std::set<std::size_t> single_replica_tasks;
807 std::vector<Distrib<T>> servtproc, thinkproc, thinktproc, tputproc, callservtproc;
813 std::vector<fluid::FluidPassage> servtcdf, callservtcdf;
814 std::vector<LnCdf> entrycdfrespt;
815 std::vector<mam::AphPair<T>> entryproc;
817 std::vector<std::vector<std::vector<fluid::FluidPassage>>> cdf_repo;
818 std::vector<double> servt_prev, residt_prev, tput_prev, thinkt_prev;
819 std::vector<double> callservt_prev, callresidt_prev;
820 std::vector<T> servt_prev_v, residt_prev_v, tput_prev_v, thinkt_prev_v, callservt_prev_v;
821 Matrix<T> servtmatrix;
824 Matrix<T> il_all, il_ph1;
825 std::vector<std::vector<std::size_t>> il_common_entries, il_src_all, il_src_ph2;
831 std::vector<std::vector<std::vector<double>>> layer_interlock;
832 std::vector<double> il_num_sources;
834 std::vector<std::vector<LayerResult<T>>> results;
835 std::vector<Matrix<T>> layer_init_sol;
840 std::vector<std::vector<double>> layer_tran_init;
852 std::vector<mva::FjMmt<T>> fj_tr;
853 std::vector<std::vector<T>> fj_lambda;
854 double relax_omega = 1.0;
856 std::shared_ptr<LnStochController<T>> stoch_ctl;
857 long averagingstart = -1;
858 bool hasconverged =
false;
860 bool moment_pass_done =
false;
861 std::vector<double> maxitererr;
862 int iterations_done = 0;
863 bool did_converge =
false;
865 std::size_t NT()
const {
return lqn.
tshift + lqn.
ntasks; }
867 static T Tzero() {
return num_traits<T>::from_int(0); }
868 static T Tone() {
return num_traits<T>::from_int(1); }
869 static double dbl(
const T& x) {
return num_traits<T>::to_double(x); }
881 reject_unsupported();
884 const std::size_t N = lqn.nidx;
885 ignore.assign(N + 1,
false);
888 std::vector<long> comp(N + 1, -1);
890 for (std::size_t v = 1; v <= N; ++v) {
891 if (comp[v] >= 0)
continue;
892 std::vector<std::size_t> stack{v};
894 while (!stack.empty()) {
895 const std::size_t u = stack.back();
897 for (std::size_t w : lqn.graph.succ(u))
902 for (std::size_t w : lqn.graph.pred(u))
911 std::vector<bool> has_ref(ncomp,
false);
912 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
913 const std::size_t tidx = lqn.tshift + t;
914 if (lqn.sched[tidx] == SchedStrategy::REF) has_ref[comp[tidx]] =
true;
916 for (std::size_t v = 1; v <= N; ++v)
917 if (!has_ref[comp[v]]) ignore[v] =
true;
924 for (std::size_t i = 1; i <= N; ++i) {
925 servtproc[i] = lqn.hostdem[i];
926 thinkproc[i] = lqn.think[i];
929 for (std::size_t c = 1; c <= lqn.ncalls; ++c)
930 callservtproc[c] = lqn.hostdem[lqn.callpair_dst[c]];
932 njobs = Matrix<double>(NT() + 1, NT() + 1, 0.0);
936 std::vector<bool> rr(NT() + 1,
false), sr(NT() + 1,
false);
937 for (
const RouteRow& r : route_map) rr[r.idx] =
true;
938 for (
const UpdRow& r : thinkt_map) sr[r.idx] =
true;
939 for (
const UpdRow& r : call_map) sr[r.idx] =
true;
942 for (
const UpdRow& r : arv_call_map) rr[r.idx] =
true;
943 for (std::size_t i = 1; i <= NT(); ++i) {
944 if (rr[i] && idxhash[i] >= 0) route_reset.push_back(std::size_t(idxhash[i]));
945 if (sr[i] && idxhash[i] >= 0) svc_reset.push_back(std::size_t(idxhash[i]));
947 std::sort(route_reset.begin(), route_reset.end());
948 route_reset.erase(std::unique(route_reset.begin(), route_reset.end()), route_reset.end());
949 std::sort(svc_reset.begin(), svc_reset.end());
950 svc_reset.erase(std::unique(svc_reset.begin(), svc_reset.end()), svc_reset.end());
958 void detect_phase2() {
960 for (std::size_t a = 1; a <= lqn.nacts; ++a)
961 if (lqn.actphase[a] > 1) has_phase2 =
true;
973 bool is_fwd_target(std::size_t eidx)
const {
974 for (std::size_t c = 1; c <= lqn.ncalls; ++c)
975 if (lqn.calltype[c] == CallType::FWD && lqn.callpair_dst[c] == eidx)
return true;
992 double open_arrival_rate_of(std::size_t tidx)
const {
993 if (lqn.isref[tidx])
return 0.0;
994 for (std::size_t e : lqn.entriesof[tidx])
995 if (lqn.issynccaller.any_col(e) || lqn.isasynccaller.any_col(e) || is_fwd_target(e))
998 for (std::size_t e : lqn.entriesof[tidx]) {
999 if (!lqn.has_arrival[e])
continue;
1000 const double m = dbl(lqn.arrival[e].mean);
1007 bool has_cache_node(std::size_t e)
const {
1008 for (std::size_t n = 0; n < ensemble[e].nodes.size(); ++n)
1009 if (ensemble[e].nodes[n].nodetype == NodeType::Cache)
return true;
1014 bool async_only_activity(std::size_t aidx)
const {
1015 if (aidx <= lqn.ashift || aidx > lqn.ashift + lqn.nacts)
return false;
1016 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
1017 const std::size_t eidx = lqn.eshift + e;
1018 if (lqn.graph.get(eidx, aidx) == Tzero())
continue;
1019 return lqn.isasynccaller.any_col(eidx) && !lqn.issynccaller.any_col(eidx);
1034 void reject_unsupported()
const {}
1039 void build_layers() {
1046 const std::string requested = ln_requested_method(opt.method);
1047 assert_call_groups(requested ==
"flat.cs");
1048 if (requested ==
"flat.cs") {
1049 lnmethod =
"flat.cs";
1053 if (requested ==
"flat.ph") {
1058 ph_laws_ready =
false;
1059 lnmethod =
"flat.ph";
1060 build_layers_ph(
true);
1063 if (requested ==
"srvn.ph" || requested ==
"srvn") {
1064 ph_laws_ready =
false;
1065 if (requested ==
"srvn.ph" || probe_srvn_ph()) {
1066 lnmethod =
"srvn.ph";
1071 lnmethod = (requested ==
"moment3") ?
"moment3" :
"srvn.cs";
1072 std::vector<qn::Layer<T>> raw(NT() + 1);
1073 std::vector<bool> present(NT() + 1,
false);
1075 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx) {
1076 if (ignore[hidx])
continue;
1077 build_layer(raw[hidx], {hidx}, lqn.tasksof[hidx],
true,
false);
1078 present[hidx] =
true;
1080 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
1081 const std::size_t tidx = lqn.tshift + t;
1082 if (ignore[tidx] || lqn.isref[tidx])
continue;
1083 bool any_caller = lqn.iscaller.any_row(tidx) || lqn.iscaller.any_col(tidx);
1084 if (!any_caller)
continue;
1086 std::vector<std::size_t> callers;
1087 for (std::size_t ct = 1; ct <= lqn.ntasks; ++ct) {
1088 const std::size_t c = lqn.tshift + ct;
1090 for (std::size_t e : lqn.entriesof[tidx])
1091 if (lqn.iscaller.get(c, e)) calls =
true;
1092 if (calls) callers.push_back(c);
1094 if (callers.empty())
continue;
1095 build_layer(raw[tidx], {tidx}, callers,
false,
false);
1096 present[tidx] =
true;
1099 idxhash.assign(lqn.nidx + 1, -1);
1101 for (std::size_t i = 1; i <= NT(); ++i)
1103 idxhash[i] = next++;
1104 ensemble.push_back(std::move(raw[i]));
1106 layer_init_sol.assign(ensemble.size(), Matrix<T>());
1133 void assert_call_groups(
bool flat)
const {
1134 if (lqn.callgroups.empty())
return;
1137 "Call groups routed by a routing strategy require the squashed layering; use "
1138 "method='flat'. Under srvn the targets never share a submodel, so the dispatch "
1139 "order cannot be represented.");
1140 if (opt.layer_solver !=
"ssa")
1142 "Routed call groups need a layer solver that resolves the strategy from the "
1143 "state; set layer_solver='ssa'. MVA, NC and FLD read the routing matrix, into "
1144 "which refresh_routing has expanded the strategy as a uniform split, and would "
1145 "return that split under a round-robin or JSQ label.");
1148 std::vector<std::size_t> flat_server_set()
const {
1149 std::vector<std::size_t> servers;
1150 for (std::size_t i = 1; i <= NT(); ++i) {
1151 if (lqn.repl[i] > 1.0)
1153 "Flat layering does not support replicated processors or tasks, use the "
1154 "default 'srvn' layering.");
1157 "Flat layering does not support cache tasks, use the default 'srvn' "
1159 if (lqn.hassetup[i])
1161 "Flat layering does not support setup tasks, use the default 'srvn' "
1164 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx)
1165 if (!ignore[hidx] && !lqn.tasksof[hidx].empty()) servers.push_back(hidx);
1166 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
1167 const std::size_t tidx = lqn.tshift + t;
1168 if (ignore[tidx] || lqn.isref[tidx])
continue;
1169 if (!lqn.iscaller.any_row(tidx) && !lqn.iscaller.any_col(tidx))
continue;
1170 bool has_task_caller =
false;
1171 for (std::size_t eidx : lqn.entriesof[tidx])
1172 for (std::size_t c = 1; c <= lqn.ntasks; ++c)
1173 if (lqn.iscaller.get(lqn.tshift + c, eidx)) has_task_caller =
true;
1174 if (has_task_caller) servers.push_back(tidx);
1176 if (servers.empty())
1178 "Flat layering found no server: the model has no processor with tasks.");
1190 void build_flat_layer() {
1191 const std::vector<std::size_t> flat_servers = flat_server_set();
1192 std::vector<std::size_t> flat_callers;
1193 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
1194 const std::size_t tidx = lqn.tshift + t;
1195 if (!ignore[tidx]) flat_callers.push_back(tidx);
1198 build_layer(layer, flat_servers, flat_callers,
false,
true);
1200 ensemble.push_back(std::move(layer));
1201 idxhash.assign(lqn.nidx + 1, -1);
1202 for (std::size_t s : flat_servers) idxhash[s] = 0;
1203 layer_init_sol.assign(ensemble.size(), Matrix<T>());
1211 void build_fork_views() {
1212 fj_tr.assign(ensemble.size(), mva::FjMmt<T>());
1213 fj_lambda.assign(ensemble.size(), {});
1214 for (std::size_t e = 0; e < ensemble.size(); ++e) {
1215 if (!ensemble[e].has_fork())
continue;
1222 if (ensemble[e].sourceIdx != 0)
1224 "SolverLN: layer '" + ensemble[e].name +
1225 "' carries both an AND fork and an open stream (an async call or an entry "
1226 "arrival); the fork-join transform needs a Source of its own");
1228 fj_lambda[e].assign(fj_tr[e].V.classes.size() + 1,
1249 const std::shared_ptr<std::vector<std::vector<std::size_t>>>& chaincols, std::size_t R) {
1250 return [f, cols, chaincols, R](
const std::vector<T>& n) -> std::vector<T> {
1251 const std::size_t L = n.size();
1252 const std::vector<std::vector<std::size_t>>& use =
1253 (L == R || chaincols->empty()) ? cols : *chaincols;
1254 std::vector<T> nop(use.size(), num_traits<T>::from_int(0));
1255 for (std::size_t j = 0; j < use.size(); ++j)
1256 for (std::size_t k = 0; k < use[j].size(); ++k)
1257 if (use[j][k] < L) nop[j] = T(nop[j] + n[use[j][k]]);
1258 const std::vector<T> w = f(nop);
1259 std::vector<T> v(L, num_traits<T>::from_int(1));
1260 if (w.empty())
return v;
1261 for (std::size_t j = 0; j < use.size(); ++j) {
1262 const T wj = w[j < w.size() ? j : w.size() - 1];
1263 for (std::size_t k = 0; k < use[j].size(); ++k)
1264 if (use[j][k] < L) v[use[j][k]] = wj;
1271 static std::vector<T> layer_peak(
const std::vector<T>& peakPerOperand,
1272 const std::vector<std::vector<std::size_t>>& cols,
1274 std::vector<T> peak(R, num_traits<T>::from_int(1));
1275 if (peakPerOperand.empty())
return peak;
1276 for (std::size_t j = 0; j < cols.size(); ++j) {
1277 const T pj = peakPerOperand[j < peakPerOperand.size() ? j : peakPerOperand.size() - 1];
1278 for (std::size_t k = 0; k < cols[j].size(); ++k)
1279 if (cols[j][k] < R) peak[cols[j][k]] = pj;
1293 void build_layer(qn::Layer<T>& m,
const std::vector<std::size_t>& idxSet,
1294 const std::vector<std::size_t>& callers,
bool ishostlayer,
bool flat) {
1295 const T one = Tone();
1299 const std::size_t idx = idxSet[0];
1300 m.name = flat ? lqn.hashnames[idx] +
".Flat" : lqn.hashnames[idx];
1303 const double rawrepl = lqn.repl[idx];
1304 std::size_t nreplicas = 1;
1305 if (!flat && rawrepl > 1.0 && !callers.empty()) {
1314 for (std::size_t c : callers)
1315 if (lqn.repl[c] != rawrepl) reduce =
false;
1317 for (std::size_t c : callers)
1318 if (lqn.fanout_at(c, idx) < rawrepl) reduce =
false;
1320 nreplicas = reduce ? 1 :
static_cast<std::size_t
>(std::llround(rawrepl));
1321 if (reduce && !ishostlayer) single_replica_tasks.insert(idx);
1323 const bool reduce_fanout = (nreplicas == 1 && rawrepl > 1.0 && !callers.empty());
1324 const std::vector<double>& mult = lqn.maxmult;
1328 m.clientIdx = m.add_station(qn::Station<T>{
"Clients", NodeType::Delay, SchedStrategy::INF,
1329 std::numeric_limits<double>::infinity(),
false, 0});
1330 const std::size_t clientNode = m.node_of_station(m.clientIdx);
1331 m.serverIdx = m.clientIdx + 1;
1335 std::vector<std::vector<std::size_t>> srv(lqn.nidx + 1);
1336 std::vector<std::vector<std::size_t>> srvnode(lqn.nidx + 1);
1337 m.server_idx_of.assign(lqn.nidx + 1, 0);
1338 for (std::size_t sidx : idxSet) {
1339 const bool sishost = sidx <= lqn.nhosts;
1340 srv[sidx].resize(nreplicas);
1341 for (std::size_t r = 0; r < nreplicas; ++r) {
1343 st.name = r == 0 ? lqn.hashnames[sidx]
1344 : lqn.hashnames[sidx] +
"." + std::to_string(r + 1);
1347 lqn.sched[sidx] == SchedStrategy::INF ? NodeType::Delay : NodeType::Queue;
1348 st.sched = lqn.sched[sidx];
1350 st.nservers = lqn.sched[sidx] == SchedStrategy::INF
1351 ? std::numeric_limits<double>::infinity()
1353 st.attr_ishost = flat ? sishost : ishostlayer;
1355 srv[sidx][r] = m.add_station(st);
1356 srvnode[sidx].push_back(m.node_of_station(srv[sidx][r]));
1358 m.server_idx_of[sidx] = srv[sidx][0];
1360 m.host_stations.push_back(srv[sidx][0]);
1362 m.task_stations.push_back(srv[sidx][0]);
1365 const std::vector<std::size_t>& server = srv[idx];
1366 const std::vector<std::size_t>& serverNode = srvnode[idx];
1369 auto servers_for = [&](std::size_t elem) ->
const std::vector<std::size_t>& {
1370 static const std::vector<std::size_t> none;
1371 if (elem >= 1 && elem <= lqn.nidx && !srv[elem].empty())
return srv[elem];
1374 auto server_nodes_for = [&](std::size_t elem) ->
const std::vector<std::size_t>& {
1375 static const std::vector<std::size_t> none;
1376 if (elem >= 1 && elem <= lqn.nidx && !srvnode[elem].empty())
return srvnode[elem];
1380 auto host_is_server = [&](std::size_t tidx_) {
1381 return !servers_for(lqn.parent[tidx_]).empty();
1384 auto is_layer_client = [&](std::size_t tidx_) {
1385 for (std::size_t s : idxSet)
1386 for (std::size_t e : lqn.entriesof[s])
1387 if (lqn.issynccaller.get(tidx_, e))
return true;
1391 auto set_all_servers = [&](std::size_t cl,
const Distrib<T>& d) {
1392 for (std::size_t s : idxSet)
1393 for (std::size_t st : srv[s]) m.set_service(st, cl, d);
1401 std::vector<std::size_t> group_of_call(lqn.ncalls + 1, 0);
1402 std::vector<std::vector<std::size_t>> group_members(lqn.callgroups.size() + 1);
1403 for (std::size_t g = 0; g < lqn.callgroups.size(); ++g) {
1404 const LqnCallGroup& grp = lqn.callgroups[g];
1405 for (std::size_t tgt : grp.
targets)
1406 for (std::size_t cidx : lqn.callsof[grp.
caller])
1407 if (lqn.callpair_dst[cidx] == tgt && lqn.calltype[cidx] == CallType::SYNC &&
1408 group_of_call[cidx] == 0) {
1409 group_of_call[cidx] = g + 1;
1410 group_members[g + 1].push_back(cidx);
1415 std::vector<std::size_t> grp_router(lqn.callgroups.size() + 1, 0);
1416 std::vector<std::size_t> grp_dispatch(lqn.callgroups.size() + 1, 0);
1417 std::vector<std::size_t> grp_class(lqn.callgroups.size() + 1, 0);
1420 std::vector<std::array<std::size_t, 3>> routed_group_sites;
1423 std::vector<std::size_t> acts_in_caller;
1424 for (std::size_t c : callers)
1425 for (std::size_t a : lqn.actsof[c]) acts_in_caller.push_back(a);
1426 bool hasfork =
false, hasjoin =
false;
1427 std::size_t maxfanout = 1;
1428 for (std::size_t a : acts_in_caller) {
1429 if (lqn.actposttype[a] == PrecedenceType::POST_AND) hasfork =
true;
1430 if (lqn.actpretype[a] == PrecedenceType::PRE_AND) hasjoin =
true;
1431 std::size_t nand = 0;
1432 for (std::size_t sx : lqn.graph.succ(a))
1433 if (lqn.actposttype[sx] == PrecedenceType::POST_AND) ++nand;
1434 if (nand > maxfanout) maxfanout = nand;
1436 std::size_t forkNode = 0, joinNode = 0, joinStation = 0;
1437 std::vector<std::size_t> forkRouter;
1439 forkNode = m.add_node(
"Fork_PostAnd", NodeType::Fork,
false);
1440 for (std::size_t f = 1; f <= maxfanout; ++f)
1441 forkRouter.push_back(
1442 m.add_node(
"Fork_PostAnd_" + std::to_string(f), NodeType::Router,
true));
1446 js.name =
"Join_PreAnd";
1447 js.nodetype = NodeType::Join;
1448 js.sched = SchedStrategy::INF;
1449 js.nservers = std::numeric_limits<double>::infinity();
1450 js.attr_ishost =
false;
1452 joinStation = m.add_station(js);
1453 joinNode = m.node_of_station(joinStation);
1454 if (forkNode) m.fj.emplace_back(forkNode, joinNode);
1461 bool iscachelayer = !flat && ishostlayer && !callers.empty();
1462 for (std::size_t c : callers)
1463 if (!lqn.iscache[c]) iscachelayer =
false;
1472 std::size_t cacheNode = 0;
1473 qn::CacheParam<T> cachepar;
1475 const std::size_t ct = callers[0];
1476 cachepar.nitems = lqn.nitems[ct];
1477 cachepar.itemcap = lqn.itemcap[ct];
1478 cachepar.replacestrat = lqn.replacestrat[ct];
1479 cacheNode = m.add_node(lqn.hashnames[ct], NodeType::Cache,
true);
1489 std::vector<std::size_t> async_here;
1490 for (std::size_t c : callers)
1491 for (std::size_t aidx : lqn.actsof[c])
1492 for (std::size_t cidx : lqn.callsof[aidx])
1493 if (lqn.calltype[cidx] == CallType::ASYNC &&
1494 !servers_for(lqn.parent[lqn.callpair_dst[cidx]]).empty())
1495 async_here.push_back(cidx);
1496 std::vector<std::size_t> open_entries;
1497 for (std::size_t c : callers) {
1505 for (std::size_t eidx : lqn.entriesof[c])
1506 if (lqn.has_arrival[eidx]) open_entries.push_back(eidx);
1509 std::size_t sourceStation = 0, sourceNode = 0, sinkNode = 0;
1510 if (!async_here.empty() || !open_entries.empty()) {
1512 src.name =
"Source";
1513 src.nodetype = NodeType::Source;
1515 src.sched = SchedStrategy::EXT;
1517 src.attr_ishost =
false;
1519 sourceStation = m.add_station(src);
1520 sourceNode = m.node_of_station(sourceStation);
1521 sinkNode = m.add_node(
"Sink", NodeType::Sink,
false);
1522 m.sourceIdx = sourceStation;
1523 m.sinkNode = sinkNode;
1527 std::vector<std::size_t> forkClassStack;
1530 std::vector<std::size_t> cls(lqn.nidx + 1, 0);
1531 std::vector<std::size_t> callcls(lqn.ncalls + 1, 0);
1534 std::vector<std::size_t> auxcallcls(lqn.ncalls + 1, 0);
1542 auto caller_needs_class = [&](std::size_t tidx_caller) {
1543 if (host_is_server(tidx_caller)) {
1544 if (lqn.isref[tidx_caller])
return true;
1545 for (std::size_t e : lqn.entriesof[tidx_caller])
1546 if (lqn.issynccaller.any_col(e) || lqn.isasynccaller.any_col(e) ||
1547 is_fwd_target(e) || lqn.has_arrival[e])
1550 return is_layer_client(tidx_caller);
1552 auto caller_acts_visible = [&](std::size_t tidx_caller) {
1553 return host_is_server(tidx_caller) || is_layer_client(tidx_caller);
1557 for (std::size_t tidx_caller : callers) {
1558 if (caller_needs_class(tidx_caller)) {
1559 double nj = njobs(tidx_caller, idx);
1562 const bool caller_single_replica =
1563 reduce_fanout || single_replica_tasks.count(tidx_caller) > 0;
1564 nj = caller_single_replica ? mult[tidx_caller]
1565 : mult[tidx_caller] * lqn.repl[tidx_caller];
1566 if (std::isinf(nj)) {
1568 for (std::size_t c = 1; c <= NT(); ++c)
1569 if (lqn.taskgraph.get(c, tidx_caller) != Tzero()) s += mult[c];
1571 if (std::isinf(nj)) {
1573 for (std::size_t c = 1; c <= NT(); ++c)
1574 if (std::isfinite(mult[c])) s2 += mult[c] * lqn.repl[c];
1575 nj = std::min(s2, 1000.0);
1578 njobs(tidx_caller, idx) = nj;
1581 jc.name = lqn.hashnames[tidx_caller];
1582 jc.type = JobClassType::CLOSED;
1584 jc.refstat = m.clientIdx;
1585 jc.completes =
false;
1586 jc.is_ref_class =
true;
1587 jc.attr_kind = int(LqnElement::TASK);
1588 jc.attr_idx = tidx_caller;
1589 cls[tidx_caller] = m.add_class(jc);
1590 m.attr_tasks.emplace_back(cls[tidx_caller], tidx_caller);
1591 if (lqn.isref[tidx_caller]) {
1592 m.set_service(m.clientIdx, cls[tidx_caller], thinkproc[tidx_caller]);
1598 thinkt_map.push_back({idx, tidx_caller, m.clientIdx, cls[tidx_caller]});
1601 for (std::size_t eidx : lqn.entriesof[tidx_caller]) {
1603 ec.name = lqn.hashnames[eidx];
1604 ec.type = JobClassType::CLOSED;
1605 ec.population = 0.0;
1606 ec.refstat = m.clientIdx;
1607 ec.completes =
false;
1608 ec.attr_kind = int(LqnElement::ENTRY);
1610 cls[eidx] = m.add_class(ec);
1611 m.attr_entries.emplace_back(cls[eidx], eidx);
1616 for (std::size_t aidx : lqn.actsof[tidx_caller]) {
1617 if (caller_acts_visible(tidx_caller)) {
1619 ac.name = lqn.hashnames[aidx];
1620 ac.type = JobClassType::CLOSED;
1621 ac.population = 0.0;
1622 ac.refstat = m.clientIdx;
1623 ac.completes =
false;
1624 ac.attr_kind = int(LqnElement::ACTIVITY);
1626 cls[aidx] = m.add_class(ac);
1627 m.attr_activities.emplace_back(cls[aidx], aidx);
1631 const std::size_t hidx = lqn.parent[lqn.parent[aidx]];
1632 if (servers_for(hidx).empty())
1633 m.set_service(m.clientIdx, cls[aidx], servtproc[aidx]);
1635 for (std::size_t cidx : lqn.callsof[aidx]) {
1636 if (lqn.calltype[cidx] == CallType::ASYNC) {
1640 const std::size_t adst = lqn.parent[lqn.callpair_dst[cidx]];
1641 if (servers_for(adst).empty())
continue;
1643 oc.name = lqn.callhashnames[cidx];
1644 oc.type = JobClassType::OPEN;
1645 oc.population = std::numeric_limits<double>::infinity();
1646 oc.refstat = sourceStation;
1647 oc.completes =
false;
1648 oc.is_ref_class =
false;
1649 oc.attr_kind = int(LqnElement::CALL);
1651 callcls[cidx] = m.add_class(oc);
1652 m.attr_calls.push_back({callcls[cidx], cidx, lqn.callpair_src[cidx],
1653 lqn.callpair_dst[cidx]});
1661 m.set_service(sourceStation, callcls[cidx],
1664 T minRespTA = Tzero();
1665 for (std::size_t ta : lqn.actsof[adst]) minRespTA += lqn.hostdem[ta].mean;
1666 for (std::size_t st : servers_for(adst)) {
1668 call_map.push_back({idx, cidx, st, callcls[cidx]});
1670 arv_call_map.push_back({idx, cidx, sourceStation, callcls[cidx]});
1673 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
1674 const std::size_t gid = group_of_call[cidx];
1683 if (grp_dispatch[gid] == 0) {
1689 const std::string tag =
1690 lqn.hashnames[aidx] +
".Dispatch" + std::to_string(gid);
1691 grp_router[gid] = m.add_node(tag +
".Router", NodeType::Router,
true);
1694 dc.type = JobClassType::CLOSED;
1695 dc.population = 0.0;
1696 dc.refstat = m.clientIdx;
1697 dc.completes =
false;
1698 dc.attr_kind = int(LqnElement::CALL);
1700 grp_dispatch[gid] = m.add_class(dc);
1703 gc.name = lqn.callhashnames[cidx] +
".Group" + std::to_string(gid);
1704 gc.type = JobClassType::CLOSED;
1705 gc.population = 0.0;
1706 gc.refstat = m.clientIdx;
1707 gc.completes =
false;
1708 gc.attr_kind = int(LqnElement::CALL);
1710 grp_class[gid] = m.add_class(gc);
1712 routed_group_sites.push_back(
1713 {grp_router[gid], grp_dispatch[gid], gid});
1715 callcls[cidx] = grp_dispatch[gid];
1716 m.attr_calls.push_back({callcls[cidx], cidx, lqn.callpair_src[cidx],
1717 lqn.callpair_dst[cidx]});
1718 for (std::size_t st2 : servers_for(lqn.parent[lqn.callpair_dst[cidx]]))
1719 m.set_service(st2, callcls[cidx], callservtproc[cidx]);
1723 cc.name = lqn.callhashnames[cidx];
1724 cc.type = JobClassType::CLOSED;
1725 cc.population = 0.0;
1726 cc.refstat = m.clientIdx;
1727 cc.completes =
false;
1728 cc.attr_kind = int(LqnElement::CALL);
1730 callcls[cidx] = m.add_class(cc);
1731 m.attr_calls.push_back({callcls[cidx], cidx, lqn.callpair_src[cidx],
1732 lqn.callpair_dst[cidx]});
1740 const std::size_t seedidx = flat ? lqn.parent[lqn.callpair_dst[cidx]] : idx;
1741 T minRespT = Tzero();
1742 for (std::size_t ta : lqn.actsof[seedidx]) minRespT += lqn.hostdem[ta].mean;
1743 for (std::size_t st : servers_for(seedidx))
1746 if (lqn.callproc_mean[cidx] !=
1747 num_traits<T>::from_double(
double(nreplicas))) {
1749 xc.name = lqn.callhashnames[cidx] +
".Aux";
1750 xc.type = JobClassType::CLOSED;
1751 xc.population = 0.0;
1752 xc.refstat = m.clientIdx;
1753 xc.completes =
false;
1754 xc.attr_kind = int(LqnElement::CALL);
1756 auxcallcls[cidx] = m.add_class(xc);
1765 std::size_t curclass;
1771 std::vector<std::size_t> curnodes;
1773 const int atClient = 1, atServer = 2, atCache = 3;
1774 std::vector<int> jobposkey(lqn.nidx + 1, atClient);
1775 std::vector<std::size_t> curclasskey(lqn.nidx + 1, 0);
1776 std::vector<std::vector<std::size_t>> curnodeskey(lqn.nidx + 1);
1778 std::function<Ctx(std::size_t, std::size_t, Ctx)> recur =
1779 [&](std::size_t tidx_caller, std::size_t aidx, Ctx st) -> Ctx {
1780 jobposkey[aidx] = st.jobpos;
1781 curclasskey[aidx] = st.curclass;
1782 curnodeskey[aidx] = st.curnodes;
1783 const std::vector<std::size_t> nexts = lqn.graph.succ(aidx);
1784 std::size_t lastEntryClass = st.curclass;
1786 bool next_is_fork =
false;
1787 for (std::size_t sx : nexts)
1788 if (lqn.actposttype[sx] == PrecedenceType::POST_AND) next_is_fork =
true;
1789 const Ctx preFork = st;
1790 std::vector<std::size_t> andSuccs;
1791 for (std::size_t sx : nexts)
1792 if (lqn.actposttype[sx] == PrecedenceType::POST_AND) andSuccs.push_back(sx);
1794 for (std::size_t k = 0; k < nexts.size(); ++k) {
1795 const std::size_t nextaidx = nexts[k];
1796 if (next_is_fork) st = preFork;
1797 bool isLoop = lqn.graph.get(aidx, nextaidx) != lqn.dag.get(aidx, nextaidx);
1798 if (lqn.parent[aidx] != lqn.parent[nextaidx]) {
1800 std::size_t cidx = 0;
1801 for (std::size_t c : lqn.callsof[aidx])
1802 if (lqn.callpair_dst[c] == nextaidx) cidx = c;
1803 if (cidx == 0)
continue;
1805 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
1806 const std::size_t gid = group_of_call[cidx];
1810 if (group_members[gid].empty() || group_members[gid][0] != cidx)
continue;
1817 std::vector<std::size_t> tnode2, tstat2, tcall2;
1818 for (std::size_t mc : group_members[gid]) {
1819 const std::size_t tt = lqn.parent[lqn.callpair_dst[mc]];
1820 if (servers_for(tt).empty())
continue;
1821 tnode2.push_back(server_nodes_for(tt)[0]);
1822 tstat2.push_back(servers_for(tt)[0]);
1823 tcall2.push_back(mc);
1825 if (tnode2.size() < 2)
continue;
1826 const std::size_t fromNode2 =
1827 st.jobpos == atClient ? clientNode : st.curnodes[0];
1828 const std::size_t dc = grp_dispatch[gid], gc = grp_class[gid];
1829 m.set_route(st.curclass, dc, fromNode2, grp_router[gid], Tone());
1831 T(Tone() / num_traits<T>::from_int(
int(tnode2.size())));
1832 for (std::size_t d = 0; d < tnode2.size(); ++d) {
1833 m.set_route(dc, dc, grp_router[gid], tnode2[d], share2);
1834 m.set_route(dc, gc, tnode2[d], clientNode, Tone());
1835 m.set_service(tstat2[d], dc, callservtproc[tcall2[d]]);
1836 call_map.push_back({idx, tcall2[d], tstat2[d], dc});
1839 st.jobpos = atClient;
1840 st.curnodes.clear();
1843 const std::size_t ctgt = lqn.parent[lqn.callpair_dst[cidx]];
1844 st = route_sync_call(m, idx, cidx, st, server_nodes_for(ctgt),
1845 servers_for(ctgt), callcls, auxcallcls, clientNode,
1846 atClient, atServer, flat);
1850 bool any_entry_succ =
false;
1851 for (std::size_t sx : nexts)
1852 if (sx > lqn.eshift && sx <= lqn.eshift + lqn.nentries) any_entry_succ =
true;
1853 if (!any_entry_succ) {
1854 st.jobpos = jobposkey[aidx];
1855 st.curclass = curclasskey[aidx];
1856 st.curnodes = curnodeskey[aidx];
1858 if (k > 0 && nexts[k - 1] > lqn.eshift && nexts[k - 1] <= lqn.eshift + lqn.nentries)
1859 lastEntryClass = st.curclass;
1860 st.jobpos = atClient;
1861 st.curclass = lastEntryClass;
1862 st.curnodes.clear();
1864 const T w = lqn.graph.get(aidx, nextaidx);
1872 if (iscachelayer && lqn.nitems[aidx] > 0 && cls[nextaidx] != 0) {
1873 const std::size_t readcls = cls[nextaidx];
1874 m.set_route(st.curclass, readcls, clientNode, cacheNode, w);
1875 if (cachepar.pread.size() < m.classes.size())
1876 cachepar.pread.resize(m.classes.size());
1877 cachepar.pread[readcls - 1] = lqn.itemproc[aidx];
1878 const std::vector<std::size_t> hm = lqn.graph.succ(nextaidx);
1880 throw InputError(
"SolverLN: the cache read '" + lqn.names[nextaidx] +
1881 "' needs exactly one hit and one miss successor");
1882 if (cachepar.hitclass.size() < m.classes.size()) {
1883 cachepar.hitclass.resize(m.classes.size(), 0);
1884 cachepar.missclass.resize(m.classes.size(), 0);
1886 cachepar.hitclass[readcls - 1] = cls[hm[0]];
1887 cachepar.missclass[readcls - 1] = cls[hm[1]];
1888 st.jobpos = atCache;
1889 st.curclass = readcls;
1890 st.curnodes.clear();
1891 st = recur(tidx_caller, nextaidx, st);
1895 const bool is_and_join_tail = lqn.actpretype[aidx] == PrecedenceType::PRE_AND;
1897 std::size_t fbranch = 0;
1899 for (std::size_t q = 0; q < andSuccs.size(); ++q)
1900 if (andSuccs[q] == nextaidx) fbranch = q + 1;
1907 const std::size_t hidxOf = lqn.parent[lqn.parent[nextaidx]];
1908 const std::vector<std::size_t>& actStations = servers_for(hidxOf);
1909 const std::vector<std::size_t>& actNodes = server_nodes_for(hidxOf);
1910 const bool actAtServer = !actStations.empty();
1911 const std::size_t from = st.jobpos == atClient ? clientNode
1912 : st.jobpos == atCache ? cacheNode
1915 for (std::size_t r = 0; r < nreplicas; ++r) {
1916 const std::size_t fromNode =
1917 st.jobpos == atClient ? clientNode
1918 : st.jobpos == atCache ? cacheNode
1919 : st.curnodes[std::min(r, st.curnodes.size() - 1)];
1920 const std::size_t toNode = actAtServer ? actNodes[r] : clientNode;
1921 if (next_is_fork && fbranch > 0) {
1922 m.set_route(st.curclass, st.curclass, fromNode, forkNode, Tone());
1923 if (r == 0) forkClassStack.push_back(st.curclass);
1924 m.set_route(st.curclass, st.curclass, forkNode, forkRouter[fbranch - 1],
1926 m.set_route(st.curclass, cls[nextaidx], forkRouter[fbranch - 1], toNode,
1928 }
else if (is_and_join_tail) {
1931 if (forkClassStack.empty())
1932 throw InputError(
"SolverLN: an AND join has no matching fork in '" +
1933 lqn.names[aidx] +
"'");
1934 const std::size_t forkClass = forkClassStack.back();
1935 if (r + 1 == nreplicas) forkClassStack.pop_back();
1936 m.set_route(st.curclass, forkClass, fromNode, joinNode, Tone());
1937 m.set_route(forkClass, cls[nextaidx], joinNode, toNode, Tone());
1939 m.set_route(st.curclass, cls[nextaidx], fromNode, toNode, w);
1942 m.set_service(actStations[r], cls[nextaidx], lqn.hostdem[nextaidx]);
1946 st.jobpos = atServer;
1947 st.curclass = cls[nextaidx];
1948 st.curnodes = actNodes;
1949 servt_map.push_back({idx, nextaidx, actStations[0], cls[nextaidx]});
1951 st.jobpos = atClient;
1952 st.curclass = cls[nextaidx];
1953 st.curnodes.clear();
1954 m.set_service(m.clientIdx, cls[nextaidx], servtproc[nextaidx]);
1955 thinkt_map.push_back({idx, nextaidx, m.clientIdx, cls[nextaidx]});
1957 if (aidx != nextaidx && !isLoop) {
1958 st = recur(tidx_caller, nextaidx, st);
1960 if (st.jobpos == atClient) {
1961 m.set_route(st.curclass, cls[tidx_caller], clientNode, clientNode, Tone());
1963 for (std::size_t nd : st.curnodes)
1964 m.set_route(st.curclass, cls[tidx_caller], nd, clientNode, Tone());
1967 if (!is_aux_class(m.classes[st.curclass - 1].name))
1968 m.classes[st.curclass - 1].completes =
true;
1974 for (std::size_t tidx_caller : callers) {
1975 if (!caller_needs_class(tidx_caller))
continue;
1976 const std::vector<std::size_t>& ents = lqn.entriesof[tidx_caller];
1977 const T share = T(Tone() / num_traits<T>::from_int(
int(ents.size())));
1978 for (std::size_t eidx : ents) {
1979 m.set_route(cls[tidx_caller], cls[eidx], clientNode, clientNode, share);
1980 if (ents.size() > 1)
1981 route_map.push_back({idx, tidx_caller, eidx, m.clientIdx, m.clientIdx,
1982 cls[tidx_caller], cls[eidx]});
1983 Ctx st{cls[eidx], atClient, {}};
1984 recur(tidx_caller, eidx, st);
1994 if (sourceStation != 0) {
1995 const T nrep = num_traits<T>::from_int(
int(nreplicas));
1996 for (std::size_t cidx : async_here) {
1997 const std::size_t oc = callcls[cidx];
1998 if (oc == 0)
continue;
2001 const std::vector<std::size_t>& anode =
2002 server_nodes_for(lqn.parent[lqn.callpair_dst[cidx]]);
2003 const T callmean = lqn.callproc_mean[cidx];
2004 if (callmean < Tone()) {
2006 m.set_route(oc, oc, sourceNode, sinkNode, T(Tone() - callmean));
2007 for (std::size_t r = 0; r < anode.size(); ++r) {
2008 m.set_route(oc, oc, sourceNode, anode[r], T(callmean / nrep));
2009 m.set_route(oc, oc, anode[r], sinkNode, Tone());
2013 const T p = T(Tone() / callmean);
2014 for (std::size_t r = 0; r < anode.size(); ++r) {
2015 m.set_route(oc, oc, sourceNode, anode[r], T(Tone() / nrep));
2016 for (std::size_t q = 0; q < anode.size(); ++q)
2017 m.set_route(oc, oc, anode[r], anode[q], T((Tone() - p) / nrep));
2018 m.set_route(oc, oc, anode[r], sinkNode, p);
2022 for (std::size_t eidx : open_entries) {
2024 eo.name = lqn.hashnames[eidx] +
"_Open";
2025 eo.type = JobClassType::OPEN;
2026 eo.population = std::numeric_limits<double>::infinity();
2027 eo.refstat = sourceStation;
2028 eo.completes =
false;
2029 eo.is_ref_class =
false;
2030 eo.attr_kind = int(LqnElement::ENTRY);
2032 const std::size_t ec = m.add_class(eo);
2034 m.set_service(sourceStation, ec, lqn.arrival[eidx]);
2036 std::size_t bound = 0;
2037 for (std::size_t sx : lqn.graph.succ(eidx))
2038 if (bound == 0 && sx > lqn.ashift) bound = sx;
2039 const Distrib<T>& svc = bound != 0 ? servtproc[bound] : servtproc[eidx];
2042 const std::vector<std::size_t>& ostat =
2043 flat ? servers_for(lqn.parent[eidx]) : server;
2044 const std::vector<std::size_t>& onode =
2045 flat ? server_nodes_for(lqn.parent[eidx]) : serverNode;
2046 for (std::size_t r = 0; r < ostat.size(); ++r) {
2047 m.set_service(ostat[r], ec, svc);
2048 m.set_route(ec, ec, sourceNode, onode[r], T(Tone() / nrep));
2049 m.set_route(ec, ec, onode[r], sinkNode, Tone());
2061 if (idx < lqn.lincon_A.size() && lqn.lincon_A[idx].rows() > 0) {
2062 const Matrix<T>& Aelem = lqn.lincon_A[idx];
2063 Matrix<T> Alayer(Aelem.rows(), m.classes.size(), Tzero());
2064 const std::vector<std::size_t>& constrained =
2065 ishostlayer ? lqn.tasksof[idx] : lqn.entriesof[idx];
2067 for (std::size_t j = 0; j < constrained.size() && j < Aelem.cols(); ++j) {
2068 std::vector<std::size_t> layerClasses;
2070 for (std::size_t a : lqn.actsof[constrained[j]])
2071 if (cls[a] != 0) layerClasses.push_back(cls[a]);
2073 for (std::size_t c = 1; c <= lqn.ncalls; ++c)
2074 if (lqn.callpair_dst[c] == constrained[j] && callcls[c] != 0)
2075 layerClasses.push_back(callcls[c]);
2077 for (std::size_t k = 0; k < layerClasses.size(); ++k)
2078 for (std::size_t rr = 0; rr < Aelem.rows(); ++rr) {
2079 Alayer(rr, layerClasses[k] - 1) =
2080 T(Alayer(rr, layerClasses[k] - 1) + Aelem(rr, j));
2081 if (Aelem(rr, j) != Tzero()) any =
true;
2088 typename qn::NetworkStruct<T>::Region rg;
2089 const std::size_t M = m.stations.size(), K = m.classes.size();
2090 rg.cap.assign(M, std::vector<double>(K + 1, -1.0));
2091 rg.maxmem.assign(M, -1.0);
2092 rg.members.assign(M,
false);
2094 rg.weight.assign(K, Tone());
2095 rg.size.assign(K, Tone());
2096 for (std::size_t r = 0; r < nreplicas; ++r) rg.members[server[r] - 1] =
true;
2097 rg.lincon_A = Alayer;
2098 rg.lincon_b = lqn.lincon_b[idx];
2099 m.regions.push_back(rg);
2106 if (cacheNode != 0) {
2107 cachepar.pread.resize(m.classes.size());
2108 cachepar.hitclass.resize(m.classes.size(), 0);
2109 cachepar.missclass.resize(m.classes.size(), 0);
2110 m.nodeparam[cacheNode] = cachepar;
2118 for (
const std::array<std::size_t, 3>& site : routed_group_sites) {
2119 qn::NodeDef& nd = m.nodes[site[0] - 1];
2120 if (nd.routing.size() < m.classes.size())
2121 nd.routing.resize(m.classes.size(), RoutingStrategy::PROB);
2122 nd.routing[site[1] - 1] = lqn.callgroups[site[2] - 1].strategy;
2134 std::vector<std::pair<std::shared_ptr<std::vector<std::vector<std::size_t>>>,
2135 std::vector<std::vector<std::size_t>>>> deferred_chaincols;
2136 for (std::size_t sidx : idxSet) {
2137 const bool hasld = sidx < lqn.lldscaling.size() && !lqn.lldscaling[sidx].empty();
2138 const bool hascd = sidx < lqn.cdscaling.size() && bool(lqn.cdscaling[sidx]);
2139 const bool hasjd = sidx < lqn.jdscaling.size() && bool(lqn.jdscaling[sidx]);
2140 const bool haspools = sidx < lqn.pools.size() && !lqn.pools[sidx].empty();
2141 if (!(hasld || hascd || hasjd || haspools))
continue;
2142 const bool sishost = sidx <= lqn.nhosts;
2143 const std::vector<std::size_t>& operandIdx =
2144 sishost ? lqn.tasksof[sidx] : lqn.entriesof[sidx];
2145 std::vector<std::vector<std::size_t>> cols(operandIdx.size());
2146 for (std::size_t j = 0; j < operandIdx.size(); ++j) {
2148 for (std::size_t a : lqn.actsof[operandIdx[j]])
2149 if (cls[a] != 0) cols[j].push_back(cls[a] - 1);
2151 for (std::size_t c = 1; c <= lqn.ncalls; ++c)
2152 if (lqn.callpair_dst[c] == operandIdx[j] && callcls[c] != 0)
2153 cols[j].push_back(callcls[c] - 1);
2156 const std::size_t R = m.classes.size();
2157 std::shared_ptr<std::vector<std::vector<std::size_t>>> chaincols =
2158 std::make_shared<std::vector<std::vector<std::size_t>>>();
2159 deferred_chaincols.push_back(std::make_pair(chaincols, cols));
2160 for (std::size_t r = 0; r < nreplicas; ++r) {
2165 qn::Station<T>& stn = m.stations[srv[sidx][r] - 1];
2166 if (hasld) stn.lldscaling = lqn.lldscaling[sidx];
2172 bool one_class_each =
true;
2173 for (std::size_t j = 0; j < cols.size(); ++j)
2174 if (cols[j].size() > 1) one_class_each =
false;
2176 layer_dep_handle(lqn.cdscaling[sidx], cols, chaincols, R);
2177 const std::vector<T> pk = layer_peak(lqn.cdscalingpeak[sidx], cols, R);
2178 if (one_class_each) {
2180 stn.cdscalingpeak = pk;
2183 stn.jdscalingpeak = pk;
2187 stn.jdscaling = layer_dep_handle(lqn.jdscaling[sidx], cols, chaincols, R);
2188 stn.jdscalingpeak = layer_peak(lqn.jdscalingpeak[sidx], cols, R);
2199 const lqn::ServerPools<T>& pl = lqn.pools[sidx];
2200 const Matrix<T> compat = pl.compat;
2201 const std::vector<double> counts = pl.counts;
2202 const std::vector<T> rates = pl.rates;
2204 rates](
const std::vector<T>& nop) {
2205 return std::vector<T>(
2208 stn.jdscaling = layer_dep_handle(etaPool, cols, chaincols, R);
2209 stn.jdscalingpeak = layer_peak(
2210 std::vector<T>(cols.size(), num_traits<T>::from_int(1)), cols, R);
2219 for (std::size_t d = 0; d < deferred_chaincols.size(); ++d) {
2220 const std::vector<std::vector<std::size_t>>& cols = deferred_chaincols[d].second;
2221 std::vector<std::vector<std::size_t>>& out = *deferred_chaincols[d].first;
2222 out.assign(cols.size(), std::vector<std::size_t>());
2223 for (std::size_t j = 0; j < cols.size(); ++j) {
2224 for (std::size_t k = 0; k < cols[j].size(); ++k)
2225 for (std::size_t cc = 0; cc < m.chains.size(); ++cc)
2226 if (cols[j][k] < m.chains[cc].size() && m.chains[cc][cols[j][k]]) {
2228 for (std::size_t q = 0; q < out[j].size(); ++q)
2229 if (out[j][q] == cc) seen =
true;
2230 if (!seen) out[j].push_back(cc);
2274 std::size_t curclass;
2279 static bool is_aux_class(
const std::string& name) {
2280 return name.size() >= 4 && name.compare(name.size() - 4, 4,
".Aux") == 0;
2282 template <
class Ctx>
2283 Ctx route_sync_call(qn::Layer<T>& m, std::size_t idx, std::size_t cidx, Ctx st,
2284 const std::vector<std::size_t>& tnode,
2285 const std::vector<std::size_t>& tstat,
2286 const std::vector<std::size_t>& callcls,
2287 const std::vector<std::size_t>& auxcallcls, std::size_t clientNode,
2288 int atClient,
int atServer,
bool flat) {
2289 const T one = Tone();
2293 const std::size_t ntgt = tnode.size();
2294 const bool to_this_server = ntgt > 0;
2295 const T nrep = num_traits<T>::from_int(
int(ntgt > 0 ? ntgt : 1));
2296 const T share = T(one / nrep);
2297 const T callmean = lqn.callproc_mean[cidx];
2298 const std::size_t cc = callcls[cidx];
2299 const std::size_t ax = auxcallcls[cidx];
2300 const bool below = callmean < one;
2301 const bool above = callmean > one;
2302 if (st.jobpos == atClient) {
2303 if (to_this_server) {
2305 m.set_route(st.curclass, ax, clientNode, clientNode, T(one - callmean));
2306 for (std::size_t r = 0; r < ntgt; ++r) {
2307 m.set_route(st.curclass, cc, clientNode, tnode[r], T(callmean / nrep));
2308 m.set_route(cc, cc, tnode[r], clientNode, one);
2311 m.set_route(ax, cc, clientNode, clientNode, one);
2313 for (std::size_t r = 0; r < ntgt; ++r) {
2314 m.set_route(st.curclass, cc, clientNode, tnode[r], share);
2315 m.set_route(cc, ax, tnode[r], clientNode, one);
2316 m.set_route(ax, cc, clientNode, tnode[r], T((one - one / callmean) / nrep));
2318 m.set_route(ax, cc, clientNode, clientNode, T(one / callmean));
2320 for (std::size_t r = 0; r < ntgt; ++r) {
2321 m.set_route(st.curclass, cc, clientNode, tnode[r], share);
2322 m.set_route(cc, cc, tnode[r], clientNode, one);
2325 for (std::size_t r = 0; r < ntgt; ++r) {
2326 m.set_service(tstat[r], cc, callservtproc[cidx]);
2327 call_map.push_back({idx, cidx, tstat[r], cc});
2330 st.jobpos = atClient;
2331 st.curnodes.clear();
2334 m.set_route(st.curclass, cc, clientNode, clientNode, one);
2335 if (below || above) {
2336 m.set_route(cc, ax, clientNode, clientNode, one);
2341 st.jobpos = atClient;
2342 st.curnodes.clear();
2343 m.set_service(m.clientIdx, cc, callservtproc[cidx]);
2344 call_map.push_back({idx, cidx, m.clientIdx, cc});
2348 auto fromNode = [&](std::size_t r) {
2349 return st.curnodes[std::min(r, st.curnodes.size() - 1)];
2351 if (to_this_server) {
2360 for (std::size_t r = 0; r < ntgt; ++r) {
2361 m.set_route(st.curclass, ax, fromNode(r), clientNode, T(one - callmean));
2362 m.set_route(st.curclass, cc, fromNode(r), tnode[r], T(callmean / nrep));
2363 m.set_route(cc, cc, tnode[r], clientNode, one);
2365 m.set_route(ax, cc, clientNode, clientNode, one);
2367 st.jobpos = atClient;
2368 st.curnodes.clear();
2374 for (std::size_t r = 0; r < ntgt; ++r) {
2375 m.set_route(st.curclass, cc, fromNode(r), tnode[r], one);
2376 m.set_route(cc, ax, tnode[r], clientNode, one);
2377 m.set_route(ax, cc, clientNode, tnode[r],
2378 T((one - one / callmean) / nrep));
2380 m.set_route(ax, cc, clientNode, clientNode, T(one / callmean));
2384 for (std::size_t r = 0; r < ntgt; ++r) {
2385 m.set_route(st.curclass, cc, fromNode(r), tnode[r], one);
2386 m.set_route(cc, cc, tnode[r], tnode[r], T(one - one / callmean));
2387 m.set_route(cc, ax, tnode[r], clientNode, T(one / callmean));
2391 st.jobpos = atClient;
2392 st.curnodes.clear();
2394 for (std::size_t r = 0; r < ntgt; ++r)
2395 m.set_route(st.curclass, cc, fromNode(r), tnode[r], one);
2399 for (std::size_t r = 0; r < ntgt; ++r)
2400 m.set_route(cc, cc, tnode[r], clientNode, one);
2402 st.jobpos = atClient;
2403 st.curnodes.clear();
2405 st.jobpos = atServer;
2406 st.curnodes = tnode;
2410 for (std::size_t r = 0; r < ntgt; ++r) {
2411 m.set_service(tstat[r], cc, callservtproc[cidx]);
2412 call_map.push_back({idx, cidx, tstat[r], cc});
2415 for (std::size_t nd : st.curnodes)
2416 m.set_route(st.curclass, cc, nd, clientNode, one);
2417 if (below || above) {
2418 m.set_route(cc, ax, clientNode, clientNode, one);
2423 st.jobpos = atClient;
2424 st.curnodes.clear();
2425 m.set_service(m.clientIdx, cc, callservtproc[cidx]);
2426 call_map.push_back({idx, cidx, m.clientIdx, cc});
2441 void build_entry_service_matrix() {
2442 const std::size_t dim = lqn.nidx + lqn.ncalls;
2443 servtmatrix =
Matrix<T>(dim + 1, dim + 1, Tzero());
2444 std::function<void(std::size_t, std::size_t)> rec = [&](std::size_t aidx, std::size_t eidx) {
2445 for (std::size_t nextaidx : lqn.graph.succ(aidx)) {
2446 const bool isLoop = lqn.graph.get(aidx, nextaidx) != lqn.dag.get(aidx, nextaidx);
2447 if (lqn.parent[aidx] != lqn.parent[nextaidx]) {
2448 for (std::size_t cidx : lqn.callsof[aidx])
2449 if (lqn.calltype[cidx] == CallType::SYNC)
2450 servtmatrix(eidx, lqn.nidx + cidx) = Tone();
2451 }
else if (nextaidx != aidx && !isLoop) {
2452 servtmatrix(eidx, nextaidx) = Tone();
2453 rec(nextaidx, eidx);
2457 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
2458 const std::size_t eidx = lqn.eshift + e;
2468 void init_interlock() {
2469 const std::size_t NE = lqn.nentries;
2470 il_all =
Matrix<T>(NE + 1, NE + 1, Tzero());
2471 il_ph1 =
Matrix<T>(NE + 1, NE + 1, Tzero());
2473 std::function<void(std::size_t, std::size_t, T, T, std::vector<bool>&,
int)> trace =
2474 [&](std::size_t eidx, std::size_t root_e, T pall, T pph1, std::vector<bool>& visited,
2476 if (eidx <= lqn.eshift || eidx > lqn.eshift + NE)
return;
2477 const std::size_t e = eidx - lqn.eshift;
2478 if (visited[e])
return;
2480 il_all(root_e, e) = T(il_all(root_e, e) + pall);
2481 il_ph1(root_e, e) = T(il_ph1(root_e, e) + pph1);
2482 for (std::size_t aidx : lqn.actsof[eidx]) {
2483 if (aidx <= lqn.ashift || aidx > lqn.ashift + lqn.nacts)
continue;
2484 const std::size_t a = aidx - lqn.ashift;
2485 if (depth > 0 && lqn.actphase[a] > 1)
continue;
2486 const bool is_ph1 = lqn.actphase[a] <= 1;
2487 for (std::size_t cidx : lqn.callsof[aidx]) {
2488 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
2489 if (!(lqn.callproc_mean[cidx] > Tzero()))
continue;
2490 const std::size_t dst = lqn.callpair_dst[cidx];
2491 if (dst <= lqn.eshift || dst > lqn.eshift + NE)
continue;
2492 trace(dst, root_e, T(pall * lqn.callproc_mean[cidx]),
2493 is_ph1 ? T(pph1 * lqn.callproc_mean[cidx]) : Tzero(), visited,
2499 for (std::size_t e = 1; e <= NE; ++e) {
2500 std::vector<bool> visited(NE + 1,
false);
2501 trace(lqn.eshift + e, e, Tone(), Tone(), visited, 0);
2504 il_common_entries.assign(NT() + 1, {});
2505 il_src_all.assign(NT() + 1, {});
2506 il_src_ph2.assign(NT() + 1, {});
2507 il_num_sources.assign(NT() + 1, 0.0);
2509 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
2510 const std::size_t tidx = lqn.tshift + t;
2511 if (lqn.isref[tidx] || lqn.sched[tidx] == SchedStrategy::INF)
continue;
2512 interlock_for_server(tidx);
2514 for (std::size_t h = 1; h <= lqn.nhosts; ++h) {
2515 if (lqn.sched[h] == SchedStrategy::INF)
continue;
2516 interlock_for_server(h);
2520 std::vector<std::size_t> server_entry_nums(std::size_t serverIdx)
const {
2521 std::vector<std::size_t> out;
2522 if (serverIdx <= lqn.nhosts) {
2523 for (std::size_t tidx : lqn.tasksof[serverIdx])
2524 for (std::size_t se : lqn.entriesof[tidx]) out.push_back(se - lqn.eshift);
2526 for (std::size_t se : lqn.entriesof[serverIdx]) out.push_back(se - lqn.eshift);
2531 std::vector<std::size_t> client_tasks(std::size_t serverIdx)
const {
2532 std::vector<std::size_t> out;
2533 if (serverIdx <= lqn.nhosts)
return lqn.tasksof[serverIdx];
2534 for (std::size_t se : lqn.entriesof[serverIdx])
2535 for (std::size_t ci : lqn.iscaller.col(se))
2536 if (ci > lqn.tshift && ci <= NT()) out.push_back(ci);
2537 std::sort(out.begin(), out.end());
2538 out.erase(std::unique(out.begin(), out.end()), out.end());
2542 std::vector<std::size_t> call_dst_tasks(std::size_t src_eidx, std::size_t target_e)
const {
2543 std::vector<std::size_t> out;
2544 for (std::size_t aidx : lqn.actsof[src_eidx]) {
2545 if (aidx <= lqn.ashift || aidx > lqn.ashift + lqn.nacts)
continue;
2546 for (std::size_t cidx : lqn.callsof[aidx]) {
2547 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
2548 const std::size_t dst = lqn.callpair_dst[cidx];
2549 const std::size_t de = dst - lqn.eshift;
2550 if (de >= 1 && de <= lqn.nentries && il_all(de, target_e) > Tzero())
2551 out.push_back(lqn.parent[dst]);
2554 std::sort(out.begin(), out.end());
2555 out.erase(std::unique(out.begin(), out.end()), out.end());
2559 bool is_branch_point(std::size_t srcX, std::size_t entryA, std::size_t srcY,
2560 std::size_t entryB)
const {
2561 const std::size_t taskA = lqn.parent[entryA], taskB = lqn.parent[entryB];
2562 const std::size_t taskX = lqn.parent[srcX];
2563 if (taskX == taskA && taskX == taskB)
return false;
2564 if (srcX == entryA || srcY == entryB)
return true;
2565 const std::vector<std::size_t> dx = call_dst_tasks(srcX, entryA - lqn.eshift);
2566 const std::vector<std::size_t> dy = call_dst_tasks(srcY, entryB - lqn.eshift);
2567 for (std::size_t a : dx)
2568 for (std::size_t b : dy)
2569 if (a != b)
return true;
2573 void trace_to_server(std::size_t eidx, std::size_t serverIdx, std::vector<bool>& visited,
2574 std::vector<std::size_t>& itasks,
bool isHead)
const {
2575 if (eidx <= lqn.eshift || eidx > lqn.eshift + lqn.nentries)
return;
2576 const std::size_t e = eidx - lqn.eshift;
2577 if (visited[e])
return;
2578 const std::size_t owner = lqn.parent[eidx];
2579 if (owner == serverIdx)
return;
2580 if (serverIdx <= lqn.nhosts && lqn.parent[owner] == serverIdx)
return;
2583 const std::vector<std::size_t> sen = server_entry_nums(serverIdx);
2584 for (std::size_t aidx : lqn.actsof[eidx]) {
2585 if (aidx <= lqn.ashift || aidx > lqn.ashift + lqn.nacts)
continue;
2586 for (std::size_t cidx : lqn.callsof[aidx]) {
2587 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
2588 const std::size_t dst = lqn.callpair_dst[cidx];
2589 const std::size_t dtask = lqn.parent[dst];
2590 bool reaches = dtask == serverIdx ||
2591 (serverIdx <= lqn.nhosts && lqn.parent[dtask] == serverIdx);
2593 const std::size_t de = dst - lqn.eshift;
2594 for (std::size_t sn : sen)
2595 if (il_all(de, sn) > Tzero()) reaches =
true;
2598 trace_to_server(dst, serverIdx, visited, itasks,
false);
2603 if (found && !isHead) {
2604 itasks.push_back(owner);
2605 std::sort(itasks.begin(), itasks.end());
2606 itasks.erase(std::unique(itasks.begin(), itasks.end()), itasks.end());
2611 void interlock_for_server(std::size_t serverIdx) {
2612 const std::vector<std::size_t> sen = server_entry_nums(serverIdx);
2613 if (sen.empty())
return;
2614 const std::vector<std::size_t> cts = client_tasks(serverIdx);
2615 if (cts.empty())
return;
2617 std::vector<std::pair<std::size_t, std::size_t>> pairs;
2618 for (std::size_t ct : cts)
2619 for (std::size_t ce : lqn.entriesof[ct]) {
2620 const std::size_t cen = ce - lqn.eshift;
2621 if (cen < 1 || cen > lqn.nentries)
continue;
2622 for (std::size_t se : sen)
2623 if (il_all(cen, se) > Tzero()) {
2624 pairs.emplace_back(ct, cen);
2628 if (pairs.size() < 2)
return;
2630 std::vector<std::size_t> common;
2631 for (std::size_t i = 0; i < pairs.size(); ++i)
2632 for (std::size_t j = i + 1; j < pairs.size(); ++j) {
2633 if (pairs[i].first == pairs[j].first)
continue;
2634 const std::size_t eA = pairs[i].second, eC = pairs[j].second;
2635 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
2636 const std::size_t tidx = lqn.tshift + t;
2637 for (std::size_t ex : lqn.entriesof[tidx])
2638 for (std::size_t ey : lqn.entriesof[tidx]) {
2639 const std::size_t xn = ex - lqn.eshift, yn = ey - lqn.eshift;
2640 if (xn < 1 || yn < 1 || xn > lqn.nentries || yn > lqn.nentries)
continue;
2641 if (il_all(xn, eA) > Tzero() && il_all(yn, eC) > Tzero() &&
2642 is_branch_point(ex, eA + lqn.eshift, ey, eC + lqn.eshift))
2643 common.push_back(ex);
2647 std::sort(common.begin(), common.end());
2648 common.erase(std::unique(common.begin(), common.end()), common.end());
2649 if (common.empty())
return;
2651 std::vector<std::size_t> interlocked;
2652 for (std::size_t ce : common) {
2653 std::vector<bool> visited(lqn.nentries + 1,
false);
2654 std::vector<std::size_t> it;
2655 trace_to_server(ce, serverIdx, visited, it,
true);
2656 for (std::size_t x : it) interlocked.push_back(x);
2658 std::sort(interlocked.begin(), interlocked.end());
2659 interlocked.erase(std::unique(interlocked.begin(), interlocked.end()), interlocked.end());
2661 std::vector<std::size_t> src_all;
2662 for (std::size_t ce : common) src_all.push_back(lqn.parent[ce]);
2663 std::sort(src_all.begin(), src_all.end());
2664 src_all.erase(std::unique(src_all.begin(), src_all.end()), src_all.end());
2666 std::vector<std::size_t> diff;
2667 for (std::size_t x : src_all)
2668 if (!std::binary_search(interlocked.begin(), interlocked.end(), x)) diff.push_back(x);
2676 std::vector<std::size_t> src_ph2;
2677 for (std::size_t it : interlocked)
2678 for (std::size_t ie : lqn.entriesof[it])
2679 for (std::size_t ci : lqn.iscaller.col(ie))
2680 if (ci > lqn.tshift && ci <= NT() &&
2681 !std::binary_search(interlocked.begin(), interlocked.end(), ci))
2682 src_all.push_back(ci);
2683 std::sort(src_all.begin(), src_all.end());
2684 src_all.erase(std::unique(src_all.begin(), src_all.end()), src_all.end());
2687 for (std::size_t st : src_all) nsrc += lqn.mult[st];
2689 il_common_entries[serverIdx] = common;
2690 il_src_all[serverIdx] = src_all;
2691 il_src_ph2[serverIdx] = src_ph2;
2692 il_num_sources[serverIdx] = nsrc;
2699 const std::size_t N = lqn.nidx;
2700 tput.assign(N + 1, Tzero());
2701 util.assign(N + 1, Tzero());
2702 servt.assign(N + 1, Tzero());
2703 residt.assign(N + 1, Tzero());
2704 thinkt.assign(N + 1, Tzero());
2705 callservt.assign(lqn.ncalls + 1, Tzero());
2706 callresidt.assign(lqn.ncalls + 1, Tzero());
2708 servt_ph1.assign(N + 1, Tzero());
2709 servt_ph2.assign(N + 1, Tzero());
2710 prOvertake.assign(lqn.nentries + 1, Tzero());
2711 build_entry_service_matrix();
2713 relax_omega = (opt.relax ==
"fixed" || opt.relax ==
"adaptive") ? opt.relax_factor : 1.0;
2715 servt_prev.assign(N + 1, std::numeric_limits<double>::quiet_NaN());
2716 residt_prev.assign(N + 1, std::numeric_limits<double>::quiet_NaN());
2717 tput_prev.assign(N + 1, std::numeric_limits<double>::quiet_NaN());
2718 thinkt_prev.assign(N + 1, std::numeric_limits<double>::quiet_NaN());
2719 callservt_prev.assign(lqn.ncalls + 1, std::numeric_limits<double>::quiet_NaN());
2720 callresidt_prev.assign(lqn.ncalls + 1, std::numeric_limits<double>::quiet_NaN());
2721 servt_prev_v.assign(N + 1, Tzero());
2722 residt_prev_v.assign(N + 1, Tzero());
2723 tput_prev_v.assign(N + 1, Tzero());
2724 thinkt_prev_v.assign(N + 1, Tzero());
2725 callservt_prev_v.assign(lqn.ncalls + 1, Tzero());
2727 unique_route_idx.clear();
2728 for (
const RouteRow& r : route_map) unique_route_idx.push_back(r.idx);
2729 std::sort(unique_route_idx.begin(), unique_route_idx.end());
2730 unique_route_idx.erase(std::unique(unique_route_idx.begin(), unique_route_idx.end()),
2731 unique_route_idx.end());
2733 if (opt.interlocking) init_interlock();
2735 maxitererr.assign(opt.iter_max + 2, 0.0);
2736 averagingstart = -1;
2737 hasconverged =
false;
2738 moment_pass_done =
false;
2741 servtcdf.assign(N + 1, fluid::FluidPassage());
2742 callservtcdf.assign(lqn.ncalls + 1, fluid::FluidPassage());
2743 entrycdfrespt.assign(lqn.nentries + 1, LnCdf());
2744 entryproc.assign(lqn.nentries + 1, mam::AphPair<T>());
2745 cdf_repo.assign(ensemble.size(), std::vector<std::vector<fluid::FluidPassage>>());
2759 const bool stoch = opt.layer_solver ==
"ssa";
2762 cfg.iter_tol = opt.iter_tol;
2763 cfg.relax_burnin = relax_omega;
2764 stoch_ctl = std::make_shared<LnStochController<T>>(cfg);
2767 while (it < opt.iter_max) {
2768 if (!stoch && converged(it))
break;
2770 results.emplace_back(ensemble.size());
2771 for (std::size_t e = 0; e < ensemble.size(); ++e) analyze(it, e);
2774 std::vector<double> jobs(ensemble.size(), 0.0);
2775 for (std::size_t e = 0; e < ensemble.size(); ++e)
2776 jobs[e] = ensemble[e].total_jobs();
2777 const bool stop = stoch_ctl->update(it, results.back(), jobs, servt, residt);
2780 relax_omega = stoch_ctl->relax_omega();
2782 did_converge =
true;
2783 hasconverged =
true;
2788 iterations_done = it;
2791 if (stoch && stoch_ctl && stoch_ctl->averaging_count() > 0) {
2792 const std::vector<LayerResult<T>>& avg = stoch_ctl->averaged_results();
2793 for (std::size_t e = 0; e < results.back().size() && e < avg.size(); ++e)
2794 results.back()[e] = avg[e];
2795 servt = stoch_ctl->averaged_servt();
2796 residt = stoch_ctl->averaged_residt();
2820 Matrix<T> filter_metric(
const qn::Layer<T>& L,
const Matrix<T>& metric,
2821 const std::vector<std::vector<bool>>* zero_mask)
const {
2822 const std::size_t M = L.nstations, K = L.nclasses;
2823 Matrix<T> out(M, K, Tzero());
2824 for (std::size_t i = 0; i < M; ++i)
2825 for (std::size_t k = 0; k < K; ++k)
2826 if (!L.disabled[i][k]) out(i, k) = metric(i, k);
2828 for (std::size_t i = 0; i < M; ++i)
2829 for (std::size_t k = 0; k < K; ++k)
2830 if ((*zero_mask)[i][k]) out(i, k) = Tzero();
2831 for (std::size_t i = 0; i < M; ++i)
2832 for (std::size_t k = 0; k < K; ++k)
2834 for (std::size_t k = 0; k < K; ++k) {
2835 std::size_t c = L.nchains;
2836 for (std::size_t cc = 0; cc < L.nchains; ++cc)
2837 if (L.chains[cc][k]) c = cc;
2838 if (c == L.nchains)
continue;
2839 for (std::size_t i = 0; i < M; ++i)
2840 if (L.visits[c](L.stateful_of_station(i + 1) - 1, k) == Tzero())
2841 out(i, k) = Tzero();
2888 mva::MvaSolution<T> solve_layer_ctmc(std::size_t e) {
2889 qn::Layer<T>& L = ensemble[e];
2894 ctmc::CtmcOptions co;
2899 mva::MvaSolution<T> out;
2900 out.method =
"ctmc";
2919 mva::MvaSolution<T> solve_layer_mam(std::size_t e) {
2920 qn::Layer<T>& L = ensemble[e];
2922 L.refresh_capacity();
2924 mo.method =
"dec.poisson";
2940 mva::MvaSolution<T> solve_layer_nc(std::size_t e) {
2941 if (fj_tr[e].active())
2943 "SolverLN: layer '" + ensemble[e].name +
2944 "' carries a fork, whose fixed point is driven by MVA in this port; solve this "
2945 "model with layer_solver 'mva'");
2946 qn::Layer<T>& L = ensemble[e];
2976 mva::MvaSolution<T> solve_layer_ssa(std::size_t e) {
2977 if (fj_tr[e].active())
2979 "SolverLN: layer '" + ensemble[e].name +
2980 "' carries a fork, whose fixed point is driven by MVA in this port; solve this "
2981 "model with layer_solver 'mva'");
2982 qn::Layer<T>& L = ensemble[e];
2985 mva::MvaSolution<T> out;
2988 out.Q =
Matrix<T>(L.nstations, L.nclasses, Tzero());
2989 out.U =
Matrix<T>(L.nstations, L.nclasses, Tzero());
2990 out.R =
Matrix<T>(L.nstations, L.nclasses, Tzero());
2991 out.Tp =
Matrix<T>(L.nstations, L.nclasses, Tzero());
2992 for (std::size_t i = 0; i < L.nstations; ++i)
2993 for (std::size_t r = 0; r < L.nclasses; ++r) {
2994 out.Q(i, r) = num_traits<T>::from_double(s.QN(i, r));
2995 out.U(i, r) = num_traits<T>::from_double(s.UN(i, r));
2996 out.R(i, r) = num_traits<T>::from_double(s.RN(i, r));
2997 out.Tp(i, r) = num_traits<T>::from_double(s.TN(i, r));
2999 out.C.assign(L.nclasses, Tzero());
3000 out.X.assign(L.nclasses, Tzero());
3001 for (std::size_t r = 0; r < L.nclasses && r < s.XN.size(); ++r) {
3002 out.X[r] = num_traits<T>::from_double(s.XN[r]);
3003 out.C[r] = num_traits<T>::from_double(s.CN[r]);
3008 mva::MvaSolution<T> solve_layer_fluid(std::size_t e) {
3009 if (fj_tr[e].active())
3011 "SolverLN: layer '" + std::to_string(e) +
3012 "' carries a fork, whose fixed point is driven by MVA in this port; solve this "
3013 "model with layer_solver 'mva'");
3014 qn::Layer<T>& L = ensemble[e];
3024 mva::MvaSolution<T> out;
3025 out.method =
"fluid";
3027 out.Q =
Matrix<T>(L.nstations, L.nclasses, Tzero());
3028 out.U =
Matrix<T>(L.nstations, L.nclasses, Tzero());
3029 out.R =
Matrix<T>(L.nstations, L.nclasses, Tzero());
3030 out.Tp =
Matrix<T>(L.nstations, L.nclasses, Tzero());
3031 out.C.assign(L.nclasses, Tzero());
3032 out.X.assign(L.nclasses, Tzero());
3033 detail::ln_fluid_solve(L, opt.layer_fluid, out);
3037 mva::MvaSolution<T> solve_layer(std::size_t e) {
3047 ensemble[e].refresh_capacity();
3055 if (!ensemble[e].regions.empty())
return solve_layer_ctmc(e);
3068 if (has_cache_node(e)) ensemble[e].refresh_rt();
3069 if (opt.layer_solver ==
"fluid")
return solve_layer_fluid(e);
3070 if (opt.layer_solver ==
"nc")
return solve_layer_nc(e);
3071 if (opt.layer_solver ==
"ssa")
return solve_layer_ssa(e);
3072 if (opt.layer_solver !=
"mva")
3074 "' is not available; use 'mva', 'nc', 'fluid' or 'ssa'");
3081 mva::MvaOptions eopt = opt.layer;
3082 if (e < layer_interlock.size()) eopt.interlock = layer_interlock[e];
3083 if (!fj_tr[e].active())
3085 mva::MvaOptions lopt = eopt;
3086 if (lopt.method ==
"default") lopt.method =
"amva";
3088 ensemble[e], fj_tr[e], fj_lambda[e], eopt,
3089 [&lopt](qn::NetworkStruct<T>& V) {
3094 void analyze(
int it, std::size_t e) {
3095 qn::Layer<T>& L = ensemble[e];
3096 const mva::MvaSolution<T> s = solve_layer(e);
3097 LayerResult<T>& r = results[it - 1][e];
3098 r.RN = filter_metric(L, s.R,
nullptr);
3099 std::vector<std::vector<bool>> zmask(L.nstations, std::vector<bool>(L.nclasses,
false));
3100 for (std::size_t i = 0; i < L.nstations; ++i)
3101 for (std::size_t k = 0; k < L.nclasses; ++k)
3103 r.QN = filter_metric(L, s.Q, &zmask);
3104 r.UN = filter_metric(L, s.U, &zmask);
3105 r.TN = filter_metric(L, s.Tp,
nullptr);
3106 r.WN = residence_from_response(L, r.RN);
3116 if (!fj_tr[e].active() && opt.layer_solver ==
"mva") {
3117 Matrix<T> Qch(L.nstations, L.nchains, Tzero());
3118 for (std::size_t c = 0; c < L.nchains; ++c)
3119 for (std::size_t i = 0; i < L.nstations; ++i) {
3121 for (std::size_t k : L.inchain[c]) s2 += s.Q(i, k - 1);
3124 layer_init_sol[e] = Qch;
3129 Matrix<T> residence_from_response(
const qn::Layer<T>& L,
const Matrix<T>& RN)
const {
3130 Matrix<T> V(L.nstations, L.nclasses, Tzero());
3131 for (std::size_t c = 0; c < L.nchains; ++c)
3132 for (std::size_t i = 0; i < L.nstations; ++i) {
3133 const std::size_t sf = L.stateful_of_station(i + 1) - 1;
3134 for (std::size_t k = 0; k < L.nclasses; ++k)
3135 V(i, k) = T(V(i, k) + L.visits[c](sf, k));
3137 Matrix<T> WN(L.nstations, L.nclasses, Tzero());
3138 for (std::size_t i = 0; i < L.nstations; ++i)
3139 for (std::size_t k = 0; k < L.nclasses; ++k) {
3140 if (L.disabled[i][k])
continue;
3141 if (!(RN(i, k) > Tzero()))
continue;
3143 WN(i, k) = RN(i, k);
3147 for (std::size_t cc = 0; cc < L.nchains; ++cc)
3148 if (L.chains[cc][k]) c = cc;
3149 const std::size_t rstat = L.classes[k].refstat;
3151 if (L.refclass[c] > 0) {
3152 den = V(rstat - 1, L.refclass[c] - 1);
3154 for (std::size_t kk : L.inchain[c]) den += V(rstat - 1, kk - 1);
3156 if (den == Tzero())
continue;
3157 WN(i, k) = T(RN(i, k) * V(i, k) / den);
3159 for (std::size_t i = 0; i < L.nstations; ++i)
3160 for (std::size_t k = 0; k < L.nclasses; ++k)
3167 if (opt.interlocking) update_populations(it);
3168 update_think_times(it);
3170 update_routing_probabilities(it);
3171 for (std::size_t e : route_reset) {
3172 ensemble[e].refresh_chains();
3178 for (std::size_t e : svc_reset) ensemble[e].refresh_rates();
3182 bool converged(
int it) {
3183 const std::size_t E = ensemble.size();
3184 const int iter_min = std::max<int>(2 *
int(E),
int(std::ceil(opt.iter_max / 4.0)));
3185 const int wnd_size = std::max(5,
int(std::ceil(iter_min / 5.0)));
3187 if (it >= iter_min &&
int(results.size()) >= wnd_size) {
3188 const T w = T(Tone() / num_traits<T>::from_int(wnd_size));
3189 for (std::size_t e = 0; e < E; ++e) {
3190 LayerResult<T>& cur = results[results.size() - 1][e];
3191 auto scale = [&](Matrix<T>& A) {
3192 for (std::size_t i = 0; i < A.rows(); ++i)
3193 for (std::size_t j = 0; j < A.cols(); ++j) A(i, j) = T(A(i, j) * w);
3195 Matrix<T> Q = cur.QN, U = cur.UN, R = cur.RN, Tp = cur.TN, W = cur.WN;
3196 scale(Q); scale(U); scale(R); scale(Tp); scale(W);
3197 for (
int k = 1; k < wnd_size; ++k) {
3198 const LayerResult<T>& old = results[results.size() - 1 - k][e];
3199 auto add = [&](Matrix<T>& A,
const Matrix<T>& B) {
3200 for (std::size_t i = 0; i < A.rows(); ++i)
3201 for (std::size_t j = 0; j < A.cols(); ++j) A(i, j) = T(A(i, j) + B(i, j) * w);
3203 add(Q, old.QN); add(U, old.UN); add(R, old.RN); add(Tp, old.TN); add(W, old.WN);
3205 cur.QN = Q; cur.UN = U; cur.RN = R; cur.TN = Tp; cur.WN = W;
3211 for (std::size_t e = 0; e < E; ++e) {
3212 const Matrix<T>& Q = results[results.size() - 1][e].QN;
3213 const Matrix<T>& Q1 = results[results.size() - 2][e].QN;
3214 const double Njobs = ensemble[e].total_jobs();
3215 if (!(Njobs > 0.0))
continue;
3217 for (std::size_t i = 0; i < Q.rows(); ++i)
3218 for (std::size_t j = 0; j < Q.cols(); ++j)
3219 mx = std::max(mx, std::fabs(dbl(Q(i, j)) - dbl(Q1(i, j))));
3222 maxitererr[it] = err;
3224 static_cast<long>(it),
3225 "layer iteration %zu: max queue-length change %.3e (tolerance %.3e)",
3226 static_cast<std::size_t
>(it), err, opt.iter_tol);
3227 if (it == iter_min) {
3229 averagingstart = it;
3233 if (it > iter_min && maxitererr[it] < opt.iter_tol && maxitererr[it - 1] < opt.iter_tol &&
3234 maxitererr[it - 2] < opt.iter_tol) {
3235 if (!hasconverged) {
3236 hasconverged =
true;
3238 did_converge =
true;
3242 hasconverged =
false;
3277 std::string lnmethod;
3279 bool ph_laws_ready =
false;
3296 static std::string ln_requested_method(
const std::string& method) {
3298 for (
char c : method) m +=
static_cast<char>(std::tolower(
static_cast<unsigned char>(c)));
3299 if (m.empty() || m ==
"srvn" || m ==
"default" || m ==
"auto")
return "srvn";
3300 if (m ==
"srvn.ph" || m ==
"ph")
return "srvn.ph";
3301 if (m ==
"srvn.cs" || m ==
"srvncs" || m ==
"cs")
return "srvn.cs";
3302 if (m ==
"flat.cs" || m ==
"flatcs" || m ==
"flat" || m ==
"squashed")
return "flat.cs";
3303 if (m ==
"flat.ph" || m ==
"flatph" || m ==
"squashed.ph")
return "flat.ph";
3304 if (m ==
"moment3")
return "moment3";
3311 bool is_srvn_ph()
const {
return lnmethod ==
"srvn.ph"; }
3319 bool is_ph_encoding()
const {
return lnmethod ==
"srvn.ph" || lnmethod ==
"flat.ph"; }
3327 static std::size_t station_idx_of(
const qn::Layer<T>& L, std::size_t elem) {
3328 if (elem >= 1 && elem < L.server_idx_of.size() && L.server_idx_of[elem] != 0)
3329 return L.server_idx_of[elem];
3337 std::size_t station_idx_of_class(
const qn::Layer<T>& L, std::size_t k)
const {
3338 const int kind = L.classes[k].attr_kind;
3339 const std::size_t a = L.classes[k].attr_idx;
3340 if (kind ==
int(LqnElement::ACTIVITY))
3341 return station_idx_of(L, lqn.parent[lqn.parent[a]]);
3342 if (kind ==
int(LqnElement::CALL))
3343 return station_idx_of(L, lqn.parent[lqn.callpair_dst[a]]);
3355 bool probe_srvn_ph() {
3357 ph_assert_supported();
3359 ph_laws_ready =
true;
3361 }
catch (
const std::exception&) {
3362 ph_laws_ready =
false;
3369 std::size_t idx = 0;
3370 bool ishost =
false;
3371 std::vector<std::size_t> callers;
3373 std::vector<std::size_t> class_of_caller;
3374 std::size_t nreplicas = 1;
3376 std::vector<std::size_t> qstations;
3378 std::vector<T> svcmean_by_class;
3380 std::vector<std::pair<std::size_t, long>> open_arrivals;
3390 std::vector<workflow::Workflow<T>> ph_wf, ph_wfhost;
3391 std::vector<std::unordered_map<std::size_t, T>> ph_execs;
3392 std::vector<bool> ph_has_wf;
3393 std::vector<workflow::PhLaw<T>> ph_hostlaw, ph_entrylaw;
3394 std::vector<T> ph_hostmean, ph_entrymean, ph_entryscv;
3395 std::vector<T> ph_share, ph_overlap, ph_setupshare, ph_xdemand;
3396 Matrix<T> ph_ncalls, ph_calltime;
3397 std::vector<T> ph_procresid, ph_actthinkt, ph_calltotal;
3398 std::vector<PHLayer> ph_layers;
3399 std::vector<bool> ph_has_layer;
3402 T act_think_time(std::size_t aidx)
const {
3403 if (aidx >= lqn.actthink.size() || lqn.actthink[aidx].disabled)
return Tzero();
3404 const double m = dbl(lqn.actthink[aidx].mean);
3406 return lqn.actthink[aidx].mean;
3427 double setup_charge(std::size_t tidx)
const {
3428 if (tidx >= lqn.hassetup.size() || !lqn.hassetup[tidx])
return 0.0;
3429 const double s = tidx < lqn.setuptime.size() && !lqn.setuptime[tidx].disabled
3430 ? dbl(lqn.setuptime[tidx].mean) : 0.0;
3431 const double d = tidx < lqn.delayofftime.size() && !lqn.delayofftime[tidx].disabled
3432 ? dbl(lqn.delayofftime[tidx].mean) : 0.0;
3434 const double mult = lqn.mult[tidx];
3435 if (!std::isfinite(mult) || mult <= 0.0)
return 0.0;
3436 if (tidx >= tput.size() || tidx >= util.size())
return s;
3437 const double X = dbl(tput[tidx]);
3439 double rho = dbl(util[tidx]);
3440 if (!std::isfinite(rho) || rho < 0.0) rho = 0.0;
3442 const double b = rho * mult;
3443 const double EI = (std::max(1.0, b) - b) / X;
3444 return s * EI / (EI + d);
3453 double ph_setup_prob(std::size_t eidx)
const {
3454 const std::size_t tidx = lqn.parent[eidx];
3455 if (tidx >= lqn.hassetup.size() || !lqn.hassetup[tidx])
return 0.0;
3456 const double s = tidx < lqn.setuptime.size() && !lqn.setuptime[tidx].disabled
3457 ? dbl(lqn.setuptime[tidx].mean) : 0.0;
3458 const double d = tidx < lqn.delayofftime.size() && !lqn.delayofftime[tidx].disabled
3459 ? dbl(lqn.delayofftime[tidx].mean) : 0.0;
3461 return std::min(1.0, std::max(0.0, setup_charge(tidx) / s));
3465 double ph_host_servers(std::size_t hidx)
const {
3466 if (lqn.sched[hidx] == SchedStrategy::INF)
return 1.0;
3467 const double m = lqn.maxmult[hidx];
3468 return (std::isfinite(m) && m > 0.0) ? m : 1.0;
3472 bool ph_any_caller_of(std::size_t eidx)
const {
3473 return lqn.issynccaller.any_col(eidx) || lqn.isasynccaller.any_col(eidx);
3481 bool ph_open_arrival_only(std::size_t tidx)
const {
3482 if (lqn.isref[tidx])
return false;
3483 for (std::size_t eidx : lqn.entriesof[tidx])
3484 if (ph_any_caller_of(eidx))
return false;
3485 for (std::size_t eidx : lqn.entriesof[tidx])
3486 if (lqn.has_arrival[eidx])
return true;
3491 std::vector<std::size_t> ph_async_calls_into(std::size_t tidx)
const {
3492 std::vector<std::size_t> out;
3493 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
3494 if (lqn.calltype[cidx] != CallType::ASYNC)
continue;
3495 for (std::size_t e : lqn.entriesof[tidx])
3496 if (lqn.callpair_dst[cidx] == e) { out.push_back(cidx);
break; }
3505 void ph_assert_supported(
bool flat =
false)
const {
3508 const std::string mname = flat ?
"flat.ph" :
"srvn.ph";
3514 "method='" + mname +
"' does not support second-phase activities: the composed "
3515 "entry law has no reply point. Use method='default'.");
3521 const std::string why = ph_method_refusal(mname);
3540 std::string ph_method_refusal(
const std::string& method)
const {
3541 if (method !=
"srvn.ph" && method !=
"flat.ph")
return std::string();
3544 if (method ==
"flat.ph") {
3545 for (std::size_t i = 1; i <= NT(); ++i) {
3546 if (lqn.repl[i] > 1.0)
3547 return "method='flat.ph' does not support replicated processors or tasks, "
3548 "whose replicas need a submodel each. Use method='srvn.ph'.";
3549 if (lqn.hassetup[i])
3550 return "method='flat.ph' does not support setup tasks, whose powered-down "
3551 "threads are per-layer state. Use method='srvn.ph'.";
3554 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx)
3555 if (lqn.calltype[cidx] == CallType::FWD)
3556 return "method='" + method +
3557 "' does not support forwarding calls, whose target is not part of the "
3558 "caller's activity graph. Use method='default'.";
3559 for (std::size_t i = 1; i < lqn.iscache.size(); ++i)
3561 return "method='" + method +
3562 "' does not support cache tasks. Use method='default'.";
3563 for (std::size_t i = 1; i < lqn.hassetup.size(); ++i) {
3564 if (!lqn.hassetup[i])
continue;
3565 if (lqn.sched[i] == SchedStrategy::INF || !std::isfinite(lqn.mult[i]))
3566 return "method='" + method +
"': task '" + lqn.names[i] +
3567 "' declares a setup time on an infinite-server task, which holds no "
3568 "thread to power down; give it a finite multiplicity.";
3570 for (std::size_t i = 0; i < lqn.lincon_A.size(); ++i)
3571 if (lqn.lincon_A[i].rows() > 0)
3572 return "method='" + method +
3573 "' does not support admission constraints on a layer station. "
3574 "Use method='default'.";
3575 for (std::size_t i = 1; i < lqn.lldscaling.size(); ++i) {
3576 const char* fname =
nullptr;
3577 if (!lqn.lldscaling[i].empty()) fname =
"a load dependence";
3578 else if (lqn.cdscaling[i]) fname =
"a class dependence";
3579 else if (lqn.jdscaling[i]) fname =
"a joint dependence";
3580 else if (!lqn.pools[i].empty()) fname =
"server pools";
3581 if (fname !=
nullptr)
3582 return "method='" + method +
3583 "' does not support queue-dependent service rates on a layer station ('" +
3584 lqn.names[i] +
"' declares " + fname +
"). Use method='srvn.cs'.";
3586 return std::string();
3590 void ph_init_laws() {
3591 const std::size_t N = lqn.nidx;
3592 ph_wf.assign(N + 1, workflow::Workflow<T>(
"empty"));
3593 ph_wfhost.assign(N + 1, workflow::Workflow<T>(
"empty"));
3594 ph_execs.assign(N + 1, {});
3595 ph_has_wf.assign(N + 1,
false);
3596 ph_hostlaw.assign(N + 1, workflow::PhLaw<T>());
3597 ph_entrylaw.assign(N + 1, workflow::PhLaw<T>());
3598 ph_hostmean.assign(N + 1, Tzero());
3599 ph_entrymean.assign(N + 1, Tzero());
3600 ph_entryscv.assign(N + 1, Tone());
3601 ph_share.assign(N + 1, Tzero());
3602 ph_overlap.assign(N + 1, Tone());
3603 ph_setupshare.assign(N + 1, Tzero());
3604 ph_xdemand.assign(N + 1, Tzero());
3605 ph_ncalls =
Matrix<T>(N + 1, N + 1, Tzero());
3606 ph_calltime =
Matrix<T>(N + 1, N + 1, Tzero());
3607 ph_procresid.assign(N + 1, Tzero());
3608 ph_actthinkt.assign(N + 1, Tzero());
3609 ph_calltotal.assign(N + 1, Tzero());
3610 ph_layers.assign(NT() + 1, PHLayer());
3611 ph_has_layer.assign(NT() + 1,
false);
3613 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
3614 const std::size_t eidx = lqn.eshift + e;
3615 const std::size_t tidx = lqn.parent[eidx];
3616 if (ignore[tidx])
continue;
3618 ph_wf[eidx] = std::move(ew.wf);
3619 ph_execs[eidx] = ew.execs;
3620 ph_has_wf[eidx] =
true;
3622 ph_wfhost[eidx] = std::move(eh.wf);
3632 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
3633 const std::size_t tidx = lqn.tshift + t;
3634 const std::size_t n = lqn.entriesof[tidx].size();
3635 if (n == 0)
continue;
3636 for (std::size_t eidx : lqn.entriesof[tidx])
3637 ph_share[eidx] = num_traits<T>::from_double(1.0 /
double(n));
3648 Distrib<T> ph_call_burst_law(std::size_t cidx)
const {
3649 const double m = dbl(lqn.callproc_mean[cidx]);
3650 const std::size_t eidx = lqn.callpair_dst[cidx];
3652 const T R = T(callservt[cidx] / lqn.callproc_mean[cidx]);
3653 double scvd = dbl(ph_entryscv[eidx]);
3659 const workflow::PhLaw<T> loop =
3673 Distrib<T> ph_station_law(
const workflow::PhLaw<T>& law)
const {
3681 static workflow::PhLaw<T> ph_immediate_law() {
3682 workflow::PhLaw<T> out;
3683 out.alpha.assign(1, Tone());
3697 void ph_compose_entry_laws() {
3698 const std::size_t N = lqn.nidx;
3699 std::vector<T> entry_setup_share(N + 1, Tzero());
3700 ph_overlap.assign(N + 1, Tone());
3702 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
3703 const std::size_t eidx = lqn.eshift + e;
3704 if (!ph_has_wf[eidx])
continue;
3705 workflow::Workflow<T>& w = ph_wf[eidx];
3706 const std::unordered_map<std::size_t, T>& ex = ph_execs[eidx];
3707 T entrysum = Tzero(), procsum = Tzero();
3708 for (std::size_t aidx : lqn.actsof[eidx]) {
3709 T m = T(residt[aidx] + act_think_time(aidx));
3710 const T xa = ex.at(aidx);
3711 procsum += T(xa * m);
3712 const double md = dbl(m);
3713 w.set_activity_demand_mean(
3717 for (std::size_t cidx : lqn.callsof[aidx]) {
3718 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
3719 w.set_activity_demand(lqn.callhashnames[cidx], ph_call_burst_law(cidx));
3720 m = T(m + callservt[cidx]);
3722 entrysum += T(xa * m);
3726 T m1 = mm.first, scv = mm.second;
3733 const T f = T(m1 / procsum);
3734 for (std::size_t i = 0; i < law.S.rows(); ++i)
3735 for (std::size_t j = 0; j < law.S.cols(); ++j) law.S(i, j) = T(law.S(i, j) * f);
3744 const double p = ph_setup_prob(eidx);
3746 const std::size_t tidx = lqn.parent[eidx];
3747 const double sm = dbl(lqn.setuptime[tidx].mean);
3748 double sscv = dbl(lqn.setuptime[tidx].scv);
3752 num_traits<T>::from_double(sm), num_traits<T>::from_double(sscv)));
3753 std::vector<workflow::PhLaw<T>> mix;
3756 std::vector<T> probs;
3757 probs.push_back(num_traits<T>::from_double(p));
3758 probs.push_back(num_traits<T>::from_double(1.0 - p));
3768 entry_setup_share[eidx] = num_traits<T>::from_double(p * sm / denom);
3771 ph_entrylaw[eidx] = law;
3772 ph_entrymean[eidx] = m1;
3773 ph_entryscv[eidx] = scv;
3775 const double r = std::min(1.0, dbl(m1) / dbl(entrysum));
3776 ph_overlap[eidx] = num_traits<T>::from_double(r);
3782 ph_setupshare.assign(N + 1, Tzero());
3783 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
3784 const std::size_t tidx = lqn.tshift + t;
3785 if (ignore[tidx])
continue;
3786 for (std::size_t eidx : lqn.entriesof[tidx])
3787 ph_setupshare[tidx] = T(ph_setupshare[tidx] + ph_share[eidx] * entry_setup_share[eidx]);
3791 ph_ncalls =
Matrix<T>(N + 1, N + 1, Tzero());
3792 ph_calltime =
Matrix<T>(N + 1, N + 1, Tzero());
3793 ph_procresid.assign(N + 1, Tzero());
3794 ph_actthinkt.assign(N + 1, Tzero());
3795 ph_calltotal.assign(N + 1, Tzero());
3796 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
3797 const std::size_t tidx = lqn.tshift + t;
3798 if (ignore[tidx])
continue;
3799 for (std::size_t eidx : lqn.entriesof[tidx]) {
3800 if (!ph_has_wf[eidx])
continue;
3801 const T w = ph_share[eidx];
3802 if (!(dbl(w) > 0.0))
continue;
3803 const std::unordered_map<std::size_t, T>& ex = ph_execs[eidx];
3804 const T r = ph_overlap[eidx];
3805 for (std::size_t aidx : lqn.actsof[eidx]) {
3806 const T xa = ex.at(aidx);
3807 ph_procresid[tidx] = T(ph_procresid[tidx] + w * r * xa * residt[aidx]);
3808 ph_actthinkt[tidx] = T(ph_actthinkt[tidx] + w * r * xa * act_think_time(aidx));
3809 for (std::size_t cidx : lqn.callsof[aidx]) {
3810 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
3811 const std::size_t tgte = lqn.callpair_dst[cidx];
3812 const std::size_t tgtt = lqn.parent[tgte];
3815 ph_ncalls(tidx, tgte) =
3816 T(ph_ncalls(tidx, tgte) + w * xa * lqn.callproc_mean[cidx]);
3817 ph_calltime(tidx, tgtt) =
3818 T(ph_calltime(tidx, tgtt) + w * r * xa * callservt[cidx]);
3819 ph_calltotal[tidx] = T(ph_calltotal[tidx] + w * r * xa * callservt[cidx]);
3827 void build_layers_ph(
bool flat =
false) {
3828 if (!ph_laws_ready) ph_assert_supported(flat);
3832 opt.interlocking =
false;
3836 if (!ph_laws_ready) ph_init_laws();
3839 const std::size_t N = lqn.nidx;
3840 residt.assign(N + 1, Tzero());
3841 servt.assign(N + 1, Tzero());
3842 callservt.assign(lqn.ncalls + 1, Tzero());
3843 callresidt.assign(lqn.ncalls + 1, Tzero());
3844 for (std::size_t aidx = lqn.ashift + 1; aidx <= lqn.ashift + lqn.nacts; ++aidx)
3845 residt[aidx] = lqn.hostdem[aidx].disabled ? Tzero() : lqn.hostdem[aidx].mean;
3846 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
3847 if (lqn.calltype[cidx] != CallType::SYNC && lqn.calltype[cidx] != CallType::ASYNC)
3849 const std::size_t eidx = lqn.callpair_dst[cidx];
3850 callservt[cidx] = T(lqn.callproc_mean[cidx] * ph_hostmean[eidx]);
3851 callresidt[cidx] = callservt[cidx];
3853 ph_compose_entry_laws();
3857 build_ph_flat_layer();
3858 tput.assign(N + 1, Tzero());
3859 util.assign(N + 1, Tzero());
3860 thinkt.assign(N + 1, Tzero());
3861 update_layers_ph(0);
3865 std::vector<qn::Layer<T>> raw(NT() + 1);
3866 std::vector<bool> present(NT() + 1,
false);
3868 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx) {
3869 if (ignore[hidx])
continue;
3870 std::vector<std::size_t> callers;
3871 for (std::size_t tidx : lqn.tasksof[hidx]) {
3872 if (ignore[tidx])
continue;
3873 if (lqn.isref[tidx]) { callers.push_back(tidx);
continue; }
3874 for (std::size_t eidx : lqn.entriesof[tidx])
3875 if (ph_any_caller_of(eidx) || lqn.has_arrival[eidx]) {
3876 callers.push_back(tidx);
3880 if (callers.empty())
continue;
3881 build_layer_ph(raw[hidx], hidx, callers,
true);
3882 present[hidx] =
true;
3884 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
3885 const std::size_t tidx = lqn.tshift + t;
3886 if (ignore[tidx] || lqn.isref[tidx])
continue;
3887 std::vector<std::size_t> callers;
3888 for (std::size_t ct = 1; ct <= lqn.ntasks; ++ct) {
3889 const std::size_t c = lqn.tshift + ct;
3890 if (c == tidx || ignore[c])
continue;
3891 for (std::size_t e : lqn.entriesof[tidx])
3892 if (lqn.issynccaller.get(c, e)) { callers.push_back(c);
break; }
3894 if (callers.empty() && ph_async_calls_into(tidx).empty())
continue;
3895 build_layer_ph(raw[tidx], tidx, callers,
false);
3896 present[tidx] =
true;
3899 idxhash.assign(lqn.nidx + 1, -1);
3901 for (std::size_t i = 1; i <= NT(); ++i)
3903 idxhash[i] = next++;
3904 ensemble.push_back(std::move(raw[i]));
3906 layer_init_sol.assign(ensemble.size(), Matrix<T>());
3911 tput.assign(N + 1, Tzero());
3912 util.assign(N + 1, Tzero());
3913 thinkt.assign(N + 1, Tzero());
3914 update_layers_ph(0);
3930 void build_ph_flat_layer() {
3931 const std::vector<std::size_t> servers = ph_flat_server_set();
3932 const std::size_t nsrv = servers.size();
3936 m.clientIdx = m.add_station(qn::Station<T>{
3937 "Clients", NodeType::Delay, SchedStrategy::INF,
3938 std::numeric_limits<double>::infinity(),
false, 0});
3939 const std::size_t clientNode = m.node_of_station(m.clientIdx);
3940 m.server_idx_of.assign(lqn.nidx + 1, 0);
3942 std::vector<std::size_t> station_of(lqn.nidx + 1, 0);
3943 std::vector<std::size_t> serverNode(nsrv, 0);
3944 std::vector<std::size_t> srvStation(nsrv, 0);
3945 for (std::size_t si = 0; si < nsrv; ++si) {
3946 const std::size_t idx = servers[si];
3947 const bool ishost = idx <= lqn.nhosts;
3949 st.name = lqn.hashnames[idx];
3950 st.nodetype = lqn.sched[idx] == SchedStrategy::INF ? NodeType::Delay : NodeType::Queue;
3951 st.sched = lqn.sched[idx];
3952 st.nservers = lqn.sched[idx] == SchedStrategy::INF
3953 ? std::numeric_limits<double>::infinity()
3955 st.attr_ishost = ishost;
3957 const std::size_t s = m.add_station(st);
3959 station_of[idx] = s;
3960 serverNode[si] = m.node_of_station(s);
3961 m.server_idx_of[idx] = s;
3963 m.host_stations.push_back(s);
3965 m.task_stations.push_back(s);
3968 m.serverIdx = station_of[servers[0]];
3972 std::vector<std::vector<std::size_t>> callers_of(lqn.nidx + 1);
3973 std::vector<std::size_t> all_callers;
3974 for (std::size_t si = 0; si < nsrv; ++si) {
3975 const std::size_t idx = servers[si];
3976 std::vector<std::size_t> cs;
3977 if (idx <= lqn.nhosts) {
3978 for (std::size_t tidx : lqn.tasksof[idx]) {
3979 if (ignore[tidx])
continue;
3980 if (lqn.isref[tidx]) { cs.push_back(tidx);
continue; }
3981 for (std::size_t eidx : lqn.entriesof[tidx])
3982 if (ph_any_caller_of(eidx) || lqn.has_arrival[eidx]) {
3988 for (std::size_t ct = 1; ct <= lqn.ntasks; ++ct) {
3989 const std::size_t c = lqn.tshift + ct;
3990 if (c == idx || ignore[c])
continue;
3991 for (std::size_t e : lqn.entriesof[idx])
3992 if (lqn.issynccaller.get(c, e)) { cs.push_back(c);
break; }
3995 callers_of[idx] = cs;
3996 for (std::size_t c : cs)
3997 if (std::find(all_callers.begin(), all_callers.end(), c) == all_callers.end())
3998 all_callers.push_back(c);
4000 std::sort(all_callers.begin(), all_callers.end());
4003 std::vector<std::size_t> class_of_caller(lqn.nidx + 1, 0);
4005 for (std::size_t c : all_callers) {
4008 double nj = lqn.maxmult[c];
4009 if (std::isinf(nj)) {
4011 for (std::size_t k = 1; k <= NT(); ++k)
4012 if (lqn.taskgraph.get(k, c) != Tzero()) sacc += lqn.maxmult[k];
4014 if (std::isinf(nj) || nj == 0.0) {
4016 for (std::size_t k = 1; k <= NT(); ++k)
4017 if (std::isfinite(lqn.maxmult[k])) s2 += lqn.maxmult[k];
4018 nj = std::min(s2, 1000.0);
4022 jc.name = lqn.hashnames[c];
4023 jc.type = JobClassType::CLOSED;
4025 jc.refstat = m.clientIdx;
4026 jc.is_ref_class =
true;
4027 jc.attr_kind = int(LqnElement::TASK);
4029 const std::size_t k = m.add_class(jc);
4030 class_of_caller[c] = k;
4031 m.attr_tasks.emplace_back(k, c);
4033 const double zt = dbl(ref_think_time(c));
4034 m.set_service(m.clientIdx, k,
4043 for (std::size_t si = 0; si < nsrv; ++si)
4045 for (std::size_t si = 0; si < nsrv; ++si) {
4046 const std::size_t idx = servers[si];
4047 const std::vector<std::size_t>& cs = callers_of[idx];
4048 if (std::find(cs.begin(), cs.end(), c) == cs.end())
continue;
4049 m.set_service(srvStation[si], k,
4053 thinkt_map.push_back({idx, c, m.clientIdx, k});
4054 servt_map.push_back({idx, c, station_of[idx], k});
4059 std::vector<std::vector<std::pair<std::size_t, long>>> open_of(lqn.nidx + 1);
4060 std::size_t sourceStation = 0, sinkNode = 0;
4061 for (std::size_t si = 0; si < nsrv; ++si) {
4062 const std::size_t hidx = servers[si];
4063 if (hidx > lqn.nhosts)
continue;
4064 for (std::size_t c : callers_of[hidx]) {
4070 if (ph_open_arrival_only(c))
continue;
4071 for (std::size_t eidx : lqn.entriesof[c]) {
4072 if (!lqn.has_arrival[eidx])
continue;
4073 if (sourceStation == 0) {
4074 sourceStation = m.add_station(qn::Station<T>{
4075 "Source", NodeType::Source, SchedStrategy::EXT,
4076 std::numeric_limits<double>::infinity(),
false, 0});
4077 m.sourceIdx = sourceStation;
4078 sinkNode = m.add_node(
"Sink", NodeType::Sink,
false);
4079 m.sinkNode = sinkNode;
4082 oc.name = lqn.hashnames[eidx] +
".Open";
4083 oc.type = JobClassType::OPEN;
4084 oc.population = std::numeric_limits<double>::infinity();
4085 oc.refstat = sourceStation;
4086 oc.attr_kind = int(LqnElement::ENTRY);
4088 const std::size_t k = m.add_class(oc);
4089 m.set_service(sourceStation, k, lqn.arrival[eidx]);
4091 for (std::size_t s2 = 0; s2 < nsrv; ++s2)
4094 m.set_service(srvStation[si], k,
4096 open_of[hidx].emplace_back(k,
long(eidx));
4097 m.attr_entries.emplace_back(k, eidx);
4101 for (std::size_t si = 0; si < nsrv; ++si) {
4102 const std::size_t tidx = servers[si];
4103 if (tidx <= lqn.nhosts)
continue;
4104 for (std::size_t cidx : ph_async_calls_into(tidx)) {
4105 if (sourceStation == 0) {
4106 sourceStation = m.add_station(qn::Station<T>{
4107 "Source", NodeType::Source, SchedStrategy::EXT,
4108 std::numeric_limits<double>::infinity(),
false, 0});
4109 m.sourceIdx = sourceStation;
4110 sinkNode = m.add_node(
"Sink", NodeType::Sink,
false);
4111 m.sinkNode = sinkNode;
4113 const std::size_t eidx = lqn.callpair_dst[cidx];
4115 oc.name = lqn.callhashnames[cidx];
4116 oc.type = JobClassType::OPEN;
4117 oc.population = std::numeric_limits<double>::infinity();
4118 oc.refstat = sourceStation;
4119 oc.attr_kind = int(LqnElement::CALL);
4121 const std::size_t k = m.add_class(oc);
4124 for (std::size_t s2 = 0; s2 < nsrv; ++s2)
4127 m.set_service(srvStation[si], k,
4129 open_of[tidx].emplace_back(k, -
long(cidx));
4130 m.attr_calls.push_back({k, cidx, lqn.callpair_src[cidx], eidx});
4131 arv_call_map.push_back({tidx, cidx, sourceStation, k});
4132 call_map.push_back({tidx, cidx, station_of[tidx], k});
4139 for (std::size_t c : all_callers) {
4140 const std::size_t k = class_of_caller[c];
4141 std::size_t prev = clientNode;
4142 bool visited =
false;
4143 for (std::size_t si = 0; si < nsrv; ++si) {
4144 const std::vector<std::size_t>& cs = callers_of[servers[si]];
4145 if (std::find(cs.begin(), cs.end(), c) == cs.end())
continue;
4146 m.set_route(k, k, prev, serverNode[si], Tone());
4147 prev = serverNode[si];
4150 if (visited) m.set_route(k, k, prev, clientNode, Tone());
4152 if (sourceStation != 0) {
4153 const std::size_t srcNode = m.node_of_station(sourceStation);
4154 for (std::size_t si = 0; si < nsrv; ++si)
4155 for (
const auto& oa : open_of[servers[si]]) {
4156 m.set_route(oa.first, oa.first, srcNode, serverNode[si], Tone());
4157 m.set_route(oa.first, oa.first, serverNode[si], sinkNode, Tone());
4162 idxhash.assign(lqn.nidx + 1, -1);
4163 for (std::size_t idx : servers) idxhash[idx] = 0;
4165 ensemble.push_back(std::move(m));
4166 layer_init_sol.assign(ensemble.size(), Matrix<T>());
4169 const std::size_t nclasses = ensemble[0].classes.size();
4170 for (std::size_t idx : servers) {
4173 L.ishost = idx <= lqn.nhosts;
4174 L.callers = callers_of[idx];
4175 L.class_of_caller = class_of_caller;
4177 L.qstations.assign(1, station_of[idx]);
4178 L.svcmean_by_class.assign(nclasses + 1, Tzero());
4179 L.open_arrivals = open_of[idx];
4180 L.npop = npop < 1.0 ? 1.0 : npop;
4182 ph_has_layer[idx] =
true;
4195 std::vector<std::size_t> ph_flat_server_set()
const {
4196 for (std::size_t i = 1; i <= NT(); ++i) {
4197 if (lqn.repl[i] > 1.0)
4199 "method='flat.ph' does not support replicated processors or tasks, whose "
4200 "replicas need a submodel each. Use method='srvn.ph'.");
4203 "method='flat.ph' does not support cache tasks. Use method='default'.");
4204 if (lqn.hassetup[i])
4206 "method='flat.ph' does not support setup tasks, whose powered-down threads "
4207 "are per-layer state. Use method='srvn.ph'.");
4209 std::vector<std::size_t> servers;
4210 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx) {
4211 if (ignore[hidx] || lqn.tasksof[hidx].empty())
continue;
4213 for (std::size_t tidx : lqn.tasksof[hidx]) {
4214 if (ignore[tidx])
continue;
4215 if (lqn.isref[tidx]) { any =
true;
break; }
4216 for (std::size_t eidx : lqn.entriesof[tidx])
4217 if (ph_any_caller_of(eidx) || lqn.has_arrival[eidx]) { any =
true;
break; }
4220 if (any) servers.push_back(hidx);
4222 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
4223 const std::size_t tidx = lqn.tshift + t;
4224 if (ignore[tidx] || lqn.isref[tidx])
continue;
4225 bool has_caller =
false;
4226 for (std::size_t ct = 1; ct <= lqn.ntasks && !has_caller; ++ct) {
4227 const std::size_t c = lqn.tshift + ct;
4228 if (c == tidx || ignore[c])
continue;
4229 for (std::size_t e : lqn.entriesof[tidx])
4230 if (lqn.issynccaller.get(c, e)) { has_caller =
true;
break; }
4232 if (!has_caller && ph_async_calls_into(tidx).empty())
continue;
4233 servers.push_back(tidx);
4235 if (servers.empty())
4237 "method='flat.ph' found no server: the model has no processor with tasks.");
4242 void build_layer_ph(qn::Layer<T>& m, std::size_t idx,
4243 const std::vector<std::size_t>& callers,
bool ishost) {
4244 m.name = lqn.hashnames[idx];
4249 const double rawrepl = lqn.repl[idx];
4250 std::size_t nreplicas = 1;
4251 if (rawrepl > 1.0 && !callers.empty()) {
4254 for (std::size_t c : callers)
4255 if (lqn.repl[c] != rawrepl) reduce =
false;
4257 for (std::size_t c : callers)
4258 if (lqn.fanout_at(c, idx) < rawrepl) reduce =
false;
4260 nreplicas = reduce ? 1 :
static_cast<std::size_t
>(std::llround(rawrepl));
4261 if (reduce && !ishost) single_replica_tasks.insert(idx);
4263 const bool reduce_fanout = (nreplicas == 1 && rawrepl > 1.0 && !callers.empty());
4265 m.clientIdx = m.add_station(qn::Station<T>{
4266 "Clients", NodeType::Delay, SchedStrategy::INF,
4267 std::numeric_limits<double>::infinity(),
false, 0});
4268 const std::size_t clientNode = m.node_of_station(m.clientIdx);
4269 m.serverIdx = m.clientIdx + 1;
4273 L.callers = callers;
4274 L.nreplicas = nreplicas;
4275 L.class_of_caller.assign(lqn.nidx + 1, 0);
4276 std::vector<std::size_t> serverNode(nreplicas);
4277 for (std::size_t r = 0; r < nreplicas; ++r) {
4279 st.name = r == 0 ? lqn.hashnames[idx] : lqn.hashnames[idx] +
"." + std::to_string(r + 1);
4280 st.nodetype = lqn.sched[idx] == SchedStrategy::INF ? NodeType::Delay : NodeType::Queue;
4281 st.sched = lqn.sched[idx];
4282 st.nservers = lqn.sched[idx] == SchedStrategy::INF
4283 ? std::numeric_limits<double>::infinity()
4285 st.attr_ishost = ishost;
4287 const std::size_t s = m.add_station(st);
4288 L.qstations.push_back(s);
4289 serverNode[r] = m.node_of_station(s);
4293 for (std::size_t c : callers) {
4294 double nj = njobs(c, idx);
4296 const bool caller_single_replica =
4297 reduce_fanout || single_replica_tasks.count(c) > 0;
4298 nj = caller_single_replica ? lqn.maxmult[c] : lqn.maxmult[c] * lqn.repl[c];
4299 if (std::isinf(nj)) {
4301 for (std::size_t k = 1; k <= NT(); ++k)
4302 if (lqn.taskgraph.get(k, c) != Tzero()) s += lqn.maxmult[k];
4304 if (std::isinf(nj) || nj == 0.0) {
4306 for (std::size_t k = 1; k <= NT(); ++k)
4307 if (std::isfinite(lqn.maxmult[k])) s2 += lqn.maxmult[k] * lqn.repl[k];
4308 nj = std::min(s2, 1000.0);
4314 jc.name = lqn.hashnames[c];
4315 jc.type = JobClassType::CLOSED;
4317 jc.refstat = m.clientIdx;
4318 jc.is_ref_class =
true;
4319 jc.attr_kind = int(LqnElement::TASK);
4321 const std::size_t k = m.add_class(jc);
4322 L.class_of_caller[c] = k;
4323 m.attr_tasks.emplace_back(k, c);
4324 const double zt = dbl(ref_think_time(c));
4325 m.set_service(m.clientIdx, k,
4328 for (std::size_t s : L.qstations)
4333 thinkt_map.push_back({idx, c, m.clientIdx, k});
4334 servt_map.push_back({idx, c, L.qstations[0], k});
4338 std::size_t sourceStation = 0, sinkNode = 0;
4340 for (std::size_t c : callers) {
4347 if (ph_open_arrival_only(c))
continue;
4348 for (std::size_t eidx : lqn.entriesof[c]) {
4349 if (!lqn.has_arrival[eidx])
continue;
4350 if (sourceStation == 0) {
4351 sourceStation = m.add_station(qn::Station<T>{
4352 "Source", NodeType::Source, SchedStrategy::EXT,
4353 std::numeric_limits<double>::infinity(),
false, 0});
4354 m.sourceIdx = sourceStation;
4355 sinkNode = m.add_node(
"Sink", NodeType::Sink,
false);
4356 m.sinkNode = sinkNode;
4359 oc.name = lqn.hashnames[eidx] +
".Open";
4360 oc.type = JobClassType::OPEN;
4361 oc.population = std::numeric_limits<double>::infinity();
4362 oc.refstat = sourceStation;
4363 oc.attr_kind = int(LqnElement::ENTRY);
4365 const std::size_t k = m.add_class(oc);
4366 m.set_service(sourceStation, k, lqn.arrival[eidx]);
4367 for (std::size_t s : L.qstations) {
4371 L.open_arrivals.emplace_back(k,
long(eidx));
4372 m.attr_entries.emplace_back(k, eidx);
4376 for (std::size_t cidx : ph_async_calls_into(idx)) {
4377 if (sourceStation == 0) {
4378 sourceStation = m.add_station(qn::Station<T>{
4379 "Source", NodeType::Source, SchedStrategy::EXT,
4380 std::numeric_limits<double>::infinity(),
false, 0});
4381 m.sourceIdx = sourceStation;
4382 sinkNode = m.add_node(
"Sink", NodeType::Sink,
false);
4383 m.sinkNode = sinkNode;
4385 const std::size_t eidx = lqn.callpair_dst[cidx];
4387 oc.name = lqn.callhashnames[cidx];
4388 oc.type = JobClassType::OPEN;
4389 oc.population = std::numeric_limits<double>::infinity();
4390 oc.refstat = sourceStation;
4391 oc.attr_kind = int(LqnElement::CALL);
4393 const std::size_t k = m.add_class(oc);
4395 for (std::size_t s : L.qstations) {
4399 L.open_arrivals.emplace_back(k, -
long(cidx));
4400 m.attr_calls.push_back({k, cidx, lqn.callpair_src[cidx], eidx});
4401 arv_call_map.push_back({idx, cidx, sourceStation, k});
4402 call_map.push_back({idx, cidx, L.qstations[0], k});
4408 const T share = num_traits<T>::from_double(1.0 /
double(nreplicas));
4409 for (std::size_t c : callers) {
4410 const std::size_t k = L.class_of_caller[c];
4411 for (std::size_t r = 0; r < nreplicas; ++r) {
4412 m.set_route(k, k, clientNode, serverNode[r], share);
4413 m.set_route(k, k, serverNode[r], clientNode, Tone());
4416 for (
const auto& oa : L.open_arrivals) {
4417 const std::size_t k = oa.first;
4418 const std::size_t srcNode = m.node_of_station(sourceStation);
4419 for (std::size_t r = 0; r < nreplicas; ++r) {
4420 m.set_route(k, k, srcNode, serverNode[r], share);
4421 m.set_route(k, k, serverNode[r], sinkNode, Tone());
4424 L.svcmean_by_class.assign(m.classes.size() + 1, Tzero());
4426 for (std::size_t c : callers) {
4427 const double v = njobs(c, idx);
4428 if (std::isfinite(v) && v > 0.0) np += v;
4430 L.npop = np < 1.0 ? 1.0 : np;
4433 ph_has_layer[idx] =
true;
4437 workflow::PhLaw<T> ph_service_law(std::size_t idx,
bool ishost, std::size_t c)
const {
4440 std::vector<workflow::PhLaw<T>> laws;
4441 std::vector<T> probs;
4443 for (std::size_t eidx : lqn.entriesof[c]) {
4444 if (ph_hostlaw[eidx].S.rows() == 0 || !(dbl(ph_share[eidx]) > 0.0))
continue;
4445 laws.push_back(ph_hostlaw[eidx]);
4446 probs.push_back(ph_share[eidx]);
4447 tot += dbl(ph_share[eidx]);
4449 if (laws.empty())
return ph_immediate_law();
4450 for (T& p : probs) p = T(p / num_traits<T>::from_double(tot));
4456 bool started =
false;
4457 workflow::PhLaw<T> out;
4458 for (std::size_t eidx : lqn.entriesof[idx]) {
4459 const T n = ph_ncalls(c, eidx);
4461 const workflow::PhLaw<T> lp =
4466 if (!started)
return ph_immediate_law();
4474 T ph_delay_mean(std::size_t idx, std::size_t c)
const {
4484 T z = T(thinkt[c] + ref_think_time(c));
4485 if (!std::isfinite(dbl(z)) || dbl(z) < 0.0) z = Tzero();
4486 z = T(z + ph_actthinkt[c]);
4487 const std::size_t hidx = lqn.parent[c];
4488 if (!ph_served_here(idx, hidx)) z = T(z + ph_procresid[c]);
4489 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
4490 const std::size_t tidx = lqn.tshift + t;
4491 if (!ph_served_here(idx, tidx)) z = T(z + ph_calltime(c, tidx));
4493 if (!std::isfinite(dbl(z)) || dbl(z) < 0.0)
4503 bool ph_served_here(std::size_t idx, std::size_t elem)
const {
4504 if (elem < 1 || elem >= idxhash.size() || idx < 1 || idx >= idxhash.size())
4506 return idxhash[elem] >= 0 && idxhash[elem] == idxhash[idx];
4518 void update_layers_ph(
int) {
4519 for (std::size_t idx = 1; idx <= NT(); ++idx) {
4520 if (idxhash[idx] < 0 || !ph_has_layer[idx])
continue;
4521 PHLayer& L = ph_layers[idx];
4522 qn::Layer<T>& m = ensemble[std::size_t(idxhash[idx])];
4523 for (std::size_t c : L.callers) {
4524 const std::size_t k = L.class_of_caller[c];
4525 const workflow::PhLaw<T> sl = ph_service_law(idx, L.ishost, c);
4527 const Distrib<T> law = ph_station_law(sl);
4528 for (std::size_t s : L.qstations) m.set_service(s, k, law);
4529 const double zd = std::max(dbl(ph_delay_mean(idx, c)),
4531 m.set_service(m.clientIdx, k,
4534 for (
const auto& oa : L.open_arrivals) {
4535 const std::size_t k = oa.first;
4536 if (oa.second > 0) {
4538 L.svcmean_by_class[k] = ph_hostmean[std::size_t(oa.second)];
4541 const std::size_t cidx = std::size_t(-oa.second);
4542 const std::size_t eidx = lqn.callpair_dst[cidx];
4543 L.svcmean_by_class[k] = ph_entrymean[eidx];
4544 const Distrib<T> law = ph_station_law(ph_entrylaw[eidx]);
4545 for (std::size_t s : L.qstations) m.set_service(s, k, law);
4546 const std::size_t aidx = lqn.callpair_src[cidx];
4547 double rate = dbl(tput[aidx]) * dbl(lqn.callproc_mean[cidx]);
4550 m.set_service(m.sourceIdx, k,
4563 static T ph_residence(
const T& Q,
const T& X,
const T& RN) {
4564 const double q = num_traits<T>::to_double(Q), x = num_traits<T>::to_double(X);
4575 static T ph_inflation_of(
const T& R,
const T& S,
double npop) {
4577 const double s = num_traits<T>::to_double(S), r = num_traits<T>::to_double(R);
4579 if (!std::isfinite(f) || f < 1.0) f = 1.0;
4580 if (std::isfinite(npop) && npop >= 1.0 && f > npop) f = npop;
4581 return num_traits<T>::from_double(f);
4585 double ph_layer_pop(
const PHLayer& L, std::size_t idx)
const {
4588 if (L.npop >= 1.0)
return L.npop;
4590 for (std::size_t c : L.callers) {
4591 const double v = njobs(c, idx);
4592 if (std::isfinite(v) && v > 0.0) n += v;
4594 return n < 1.0 ? 1.0 : n;
4610 void update_metrics_ph(
int it) {
4611 const std::size_t N = lqn.nidx;
4612 servt.assign(N + 1, Tzero());
4613 residt.assign(N + 1, Tzero());
4614 callservt.assign(lqn.ncalls + 1, Tzero());
4615 callresidt.assign(lqn.ncalls + 1, Tzero());
4617 std::vector<T> inflNum(N + 1, Tzero()), inflDen(N + 1, Tzero());
4618 std::vector<T> taskTput(N + 1, Tzero()), openTput(N + 1, Tzero());
4621 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx) {
4622 if (idxhash[hidx] < 0 || !ph_has_layer[hidx])
continue;
4623 const PHLayer& L = ph_layers[hidx];
4624 const LayerResult<T>& res = results.back()[std::size_t(idxhash[hidx])];
4625 const double npop = ph_layer_pop(L, hidx);
4626 for (std::size_t c : L.callers) {
4627 const std::size_t k = L.class_of_caller[c];
4628 T X = Tzero(), Q = Tzero();
4629 for (std::size_t s : L.qstations) {
4630 X = T(X + res.TN(s - 1, k - 1));
4631 Q = T(Q + res.QN(s - 1, k - 1));
4633 const T R = ph_residence(Q, X, res.RN(L.qstations[0] - 1, k - 1));
4634 const T f = ph_inflation_of(R, L.svcmean_by_class[k], npop);
4635 if (!std::isfinite(dbl(X)) || dbl(X) < 0.0) X = Tzero();
4641 taskTput[c] = T(taskTput[c] + num_traits<T>::from_double(lqn.repl[c]) * X);
4642 for (std::size_t eidx : lqn.entriesof[c]) {
4643 const T sh = dbl(ph_share[eidx]) > 0.0 ? ph_share[eidx] : Tzero();
4644 const T w = T(sh * X);
4645 inflNum[eidx] = T(inflNum[eidx] + w * f);
4646 inflDen[eidx] = T(inflDen[eidx] + w);
4649 for (
const auto& oa : L.open_arrivals) {
4650 if (oa.second <= 0)
continue;
4651 const std::size_t k = oa.first;
4652 const std::size_t eidx = std::size_t(oa.second);
4653 T X = Tzero(), Q = Tzero();
4654 for (std::size_t s : L.qstations) {
4655 X = T(X + res.TN(s - 1, k - 1));
4656 Q = T(Q + res.QN(s - 1, k - 1));
4658 if (!std::isfinite(dbl(X)) || dbl(X) <= 0.0)
continue;
4659 const T f = ph_inflation_of(ph_residence(Q, X, res.RN(L.qstations[0] - 1, k - 1)),
4660 L.svcmean_by_class[k], npop);
4661 inflNum[eidx] = T(inflNum[eidx] + X * f);
4662 inflDen[eidx] = T(inflDen[eidx] + X);
4663 openTput[eidx] = T(openTput[eidx] + X);
4664 taskTput[lqn.parent[eidx]] = T(taskTput[lqn.parent[eidx]] + X);
4668 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
4669 const std::size_t eidx = lqn.eshift + e;
4672 f = dbl(inflNum[eidx]) / dbl(inflDen[eidx]);
4673 if (!std::isfinite(f) || f < 1.0) f = 1.0;
4674 const T fv = num_traits<T>::from_double(f);
4675 for (std::size_t aidx : lqn.actsof[eidx])
4676 residt[aidx] = T(fv * (lqn.hostdem[aidx].disabled ? Tzero()
4677 : lqn.hostdem[aidx].mean));
4681 std::vector<T> relw(N + 1, Tzero());
4682 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
4683 const std::size_t tidx = lqn.tshift + t;
4684 if (idxhash[tidx] < 0 || !ph_has_layer[tidx])
continue;
4685 const PHLayer& L = ph_layers[tidx];
4686 const LayerResult<T>& res = results.back()[std::size_t(idxhash[tidx])];
4687 const double npop = ph_layer_pop(L, tidx);
4688 for (std::size_t c : L.callers) {
4689 const std::size_t k = L.class_of_caller[c];
4690 T X = Tzero(), Q = Tzero();
4691 for (std::size_t s : L.qstations) {
4692 X = T(X + res.TN(s - 1, k - 1));
4693 Q = T(Q + res.QN(s - 1, k - 1));
4695 if (!std::isfinite(dbl(X)) || dbl(X) < 0.0) X = Tzero();
4696 const T g = ph_inflation_of(ph_residence(Q, X, res.RN(L.qstations[0] - 1, k - 1)),
4697 L.svcmean_by_class[k], npop);
4698 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
4699 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
4700 if (lqn.parent[lqn.callpair_src[cidx]] != c)
continue;
4701 if (lqn.parent[lqn.callpair_dst[cidx]] != tidx)
continue;
4702 const std::size_t eidx = lqn.callpair_dst[cidx];
4703 callservt[cidx] = T(lqn.callproc_mean[cidx] * g * ph_entrymean[eidx]);
4704 callresidt[cidx] = callservt[cidx];
4706 for (std::size_t eidx : lqn.entriesof[tidx])
4707 relw[eidx] = T(relw[eidx] + X * ph_ncalls(c, eidx));
4709 for (
const auto& oa : L.open_arrivals) {
4710 if (oa.second >= 0)
continue;
4711 const std::size_t k = oa.first;
4712 const std::size_t cidx = std::size_t(-oa.second);
4713 const std::size_t eidx = lqn.callpair_dst[cidx];
4715 for (std::size_t s : L.qstations) X = T(X + res.TN(s - 1, k - 1));
4716 const T R = res.RN(L.qstations[0] - 1, k - 1);
4717 if (std::isfinite(dbl(R)) && dbl(R) > 0.0) {
4718 callservt[cidx] = T(R * lqn.callproc_mean[cidx]);
4719 callresidt[cidx] = callservt[cidx];
4721 if (std::isfinite(dbl(X)) && dbl(X) > 0.0) relw[eidx] = T(relw[eidx] + X);
4726 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
4727 const std::size_t tidx = lqn.tshift + t;
4728 const std::vector<std::size_t>& entries = lqn.entriesof[tidx];
4729 if (entries.empty())
continue;
4736 for (std::size_t eidx : entries) tot = T(tot + relw[eidx] + openTput[eidx]);
4738 for (std::size_t eidx : entries)
4739 ph_share[eidx] = T((relw[eidx] + openTput[eidx]) / tot);
4741 for (std::size_t eidx : entries)
4742 ph_share[eidx] = num_traits<T>::from_double(1.0 /
double(entries.size()));
4750 for (std::size_t eidx : entries) tput[eidx] = T(tput[tidx] * ph_share[eidx]);
4756 const T nrep = num_traits<T>::from_double(std::max(1.0, lqn.repl[tidx]));
4758 : T(tput[tidx] / nrep);
4762 for (std::size_t aidx = lqn.ashift + 1; aidx <= lqn.ashift + lqn.nacts; ++aidx) {
4764 if (!std::isfinite(dbl(v)) && it > 1 && !std::isnan(residt_prev[aidx]))
4765 v = residt_prev_v[aidx];
4766 if (relax_omega < 1.0 && it > 1 && !std::isnan(residt_prev[aidx])) {
4767 const T om = num_traits<T>::from_double(relax_omega);
4768 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
4769 v = T(om * v + om1 * residt_prev_v[aidx]);
4772 residt_prev[aidx] = dbl(v);
4773 residt_prev_v[aidx] = v;
4775 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
4776 T v = callservt[cidx];
4777 if (!std::isfinite(dbl(v)))
4778 v = (it > 1 && std::isfinite(callservt_prev[cidx])) ? callservt_prev_v[cidx]
4780 if (relax_omega < 1.0 && it > 1 && !std::isnan(callservt_prev[cidx])) {
4781 const T om = num_traits<T>::from_double(relax_omega);
4782 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
4783 v = T(om * v + om1 * callservt_prev_v[cidx]);
4785 callservt[cidx] = v;
4786 callresidt[cidx] = v;
4787 callservt_prev[cidx] = dbl(v);
4788 callresidt_prev[cidx] = dbl(v);
4789 callservt_prev_v[cidx] = v;
4797 ph_compose_entry_laws();
4799 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
4800 const std::size_t eidx = lqn.eshift + e;
4801 if (!ph_has_wf[eidx])
continue;
4802 const std::unordered_map<std::size_t, T>& ex = ph_execs[eidx];
4803 for (std::size_t aidx : lqn.actsof[eidx]) {
4804 T sa = T(residt[aidx] + act_think_time(aidx));
4805 for (std::size_t cidx : lqn.callsof[aidx])
4806 if (lqn.calltype[cidx] == CallType::SYNC) sa = T(sa + callservt[cidx]);
4808 servt_prev[aidx] = dbl(sa);
4809 servt_prev_v[aidx] = sa;
4810 tput[aidx] = T(tput[eidx] * ex.at(aidx));
4811 tput_prev[aidx] = dbl(tput[aidx]);
4812 tput_prev_v[aidx] = tput[aidx];
4815 : Distrib<T>::disabled_dist();
4818 servt[eidx] = ph_entrymean[eidx];
4819 residt[eidx] = ph_entrymean[eidx];
4835 void update_think_times_ph(
int it) {
4838 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
4839 const std::size_t tidx = lqn.tshift + t;
4840 if (ignore[tidx])
continue;
4841 const T ztask = ref_think_time(tidx);
4842 if (idxhash[tidx] < 0) {
4847 const double arvrate = open_arrival_rate_of(tidx);
4849 double nja = lqn.maxmult[tidx];
4850 if (!std::isfinite(nja) || nja <= 0.0) {
4852 for (std::size_t c = 1; c <= NT(); ++c) nja = std::max(nja, njobs(tidx, c));
4855 const std::size_t hidx = lqn.parent[tidx];
4856 if (hidx >= 1 && hidx < idxhash.size() && idxhash[hidx] >= 0 &&
4857 ph_has_layer[hidx] && ph_layers[hidx].class_of_caller[tidx] > 0) {
4858 const PHLayer& HL = ph_layers[hidx];
4859 const LayerResult<T>& hr = results.back()[std::size_t(idxhash[hidx])];
4860 const T rr = hr.RN(HL.qstations[0] - 1, HL.class_of_caller[tidx] - 1);
4861 if (!std::isnan(dbl(rr))) hres = rr;
4863 T za = T(num_traits<T>::from_double(nja / arvrate) - hres - ztask);
4864 if (za < floorv) za = floorv;
4865 if (relax_omega < 1.0 && it > 1 && !std::isnan(thinkt_prev[tidx])) {
4866 const T om = num_traits<T>::from_double(relax_omega);
4867 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
4868 za = T(om * za + om1 * thinkt_prev_v[tidx]);
4870 tput[tidx] = num_traits<T>::from_double(arvrate);
4872 thinkt_prev[tidx] = dbl(za);
4873 thinkt_prev_v[tidx] = za;
4883 const PHLayer& L = ph_layers[tidx];
4884 const qn::Layer<T>& m = ensemble[std::size_t(idxhash[tidx])];
4885 const LayerResult<T>& r = results.back()[std::size_t(idxhash[tidx])];
4887 for (std::size_t k = 0; k < m.nclasses; ++k) {
4888 const T v = r.UN(L.qstations[0] - 1, k);
4889 if (!std::isnan(dbl(v))) U = T(U + v);
4897 if (dbl(ph_setupshare[tidx]) > 0.0)
4898 U = T(U * (Tone() - ph_setupshare[tidx]));
4901 T X = ph_xdemand[tidx];
4904 double nj = lqn.maxmult[tidx];
4905 if (!std::isfinite(nj) || nj <= 0.0) {
4907 for (std::size_t c = 1; c <= NT(); ++c) nj = std::max(nj, njobs(tidx, c));
4911 if (lqn.sched[tidx] == SchedStrategy::INF) {
4913 z = T((num_traits<T>::from_double(nj) - U) / X - ztask);
4915 const T om = U > Tone() ? T(U - Tone()) : T(Tone() - U);
4916 z = T(num_traits<T>::from_double(nj) * om / X - ztask);
4921 if (z < floorv) z = floorv;
4922 if (it > 1 && !std::isnan(thinkt_prev[tidx]) && !std::isfinite(dbl(z)))
4923 z = thinkt_prev_v[tidx];
4924 if (relax_omega < 1.0 && it > 1 && !std::isnan(thinkt_prev[tidx])) {
4925 const T om = num_traits<T>::from_double(relax_omega);
4926 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
4927 z = T(om * z + om1 * thinkt_prev_v[tidx]);
4930 thinkt_prev[tidx] = dbl(z);
4931 thinkt_prev_v[tidx] = z;
4945 LnSolution<T> aggregate_ph() {
4946 const std::size_t N = lqn.nidx;
4948 out.QN.assign(N + 1, Tzero());
4949 out.UN.assign(N + 1, Tzero());
4950 out.RN.assign(N + 1, Tzero());
4951 out.TN.assign(N + 1, Tzero());
4952 out.AN.assign(N + 1, Tzero());
4953 out.WN.assign(N + 1, Tzero());
4954 out.defined_Q.assign(N + 1,
false);
4955 out.defined_U.assign(N + 1,
false);
4956 out.defined_R.assign(N + 1,
false);
4957 out.defined_T.assign(N + 1,
false);
4958 out.defined_A.assign(N + 1,
false);
4959 out.defined_W.assign(N + 1,
false);
4960 out.iterations = iterations_done;
4961 out.converged = did_converge;
4963 std::vector<T> PN(N + 1, Tzero()), UT(N + 1, Tzero());
4964 std::vector<bool> hasPN(N + 1,
false), hasUT(N + 1,
false);
4966 for (std::size_t a = 1; a <= lqn.nacts; ++a) {
4967 const std::size_t aidx = lqn.ashift + a;
4968 const std::size_t tidx = lqn.parent[aidx];
4969 if (ignore[tidx])
continue;
4970 const std::size_t hidx = lqn.parent[tidx];
4971 out.TN[aidx] = tput[aidx];
4972 out.defined_T[aidx] =
true;
4973 out.RN[aidx] = servt[aidx];
4974 out.defined_R[aidx] =
true;
4975 UT[aidx] = T(tput[aidx] * servt[aidx]);
4981 const T hd = lqn.hostdem[aidx].disabled ? Tzero() : lqn.hostdem[aidx].mean;
4982 PN[aidx] = T(tput[aidx] * hd / num_traits<T>::from_double(ph_host_servers(hidx)));
4984 PN[hidx] = T(PN[hidx] + PN[aidx]);
4988 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
4989 const std::size_t eidx = lqn.eshift + e;
4990 const std::size_t tidx = lqn.parent[eidx];
4991 if (ignore[tidx])
continue;
4992 out.TN[eidx] = tput[eidx];
4993 out.defined_T[eidx] =
true;
4994 out.RN[eidx] = servt[eidx];
4995 out.defined_R[eidx] =
true;
4996 UT[eidx] = T(tput[eidx] * servt[eidx]);
4998 for (std::size_t aidx : lqn.actsof[eidx]) {
4999 PN[eidx] = T(PN[eidx] + PN[aidx]);
5006 if (ph_has_wf[eidx]) {
5007 const std::unordered_map<std::size_t, T>& ex = ph_execs[eidx];
5008 for (std::size_t aidx : lqn.actsof[eidx]) {
5009 out.WN[aidx] = T(ph_share[eidx] * ex.at(aidx) * residt[aidx]);
5010 out.defined_W[aidx] =
true;
5013 UT[tidx] = T(UT[tidx] + UT[eidx]);
5017 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
5018 const std::size_t tidx = lqn.tshift + t;
5019 if (ignore[tidx])
continue;
5020 out.TN[tidx] = tput[tidx];
5021 out.defined_T[tidx] =
true;
5024 for (std::size_t aidx : lqn.actsof[tidx]) {
5025 PN[tidx] = T(PN[tidx] + PN[aidx]);
5027 if (out.defined_W[aidx]) { w = T(w + out.WN[aidx]); anyw =
true; }
5029 if (anyw) { out.WN[tidx] = w; out.defined_W[tidx] =
true; }
5032 for (std::size_t hidx = 1; hidx <= lqn.nhosts; ++hidx)
5033 out.defined_T[hidx] =
false;
5035 for (std::size_t idx = 1; idx <= N; ++idx) {
5036 out.QN[idx] = UT[idx];
5037 out.defined_Q[idx] = hasUT[idx];
5038 out.UN[idx] = PN[idx];
5039 out.defined_U[idx] = hasPN[idx];
5045 out.UN[idx] = Tzero();
5046 out.defined_U[idx] =
true;
5047 out.defined_A[idx] =
false;
5048 const bool host = lqn.type[idx] == LqnElement::HOST;
5049 const bool task = lqn.type[idx] == LqnElement::TASK;
5050 const bool entry = lqn.type[idx] == LqnElement::ENTRY;
5051 out.QN[idx] = Tzero();
5052 out.defined_Q[idx] = !host;
5053 out.RN[idx] = Tzero();
5054 out.defined_R[idx] = !host && !task;
5055 out.WN[idx] = Tzero();
5056 out.defined_W[idx] = !host && !entry;
5057 out.TN[idx] = Tzero();
5058 out.defined_T[idx] = !host;
5064 void update_metrics(
int it) {
5065 if (is_ph_encoding()) {
5066 update_metrics_ph(it);
5069 if (lnmethod ==
"moment3") {
5070 update_metrics_moment_based(it);
5073 update_metrics_default(it);
5076 void update_metrics_default(
int it) {
5077 const std::size_t N = lqn.nidx;
5078 servt.assign(N + 1, Tzero());
5079 residt.assign(N + 1, Tzero());
5080 const int iter_min = std::min(30,
int(std::ceil(opt.iter_max / 4.0)));
5081 const bool averaging = averagingstart >= 0 && it >= iter_min;
5082 const int wnd = averaging ? int(it - averagingstart + 1) : 1;
5084 for (
const UpdRow& row : servt_map) {
5085 const std::size_t e = std::size_t(idxhash[row.idx]);
5086 const qn::Layer<T>& L = ensemble[e];
5087 const std::size_t k = row.cls - 1;
5089 for (std::size_t cc = 0; cc < L.nchains; ++cc)
5090 if (L.chains[cc][k]) c = cc;
5091 const std::size_t refclass_c = L.refclass[c];
5092 const std::size_t refstat_k = L.classes[k].refstat;
5094 T sv = Tzero(), rs = Tzero(), tp = Tzero();
5095 const T wT = T(Tone() / num_traits<T>::from_int(wnd));
5096 for (
int w = 0; w < wnd; ++w) {
5097 const LayerResult<T>& r = results[results.size() - 1 - w][e];
5098 sv += r.RN(row.node - 1, k) * wT;
5099 const T TN_ref = (refclass_c > 0 && refstat_k > 0) ? r.TN(refstat_k - 1, refclass_c - 1) : Tzero();
5101 rs += r.QN(row.node - 1, k) / TN_ref * wT;
5103 rs += r.WN(row.node - 1, k) * wT;
5104 tp += r.TN(row.node - 1, k) * wT;
5106 servt[row.aidx] = sv;
5107 residt[row.aidx] = rs;
5108 tput[row.aidx] = tp;
5111 const Distrib<T>& at = lqn.actthink[row.aidx];
5113 servt[row.aidx] = T(servt[row.aidx] + at.mean);
5114 residt[row.aidx] = T(residt[row.aidx] + at.mean);
5122 if (async_only_activity(row.aidx)) residt[row.aidx] = servt[row.aidx];
5124 if (relax_omega < 1.0 && it > 1) {
5125 const T om = num_traits<T>::from_double(relax_omega);
5126 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
5127 if (!std::isnan(servt_prev[row.aidx]))
5128 servt[row.aidx] = T(om * servt[row.aidx] + om1 * servt_prev_v[row.aidx]);
5129 if (!std::isnan(residt_prev[row.aidx]))
5130 residt[row.aidx] = T(om * residt[row.aidx] + om1 * residt_prev_v[row.aidx]);
5131 if (!std::isnan(tput_prev[row.aidx]))
5132 tput[row.aidx] = T(om * tput[row.aidx] + om1 * tput_prev_v[row.aidx]);
5134 servt_prev[row.aidx] = dbl(servt[row.aidx]);
5135 residt_prev[row.aidx] = dbl(residt[row.aidx]);
5136 tput_prev[row.aidx] = dbl(tput[row.aidx]);
5137 servt_prev_v[row.aidx] = servt[row.aidx];
5138 residt_prev_v[row.aidx] = residt[row.aidx];
5139 tput_prev_v[row.aidx] = tput[row.aidx];
5141 if (servt[row.aidx] > Tzero() && dbl(servt[row.aidx]) <= 1e10)
5151 servt_ph1.assign(N + 1, Tzero());
5152 servt_ph2.assign(N + 1, Tzero());
5153 for (std::size_t a = 1; a <= lqn.nacts; ++a) {
5154 const std::size_t aidx = lqn.ashift + a;
5155 if (lqn.actphase[a] == 1)
5156 servt_ph1[aidx] = servt[aidx];
5158 servt_ph2[aidx] = servt[aidx];
5160 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5161 const std::size_t eidx = lqn.eshift + e;
5162 for (std::size_t aidx : lqn.actsof[eidx]) {
5163 if (aidx <= lqn.ashift)
continue;
5164 const std::size_t a = aidx - lqn.ashift;
5165 if (a > lqn.nacts)
continue;
5166 if (lqn.actphase[a] == 1)
5167 servt_ph1[eidx] = T(servt_ph1[eidx] + servt_ph1[aidx]);
5169 servt_ph2[eidx] = T(servt_ph2[eidx] + servt_ph2[aidx]);
5175 for (
const UpdRow& row : thinkt_map) {
5176 if (!tputproc[row.aidx].disabled)
continue;
5177 const std::size_t e = std::size_t(idxhash[row.idx]);
5179 const T wT = T(Tone() / num_traits<T>::from_int(wnd));
5180 for (
int w = 0; w < wnd; ++w)
5181 tp += results[results.size() - 1 - w][e].TN(row.node - 1, row.cls - 1) * wT;
5182 tput[row.aidx] = tp;
5187 callservt.assign(lqn.ncalls + 1, Tzero());
5188 callresidt.assign(lqn.ncalls + 1, Tzero());
5189 for (
const UpdRow& row : call_map) {
5190 if (row.node <= 1)
continue;
5191 const std::size_t e = std::size_t(idxhash[row.idx]);
5192 const LayerResult<T>& r = results.back()[e];
5193 callservt[row.aidx] =
5194 T(r.RN(row.node - 1, row.cls - 1) * lqn.callproc_mean[row.aidx]);
5200 const qn::Layer<T>& L = ensemble[e];
5201 const std::size_t k = row.cls - 1;
5203 for (std::size_t cc = 0; cc < L.nchains; ++cc)
5204 if (L.chains[cc][k]) c = cc;
5205 const std::size_t refclass_c = L.refclass[c];
5206 const std::size_t refstat_k = L.classes[k].refstat;
5207 const T TN_ref = (refclass_c > 0 && refstat_k > 0) ? r.TN(refstat_k - 1, refclass_c - 1) : Tzero();
5209 ? T(r.QN(row.node - 1, k) / TN_ref)
5210 : r.WN(row.node - 1, k);
5212 const T rw = region_wait(e, row.cls);
5214 callservt[row.aidx] = T(callservt[row.aidx] + rw);
5215 callresidt[row.aidx] = T(callresidt[row.aidx] + rw);
5217 if (relax_omega < 1.0 && it > 1 && !std::isnan(callservt_prev[row.aidx])) {
5218 const T om = num_traits<T>::from_double(relax_omega);
5219 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
5220 callservt[row.aidx] =
5221 T(om * callservt[row.aidx] + om1 * callservt_prev_v[row.aidx]);
5223 callservt_prev[row.aidx] = dbl(callservt[row.aidx]);
5224 callservt_prev_v[row.aidx] = callservt[row.aidx];
5225 callresidt_prev[row.aidx] = dbl(callresidt[row.aidx]);
5228 resolve_entry_service();
5238 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5239 const std::size_t eidx = lqn.eshift + e;
5240 const std::size_t tidx = lqn.parent[eidx];
5242 if (lqn.isref[tidx] || !lqn.issynccaller.any_col(eidx)) {
5243 residt[eidx] = servt[eidx];
5246 T entry_tput = Tzero();
5248 entry_tput = tput[eidx];
5250 entry_tput = tput[tidx];
5256 residt[eidx] = T(servt_ph1[eidx] + prOvertake[e] * servt_ph2[eidx]);
5260 for (
const UpdRow& row : call_map) {
5261 if (row.node <= 1)
continue;
5262 const std::size_t eidx = lqn.callpair_dst[row.aidx];
5265 for (
const UpdRow& row : call_map) {
5266 if (row.node <= 1)
continue;
5267 const std::size_t eidx = lqn.callpair_dst[row.aidx];
5269 callservt[row.aidx] = servt[eidx];
5270 callservtproc[row.aidx] = servtproc[eidx];
5271 }
else if (callservt[row.aidx] > Tzero()) {
5281 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
5282 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
5283 const std::size_t target = lqn.callpair_dst[cidx];
5284 if (target <= lqn.eshift || target > lqn.eshift + lqn.nentries)
continue;
5286 const T eff = residt[target];
5287 if (!(eff > Tzero()))
continue;
5288 const T w = T(eff * lqn.callproc_mean[cidx]);
5289 callservt[cidx] = w;
5290 callresidt[cidx] = w;
5310 static void cdf_moments(
const fluid::FluidPassage& c,
double& m1,
double& m2,
double& m3) {
5312 for (std::size_t i = 0; i + 1 < c.t.size(); ++i) {
5313 const double x = 0.5 * (c.t[i + 1] + c.t[i]);
5314 const double w = c.cdf[i + 1] - c.cdf[i];
5317 m3 += x * x * x * w;
5322 const std::vector<std::vector<fluid::FluidPassage>>& layer_cdf(std::size_t e) {
5323 if (cdf_repo[e].empty()) {
5330 ensemble[e].refresh_rt();
5332 cdf_repo[e] = detail::ln_fluid_cdf_respt(ensemble[e], opt.layer_fluid);
5333 }
catch (
const std::exception&) {
5340 cdf_repo[e] = layer_cdf_exp_fallback(e);
5348 std::vector<std::vector<fluid::FluidPassage>> layer_cdf_exp_fallback(std::size_t e) {
5349 const qn::NetworkStruct<T>& L = ensemble[e];
5350 std::vector<std::vector<fluid::FluidPassage>> out(
5351 L.nstations, std::vector<fluid::FluidPassage>(L.nclasses));
5352 if (results.empty() || e >= results.back().size())
return out;
5353 const Matrix<T>& RN = results.back()[e].RN;
5354 const std::size_t npts = 100;
5355 for (std::size_t i = 0; i < L.nstations; ++i) {
5356 if (L.stations[i].nodetype == qn::NodeType::Source)
continue;
5357 for (std::size_t r = 0; r < L.nclasses; ++r) {
5358 if (L.disabled[i][r])
continue;
5359 const double rn = (i < static_cast<std::size_t>(RN.rows()) &&
5360 r < static_cast<std::size_t>(RN.cols()))
5363 if (!(std::isfinite(rn) && rn > 0.0))
continue;
5364 fluid::FluidPassage& cell = out[i][r];
5365 cell.t.reserve(npts);
5366 cell.cdf.reserve(npts);
5367 for (std::size_t j = 0; j < npts; ++j) {
5369 0.001 + (0.999 - 0.001) *
static_cast<double>(j) / (npts - 1);
5370 cell.cdf.push_back(q);
5371 cell.t.push_back(-std::log(1.0 - q) * rn);
5392 int entry_visit_ratio(std::size_t eidx, T& ratio)
const {
5393 const std::size_t tidx = lqn.parent[eidx];
5394 const std::size_t hidx = lqn.parent[tidx];
5395 if (ignore[tidx] || ignore[hidx] || idxhash[hidx] < 0)
return 0;
5396 if (!lqn.issynccaller.any_col(eidx))
return 1;
5397 const std::size_t hl = std::size_t(idxhash[hidx]);
5398 const qn::Layer<T>& L = ensemble[hl];
5399 const LayerResult<T>& r = results.back()[hl];
5400 T task_tput = Tzero(), entry_tput = Tzero();
5401 for (
const auto& kv : L.attr_tasks)
5402 if (kv.second == tidx) task_tput += r.TN(L.clientIdx - 1, kv.first - 1);
5403 for (
const auto& kv : L.attr_entries)
5404 if (kv.second == eidx) entry_tput += r.TN(L.clientIdx - 1, kv.first - 1);
5406 ratio = T(task_tput / (entry_tput > floor ? entry_tput : floor));
5411 Matrix<T> entry_service_resolvent()
const {
5412 const std::size_t dim = lqn.nidx + lqn.ncalls;
5413 Matrix<T> A(dim + 1, dim + 1, Tzero());
5414 for (std::size_t i = 0; i <= dim; ++i) {
5416 for (std::size_t j = 0; j <= dim; ++j) A(i, j) = T(A(i, j) - servtmatrix(i, j));
5418 return ::line::inverse(A);
5444 void update_metrics_moment_based(
int it) {
5445 const std::size_t N = lqn.nidx;
5446 servt.assign(N + 1, Tzero());
5447 residt.assign(N + 1, Tzero());
5448 callservt.assign(lqn.ncalls + 1, Tzero());
5449 callresidt.assign(lqn.ncalls + 1, Tzero());
5452 for (
const UpdRow& row : servt_map) {
5453 const std::size_t e = std::size_t(idxhash[row.idx]);
5454 const qn::Layer<T>& L = ensemble[e];
5455 const LayerResult<T>& r = results.back()[e];
5456 const std::size_t k = row.cls - 1;
5458 for (std::size_t cc = 0; cc < L.nchains; ++cc)
5459 if (L.chains[cc][k]) c = cc;
5460 const std::size_t refclass_c = L.refclass[c];
5461 const std::size_t refstat_k = L.classes[k].refstat;
5462 const T TN_ref = (refclass_c > 0 && refstat_k > 0) ? r.TN(refstat_k - 1, refclass_c - 1) : Tzero();
5464 tput[row.aidx] = r.TN(row.node - 1, k);
5466 ? T(r.QN(row.node - 1, k) / TN_ref)
5467 : r.WN(row.node - 1, k);
5468 if (!hasconverged) {
5469 servt[row.aidx] = r.RN(row.node - 1, k);
5471 const Distrib<T>& at = lqn.actthink[row.aidx];
5473 servt[row.aidx] = T(servt[row.aidx] + at.mean);
5474 residt[row.aidx] = T(residt[row.aidx] + at.mean);
5479 if (async_only_activity(row.aidx)) residt[row.aidx] = servt[row.aidx];
5481 servtcdf[row.aidx] = layer_cdf(e)[row.node - 1][k];
5485 for (
const UpdRow& row : call_map) {
5486 if (row.node <= 1)
continue;
5487 const std::size_t e = std::size_t(idxhash[row.idx]);
5488 const LayerResult<T>& r = results.back()[e];
5489 callresidt[row.aidx] = r.WN(row.node - 1, row.cls - 1);
5491 callservt[row.aidx] =
5492 T(r.RN(row.node - 1, row.cls - 1) * lqn.callproc_mean[row.aidx]);
5494 callservtcdf[row.aidx] = layer_cdf(e)[row.node - 1][row.cls - 1];
5498 moment3_means_pass(it);
5500 moment3_distribution_pass(it);
5504 void moment3_means_pass(
int it) {
5505 const std::size_t dim = lqn.nidx + lqn.ncalls;
5506 std::vector<T> x(dim + 1, Tzero());
5507 for (std::size_t i = 1; i <= lqn.nidx; ++i) x[i] = residt[i];
5508 for (std::size_t c = 1; c <= lqn.ncalls; ++c) x[lqn.nidx + c] = callresidt[c];
5511 const Matrix<T> Rinv = entry_service_resolvent();
5512 std::vector<T> entry_servt(dim + 1, Tzero());
5513 for (std::size_t i = 1; i <= dim; ++i) {
5515 for (std::size_t j = 1; j <= dim; ++j)
5516 if (Rinv(i, j) != Tzero()) s += Rinv(i, j) * x[j];
5519 for (std::size_t i = 1; i <= lqn.eshift; ++i) entry_servt[i] = Tzero();
5530 for (std::size_t e = 1; e <= lqn.nentries; ++e)
5531 servt[lqn.eshift + e] = entry_servt[lqn.eshift + e];
5534 std::vector<T> entry_residt(dim + 1, Tzero());
5535 for (std::size_t i = lqn.eshift + 1; i <= lqn.eshift + lqn.nentries; ++i) {
5537 for (std::size_t j = 1; j <= dim; ++j)
5538 if (servtmatrix(i, j) != Tzero()) s += servtmatrix(i, j) * x[j];
5539 entry_residt[i] = s;
5542 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5543 const std::size_t eidx = lqn.eshift + e;
5545 const int state = entry_visit_ratio(eidx, ratio);
5546 if (state == 0)
continue;
5548 residt[eidx] = entry_residt[eidx];
5551 servt[eidx] = T(entry_servt[eidx] * ratio);
5552 residt[eidx] = T(entry_residt[eidx] * ratio);
5555 for (
const UpdRow& row : call_map) {
5556 if (row.node <= 1)
continue;
5557 const std::size_t eidx = lqn.callpair_dst[row.aidx];
5560 for (
const UpdRow& row : call_map) {
5561 if (row.node <= 1)
continue;
5562 const std::size_t eidx = lqn.callpair_dst[row.aidx];
5565 callservt[row.aidx] = servt[eidx];
5566 callservtproc[row.aidx] = servtproc[eidx];
5567 }
else if (callservt[row.aidx] > Tzero()) {
5583 void moment3_distribution_pass(
int it) {
5584 if constexpr (!num_traits<T>::has_transcendental) {
5586 "SolverLN: the 'moment3' method fits an APH to a response-time CDF, which needs "
5587 "square roots and a matrix exponential; rerun with --arith double or real");
5589 const std::size_t dim = lqn.nidx + lqn.ncalls;
5590 const Matrix<T> Rinv = entry_service_resolvent();
5593 mam::AphPair<T> zero_law;
5594 zero_law.alpha.push_back(Tone());
5597 for (std::size_t en = 1; en <= lqn.nentries; ++en) {
5598 const std::size_t eidx = lqn.eshift + en;
5599 std::vector<mam::AphPair<T>> seq;
5600 for (std::size_t fitidx = 1; fitidx <= dim; ++fitidx) {
5601 if (!(Rinv(eidx, fitidx) > Tzero()))
continue;
5605 if (fitidx <= lqn.eshift + lqn.nentries)
continue;
5606 const bool is_call = fitidx > lqn.nidx;
5607 const fluid::FluidPassage& curve =
5608 is_call ? callservtcdf[fitidx - lqn.nidx] : servtcdf[fitidx];
5609 double m1 = 0.0, m2 = 0.0, m3 = 0.0;
5610 cdf_moments(curve, m1, m2, m3);
5614 if (!is_call && !lqn.actthink[fitidx].disabled &&
5616 const Distrib<T>& zt = lqn.actthink[fitidx];
5620 m3 = m3 + 3.0 * m2 * t1 + 3.0 * m1 * t2 + t3;
5621 m2 = m2 + 2.0 * m1 * t1 + t2;
5632 num_traits<T>::from_double(m1), num_traits<T>::from_double(m2),
5633 num_traits<T>::from_double(m3), 10u,
5635 const mam::AphPair<T> law = aph_pair_of_map(fit.aph);
5637 double reps = dbl(Rinv(eidx, fitidx));
5640 if (is_call) reps *= dbl(lqn.callproc_mean[fitidx - lqn.nidx]);
5641 const long whole =
static_cast<long>(std::floor(reps));
5642 const double frac = reps -
static_cast<double>(whole);
5643 for (
long q = 0; q < whole; ++q) seq.push_back(law);
5646 num_traits<T>::from_double(frac),
5647 num_traits<T>::from_double(1.0 - frac),
5651 const std::size_t cidx = fitidx - lqn.nidx;
5652 callservt[cidx] = num_traits<T>::from_double(m1);
5655 servt[fitidx] = num_traits<T>::from_double(m1);
5661 servt[eidx] = Tzero();
5665 entryproc[en] = entry_law;
5668 entrycdfrespt[en] = aph_eval_cdf(entry_law);
5676 std::vector<T> x(dim + 1, Tzero());
5677 for (std::size_t i = 1; i <= lqn.nidx; ++i) x[i] = residt[i];
5678 for (std::size_t c = 1; c <= lqn.ncalls; ++c) x[lqn.nidx + c] = callresidt[c];
5679 for (std::size_t en = 1; en <= lqn.nentries; ++en) {
5680 const std::size_t eidx = lqn.eshift + en;
5682 for (std::size_t j = 1; j <= dim; ++j)
5683 if (servtmatrix(eidx, j) != Tzero()) s += servtmatrix(eidx, j) * x[j];
5685 const int state = entry_visit_ratio(eidx, ratio);
5686 if (state == 0)
continue;
5687 residt[eidx] = state == 2 ? T(s * ratio) : s;
5691 for (
const UpdRow& row : call_map) {
5692 if (row.node <= 1)
continue;
5693 const std::size_t eidx = lqn.callpair_dst[row.aidx];
5694 callservt[row.aidx] = servt[eidx];
5701 moment_pass_done =
true;
5713 static mam::AphPair<T> aph_pair_of_map(
const mam::Map<T>& m) {
5714 mam::AphPair<T> out;
5715 const std::size_t n = m.D0.rows();
5717 out.alpha.assign(n, Tzero());
5718 for (std::size_t i = 0; i < n; ++i) {
5720 for (std::size_t j = 0; j < n; ++j) rowsum += m.D1(i, j);
5722 for (std::size_t j = 0; j < n; ++j) out.alpha[j] = T(m.D1(i, j) / rowsum);
5727 out.alpha[0] = Tone();
5735 static LnCdf aph_eval_cdf(
const mam::AphPair<T>& law) {
5739 const double var = m2 - m1 * m1;
5740 const double sigma = var > 0.0 ? std::sqrt(var) : 0.0;
5741 const double tmax = m1 + 10.0 * sigma;
5742 const std::size_t P = 500;
5745 for (std::size_t k = 0; k < P; ++k) {
5746 const double t = tmax *
static_cast<double>(k) /
static_cast<double>(P - 1);
5748 Matrix<T> St(law.S.rows(), law.S.cols(), Tzero());
5749 for (std::size_t i = 0; i < law.S.rows(); ++i)
5750 for (std::size_t j = 0; j < law.S.cols(); ++j)
5751 St(i, j) = T(law.S(i, j) * num_traits<T>::from_double(t));
5754 for (std::size_t i = 0; i < E.rows(); ++i)
5755 for (std::size_t j = 0; j < E.cols(); ++j) surv += law.alpha[i] * E(i, j);
5756 out.cdf[k] = 1.0 - dbl(surv);
5775 std::vector<T> join_excess()
const {
5776 std::vector<T> excess(lqn.nidx + 1, Tzero());
5777 std::vector<std::size_t> joined;
5778 for (std::size_t tail = 1; tail <= lqn.nidx; ++tail) {
5779 if (lqn.actpretype[tail] != PrecedenceType::PRE_AND)
continue;
5780 for (std::size_t sx : lqn.graph.succ(tail))
5781 if (sx > lqn.ashift && sx <= lqn.ashift + lqn.nacts) joined.push_back(sx);
5783 std::sort(joined.begin(), joined.end());
5784 joined.erase(std::unique(joined.begin(), joined.end()), joined.end());
5785 if (joined.empty())
return excess;
5787 fj::LqnBranchView<T> view;
5788 view.graph =
Matrix<T>(lqn.nidx, lqn.nidx, Tzero());
5789 for (std::size_t i = 1; i <= lqn.nidx; ++i)
5790 for (
const auto& e : lqn.graph.row[i]) view.graph(i - 1, e.first - 1) = e.second;
5791 view.ashift = lqn.ashift;
5792 view.nacts = lqn.nacts;
5793 view.actposttype.assign(lqn.nidx, 0);
5794 for (std::size_t i = 1; i <= lqn.nidx; ++i)
5795 view.actposttype[i - 1] =
static_cast<int>(lqn.actposttype[i]);
5797 for (std::size_t aidx : joined) {
5798 const std::vector<std::vector<std::size_t>> branches =
5800 if (branches.empty())
continue;
5801 std::vector<T> means;
5802 for (
const auto& br : branches) {
5804 for (std::size_t a : br) {
5808 if (a < lqn.callsof.size())
5809 for (std::size_t cidx : lqn.callsof[a])
5810 if (cidx < lqn.calltype.size() && lqn.calltype[cidx] == CallType::SYNC
5811 && cidx < callresidt.size())
5812 s += callresidt[cidx];
5816 if (means.size() == 1)
continue;
5817 std::size_t quorum = means.size();
5818 if (lqn.actquorum[aidx] >= 1 && lqn.actquorum[aidx] <= means.size())
5819 quorum = lqn.actquorum[aidx];
5820 std::vector<T> vars;
5821 for (
const T& m2 : means) vars.push_back(T(m2 * m2));
5823 for (
const T& m2 : means) serial += m2;
5824 if constexpr (num_traits<T>::has_transcendental) {
5826 excess[aidx] = T(q.m - serial);
5829 "SolverLN: the AND-join completion time is a k-th order statistic fitted "
5830 "through a three-point distribution, which needs a square root; use the "
5831 "double or real backend for a model with an AND join");
5838 void resolve_entry_service() {
5839 const std::size_t dim = lqn.nidx + lqn.ncalls;
5840 std::vector<T> x(dim + 1, Tzero());
5841 for (std::size_t i = 1; i <= lqn.nidx; ++i) x[i] = residt[i];
5842 for (std::size_t c = 1; c <= lqn.ncalls; ++c) x[lqn.nidx + c] = callresidt[c];
5843 std::vector<T> entry_servt(dim + 1, Tzero());
5844 for (std::size_t i = 1; i <= dim; ++i) {
5846 for (std::size_t j = 1; j <= dim; ++j)
5847 if (servtmatrix(i, j) != Tzero()) s += servtmatrix(i, j) * x[j];
5853 const std::vector<T> excess = join_excess();
5854 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5855 const std::size_t eidx = lqn.eshift + e;
5856 for (std::size_t aidx = 1; aidx <= lqn.nidx; ++aidx) {
5857 if (excess[aidx] == Tzero())
continue;
5858 if (servtmatrix(eidx, aidx) > Tzero())
5859 entry_servt[eidx] = T(entry_servt[eidx] + excess[aidx]);
5861 if (entry_servt[eidx] < Tzero()) entry_servt[eidx] = Tzero();
5871 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5872 const std::size_t eidx = lqn.eshift + e;
5873 const double c = setup_charge(lqn.parent[eidx]);
5875 entry_servt[eidx] = T(entry_servt[eidx] + num_traits<T>::from_double(c));
5880 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
5881 const std::size_t eidx = lqn.eshift + e;
5882 const std::size_t tidx = lqn.parent[eidx];
5883 const std::size_t hidx = lqn.parent[tidx];
5884 if (ignore[tidx] || ignore[hidx])
continue;
5885 if (idxhash[hidx] < 0)
continue;
5886 const bool has_sync = lqn.issynccaller.any_col(eidx);
5888 servt[eidx] = entry_servt[eidx];
5889 residt[eidx] = entry_servt[eidx];
5892 const std::size_t hl = std::size_t(idxhash[hidx]);
5893 const qn::Layer<T>& L = ensemble[hl];
5894 const LayerResult<T>& r = results.back()[hl];
5895 T task_tput = Tzero(), entry_tput = Tzero();
5896 for (
const auto& kv : L.attr_tasks)
5897 if (kv.second == tidx) task_tput += r.TN(L.clientIdx - 1, kv.first - 1);
5898 for (
const auto& kv : L.attr_entries)
5899 if (kv.second == eidx) entry_tput += r.TN(L.clientIdx - 1, kv.first - 1);
5901 servt[eidx] = T(entry_servt[eidx] * task_tput / entry_tput);
5902 residt[eidx] = servt[eidx];
5904 servt[eidx] = entry_servt[eidx];
5905 residt[eidx] = entry_servt[eidx];
5921 bool layer_takes_interlock(std::size_t e)
const {
5922 return opt.layer_solver ==
"mva" && e < ensemble.size() &&
5938 std::vector<std::vector<double>> build_layer_interlock(std::size_t e,
5939 const std::vector<std::size_t>& hts,
5940 const std::vector<T>& prIL,
5941 const std::vector<T>& PrIL)
const {
5942 const qn::Layer<T>& L = ensemble[e];
5943 const std::size_t R = L.nclasses;
5944 std::vector<double> cls_ir(R, 0.0);
5945 std::vector<double> cls_pr(R, 0.0);
5946 auto stamp = [&](std::size_t classIdx, std::size_t tidx) {
5947 if (classIdx < 1 || classIdx > R)
return;
5948 for (std::size_t i = 0; i < hts.size(); ++i)
5949 if (hts[i] == tidx) {
5950 cls_ir[classIdx - 1] = dbl(prIL[i]);
5951 cls_pr[classIdx - 1] = dbl(PrIL[i]);
5954 for (
const auto& a : L.attr_tasks) stamp(a.first, a.second);
5955 for (
const auto& a : L.attr_entries) stamp(a.first, lqn.parent[a.second]);
5956 for (
const auto& a : L.attr_activities) stamp(a.first, lqn.parent[a.second]);
5957 for (
const auto& a : L.attr_calls) stamp(a[0], lqn.parent[a[2]]);
5959 std::vector<std::vector<double>> IL(R, std::vector<double>(R, 0.0));
5961 for (std::size_t r = 0; r < R; ++r) {
5963 for (std::size_t sIl = 0; sIl < R; ++sIl) {
5965 IL[r][sIl] = cls_pr[sIl] * cls_ir[sIl] * cls_ir[r];
5969 return any ? IL : std::vector<std::vector<double>>();
5982 std::pair<T, T> interlock_prob(std::size_t client_tidx, std::size_t server_idx,
5983 bool isProcessorHost)
const {
5984 const std::vector<std::size_t>& common = il_common_entries[server_idx];
5985 const double nsrc = il_num_sources[server_idx];
5986 if (nsrc == 0.0 || common.empty())
return std::make_pair(Tzero(), Tzero());
5987 const std::vector<std::size_t>& allsrc = il_src_all[server_idx];
5988 T sum_flow = Tzero();
5989 T sum_pril = Tzero();
5990 for (std::size_t ce : common) {
5991 const std::size_t srcTask = lqn.parent[ce];
5992 const std::size_t cen = ce - lqn.eshift;
5994 double m_src = lqn.mult[srcTask];
5995 if (!std::isfinite(m_src) || m_src < 1.0) m_src = 1.0;
5996 double m_eff = m_src;
5997 if (isProcessorHost) m_eff = (m_src > 3.0) ? (m_src + m_src) : (m_src * m_src);
5998 for (std::size_t dst : lqn.entriesof[client_tidx]) {
5999 const std::size_t dn = dst - lqn.eshift;
6000 if (dn < 1 || dn > lqn.nentries)
continue;
6001 if (!(il_all(cen, dn) > Tzero()))
continue;
6002 T ce_tput = tput[ce];
6004 ce_tput = tput[lqn.actsof[ce][0]];
6008 if (std::find(allsrc.begin(), allsrc.end(), srcTask) != allsrc.end()) {
6009 const T contrib = T(ce_tput * il_all(cen, dn));
6010 sum_flow += contrib;
6011 sum_pril += T(contrib / num_traits<T>::from_double(m_eff));
6015 T client_tput = tput[client_tidx];
6017 for (std::size_t e : lqn.entriesof[client_tidx]) {
6020 et = tput[lqn.actsof[e][0]];
6021 client_tput = T(client_tput + et);
6025 return std::make_pair(Tzero(), Tzero());
6026 const T flow = sum_flow < client_tput ? sum_flow : client_tput;
6027 T IR = T(flow / client_tput);
6028 if (IR > Tone()) IR = Tone();
6029 if (dbl(IR) < 0.0) IR = Tzero();
6032 if (pr > Tone()) pr = Tone();
6033 if (dbl(pr) < 0.0) pr = Tzero();
6034 return std::make_pair(IR, pr);
6037 void update_populations(
int) {
6038 const std::vector<T> residt_orig = residt;
6039 const std::vector<T> callresidt_orig = callresidt;
6040 bool adjusted =
false;
6042 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
6043 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
6044 const std::size_t dst = lqn.callpair_dst[cidx];
6045 const std::size_t server_tidx = lqn.parent[dst];
6046 std::size_t server = 0;
6047 if (server_tidx <= NT() && !il_common_entries[server_tidx].empty()) {
6048 server = server_tidx;
6049 }
else if (server_tidx > lqn.tshift) {
6050 const std::size_t h = lqn.parent[server_tidx];
6051 if (h >= 1 && h <= NT() && !il_common_entries[h].empty()) server = h;
6053 if (server == 0)
continue;
6054 const std::size_t client_tidx = lqn.parent[lqn.callpair_src[cidx]];
6058 const std::pair<T, T> ilp = interlock_prob(client_tidx, server,
false);
6059 const T prIL = T(ilp.first * ilp.second);
6061 const T S = servt[dst];
6062 const T cm = lqn.callproc_mean[cidx];
6063 if (!(cm > Tzero()) || !(callservt[cidx] > Tzero()))
continue;
6064 const T RN = T(callservt[cidx] / cm);
6065 const T W = RN > S ? T(RN - S) : Tzero();
6067 const T RN_adj = T(S + (Tone() - prIL) * W);
6068 const T scale = T(RN_adj / RN);
6069 callservt[cidx] = T(callservt[cidx] * scale);
6070 callresidt[cidx] = T(callresidt[cidx] * scale);
6071 if (callservt[cidx] > Tzero())
6078 layer_interlock.assign(ensemble.size(), {});
6079 for (std::size_t h = 1; h <= lqn.nhosts; ++h) {
6080 if (il_common_entries[h].empty())
continue;
6081 const std::vector<std::size_t>& hts = lqn.tasksof[h];
6082 std::vector<T> prIL(hts.size(), Tzero()), PrIL(hts.size(), Tzero()),
6083 tutil(hts.size(), Tzero());
6084 for (std::size_t i = 0; i < hts.size(); ++i) {
6086 const std::pair<T, T> ilp = interlock_prob(hts[i], h,
true);
6087 prIL[i] = ilp.first;
6088 PrIL[i] = ilp.second;
6089 for (std::size_t e : lqn.entriesof[hts[i]])
6090 for (std::size_t a : lqn.actsof[e]) tutil[i] += tput[a] * lqn.hostdem[a].mean;
6092 T Utot = Tzero(), Uil = Tzero();
6093 for (std::size_t i = 0; i < hts.size(); ++i) {
6099 const T frac = T(Uil / Utot);
6107 if (idxhash[h] >= 0 && layer_takes_interlock(std::size_t(idxhash[h]))) {
6108 const std::size_t e = std::size_t(idxhash[h]);
6109 layer_interlock[e] = build_layer_interlock(e, hts, prIL, PrIL);
6113 for (std::size_t i = 0; i < hts.size(); ++i) {
6121 const T eff = T(prIL[i] * PrIL[i] * frac);
6122 for (std::size_t e : lqn.entriesof[hts[i]])
6123 for (std::size_t a : lqn.actsof[e]) {
6124 const T D = lqn.hostdem[a].mean;
6126 residt[a] = T(D + (Tone() - eff) * (residt[a] - D));
6133 if (!adjusted)
return;
6134 const std::size_t dim = lqn.nidx + lqn.ncalls;
6135 auto esum = [&](
const std::vector<T>& rs,
const std::vector<T>& cr, std::size_t i) {
6137 for (std::size_t j = 1; j <= dim; ++j) {
6138 if (servtmatrix(i, j) == Tzero())
continue;
6139 s += servtmatrix(i, j) * (j <= lqn.nidx ? rs[j] : cr[j - lqn.nidx]);
6143 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
6144 const std::size_t eidx = lqn.eshift + e;
6145 const T oldv = esum(residt_orig, callresidt_orig, eidx);
6147 const T ratio = T(esum(residt, callresidt, eidx) / oldv);
6156 if (!moment_pass_done) {
6157 servt[eidx] = T(servt[eidx] * ratio);
6160 residt[eidx] = T(residt[eidx] * ratio);
6176 T ref_think_time(std::size_t tidx)
const {
6177 if (tidx >= lqn.isref.size() || !lqn.isref[tidx])
return Tzero();
6178 if (lqn.think[tidx].disabled)
return Tzero();
6179 return lqn.think[tidx].mean;
6184 void update_think_times(
int it) {
6187 if (is_ph_encoding()) {
6188 update_think_times_ph(it);
6192 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
6193 const std::size_t tidx = lqn.tshift + t;
6195 const T ztask = ref_think_time(tidx);
6196 if (idxhash[tidx] < 0) {
6200 const double arvrate = open_arrival_rate_of(tidx);
6203 for (std::size_t c = 1; c <= NT(); ++c) nja = std::max(nja, njobs(tidx, c));
6204 if (!(nja > 0.0)) nja = lqn.maxmult[tidx];
6206 for (std::size_t eidx : lqn.entriesof[tidx])
6207 for (std::size_t aidx : lqn.actsof[eidx])
6208 if (!std::isnan(dbl(residt[aidx]))) hres += residt[aidx];
6210 T za = T(num_traits<T>::from_double(nja / arvrate) - hres - ztask);
6211 if (za < floora) za = floora;
6212 if (relax_omega < 1.0 && it > 1 && !std::isnan(thinkt_prev[tidx])) {
6213 const T om = num_traits<T>::from_double(relax_omega);
6214 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
6215 za = T(om * za + om1 * thinkt_prev_v[tidx]);
6217 tput[tidx] = num_traits<T>::from_double(arvrate);
6219 thinkt_prev[tidx] = dbl(za);
6220 thinkt_prev_v[tidx] = za;
6228 const std::size_t e = std::size_t(idxhash[tidx]);
6229 const qn::Layer<T>& L = ensemble[e];
6230 const LayerResult<T>& r = results.back()[e];
6232 for (std::size_t c = 1; c <= NT(); ++c) nj = std::max(nj, njobs(tidx, c));
6234 const std::size_t ts = station_idx_of(L, tidx);
6235 T tp = Tzero(), ut = Tzero();
6236 for (std::size_t k = 0; k < L.nclasses; ++k) {
6237 tp += r.TN(ts - 1, k);
6238 ut += r.UN(ts - 1, k);
6240 tput[tidx] = T(num_traits<T>::from_double(lqn.repl[tidx]) * tp);
6243 if (lqn.sched[tidx] == SchedStrategy::INF) {
6245 raw = tput[tidx] == Tzero()
6247 : T((num_traits<T>::from_double(nj) - util[tidx]) / tput[tidx] - ztask);
6249 const T om = util[tidx] > Tone() ? T(util[tidx] - Tone()) : T(Tone() - util[tidx]);
6250 raw = tput[tidx] == Tzero()
6252 : T(num_traits<T>::from_double(nj) * om / tput[tidx] - ztask);
6255 thinkt[tidx] = raw < floorv ? floorv : raw;
6266 const double c = setup_charge(tidx);
6268 const T withc = T(thinkt[tidx] + num_traits<T>::from_double(c));
6269 thinkt[tidx] = withc < floorv ? floorv : withc;
6273 if (relax_omega < 1.0 && it > 1 && !std::isnan(thinkt_prev[tidx])) {
6274 const double rawd = dbl(thinkt[tidx]);
6276 thinkt_prev[tidx] = rawd;
6277 thinkt_prev_v[tidx] = thinkt[tidx];
6279 const T om = num_traits<T>::from_double(relax_omega);
6280 const T om1 = num_traits<T>::from_double(1.0 - relax_omega);
6281 thinkt[tidx] = T(om * thinkt[tidx] + om1 * thinkt_prev_v[tidx]);
6283 thinkt_prev[tidx] = dbl(thinkt[tidx]);
6284 thinkt_prev_v[tidx] = thinkt[tidx];
6292 void update_layers(
int it) {
6295 if (is_ph_encoding()) {
6296 update_layers_ph(it);
6299 const bool elevator = (it % 2) == 1;
6300 const std::size_t nt = thinkt_map.size();
6301 for (std::size_t r = 0; r < nt; ++r) {
6302 const UpdRow& row = thinkt_map[elevator ? nt - 1 - r : r];
6303 qn::Layer<T>& L = ensemble[std::size_t(idxhash[row.idx])];
6304 if (row.aidx <= NT() && opt.interlocking &&
6305 L.classes[row.cls - 1].type == JobClassType::CLOSED)
6306 L.classes[row.cls - 1].population = njobs(row.aidx, row.idx);
6307 if (row.node == L.clientIdx) {
6308 if (lqn.type[row.aidx] == LqnElement::TASK) {
6309 if (lqn.sched[row.aidx] != SchedStrategy::REF) {
6310 if (!thinktproc[row.aidx].disabled)
6311 L.set_service(row.node, row.cls, thinktproc[row.aidx]);
6313 L.set_service(row.node, row.cls, servtproc[row.aidx]);
6316 L.set_service(row.node, row.cls, servtproc[row.aidx]);
6319 L.set_service(row.node, row.cls, servtproc[row.aidx]);
6322 const std::size_t nc = call_map.size();
6323 for (std::size_t r = 0; r < nc; ++r) {
6324 const UpdRow& row = call_map[elevator ? nc - 1 - r : r];
6325 qn::Layer<T>& L = ensemble[std::size_t(idxhash[row.idx])];
6326 if (row.node == L.clientIdx)
6327 L.set_service(row.node, row.cls, callservtproc[row.aidx]);
6329 L.set_service(row.node, row.cls, servtproc[lqn.callpair_dst[row.aidx]]);
6336 for (
const UpdRow& row : arv_call_map) {
6337 qn::Layer<T>& L = ensemble[std::size_t(idxhash[row.idx])];
6338 const Distrib<T>& d = tputproc[lqn.callpair_src[row.aidx]];
6339 if (!d.disabled) L.set_service(row.node, row.cls, d);
6356 T region_wait(std::size_t e, std::size_t cls)
const {
6357 const qn::Layer<T>& L = ensemble[e];
6358 if (L.regions.empty() || results.empty())
return Tzero();
6359 const LayerResult<T>& r = results.back()[e];
6360 const std::size_t k = cls - 1;
6361 std::size_t c = L.nchains;
6362 for (std::size_t cc = 0; cc < L.nchains; ++cc)
6363 if (L.chains[cc][k]) c = cc;
6364 if (c == L.nchains)
return Tzero();
6368 std::size_t rs = L.serverIdx;
6369 for (std::size_t i = 0; i < L.regions[0].members.size(); ++i)
6370 if (L.regions[0].members[i]) { rs = i + 1;
break; }
6373 T inside = Tzero(), tput_srv = Tzero();
6374 for (std::size_t j = 0; j < L.nclasses; ++j) {
6375 if (!L.chains[c][j])
continue;
6376 const double p = L.classes[j].population;
6377 if (!std::isfinite(p))
return Tzero();
6379 for (std::size_t i = 0; i < L.nstations; ++i) inside += r.QN(i, j);
6380 tput_srv += r.TN(rs - 1, j);
6382 const T deficit = T(num_traits<T>::from_double(pop) - inside);
6383 if (!(deficit > Tzero()) || !(tput_srv > Tzero()))
return Tzero();
6384 return T(deficit / tput_srv);
6388 void update_routing_probabilities(
int) {
6389 for (std::size_t u = 0; u < unique_route_idx.size(); ++u) {
6392 const std::size_t idx = unique_route_idx[unique_route_idx.size() - 1 - u];
6393 qn::Layer<T>& L = ensemble[std::size_t(idxhash[idx])];
6394 bool updated =
false;
6395 for (
const RouteRow& r : route_map) {
6396 if (r.idx != idx)
continue;
6397 if (idxhash[r.tidx_caller] < 0)
continue;
6398 const std::size_t cl = std::size_t(idxhash[r.tidx_caller]);
6399 const qn::Layer<T>& CL = ensemble[cl];
6400 const LayerResult<T>& cr = results.back()[cl];
6403 const std::size_t cs = station_idx_of(CL, r.tidx_caller);
6405 for (std::size_t k = 0; k < CL.nclasses; ++k) Xtot += cr.TN(cs - 1, k);
6406 if (!(Xtot > Tzero()))
continue;
6407 T entry_tput = Tzero();
6408 for (
const auto& a : CL.attr_calls)
6410 entry_tput += cr.TN(station_idx_of(CL, lqn.parent[a[3]]) - 1, a[0] - 1);
6411 L.set_route(r.cfrom, r.cto, r.nodefrom, r.nodeto, T(entry_tput / Xtot));
6414 if (updated) L.refresh_chains();
6433 void require_transient_ready()
const {
6434 if (opt.layer_solver !=
"fluid")
6436 "SolverLN: the layered transient is the transient OF EACH LAYER, which only the "
6437 "fluid layer solver has; rebuild the ensemble with layer_solver 'fluid'");
6441 bool has_finite_horizon()
const {
6442 return std::isfinite(opt.timespan_end) && opt.timespan_end > 0.0;
6447 std::vector<std::vector<std::vector<double>>> Q, U, Tp, R;
6454 void run_layer_transient(std::size_t e,
const std::vector<FluidRateSched>& sched,
6455 LnTranLayer& block, TranTraj& traj) {
6456 qn::Layer<T>& L = ensemble[e];
6458 fluid::FluidOptions fo = opt.layer_fluid;
6459 fo.rate_sched = sched;
6465 if (e < layer_tran_init.size() && !layer_tran_init[e].empty())
6466 fo.init_sol = layer_tran_init[e];
6471 const double t_end =
6473 const std::vector<fluid::FluidTranPoint> pts =
6474 detail::ln_fluid_transient(L, fo, t_end, opt.tran_points, opt.tran_grid);
6475 const std::size_t M = L.nstations, K = L.nclasses, P = pts.size();
6477 auto alloc = [&](std::vector<std::vector<std::vector<double>>>& A) {
6478 A.assign(M, std::vector<std::vector<double>>(K, std::vector<double>(P, 0.0)));
6487 for (std::size_t p = 0; p < P; ++p) {
6488 block.t[p] = pts[p].t;
6489 for (std::size_t i = 0; i < M; ++i)
6490 for (std::size_t r = 0; r < K; ++r) {
6491 const double q = pts[p].QN(i, r), u = pts[p].UN(i, r), x = pts[p].TN(i, r);
6492 block.QN[i][r][p] = q;
6493 block.UN[i][r][p] = u;
6494 block.TN[i][r][p] = x;
6495 traj.Q[i][r][p] = q;
6496 traj.U[i][r][p] = u;
6497 traj.Tp[i][r][p] = x;
6509 LnTranSolution tran_avg_decoupled() {
6510 require_transient_ready();
6511 if (results.empty()) iterate();
6513 out.mode =
"decoupled";
6514 out.layers.resize(ensemble.size());
6515 for (std::size_t e = 0; e < ensemble.size(); ++e) {
6517 run_layer_transient(e, std::vector<FluidRateSched>(), out.layers[e], tj);
6539 LnTranSolution tran_avg_coupled() {
6540 require_transient_ready();
6545 if (!has_finite_horizon())
return tran_avg_decoupled();
6546 if (results.empty()) iterate();
6547 const std::size_t E = ensemble.size();
6550 out.mode =
"coupled";
6551 out.layers.resize(E);
6552 std::vector<TranTraj> traj(E), prev(E);
6553 for (std::size_t e = 0; e < E; ++e)
6554 run_layer_transient(e, std::vector<FluidRateSched>(), out.layers[e], traj[e]);
6555 const std::vector<double> tgrid = out.layers.empty() ? std::vector<double>() : out.
layers[0].t;
6557 for (
long iter = 1; iter <= opt.ln_transient_iter_max; ++iter) {
6559 const std::vector<std::vector<FluidRateSched>> sched =
6560 build_rate_sched(recompute_demand(traj, tgrid), tgrid);
6561 for (std::size_t e = 0; e < E; ++e)
6562 run_layer_transient(e, sched[e], out.layers[e], traj[e]);
6564 for (std::size_t e = 0; e < E; ++e)
6565 for (std::size_t i = 0; i < traj[e].Q.size(); ++i)
6566 for (std::size_t r = 0; r < traj[e].Q[i].size(); ++r)
6567 for (std::size_t p = 0; p < traj[e].Q[i][r].size(); ++p)
6568 gap = std::max(gap, std::fabs(traj[e].Q[i][r][p] - prev[e].Q[i][r][p]));
6569 out.iterations = iter;
6571 if (gap < opt.ln_transient_tol)
break;
6578 std::map<std::size_t, std::vector<double>> thinkt;
6579 std::map<std::size_t, std::vector<double>> callservt;
6590 TranDemand recompute_demand(
const std::vector<TranTraj>& traj,
6591 const std::vector<double>& tgrid)
const {
6593 const std::size_t ng = tgrid.size();
6595 for (std::size_t t = 1; t <= lqn.ntasks; ++t) {
6596 const std::size_t tidx = lqn.tshift + t;
6597 if (idxhash[tidx] < 0 || lqn.isref[tidx])
continue;
6598 const std::size_t e = std::size_t(idxhash[tidx]);
6599 const qn::Layer<T>& L = ensemble[e];
6600 const std::size_t s = station_idx_of(L, tidx) - 1;
6602 for (std::size_t c = 1; c <= NT(); ++c) nj = std::max(nj, njobs(tidx, c));
6604 const double userthink = dbl(ref_think_time(tidx));
6605 std::vector<double> tk(ng, 0.0);
6606 for (std::size_t p = 0; p < ng; ++p) {
6607 double U = 0.0, X = 0.0;
6608 for (std::size_t r = 0; r < L.nclasses; ++r) {
6609 U += traj[e].U[s][r][p];
6610 X += traj[e].Tp[s][r][p];
6613 double v = lqn.sched[tidx] == SchedStrategy::INF
6614 ? (nj - U) / Xs - userthink
6615 : nj * std::fabs(1.0 - U) / Xs - userthink;
6617 tk[p] = v + userthink;
6619 out.thinkt[tidx] = tk;
6622 for (std::size_t cidx = 1; cidx <= lqn.ncalls; ++cidx) {
6623 if (lqn.calltype[cidx] != CallType::SYNC)
continue;
6624 const std::size_t eidx = lqn.callpair_dst[cidx];
6625 const std::size_t tidx = lqn.parent[eidx];
6626 if (tidx > NT() || idxhash[tidx] < 0)
continue;
6627 const std::size_t e = std::size_t(idxhash[tidx]);
6628 const qn::Layer<T>& L = ensemble[e];
6629 const std::size_t s = station_idx_of(L, tidx) - 1;
6630 std::vector<double> Rc(ng, 0.0);
6632 for (std::size_t r = 0; r < L.nclasses; ++r) {
6633 if (L.classes[r].attr_kind !=
int(LqnElement::ENTRY))
continue;
6634 if (L.classes[r].attr_idx != eidx)
continue;
6636 for (std::size_t p = 0; p < ng; ++p) Rc[p] += traj[e].R[s][r][p];
6641 for (std::size_t r = 0; r < L.nclasses; ++r)
6642 for (std::size_t p = 0; p < ng; ++p) Rc[p] += traj[e].R[s][r][p];
6644 const double cm = dbl(lqn.callproc_mean[cidx]);
6645 for (std::size_t p = 0; p < ng; ++p) Rc[p] *= cm;
6646 out.callservt[cidx] = Rc;
6662 std::vector<std::vector<FluidRateSched>> build_rate_sched(
6663 const TranDemand& demand,
const std::vector<double>& tgrid)
const {
6664 std::vector<std::vector<FluidRateSched>> out(ensemble.size());
6665 const bool want_think =
6666 opt.ln_transient_channels ==
"both" || opt.ln_transient_channels ==
"thinkt";
6667 const bool want_call =
6668 opt.ln_transient_channels ==
"both" || opt.ln_transient_channels ==
"callservt";
6670 auto add = [&](std::size_t e, std::size_t station, std::size_t cls,
6671 const std::vector<double>& d) {
6672 if (d.empty())
return;
6673 const double dend = d.back();
6678 const double cap = 20.0;
6683 s.rates.resize(d.size());
6684 for (std::size_t p = 0; p < d.size(); ++p) {
6685 const double v = std::min(std::max(d[p], dend / cap), dend * cap);
6686 s.rates[p] = 1.0 / v;
6688 s.nominal = 1.0 / dend;
6689 out[e].push_back(s);
6693 for (
const UpdRow& row : thinkt_map) {
6694 if (idxhash[row.idx] < 0)
continue;
6695 const std::size_t e = std::size_t(idxhash[row.idx]);
6696 if (row.node != ensemble[e].clientIdx)
continue;
6697 if (lqn.type[row.aidx] != LqnElement::TASK)
continue;
6698 if (lqn.sched[row.aidx] == SchedStrategy::REF)
continue;
6699 const auto it = demand.thinkt.find(row.aidx);
6700 if (it == demand.thinkt.end())
continue;
6701 add(e, row.node, row.cls, it->second);
6704 for (
const UpdRow& row : call_map) {
6705 if (idxhash[row.idx] < 0)
continue;
6706 const std::size_t e = std::size_t(idxhash[row.idx]);
6707 if (row.node != ensemble[e].clientIdx)
continue;
6708 const auto it = demand.callservt.find(row.aidx);
6709 if (it == demand.callservt.end())
continue;
6710 add(e, row.node, row.cls, it->second);
6718 LnSolution<T> aggregate() {
6719 const std::size_t N = lqn.nidx;
6721 auto mk = [&](std::vector<T>& v, std::vector<bool>& d) {
6722 v.assign(N + 1, Tzero());
6723 d.assign(N + 1,
false);
6725 std::vector<T>
QN,
UN, RN, TN, PN, SN, WN, AN;
6726 std::vector<bool> dQ, dU, dR, dT, dP, dS, dW, dA;
6727 mk(QN, dQ); mk(UN, dU); mk(RN, dR); mk(TN, dT);
6728 mk(PN, dP); mk(SN, dS); mk(WN, dW); mk(AN, dA);
6729 std::vector<bool> wn_done(N + 1,
false);
6731 for (std::size_t e = 0; e < ensemble.size(); ++e) {
6732 const qn::Layer<T>& L = ensemble[e];
6733 const LayerResult<T>& r = results.back()[e];
6734 const std::size_t clientIdx = L.clientIdx;
6739 const bool has_host_server = !L.host_stations.empty();
6740 for (std::size_t hs : L.host_stations) {
6741 const std::size_t hidx = L.stations[hs - 1].attr_idx;
6746 for (std::size_t k = 0; k < L.nclasses; ++k) {
6747 if (L.classes[k].completes) {
6748 T t = clientIdx > 0 ? r.TN(clientIdx - 1, k) : Tzero();
6749 const T ts = r.TN(hs - 1, k);
6750 TN[hidx] = T(TN[hidx] + (t > ts ? t : ts));
6752 if (L.classes[k].attr_kind ==
int(LqnElement::ACTIVITY)) {
6754 if (station_idx_of_class(L, k) != hs)
continue;
6755 const std::size_t aidx = L.classes[k].attr_idx;
6756 const std::size_t tidx = lqn.parent[aidx];
6759 PN[aidx] = T(PN[aidx] + r.UN(hs - 1, k));
6760 PN[tidx] = T(PN[tidx] + r.UN(hs - 1, k));
6761 PN[hidx] = T(PN[hidx] + r.UN(hs - 1, k));
6767 for (std::size_t k = 0; k < L.nclasses; ++k) {
6768 const int kind = L.classes[k].attr_kind;
6772 const std::size_t serverIdx = station_idx_of_class(L, k);
6773 if (kind ==
int(LqnElement::TASK)) {
6774 const std::size_t tidx = L.classes[k].attr_idx;
6775 if (has_host_server && !dT[tidx]) {
6777 TN[tidx] = r.TN(clientIdx - 1, k);
6779 }
else if (kind ==
int(LqnElement::ENTRY)) {
6780 const std::size_t eidx = L.classes[k].attr_idx;
6787 if (has_host_server && !dT[eidx]) {
6789 TN[eidx] = r.TN(clientIdx - 1, k);
6791 }
else if (kind ==
int(LqnElement::CALL)) {
6792 const std::size_t cidx = L.classes[k].attr_idx;
6793 const std::size_t aidx = lqn.callpair_src[cidx];
6794 if (lqn.calltype[cidx] == CallType::SYNC) {
6796 SN[aidx] = T(SN[aidx] + r.RN(serverIdx - 1, k) * lqn.callproc_mean[cidx]);
6799 QN[aidx] = T(QN[aidx] + r.QN(serverIdx - 1, k));
6800 }
else if (kind ==
int(LqnElement::ACTIVITY)) {
6801 const std::size_t aidx = L.classes[k].attr_idx;
6802 const std::size_t tidx = lqn.parent[aidx];
6804 QN[tidx] = T(QN[tidx] + r.QN(serverIdx - 1, k));
6807 TN[aidx] = T(TN[aidx] + r.TN(serverIdx - 1, k));
6809 SN[aidx] = T(SN[aidx] + r.RN(serverIdx - 1, k));
6811 RN[aidx] = T(RN[aidx] + r.RN(serverIdx - 1, k));
6814 WN[aidx] = residt[aidx];
6815 if (!wn_done[aidx]) {
6816 WN[tidx] = T(WN[tidx] + residt[aidx]);
6817 wn_done[aidx] =
true;
6819 QN[aidx] = T(QN[aidx] + r.QN(serverIdx - 1, k));
6824 for (std::size_t e = 1; e <= lqn.nentries; ++e) {
6825 const std::size_t eidx = lqn.eshift + e;
6826 const std::size_t tidx = lqn.parent[eidx];
6833 ? T(TN[eidx] * (servt_ph1[eidx] + servt_ph2[eidx]))
6834 : T(TN[eidx] * SN[eidx]);
6837 for (std::size_t a : lqn.actsof[eidx])
6842 if (!lqn.actsof[eidx].empty()) {
6844 PN[eidx] = any ? ps : Tzero();
6846 for (std::size_t a : lqn.actsof[tidx]) {
6848 UN[a] = T(TN[a] * SN[a]);
6850 UN[tidx] = T(UN[tidx] + UN[eidx]);
6862 for (std::size_t i = 1; i <= N; ++i)
6867 const bool host = lqn.type[i] == LqnElement::HOST;
6868 const bool task = lqn.type[i] == LqnElement::TASK;
6869 const bool entry = lqn.type[i] == LqnElement::ENTRY;
6873 dS[i] = !host && !task;
6875 dW[i] = !host && !entry;
6881 s.QN =
UN; s.defined_Q = dU;
6882 s.UN = PN; s.defined_U = dP;
6883 s.RN = SN; s.defined_R = dS;
6884 s.TN = TN; s.defined_T = dT;
6885 s.AN = AN; s.defined_A = dA;
6886 s.WN = WN; s.defined_W = dW;
6887 s.iterations = iterations_done;
6888 s.converged = did_converge;