1009 ldes_engine_reject(
sn, o);
1011 const std::size_t M =
sn.nstations, K =
sn.nclasses;
1012 if (M == 0 || K == 0)
throw InputError(
"SolverLDES (native engine): empty model");
1014 const std::uint64_t max_events =
1016 if (max_events == 0)
1017 throw InputError(
"SolverLDES (native engine): the completion budget is zero");
1020 const std::uint64_t base =
1021 (o.
seed >= 0) ?
static_cast<std::uint64_t
>(o.
seed) : std::random_device{}();
1031 const bool slotted = o.
slotted;
1033 auto slot_snap = [&](
double v,
const char* what) ->
double {
1034 if (!slotted || v == 0.0)
return v;
1035 const double slots = v / slot_len;
1036 const double rounded = std::floor(slots + 0.5);
1037 if (rounded < 1.0 || std::fabs(slots - rounded) > 1e-9 * std::max(1.0, slots))
1038 throw InputError(std::string(
"SolverLDES (native engine): slotted mode sampled a ") +
1039 what +
" of " + std::to_string(v) +
1040 ", which is not a positive multiple of the slot length " +
1041 std::to_string(slot_len) +
1042 "; a discrete-time model needs lattice-valued interarrival and "
1043 "service times, e.g. Geometric or Det on an integral slot count");
1044 return rounded * slot_len;
1072 const long long seed_ll =
static_cast<long long>(base);
1073 const long long Kll =
static_cast<long long>(K);
1074 std::size_t num_sources = 0, num_service_nodes = 0;
1075 std::vector<std::size_t> svc_index(M, 0);
1076 for (std::size_t i = 0; i < M; ++i) {
1077 if (
sn.stations[i].nodetype == NodeType::Source) {
1080 svc_index[i] = num_service_nodes++;
1083 const long long nsrc =
static_cast<long long>(num_sources);
1084 const long long nsvc =
static_cast<long long>(num_service_nodes);
1086 std::vector<Rng> g_arr;
1088 for (std::size_t k = 0; k < K; ++k) {
1089 const long long off = (
static_cast<long long>(k)) * 10;
1090 g_arr.push_back(Rng(seed_ll, off, off + 2000));
1092 std::vector<std::vector<Rng>> g_svc;
1094 for (std::size_t i = 0; i < M; ++i) {
1095 std::vector<Rng> row;
1097 for (std::size_t k = 0; k < K; ++k) {
1098 const long long off =
1099 ((nsrc +
static_cast<long long>(svc_index[i])) * Kll +
static_cast<long long>(k)) *
1101 row.push_back(Rng(seed_ll, off));
1103 g_svc.push_back(row);
1108 std::vector<std::vector<std::vector<Rng>>> g_hsvc(M);
1109 for (std::size_t i = 0; i < M; ++i) {
1110 const std::size_t nT =
sn.stations[i].server_types.size();
1111 g_hsvc[i].resize(nT);
1112 for (std::size_t t = 0; t < nT; ++t) {
1113 g_hsvc[i][t].reserve(K);
1114 for (std::size_t k = 0; k < K; ++k) {
1115 const long long off =
1116 ((nsrc +
static_cast<long long>(svc_index[i])) * Kll +
1117 static_cast<long long>(k)) * 10 +
1118 600000 +
static_cast<long long>(t) * 137;
1119 g_hsvc[i][t].push_back(Rng(seed_ll, off));
1123 std::vector<Rng> g_aux;
1125 for (std::size_t i = 0; i < M; ++i) {
1126 const long long off = ((nsrc + nsvc +
static_cast<long long>(svc_index[i])) * Kll) * 10;
1127 g_aux.push_back(Rng(seed_ll, off));
1129 Rng g_routing(seed_ll, 900000);
1132 Rng g_spn(seed_ll, 5000);
1136 Rng g_qplace(seed_ll, 5001);
1141 Rng g_fork(seed_ll, 88888);
1146 Rng g_siro(seed_ll, 99999);
1154 auto fork_is_variable = [&](std::size_t node) ->
bool {
1156 if (fp == 0)
return false;
1157 const double tpl =
sn.nodes[node - 1].tasks_per_link;
1158 for (std::size_t k = 0; k < fp->
fan_out_link.rows(); ++k)
1159 for (std::size_t r = 0; r < fp->
fan_out_link.cols(); ++r) {
1161 if (p == 0.0)
continue;
1162 if (p != 1.0)
return true;
1172 for (std::size_t e = 0; e < d.params.size(); ++e) tot += d.params[e];
1173 const double u = g_fork.aux.next_double() * tot;
1175 for (std::size_t e = 0; e < d.params.size(); ++e) {
1176 accp += d.params[e];
1178 return static_cast<int>((d.trace.empty() ?
static_cast<double>(e + 1) : d.trace[e]) +
1181 return static_cast<int>(
1182 (d.trace.empty() ?
static_cast<double>(d.params.size()) : d.trace.back()) + 0.5);
1186 std::vector<StationState> S(M);
1187 std::size_t source_st = M;
1188 std::vector<int> classprio(K, 0);
1189 std::vector<double> classdeadline(K, std::numeric_limits<double>::infinity());
1190 for (std::size_t r = 0; r < K; ++r) {
1191 classprio[r] =
sn.classes[r].prio;
1192 classdeadline[r] =
sn.classes[r].deadline;
1195 for (std::size_t i = 0; i < M; ++i) {
1196 const auto& st =
sn.stations[i];
1197 StationState& s = S[i];
1199 s.off.assign(K,
true);
1200 s.class_mean.assign(K, 0.0);
1201 s.weight.assign(K, 1.0);
1202 for (std::size_t r = 0; r < K; ++r)
1205 if (st.nodetype == NodeType::Source) {
1206 s.role = Role::Source;
1208 }
else if (st.nodetype == NodeType::Fork || st.nodetype == NodeType::Join ||
1209 st.nodetype == NodeType::Place || st.nodetype == NodeType::Transition) {
1212 s.role = Role::Synchronization;
1218 if (st.nodetype == NodeType::Join) s.off.assign(K,
false);
1229 if (st.nodetype == NodeType::Place) {
1230 s.off.assign(K,
false);
1231 if (std::isfinite(st.nservers))
1232 s.nservers =
static_cast<std::size_t
>(st.nservers + 0.5);
1234 }
else if (st.nodetype == NodeType::Delay || st.sched == SchedStrategy::INF) {
1235 s.role = Role::Delay;
1237 s.role = Role::Queue;
1238 s.ps = is_ps_family(st.sched);
1239 s.preemptive = is_preemptive(st.sched);
1240 s.resume = is_preemptive_resume(st.sched);
1242 s.size_ordered = (st.sched == SchedStrategy::SJF ||
1243 st.sched == SchedStrategy::LJF ||
1244 st.sched == SchedStrategy::SRPT ||
1245 st.sched == SchedStrategy::SRPTPRIO ||
1246 st.sched == SchedStrategy::PSJF ||
1247 st.sched == SchedStrategy::LRPT ||
1248 st.sched == SchedStrategy::FSP);
1249 const double c = st.nservers;
1250 if (!(c >= 1.0) || !std::isfinite(c))
1251 throw InputError(
"SolverLDES (native engine): station '" + st.name +
1252 "' has a server count that is neither finite nor at least one");
1253 s.nservers =
static_cast<std::size_t
>(c + 0.5);
1255 s.classcap.assign(K, std::numeric_limits<double>::infinity());
1256 s.blocked_at.assign(K, 0.0);
1258 for (std::size_t r = 0; r < K; ++r)
1259 if (i <
sn.droprule.size() && r <
sn.droprule[i].size())
1260 s.droprule[r] =
sn.droprule[i][r];
1261 for (std::size_t r = 0; r < K; ++r)
1262 if (i <
sn.classcap.size() && r <
sn.classcap[i].size())
1263 s.classcap[r] =
sn.classcap[i][r];
1264 if (st.sched == SchedStrategy::LPS) s.lps_limit = s.nservers;
1268 s.parallelism.assign(K, 1);
1269 for (std::size_t r = 0; r < K && r < st.server_parallelism.size(); ++r)
1270 if (st.server_parallelism[r] > 1) {
1271 s.parallelism[r] = st.server_parallelism[r];
1272 s.has_parallelism =
true;
1278 s.ps_cd.assign(K, 1.0);
1279 s.gd_peak.assign(K, 1.0);
1285 s.util_peak =
static_cast<double>(s.nservers);
1286 for (
double a : s.lld) s.util_peak = std::max(s.util_peak, a);
1294 if (st.cdscaling || st.jdscaling) {
1304 double dep_peak = 1.0;
1305 const std::vector<T>* pks[2] = {&st.cdscalingpeak, &st.jdscalingpeak};
1306 const bool on[2] = {
static_cast<bool>(st.cdscaling),
static_cast<bool>(st.jdscaling)};
1307 const char* names[2] = {
"setClassDependence",
"setJointDependence"};
1308 for (std::size_t h = 0; h < 2; ++h) {
1309 if (!on[h])
continue;
1310 if (pks[h]->empty())
1312 std::string(
"SolverLDES: station '") + st.name +
1313 "' declares a dependent scaling with no declared peak rate. Utilization "
1314 "there is T*E[S]/peak, so pass the peak to " + names[h]);
1318 throw InputError(std::string(
"SolverLDES: station '") + st.name +
1319 "' declares a non-positive peak rate for " + names[h]);
1322 s.util_peak = dep_peak;
1323 s.cd_peak_class.assign(
sn.nclasses, 1.0);
1324 for (std::size_t r = 0; r <
sn.nclasses; ++r) {
1325 std::size_t cr =
sn.chains.size();
1326 for (std::size_t c = 0; c <
sn.chains.size(); ++c)
1327 if (r <
sn.chains[c].size() &&
sn.chains[c][r]) { cr = c;
break; }
1328 for (std::size_t h = 0; h < 2; ++h) {
1329 if (!on[h])
continue;
1331 for (std::size_t k = 0; k <
sn.nclasses; ++k) {
1332 const bool same = (cr <
sn.chains.size()) ? (k <
sn.chains[cr].size() &&
sn.chains[cr][k])
1334 if (!same)
continue;
1335 const std::size_t e = pks[h]->size() > 1 ? k : 0;
1338 s.cd_peak_class[r] *= (m > 0.0) ? m : 1.0;
1341 const auto& beta = st.cdscaling;
1342 const auto& eta = st.jdscaling;
1343 s.cd = [beta, eta](
const std::vector<double>& n) {
1344 std::vector<T> nt(n.size());
1346 std::vector<double> bd, ed;
1348 const std::vector<T> b = beta(nt);
1349 bd.resize(b.size());
1353 const std::vector<T> e = eta(nt);
1354 ed.resize(e.size());
1357 if (bd.empty())
return ed;
1358 if (ed.empty())
return bd;
1361 const std::size_t n2 = std::max(bd.size(), ed.size());
1362 std::vector<double> out(n2, 1.0);
1363 for (std::size_t k = 0; k < n2; ++k)
1364 out[k] = bd[bd.size() > 1 ? k : 0] * ed[ed.size() > 1 ? k : 0];
1368 if (st.sched == SchedStrategy::PAS || st.sched == SchedStrategy::OI) {
1370 auto pp =
sn.pasparam.find(i + 1);
1371 if (pp !=
sn.pasparam.end()) {
1372 if (pp->second.svc_rate_fun) {
1373 const auto f = pp->second.svc_rate_fun;
1382 s.pas_rate = [f](
const std::vector<std::size_t>& seq) {
1383 std::vector<std::size_t> tags(seq.size());
1384 for (std::size_t k = 0; k < seq.size(); ++k) tags[k] = seq[k] + 1;
1388 s.pas_swap = pp->second.swap_graph;
1391 throw InputError(
"SolverLDES (native engine): station '" + st.name +
1392 "' is a pass-and-swap station with no service rate function");
1394 if (st.sched == SchedStrategy::POLLING) {
1396 const auto pp =
sn.effective_polling(i + 1);
1397 s.poll_type = pp.ptype;
1398 s.poll_k = (pp.pk >= 1) ? pp.pk : 1;
1400 ? std::numeric_limits<std::size_t>::max()
1402 s.poll_parked =
true;
1403 s.switchover.resize(K);
1404 s.has_switchover.assign(K,
false);
1405 for (std::size_t r = 0; r < K && r < pp.switchover.size(); ++r)
1406 if (!pp.switchover[r].disabled &&
1408 s.switchover[r] = Sampler(pp.switchover[r],
1409 "the switchover of station '" + st.name +
1410 "', class '" +
sn.classes[r].name +
"'");
1411 s.has_switchover[r] =
true;
1414 s.retrial.resize(K);
1415 s.has_retrial.assign(K,
false);
1416 s.max_attempts.assign(K, 0);
1417 s.retrial_policy.assign(K, 1);
1418 s.const_pending.assign(K,
false);
1419 s.orbit.assign(K, std::deque<Job>());
1420 s.orbit_abandon.resize(K);
1421 s.has_orbit_abandon.assign(K,
false);
1422 s.orbit_size.assign(K, 0.0);
1423 s.tot_orbit.assign(K, 0.0);
1424 s.retried.assign(K, 0.0);
1425 s.retrial_lost.assign(K, 0.0);
1427 auto rp =
sn.retrialparam.find(i + 1);
1428 if (rp !=
sn.retrialparam.end())
1429 for (std::size_t r = 0; r < K && r < rp->second.retrial_proc.size(); ++r)
1430 if (!rp->second.retrial_proc[r].disabled) {
1431 s.retrial[r] = Sampler(rp->second.retrial_proc[r],
1432 "the retrial process of station '" + st.name +
1433 "', class '" +
sn.classes[r].name +
"'");
1434 s.has_retrial[r] =
true;
1435 if (r < rp->second.max_attempts.size())
1436 s.max_attempts[r] = rp->second.max_attempts[r];
1437 if (r < rp->second.policy.size())
1438 s.retrial_policy[r] = rp->second.policy[r];
1439 if (r < st.orbit_impatience.size() && !st.orbit_impatience[r].disabled) {
1440 s.orbit_abandon[r] = Sampler(st.orbit_impatience[r],
1441 "the orbit impatience of station '" + st.name +
1442 "', class '" +
sn.classes[r].name +
"'");
1443 s.has_orbit_abandon[r] =
true;
1447 s.patience.resize(K);
1448 s.has_patience.assign(K,
false);
1449 s.balk.assign(K, std::vector<StationState::BalkRule>());
1450 for (std::size_t r = 0; r < K; ++r) {
1452 r < st.patience.size() && !st.patience[r].disabled) {
1453 s.patience[r] = Sampler(st.patience[r],
"the patience of station '" + st.name +
1454 "', class '" +
sn.classes[r].name +
"'");
1455 s.has_patience[r] =
true;
1461 if (r < st.balking.size() &&
1463 for (
const auto& th : st.balking[r].thresholds) {
1464 StationState::BalkRule br;
1465 br.min_jobs = th.min_jobs;
1466 br.max_jobs = th.max_jobs;
1468 s.balk[r].push_back(br);
1474 auto sp =
sn.setupparam.find(i + 1);
1475 if (sp !=
sn.setupparam.end()) {
1477 if (sp->second.last(su, doff) && !su.
disabled) {
1479 s.setup_time = Sampler(su,
"the setup time of station '" + st.name +
"'");
1486 Sampler(doff,
"the delay-off time of station '" + st.name +
"'");
1488 s.delayoff_time = Sampler();
1493 typename std::map<std::size_t, qn::BreakdownParam<T> >::const_iterator bp =
1494 sn.breakdownparam.find(i + 1);
1495 if (bp !=
sn.breakdownparam.end()) {
1496 s.has_breakdown =
true;
1499 Sampler(bp->second.failure,
"the failure time of station '" + st.name +
"'");
1501 Sampler(bp->second.repair,
"the repair time of station '" + st.name +
"'");
1502 s.down_scale.assign(K, 0.0);
1503 s.down_rate_raw.assign(K, 0.0);
1504 for (std::size_t r = 0; r < K && r < bp->second.down_service_rates.size(); ++r)
1510 s.state_dependent = (!s.lld.empty() || s.has_cd || s.has_breakdown);
1513 if (s.role == Role::Synchronization)
continue;
1514 for (std::size_t r = 0; r < K; ++r) {
1515 s.off[r] =
sn.disabled[i][r];
1516 if (s.off[r])
continue;
1517 const std::string where =
"station '" + st.name +
"', class '" +
sn.classes[r].name +
"'";
1518 s.svc[r] = Sampler(
sn.service[i][r], where);
1519 s.class_mean[r] = s.svc[r].mean();
1520 if (!(s.class_mean[r] > 0.0))
1521 throw InputError(
"SolverLDES (native engine): " + where +
1522 " has a non-positive mean service time");
1534 if (s.role != Role::Source) {
1535 std::size_t nbatch = 0, enabled = 0, bcls = K;
1536 for (std::size_t r = 0; r < K; ++r) {
1537 if (s.off[r])
continue;
1546 const std::string at = fam +
" service (BMSP) at station '" + st.name +
"'";
1547 if (s.role == Role::Delay)
1549 " is not supported at an infinite-server station; "
1550 "batch service requires a single-server FCFS queue");
1551 if (s.nservers != 1)
1553 " requires a single server; multiserver bulk service "
1554 "is not supported");
1557 " requires FCFS scheduling");
1560 " requires exactly one job class; multiclass bulk "
1561 "service is not supported");
1562 if (s.pas || s.preemptive || s.polling)
1564 " is incompatible with PAS/preemptive/polling "
1568 " is incompatible with server setup/delayoff");
1569 if (s.state_dependent)
1571 " is not supported with load or class dependence or "
1572 "breakdowns: the rate rescaling cannot be applied to "
1573 "a firing clock armed at the instant the station "
1575 if (!st.server_types.empty())
1577 " is not supported with heterogeneous server types; a "
1578 "bulk server is a single server by construction");
1593 if (!st.server_types.empty()) {
1594 const std::size_t nT = st.server_types.size();
1596 s.hetero_policy = st.hetero_policy;
1597 s.type_count.assign(nT, 0);
1598 s.type_first.assign(nT, 0);
1599 s.type_compat.assign(nT, std::vector<bool>(K,
true));
1600 s.type_svc.assign(nT, std::vector<Sampler>(K));
1601 s.type_has_svc.assign(nT, std::vector<bool>(K,
false));
1602 s.type_rate.assign(nT, std::vector<double>(K, 0.0));
1603 std::size_t total = 0;
1604 for (std::size_t t = 0; t < nT; ++t) {
1606 const double c = pt.
count;
1608 throw InputError(
"SolverLDES (native engine): station '" + st.name +
1609 "' declares server pool '" + pt.
name +
1610 "' with fewer than one server");
1611 s.type_first[t] = total;
1612 s.type_count[t] =
static_cast<std::size_t
>(c + 0.5);
1613 total += s.type_count[t];
1614 for (std::size_t r = 0; r < K; ++r) {
1620 if (s.off[r]) s.type_compat[t][r] =
false;
1622 const std::string where =
"station '" + st.name +
"', pool '" + pt.
name +
1623 "', class '" +
sn.classes[r].name +
"'";
1624 s.type_svc[t][r] = Sampler(pt.
service[r], where);
1625 s.type_has_svc[t][r] =
true;
1626 const double mu = s.type_svc[t][r].mean();
1628 throw InputError(
"SolverLDES (native engine): " + where +
1629 " has a non-positive mean service time");
1630 s.type_rate[t][r] = 1.0 / mu;
1631 }
else if (!s.off[r]) {
1633 (s.class_mean[r] > 0.0) ? 1.0 / s.class_mean[r] : 0.0;
1637 for (std::size_t r = 0; r < K; ++r) {
1638 if (s.off[r])
continue;
1639 bool served =
false;
1640 for (std::size_t t = 0; t < nT && !served; ++t) served = s.type_compat[t][r];
1642 throw InputError(
"SolverLDES (native engine): station '" + st.name +
1643 "' declares no server pool compatible with class '" +
1644 sn.classes[r].name +
1645 "', so a job of that class would wait forever");
1651 s.util_peak =
static_cast<double>(total);
1653 s.server_type.assign(total, 0);
1654 for (std::size_t t = 0; t < nT; ++t)
1655 for (std::size_t j = 0; j < s.type_count[t]; ++j)
1656 s.server_type[s.type_first[t] + j] = t;
1657 s.type_order.resize(nT);
1658 for (std::size_t t = 0; t < nT; ++t) s.type_order[t] = t;
1661 s.alfs_order = s.type_order;
1662 std::stable_sort(s.alfs_order.begin(), s.alfs_order.end(),
1663 [&](std::size_t a, std::size_t b) {
1664 std::size_t ca = 0, cb = 0;
1665 for (std::size_t r = 0; r < K; ++r) {
1666 if (s.type_compat[a][r]) ++ca;
1667 if (s.type_compat[b][r]) ++cb;
1676 if (s.has_breakdown) {
1677 s.down_scale.assign(K, 0.0);
1678 for (std::size_t r = 0; r < K; ++r) {
1679 if (s.off[r] || !(s.class_mean[r] > 0.0))
continue;
1680 if (r < s.down_rate_raw.size() && s.down_rate_raw[r] > 0.0)
1681 s.down_scale[r] = s.down_rate_raw[r] * s.class_mean[r];
1685 s.cmp.sched = s.sched;
1686 s.cmp.class_mean = &s.class_mean;
1687 if (s.role == Role::Queue) {
1688 s.server.assign(s.nservers, Job());
1689 s.server_busy.assign(s.nservers,
false);
1690 s.server_start.assign(s.nservers, 0.0);
1691 s.server_tag.assign(s.nservers, 0);
1692 s.server_blocked.assign(s.nservers,
false);
1693 s.server_held.assign(s.nservers,
false);
1694 s.held_cls.assign(s.nservers, 0);
1695 s.blocked_job.assign(s.nservers, Job());
1696 s.blocked_dest.assign(s.nservers, 0);
1697 s.blocked_dest_cls.assign(s.nservers, 0);
1698 if (s.has_parallelism) s.held_extra.assign(s.nservers, std::vector<std::size_t>());
1714 for (std::size_t i = 0; i < M; ++i) {
1715 const auto& st =
sn.stations[i];
1716 std::size_t maxpar = 1;
1717 for (std::size_t r = 0; r < st.server_parallelism.size(); ++r)
1718 maxpar = std::max(maxpar, st.server_parallelism[r]);
1719 if (maxpar <= 1)
continue;
1720 StationState& s = S[i];
1721 const std::string where = std::string(
"SolverLDES (native engine): station '") +
1722 st.name +
"' declares server parallelism";
1723 if (s.role != Role::Queue)
1725 "wait for a slot; drop the declaration or give the "
1726 "station a finite server count");
1729 ", which shares its whole capacity rather than handing out "
1730 "slots, so a job cannot hold n of them; use a queueing "
1731 "discipline, or a class weight to slow the class down");
1732 if (s.polling || s.pas)
1734 ", whose controller picks the server itself and has no pool "
1735 "to seize n slots from");
1737 throw UnsupportedError(ErrorKind{
"ServerParallelism"}, where +
" at a bulk (BMSP) server, which has only a "
1738 "station-level firing clock and no per-job server "
1741 throw UnsupportedError(ErrorKind{
"ServerParallelism"}, where +
" together with heterogeneous server pools: a job is "
1742 "served at the law of the one pool it occupies, which "
1743 "n servers of mixed types cannot express");
1745 throw UnsupportedError(ErrorKind{
"ServerParallelism"}, where +
" together with setup/delayoff: the slots of such a "
1746 "station power up one at a time, so 'n free at once' "
1747 "is not a state it passes through");
1748 if (maxpar > s.nservers)
1749 throw InputError(where +
" of " + std::to_string(maxpar) +
" above its " +
1750 std::to_string(s.nservers) +
" servers, so the job could never "
1767 const bool has_gd =
static_cast<bool>(sn.gdscaling);
1768 std::vector<std::vector<double> > gd_cache(M, std::vector<double>(K, 1.0));
1770 if (sn.gdscalingpeak.size() < M * K)
1772 "SolverLDES (native engine): a global dependence requires an explicit peak rate, "
1773 "because utilization under it is reported as T*E[S]/peak; pass one to "
1774 "setGlobalDependence");
1775 for (std::size_t i = 0; i < M; ++i) {
1776 StationState& s = S[i];
1777 const std::string& nm = sn.stations[i].name;
1780 "SolverLDES (native engine): station '" + nm +
1781 "' is a bulk (BMSP) server under a global dependence, which holds only a "
1782 "station-level firing clock and no per-job residual work to rescale");
1785 "SolverLDES (native engine): station '" + nm +
1786 "' is a pass-and-swap station under a global dependence, whose single "
1787 "aggregate clock is armed from its own list and carries no residual work "
1789 for (std::size_t r = 0; r < K; ++r) {
1790 const double pk = num_traits<T>::to_double(sn.gdscalingpeak[i * K + r]);
1791 if (!(pk > 0.0) || !std::isfinite(pk))
1792 throw InputError(
"SolverLDES (native engine): the global dependence declares "
1793 "a peak rate that is not finite and positive at station '" +
1801 if (s.role == Role::Queue) s.state_dependent =
true;
1811 const std::size_t nnodes = sn.nof_nodes();
1812 std::vector<std::size_t> node_to_station(nnodes + 1, M);
1813 for (std::size_t i = 0; i < M; ++i) node_to_station[sn.station_to_node[i]] = i;
1826 std::vector<std::vector<std::vector<RouteDest>>> nroute(
1827 nnodes + 1, std::vector<std::vector<RouteDest>>(K));
1836 std::vector<std::vector<lang::RoutingStrategy>> node_routing(
1838 std::vector<std::vector<std::size_t>> rr_counter(nnodes + 1, std::vector<std::size_t>(K, 0));
1839 for (std::size_t inode = 1; inode <= nnodes; ++inode)
1840 for (std::size_t r = 0; r < K; ++r)
1841 if (sn.nodes[inode - 1].routing.size() > r)
1842 node_routing[inode][r] = sn.nodes[inode - 1].routing[r];
1843 for (std::size_t inode = 1; inode <= nnodes; ++inode) {
1850 if (sn.nodes[inode - 1].nodetype == NodeType::Sink)
continue;
1854 if (sn.nodes[inode - 1].nodetype == NodeType::Place ||
1855 sn.nodes[inode - 1].nodetype == NodeType::Transition)
1857 for (std::size_t r = 0; r < K; ++r) {
1860 for (std::size_t j = 1; j <= nnodes; ++j) {
1863 d.sink = (sn.nodes[j - 1].nodetype == NodeType::Sink);
1864 d.station = node_to_station[j];
1865 for (std::size_t s2 = 0; s2 < K; ++s2) {
1866 const double p = num_traits<T>::to_double(sn.route_eff(r + 1, s2 + 1, inode, j));
1867 if (!(p > 0.0))
continue;
1869 d.cls.push_back(std::make_pair(s2, p));
1871 if (d.cls.empty())
continue;
1878 if (!d.sink && d.station >= M && sn.nodes[j - 1].nodetype != NodeType::Fork &&
1879 sn.nodes[j - 1].nodetype != NodeType::Join &&
1880 sn.nodes[j - 1].nodetype != NodeType::Cache &&
1881 sn.nodes[j - 1].nodetype != NodeType::ClassSwitch &&
1882 sn.nodes[j - 1].nodetype != NodeType::Router &&
1883 sn.nodes[j - 1].nodetype != NodeType::Logger &&
1884 sn.nodes[j - 1].nodetype != NodeType::Place &&
1885 sn.nodes[j - 1].nodetype != NodeType::Transition)
1887 "SolverLDES (native engine): the routing crosses node '" +
1888 sn.nodes[j - 1].name +
1889 "', which the refresh did not fold into the station-to-station "
1891 if (d.station < M && S[d.station].role == Role::Source)
1892 throw InputError(
"SolverLDES (native engine): the routing sends a "
1893 "job back into the Source");
1894 nroute[inode][r].push_back(d);
1904 std::vector<int> join_of_fork(nnodes + 1, -1);
1905 std::vector<bool> is_fork(nnodes + 1,
false), is_join(nnodes + 1,
false);
1906 for (
const auto& fjp : sn.fj) {
1907 if (fjp.first > nnodes || fjp.second > nnodes)
continue;
1908 is_fork[fjp.first] =
true;
1909 is_join[fjp.second] =
true;
1910 join_of_fork[fjp.first] =
static_cast<int>(fjp.second);
1914 std::vector<bool> spawn_joins_at(K,
false);
1915 for (std::size_t inode = 1; inode <= nnodes; ++inode)
1916 for (std::size_t r = 0; r < K; ++r)
1917 for (
const RouteDest& d : nroute[inode][r])
1918 if (d.node >= 1 && d.node <= nnodes &&
1919 sn.nodes[d.node - 1].nodetype == NodeType::Join)
1920 for (std::size_t c = 0; c < d.cls.size(); ++c)
1921 spawn_joins_at[d.cls[c].first] =
true;
1922 for (std::size_t nd = 1; nd <= nnodes; ++nd) {
1923 if (sn.nodes[nd - 1].nodetype == NodeType::Fork && !is_fork[nd])
1924 throw InputError(
"SolverLDES (native engine): Fork node '" + sn.nodes[nd - 1].name +
1925 "' has no matching Join");
1926 if (sn.nodes[nd - 1].nodetype == NodeType::Join && !is_join[nd])
1927 throw InputError(
"SolverLDES (native engine): Join node '" + sn.nodes[nd - 1].name +
1928 "' closes no Fork");
1940 std::size_t cls = 0;
1952 std::vector<std::pair<std::size_t, double>> siblings;
1954 std::uint64_t outer = 0;
1955 double t_fork = 0.0;
1958 std::map<std::uint64_t, ForkSync> fork_sync;
1959 std::uint64_t next_parent = 0;
1964 std::vector<int> fj_pair_of_fork(nnodes + 1, -1);
1965 for (std::size_t p = 0; p < fj_pairs.size(); ++p)
1966 if (fj_pairs[p].first <= nnodes &&
1967 join_of_fork[fj_pairs[p].first] ==
static_cast<int>(fj_pairs[p].second))
1968 fj_pair_of_fork[fj_pairs[p].first] =
static_cast<int>(p);
1969 std::vector<std::vector<double>> fj_resp_sum(fj_pairs.size(), std::vector<double>(K, 0.0));
1970 std::vector<std::vector<double>> fj_resp_cnt(fj_pairs.size(), std::vector<double>(K, 0.0));
1971 std::vector<std::vector<std::vector<double>>> fj_release(fj_pairs.size(),
1972 std::vector<std::vector<double>>(K));
1980 std::map<std::size_t, CacheState> caches;
1981 for (
const auto& np : sn.nodeparam) {
1982 const std::size_t nd = np.first;
1983 if (nd == 0 || nd > nnodes)
continue;
1984 if (sn.nodes[nd - 1].nodetype != NodeType::Cache)
continue;
1985 const auto& cp = np.second;
1988 cs.nitems = cp.nitems;
1989 cs.policy = cp.replacestrat;
1990 for (
int c : cp.itemcap)
1991 if (c > 0) cs.capacity.push_back(
static_cast<std::size_t
>(c));
1992 if (cs.capacity.empty()) cs.capacity.push_back(1);
1993 cs.lists.assign(cs.capacity.size(), std::list<std::size_t>());
1994 for (
int v : cp.itemsize) cs.item_size.push_back(v);
1995 for (
int v : cp.costcap) cs.cost_cap.push_back(v);
1996 cs.hits.assign(K, 0.0);
1997 cs.misses.assign(K, 0.0);
1998 cs.hits_per_list.assign(K, std::vector<double>(cs.lists.size(), 0.0));
1999 cs.item_level_time.assign(cs.nitems, std::vector<double>(cs.lists.size(), 0.0));
2000 cs.hit_class.assign(K, -1);
2001 cs.miss_class.assign(K, -1);
2002 for (std::size_t r = 0; r < K && r < cp.hitclass.size(); ++r)
2003 if (cp.hitclass[r] > 0) cs.hit_class[r] =
static_cast<int>(cp.hitclass[r] - 1);
2004 for (std::size_t r = 0; r < K && r < cp.missclass.size(); ++r)
2005 if (cp.missclass[r] > 0) cs.miss_class[r] =
static_cast<int>(cp.missclass[r] - 1);
2008 cs.popularity.assign(K, std::vector<double>());
2009 for (std::size_t r = 0; r < K && r < cp.pread.size(); ++r) {
2011 for (
const T& v : cp.pread[r]) {
2012 acc += num_traits<T>::to_double(v);
2013 cs.popularity[r].push_back(acc);
2015 if (!cs.popularity[r].empty()) cs.popularity[r].back() = 1.0;
2019 if (cp.retrieval_capacity > 0 && !cp.retrieval_classes.empty()) {
2020 cs.has_retrieval =
true;
2021 cs.retrieval_class = cp.retrieval_classes;
2022 cs.retrieval_class_to_item.assign(K, -1);
2023 for (std::size_t it = 0; it < cs.retrieval_class.size(); ++it)
2024 for (std::size_t r = 0; r < cs.retrieval_class[it].size(); ++r) {
2025 const std::size_t rc = cs.retrieval_class[it][r];
2026 if (rc > 0 && rc <= K) cs.retrieval_class_to_item[rc - 1] =
static_cast<int>(it);
2028 cs.in_flight.assign(cs.nitems, 0);
2029 cs.fetch_start.assign(cs.nitems, 0.0);
2030 cs.held.assign(cs.nitems, std::vector<CacheState::HeldRequest>());
2031 cs.delayed.assign(K, 0.0);
2041 std::vector<std::size_t> place_of(nnodes + 1, M);
2042 std::vector<std::size_t> place_nodes;
2043 for (std::size_t nd = 1; nd <= nnodes; ++nd)
2044 if (sn.nodes[nd - 1].nodetype == NodeType::Place) {
2045 place_of[nd] = place_nodes.size();
2046 place_nodes.push_back(nd);
2048 std::vector<std::size_t> place_station(place_nodes.size(), M);
2049 for (std::size_t p = 0; p < place_nodes.size(); ++p)
2050 place_station[p] = node_to_station[place_nodes[p]];
2063 std::vector<std::vector<double>> marking(place_nodes.size(), std::vector<double>(K, 0.0));
2064 for (std::size_t p = 0; p < place_nodes.size(); ++p) {
2065 auto im = sn.initmarking.find(place_nodes[p]);
2066 if (im != sn.initmarking.end())
2067 for (std::size_t r = 0; r < K && r < im->second.size(); ++r)
2068 marking[p][r] = num_traits<T>::to_double(im->second[r]);
2069 for (std::size_t r = 0; r < K; ++r) {
2071 const double n = sn.classes[r].population;
2072 if (!std::isfinite(n) || !(n > 0.0))
continue;
2073 if (sn.classes[r].refstat != place_station[p] + 1)
continue;
2074 marking[p][r] = std::max(marking[p][r], n);
2096 std::vector<std::vector<double>> avail(place_nodes.size(), std::vector<double>(K, 0.0));
2097 std::vector<bool> qplace(place_nodes.size(),
false);
2098 std::vector<std::size_t> place_servers(place_nodes.size(), 1);
2099 std::vector<std::deque<std::size_t>> place_wait(place_nodes.size());
2100 std::vector<std::size_t> place_busy(place_nodes.size(), 0);
2101 std::vector<std::vector<double>> place_insvc(place_nodes.size(), std::vector<double>(K, 0.0));
2102 std::vector<std::vector<Sampler>> place_svc(place_nodes.size(), std::vector<Sampler>(K));
2114 std::vector<std::vector<Rng>> g_place(place_nodes.size());
2115 for (std::size_t p = 0; p < place_nodes.size(); ++p) {
2116 const std::size_t ist = place_station[p];
2117 g_place[p].reserve(K);
2118 for (std::size_t k = 0; k < K; ++k) {
2119 const long long off =
2120 ((nsrc + nsvc +
static_cast<long long>(p)) * Kll +
static_cast<long long>(k)) *
2122 g_place[p].push_back(Rng(seed_ll, off));
2124 if (ist >= M)
continue;
2125 qplace[p] = sn.is_queueing_place(ist);
2127 for (std::size_t r = 0; r < K; ++r) avail[p][r] = marking[p][r];
2130 const auto& pst = sn.stations[ist];
2136 if (pst.sched != SchedStrategy::FCFS && pst.sched != SchedStrategy::LCFS &&
2137 pst.sched != SchedStrategy::SIRO && pst.sched != SchedStrategy::INF)
2139 "SolverLDES (native engine): queueing place '" + pst.name +
2141 ", which the embedded-queue algorithm does not support "
2142 "(supported: FCFS, LCFS, SIRO, INF)");
2143 place_servers[p] = (pst.sched == SchedStrategy::INF || !std::isfinite(pst.nservers))
2144 ? std::numeric_limits<std::size_t>::max()
2145 :
static_cast<std::size_t
>(pst.nservers + 0.5);
2146 for (std::size_t r = 0; r < K; ++r) {
2147 if (sn.service[ist][r].disabled)
continue;
2149 Sampler(sn.service[ist][r],
2150 "queueing place '" + pst.name +
"', class '" + sn.classes[r].name +
"'");
2156 for (std::size_t k = 0; k < K; ++k)
2157 for (
double t = 0.0; t + 0.5 < marking[p][k]; t += 1.0) place_wait[p].push_back(k);
2160 std::vector<SpnTransition> transitions;
2162 std::vector<std::vector<Sampler>> spn_law;
2163 for (
const auto& tp : sn.transparam) {
2164 const std::size_t nd = tp.first;
2165 if (nd == 0 || nd > nnodes)
continue;
2168 tr.places = place_nodes;
2169 for (std::size_t m = 0; m < tp.second.nmodes; ++m) {
2173 for (std::size_t p = 0; p < place_nodes.size(); ++p) {
2174 const std::size_t pn = place_nodes[p];
2175 for (std::size_t r = 0; r < K; ++r) {
2176 mode.enabling.push_back(
2177 (m < tp.second.enabling.size() && pn - 1 < tp.second.enabling[m].rows() &&
2178 r < tp.second.enabling[m].cols())
2179 ? num_traits<T>::to_double(tp.second.enabling[m](pn - 1, r))
2181 mode.inhibiting.push_back(
2182 (m < tp.second.inhibiting.size() &&
2183 pn - 1 < tp.second.inhibiting[m].rows() &&
2184 r < tp.second.inhibiting[m].cols())
2185 ? num_traits<T>::to_double(tp.second.inhibiting[m](pn - 1, r))
2186 : std::numeric_limits<double>::infinity());
2187 mode.firing.push_back(
2188 (m < tp.second.firing.size() && pn - 1 < tp.second.firing[m].rows() &&
2189 r < tp.second.firing[m].cols())
2190 ? num_traits<T>::to_double(tp.second.firing[m](pn - 1, r))
2194 if (m < tp.second.firing.size())
2195 for (std::size_t nd2 = 1; nd2 <= nnodes && nd2 - 1 < tp.second.firing[m].rows(); ++nd2) {
2196 if (sn.nodes[nd2 - 1].nodetype == NodeType::Place)
continue;
2197 for (std::size_t r = 0; r < K && r < tp.second.firing[m].cols(); ++r) {
2198 const double f = num_traits<T>::to_double(tp.second.firing[m](nd2 - 1, r));
2199 if (f > 0.0) mode.exits.emplace_back(nd2, r, f);
2202 if (m < tp.second.nmodeservers.size()) mode.servers = tp.second.nmodeservers[m];
2203 if (m < tp.second.firingprio.size()) mode.priority = tp.second.firingprio[m];
2204 if (m < tp.second.fireweight.size())
2205 mode.weight = num_traits<T>::to_double(tp.second.fireweight[m]);
2206 mode.immediate = (m < tp.second.timing.size() &&
2208 if (m < tp.second.firingproc.size() && !tp.second.firingproc[m].disabled) {
2209 const double mean = num_traits<T>::to_double(tp.second.firingproc[m].mean);
2210 if (mean > 0.0) mode.rate = 1.0 / mean;
2212 if (m < tp.second.firingdep.size() && tp.second.firingdep[m]) {
2213 const auto f = tp.second.firingdep[m];
2214 const std::vector<std::size_t> pl = place_nodes;
2215 const std::size_t nn = nnodes, KK = K;
2218 mode.dep = [f, pl, nn, KK](
const std::vector<double>& tok) {
2219 std::vector<T> mk(nn * KK, num_traits<T>::from_int(0));
2220 for (std::size_t p = 0; p < pl.size(); ++p) {
2221 if (pl[p] < 1 || pl[p] - 1 >= nn)
continue;
2222 for (std::size_t r = 0; r < KK; ++r)
2223 if (p * KK + r < tok.size())
2224 mk[(pl[p] - 1) * KK + r] = num_traits<T>::from_double(tok[p * KK + r]);
2226 return num_traits<T>::to_double(f(mk));
2229 tr.modes.push_back(mode);
2231 tr.fired.assign(tr.modes.size(), 0.0);
2232 std::vector<Sampler> law(tr.modes.size());
2233 const std::string tname = sn.nodes[nd - 1].name;
2234 for (std::size_t m = 0; m < tr.modes.size(); ++m)
2235 if (!tr.modes[m].immediate && m < tp.second.firingproc.size() &&
2236 !tp.second.firingproc[m].disabled)
2237 law[m] = Sampler(tp.second.firingproc[m],
"the firing time of mode " +
2238 std::to_string(m + 1) +
" of transition '" + tname +
"'");
2239 spn_law.push_back(law);
2240 transitions.push_back(tr);
2245 for (std::size_t i = 0; i < M; ++i) acc.util_peak[i] = S[i].util_peak;
2250 for (std::size_t p = 0; p < place_nodes.size(); ++p) {
2251 if (place_station[p] >= M)
continue;
2252 for (std::size_t r = 0; r < K; ++r) acc.qlen[place_station[p]][r] = marking[p][r];
2257 std::vector<double> sys_resp_sum(K, 0.0), sys_resp_cnt(K, 0.0), sys_completed(K, 0.0);
2260 std::vector<std::vector<double>> dropped(M, std::vector<double>(K, 0.0));
2261 std::vector<std::vector<double>> balked(M, std::vector<double>(K, 0.0));
2262 std::vector<std::vector<double>> reneged(M, std::vector<double>(K, 0.0));
2265 std::vector<std::vector<double>> last_lost(M, std::vector<double>(K, -1.0));
2266 std::vector<std::vector<double>> blocked_count(M, std::vector<double>(K, 0.0));
2275 const bool want_respt = o.export_respt || o.export_trajectory;
2276 std::vector<std::vector<std::vector<double>>> resp_samples(
2277 want_respt ? M : 0, std::vector<std::vector<double>>(K));
2278 std::uint64_t job_id = 0;
2294 struct LoggerTrace {
2295 std::unique_ptr<std::ofstream> out;
2296 qn::NodeDef::LoggerParam par;
2298 double last_any = 0.0;
2299 std::vector<double> last_cls;
2301 std::map<std::size_t, LoggerTrace> logger_traces;
2302 std::string log_start_time;
2303 std::uint64_t log_walk_id = 0;
2304 for (std::size_t nd = 1; nd <= sn.nof_nodes(); ++nd) {
2305 if (sn.nodes[nd - 1].nodetype != NodeType::Logger)
continue;
2306 if (log_start_time.empty()) {
2307 const std::time_t tt = std::time(
nullptr);
2309 std::strftime(buf,
sizeof buf,
"%Y-%m-%d %H:%M:%S", std::localtime(&tt));
2310 log_start_time = buf;
2313 lt.par = sn.nodes[nd - 1].logger;
2314 lt.name = sn.nodes[nd - 1].name;
2315 lt.last_cls.assign(K, 0.0);
2316 const std::string file =
2317 lt.par.file_name.empty() ?
"logger_" + lt.name +
".csv" : lt.par.file_name;
2318 std::string dir = lt.par.file_path.empty() ? sn.log_path : lt.par.file_path;
2319 if (!dir.empty() && dir[dir.size() - 1] !=
'/') dir +=
'/';
2320 const std::string path = dir + file;
2321 lt.out.reset(
new std::ofstream(path.c_str(), std::ios::out | std::ios::trunc));
2324 std::fprintf(stderr,
"Warning: solver_ldes: failed to create logger file for %s: %s\n",
2325 lt.name.c_str(), path.c_str());
2328 *lt.out <<
"LOGGERNAME,TIMESTAMP,JOB_ID,CLASS_ID,INTERARRIVAL_SAMECLASS,"
2329 "INTERARRIVAL_ANYCLASS,SIMUL_START_TIME\n";
2330 logger_traces[nd] = std::move(lt);
2332 auto log_passage = [&](std::size_t nd, std::size_t r, std::uint64_t id) {
2333 typename std::map<std::size_t, LoggerTrace>::iterator it = logger_traces.find(nd);
2334 if (it == logger_traces.end())
return;
2335 LoggerTrace& lt = it->second;
2336 const double same = now - lt.last_cls[r], any = now - lt.last_any;
2337 lt.last_cls[r] = now;
2339 std::ostream& w = *lt.out;
2340 if (lt.par.logger_name) w << lt.name;
2344 if (lt.par.job_id) w << id;
2346 if (lt.par.job_class) w << sn.classes[r].name;
2352 if (lt.par.start_time) w << log_start_time;
2366 std::vector<bool> is_removal_signal(K,
false);
2367 bool has_removal_signal =
false;
2368 for (std::size_t r = 0; r < K; ++r) {
2369 if (r >= sn.issignal.size() || !sn.issignal[r])
continue;
2373 is_removal_signal[r] =
true;
2374 has_removal_signal =
true;
2385 const bool init_jobs_on = o.init_jobs.rows() > 0;
2386 const bool track_jobs = o.track_jobs || init_jobs_on;
2387 const bool track_delay_jobs = has_removal_signal || has_gd || track_jobs;
2389 std::vector<std::array<double, 4>> job_departure_log;
2400 std::vector<std::vector<Job>> delay_live(track_delay_jobs ? M : 0);
2411 std::vector<std::size_t> sync_reply(K, K);
2412 std::vector<bool> is_reply_signal(K,
false);
2413 bool has_sync_call =
false;
2414 for (std::size_t r = 0; r < K; ++r) {
2415 if (r < sn.syncreply.size() && sn.syncreply[r] >= 1 && sn.syncreply[r] <= K) {
2416 sync_reply[r] = sn.syncreply[r] - 1;
2417 has_sync_call =
true;
2419 if (r < sn.issignal.size() && sn.issignal[r] && r < sn.signaltype.size() &&
2421 is_reply_signal[r] =
true;
2424 struct PendingCall {
2425 std::size_t station = 0;
2426 std::size_t slot = 0;
2427 std::size_t cls = 0;
2430 std::map<std::uint64_t, PendingCall> pending_reply;
2431 std::uint64_t next_call = 0;
2434 bool has_immfeed =
false;
2435 for (std::size_t i = 0; i < M && i < sn.immfeed.size(); ++i)
2436 for (std::size_t r = 0; r < sn.immfeed[i].size(); ++r)
2437 if (sn.immfeed[i][r]) has_immfeed =
true;
2447 std::vector<std::size_t> spawn_of(K, K);
2448 bool has_spawn =
false;
2449 for (std::size_t r = 0; r < K; ++r)
2450 if (sn.classes[r].spawn >= 1 && sn.classes[r].spawn <= K) {
2451 spawn_of[r] = sn.classes[r].spawn - 1;
2456 std::priority_queue<Event, std::vector<Event>, EventLater> evq;
2457 std::uint64_t seq = 0, ps_tag = 0;
2459 auto push = [&](Event e) {
2477 std::function<void(std::size_t)> place_try_start = [&](std::size_t p) {
2478 const std::size_t ist = place_station[p];
2479 if (ist >= M)
return;
2480 while (place_busy[p] < place_servers[p] && !place_wait[p].empty()) {
2482 if (S[ist].sched == SchedStrategy::LCFS) {
2483 r = place_wait[p].back();
2484 place_wait[p].pop_back();
2485 }
else if (S[ist].sched == SchedStrategy::SIRO) {
2486 const std::size_t at = std::min(
2487 place_wait[p].size() - 1,
2488 static_cast<std::size_t
>(uniform01(g_qplace) *
2489 static_cast<double>(place_wait[p].size())));
2490 r = place_wait[p][at];
2491 place_wait[p].erase(place_wait[p].begin() +
static_cast<std::ptrdiff_t
>(at));
2493 r = place_wait[p].front();
2494 place_wait[p].pop_front();
2501 if (place_svc[p][r].disabled())
2502 throw InputError(
"SolverLDES (native engine): a token of class '" +
2503 sn.classes[r].name +
"' entered queueing place '" +
2504 sn.stations[ist].name +
2505 "', which has no service process for that class; call "
2506 "setService for every token colour of a queueing place");
2507 acc.update_busy(ist, r, now);
2509 place_insvc[p][r] += 1.0;
2510 acc.busy[ist][r] = place_insvc[p][r];
2512 e.t = now + place_svc[p][r].next(g_place[p][r]);
2521 auto place_deposit = [&](std::size_t p, std::size_t r,
double n) {
2522 const std::size_t ist = place_station[p];
2524 acc.update_qlen(ist, r, now);
2526 acc.qlen[ist][r] = marking[p][r];
2527 acc.arrived[ist][r] += n;
2535 for (
double t = 0.0; t + 0.5 < n; t += 1.0) place_wait[p].push_back(r);
2547 auto place_take = [&](std::size_t p, std::size_t r,
double n) {
2548 const std::size_t ist = place_station[p];
2552 acc.update_qlen(ist, r, now);
2553 acc.qlen[ist][r] = marking[p][r];
2554 acc.completed[ist][r] += n;
2567 auto buffer_push = [&](std::size_t i,
const Job& job) {
2568 StationState& s = S[i];
2569 s.buffer.push_back(job);
2570 if (s.sched != SchedStrategy::FSP)
2571 std::push_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
2573 auto buffer_pop = [&](std::size_t i) -> Job {
2574 StationState& s = S[i];
2575 if (s.sched != SchedStrategy::FSP) {
2576 std::pop_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
2577 Job job = s.buffer.back();
2578 s.buffer.pop_back();
2581 std::vector<double> works;
2582 for (
const Job& j : s.buffer) works.push_back(j.remaining);
2583 for (std::size_t sl = 0; sl < s.nservers; ++sl)
2584 if (s.server_busy[sl])
2585 works.push_back(std::max(0.0, s.server[sl].remaining -
2586 (now - s.server_start[sl])));
2587 std::size_t best = 0;
2588 double best_vft = std::numeric_limits<double>::infinity();
2589 for (std::size_t j = 0; j < s.buffer.size(); ++j) {
2591 static_cast<double>(s.nservers), now);
2597 Job job = s.buffer[best];
2598 s.buffer.erase(s.buffer.begin() +
static_cast<std::ptrdiff_t
>(best));
2613 auto slot_ok_for = [&](
const StationState& s, std::size_t sl, std::size_t cls) ->
bool {
2614 if (s.server_busy[sl] || s.server_blocked[sl] || s.server_held[sl])
return false;
2615 if (!s.has_pools)
return true;
2616 return s.type_compat[s.server_type[sl]][cls];
2620 auto free_slot_count = [&](std::size_t i) -> std::size_t {
2621 const StationState& s = S[i];
2623 for (std::size_t sl = 0; sl < s.server_busy.size(); ++sl)
2624 if (!s.server_busy[sl] && !s.server_blocked[sl] && !s.server_held[sl]) ++n;
2636 auto free_slot_for = [&](std::size_t i, std::size_t cls) -> std::size_t {
2637 StationState& s = S[i];
2641 if (s.has_parallelism && free_slot_count(i) < s.parallelism[cls])
return s.nservers;
2643 for (std::size_t sl = 0; sl < s.nservers; ++sl)
2644 if (slot_ok_for(s, sl, cls))
return sl;
2648 std::vector<std::size_t> cand;
2649 for (std::size_t t = 0; t < s.type_count.size(); ++t) {
2650 if (!s.type_compat[t][cls])
continue;
2651 for (std::size_t j = 0; j < s.type_count[t]; ++j)
2652 if (slot_ok_for(s, s.type_first[t] + j, cls)) {
2657 if (cand.empty())
return s.nservers;
2658 std::size_t chosen = cand[0];
2659 if (cand.size() > 1) {
2660 switch (s.hetero_policy) {
2666 for (std::size_t k = 0; k < s.type_order.size(); ++k) {
2667 const std::size_t t = s.type_order[k];
2668 if (std::find(cand.begin(), cand.end(), t) == cand.end())
continue;
2670 s.type_order.erase(s.type_order.begin() +
2671 static_cast<std::ptrdiff_t
>(k));
2672 s.type_order.push_back(t);
2680 for (std::size_t k = 0; k < s.alfs_order.size(); ++k)
2681 if (std::find(cand.begin(), cand.end(), s.alfs_order[k]) != cand.end()) {
2682 chosen = s.alfs_order[k];
2690 for (std::size_t k = 0; k < cand.size(); ++k) {
2691 const double rt = s.type_rate[cand[k]][cls];
2702 std::size_t at =
static_cast<std::size_t
>(
2703 uniform01(g_routing) *
static_cast<double>(cand.size()));
2704 if (at >= cand.size()) at = cand.size() - 1;
2716 for (std::size_t j = 0; j < s.type_count[chosen]; ++j) {
2717 const std::size_t sl = s.type_first[chosen] + j;
2718 if (slot_ok_for(s, sl, cls))
return sl;
2724 auto buffer_has_for_slot = [&](std::size_t i, std::size_t sl) ->
bool {
2725 StationState& s = S[i];
2726 if (s.buffer.empty())
return false;
2727 if (!s.has_pools)
return true;
2728 const std::vector<bool>& ok = s.type_compat[s.server_type[sl]];
2729 for (std::size_t j = 0; j < s.buffer.size(); ++j)
2730 if (ok[s.buffer[j].cls])
return true;
2745 auto buffer_pop_for_slot = [&](std::size_t i, std::size_t sl) -> Job {
2746 StationState& s = S[i];
2747 if (!s.has_pools)
return buffer_pop(i);
2748 const std::vector<bool>& ok = s.type_compat[s.server_type[sl]];
2749 std::vector<Job> skipped;
2752 while (!s.buffer.empty()) {
2753 Job cand = buffer_pop(i);
2759 skipped.push_back(cand);
2761 for (std::size_t j = 0; j < skipped.size(); ++j) buffer_push(i, skipped[j]);
2763 throw InputError(
"SolverLDES (native engine): buffer_pop_for_slot was asked for a "
2764 "job no pool of this slot can serve; guard with buffer_has_for_slot");
2787 auto lld_factor = [&](std::size_t i) ->
double {
2788 StationState& s = S[i];
2789 if (s.lld.empty() || s.ps)
return 1.0;
2791 for (std::size_t k = 0; k < K; ++k) total += acc.qlen[i][k];
2792 if (!(total > 0.0))
return 1.0;
2793 const std::size_t idx = std::min(
static_cast<std::size_t
>(total) - 1, s.lld.size() - 1);
2794 return s.lld[idx] > 0.0 ? s.lld[idx] : 1.0;
2804 auto cd_factor = [&](std::size_t i, std::size_t r) ->
double {
2805 StationState& s = S[i];
2806 if (!s.has_cd)
return 1.0;
2807 std::vector<double> nvec(K, 0.0);
2808 for (std::size_t k = 0; k < K; ++k) nvec[k] = acc.qlen[i][k];
2809 const std::vector<double> beta = s.cd(nvec);
2810 const double b = (beta.size() > 1) ? beta[r] : (beta.empty() ? 1.0 : beta[0]);
2811 return (b > 1e-10) ? b : 1.0;
2822 auto gd_refresh = [&]() {
2823 std::vector<T> nvec(M * K);
2824 for (std::size_t i = 0; i < M; ++i)
2825 for (std::size_t r = 0; r < K; ++r)
2826 nvec[i * K + r] = num_traits<T>::from_double(acc.qlen[i][r]);
2827 const std::vector<T> v = sn.gdscaling(nvec);
2829 throw InputError(
"SolverLDES (native engine): the global dependence handle returned "
2831 for (std::size_t i = 0; i < M; ++i)
2832 for (std::size_t r = 0; r < K; ++r) {
2835 const double f = num_traits<T>::to_double(
2836 (v.size() == 1) ? v[0] : ((v.size() == M) ? v[i] : v[i * K + r]));
2837 if (!std::isfinite(f) || f < 0.0)
2838 throw InputError(
"SolverLDES (native engine): the global dependence handle "
2839 "returned a scaling that is not finite and nonnegative at "
2840 "station '" + sn.stations[i].name +
"'");
2841 if (!(f > 0.0) && acc.qlen[i][r] > 0.0)
2843 "SolverLDES (native engine): the global dependence handle returned a "
2844 "zero scaling at station '" + sn.stations[i].name +
2845 "', which holds jobs of class '" + sn.classes[r].name +
2846 "'. A station that stops while it holds work has no completion left to "
2847 "schedule; use Queue.setBreakdown for a server that genuinely stops");
2852 auto rate_scaling = [&](std::size_t i, std::size_t r) ->
double {
2853 StationState& s = S[i];
2857 if (has_gd) factor *= gd_cache[i][r];
2862 if (s.has_breakdown && !s.up)
2863 factor *= (r < s.down_scale.size() ? s.down_scale[r] : 0.0);
2873 auto ps_servers = [&](std::size_t i) ->
double {
2874 StationState& s = S[i];
2875 if (s.lld.empty())
return static_cast<double>(s.nservers);
2877 for (std::size_t k = 0; k < K; ++k) total += acc.qlen[i][k];
2878 if (!(total > 0.0))
return 1.0;
2879 const std::size_t idx = std::min(
static_cast<std::size_t
>(total) - 1, s.lld.size() - 1);
2889 auto ps_advance_one = [&](std::size_t i) {
2890 StationState& s = S[i];
2891 const double dt = now - s.ps_last_update;
2892 s.ps_last_update = now;
2900 for (std::size_t r = 0; r < K; ++r) acc.last_busy[i][r] = now;
2901 if (!(dt > 0.0) || s.ps_jobs.empty())
return;
2902 std::vector<double> rates =
ps_shares(s.sched, s.ps_jobs, ps_servers(i), s.weight, K);
2904 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2905 rates[j] *= s.ps_cd[s.ps_jobs[j].cls];
2907 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2908 rates[j] *= gd_cache[i][s.ps_jobs[j].cls];
2909 if (s.has_breakdown && !s.up)
2910 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2911 rates[j] *= (s.ps_jobs[j].cls < s.down_scale.size() ? s.down_scale[s.ps_jobs[j].cls]
2913 for (std::size_t j = 0; j < s.ps_jobs.size(); ++j)
2914 if (rates[j] > 0.0) {
2915 acc.tot_busy[i][s.ps_jobs[j].cls] += rates[j] * dt;
2916 s.ps_jobs[j].remaining = std::max(0.0, s.ps_jobs[j].remaining - rates[j] * dt);
2919 auto ps_reschedule_one = [&](std::size_t i) {
2920 StationState& s = S[i];
2924 for (std::size_t r = 0; r < K; ++r) s.ps_cd[r] =
cd_factor(i, r);
2925 if (s.ps_jobs.empty())
return;
2926 std::vector<double> rates =
ps_shares(s.sched, s.ps_jobs, ps_servers(i), s.weight, K);
2928 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2929 rates[j] *= s.ps_cd[s.ps_jobs[j].cls];
2931 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2932 rates[j] *= gd_cache[i][s.ps_jobs[j].cls];
2933 if (s.has_breakdown && !s.up)
2934 for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
2935 rates[j] *= (s.ps_jobs[j].cls < s.down_scale.size() ? s.down_scale[s.ps_jobs[j].cls]
2937 for (std::size_t j = 0; j < s.ps_jobs.size(); ++j) {
2938 PsJob& pj = s.ps_jobs[j];
2941 const double rate = rates[j];
2942 if (!(rate > 0.0))
continue;
2944 e.t = now + ((pj.remaining <= 1e-12) ? 1e-12 : pj.remaining / rate);
2963 auto sd_advance_one = [&](std::size_t i) {
2964 StationState& s = S[i];
2965 if (!s.state_dependent || s.ps)
return;
2966 const double dt = now - s.sd_last_update;
2967 s.sd_last_update = now;
2968 if (!(dt > 0.0))
return;
2969 for (std::size_t sl = 0; sl < s.nservers; ++sl)
2970 if (s.server_busy[sl]) {
2971 const double scale = rate_scaling(i, s.server[sl].cls);
2972 s.server[sl].remaining = std::max(0.0, s.server[sl].remaining - dt * scale);
2973 s.server[sl].elapsed += dt * scale;
2974 s.server_start[sl] = now;
2977 auto sd_reschedule_one = [&](std::size_t i) {
2978 StationState& s = S[i];
2979 if (!s.state_dependent || s.ps)
return;
2985 for (std::size_t sl = 0; sl < s.nservers; ++sl)
2986 if (s.server_busy[sl]) {
2987 const double scale = rate_scaling(i, s.server[sl].cls);
2997 s.server_tag[sl] = ++ps_tag;
2998 if (!(scale > 0.0))
continue;
3000 e.t = now + s.server[sl].remaining / scale;
3003 e.cls = s.server[sl].cls;
3005 e.tag = s.server_tag[sl];
3020 auto gd_delay_advance = [&](std::size_t i) {
3021 StationState& s = S[i];
3022 const double dt = now - s.sd_last_update;
3023 s.sd_last_update = now;
3024 if (!(dt > 0.0))
return;
3025 for (std::size_t j = 0; j < delay_live[i].size(); ++j) {
3026 Job& d = delay_live[i][j];
3027 const double scale = gd_cache[i][d.cls];
3028 d.remaining = std::max(0.0, d.remaining - dt * scale);
3029 d.elapsed += dt * scale;
3032 auto gd_delay_reschedule = [&](std::size_t i) {
3033 for (std::size_t j = 0; j < delay_live[i].size(); ++j) {
3034 Job& d = delay_live[i][j];
3038 const double scale = gd_cache[i][d.cls];
3039 if (!(scale > 0.0))
continue;
3041 e.t = now + ((d.remaining <= 1e-12) ? 1e-12 : d.remaining / scale);
3045 e.slot = std::numeric_limits<std::size_t>::max();
3065 auto gd_advance_all = [&]() {
3066 for (std::size_t a = 0; a < M; ++a) {
3067 if (S[a].role == Role::Delay)
3068 gd_delay_advance(a);
3069 else if (S[a].role != Role::Queue)
3077 auto gd_reschedule_all = [&]() {
3079 for (std::size_t a = 0; a < M; ++a) {
3080 if (S[a].role == Role::Delay)
3081 gd_delay_reschedule(a);
3082 else if (S[a].role != Role::Queue)
3085 ps_reschedule_one(a);
3087 sd_reschedule_one(a);
3090 auto ps_advance = [&](std::size_t i) {
3091 if (has_gd) gd_advance_all();
else ps_advance_one(i);
3093 auto ps_reschedule = [&](std::size_t i) {
3094 if (has_gd) gd_reschedule_all();
else ps_reschedule_one(i);
3096 auto sd_advance = [&](std::size_t i) {
3097 if (has_gd) gd_advance_all();
else sd_advance_one(i);
3099 auto sd_reschedule = [&](std::size_t i) {
3100 if (has_gd) gd_reschedule_all();
else sd_reschedule_one(i);
3116 auto dest_has_room = [&](std::size_t j, std::size_t r) ->
bool {
3117 StationState& d = S[j];
3118 if (d.role != Role::Queue)
return true;
3120 for (std::size_t k = 0; k < K; ++k) total += acc.qlen[j][k] - d.blocked_at[k];
3121 return (total + 1.0 <= d.cap) &&
3122 (acc.qlen[j][r] - d.blocked_at[r] + 1.0 <= d.classcap[r]);
3127 return (S[j].role == Role::Queue && r < S[j].droprule.size()) ? S[j].droprule[r]
3135 std::vector<Region> regions;
3136 std::vector<int> region_of(M, -1);
3137 for (
const auto& rg : sn.regions) {
3140 R.members.assign(M,
false);
3141 R.class_cap.assign(K, -1.0);
3142 R.class_size.assign(K, 1.0);
3143 R.class_weight.assign(K, 1.0);
3145 R.jobs.assign(K, 0.0);
3146 R.blocked.assign(K, 0.0);
3147 R.dropped.assign(K, 0.0);
3148 R.tot_jobs.assign(K, 0.0);
3149 R.tot_weight.assign(K, 0.0);
3150 R.tot_mem.assign(K, 0.0);
3151 R.completed.assign(K, 0.0);
3152 R.resp_sum.assign(K, 0.0);
3153 R.resp_cnt.assign(K, 0.0);
3154 for (std::size_t i = 0; i < M && i < rg.members.size(); ++i)
3155 if (rg.members[i]) {
3156 R.members[i] =
true;
3157 if (region_of[i] >= 0)
3159 "SolverLDES (native engine): station '" + sn.stations[i].name +
3160 "' belongs to more than one finite capacity region");
3161 region_of[i] =
static_cast<int>(regions.size());
3166 for (std::size_t i = 0; i < rg.cap.size(); ++i) {
3167 if (i >= M || !R.members[i])
continue;
3168 for (std::size_t r = 0; r < K && r < rg.cap[i].size(); ++r)
3169 if (rg.cap[i][r] >= 0.0)
3170 R.class_cap[r] = (R.class_cap[r] < 0.0) ? rg.cap[i][r]
3171 : std::min(R.class_cap[r], rg.cap[i][r]);
3172 if (rg.cap[i].size() > K && rg.cap[i][K] >= 0.0)
3173 R.global_cap = (R.global_cap < 0.0) ? rg.cap[i][K]
3174 : std::min(R.global_cap, rg.cap[i][K]);
3176 for (
double mm : rg.maxmem)
3177 if (mm >= 0.0) R.max_mem = (R.max_mem < 0.0) ? mm : std::min(R.max_mem, mm);
3178 for (std::size_t r = 0; r < K && r < rg.rule.size(); ++r) R.rule[r] = rg.rule[r];
3179 for (std::size_t r = 0; r < K && r < rg.size.size(); ++r)
3180 R.class_size[r] = num_traits<T>::to_double(rg.size[r]);
3181 for (std::size_t r = 0; r < K && r < rg.weight.size(); ++r)
3182 R.class_weight[r] = num_traits<T>::to_double(rg.weight[r]);
3183 for (std::size_t c = 0; c < rg.lincon_A.rows(); ++c) {
3184 std::vector<double> row;
3185 for (std::size_t r = 0; r < rg.lincon_A.cols(); ++r)
3186 row.push_back(num_traits<T>::to_double(rg.lincon_A(c, r)));
3187 R.lincon_A.push_back(row);
3188 R.lincon_b.push_back(c < rg.lincon_b.size() ? num_traits<T>::to_double(rg.lincon_b[c])
3191 regions.push_back(R);
3199 struct RegionWaiter {
3200 std::size_t station = 0;
3203 std::vector<std::vector<RegionWaiter>> region_wait(regions.size());
3219 auto draw_at_arrival = [&](std::size_t i, std::size_t r) ->
bool {
3220 const StationState& s = S[i];
3225 if (s.bmsp)
return false;
3226 if (s.role != Role::Queue || s.ps || s.size_ordered)
return true;
3227 return !s.svc[r].time_varying();
3245 auto bmsp_arm = [&](std::size_t i) {
3246 StationState& s = S[i];
3247 const std::size_t r = s.bmsp_cls;
3248 const double gap = slot_snap(s.svc[r].next_at(g_svc[i][r], now),
"firing interval");
3249 const int b = s.svc[r].last_batch();
3250 s.bmsp_pending = (b > 0) ?
static_cast<std::size_t
>(b) : 1;
3251 if (!(gap > 0.0) && s.svc[r].time_varying())
return;
3263 auto start_service = [&](std::size_t i, std::size_t slot,
const Job& job_in) {
3264 StationState& s = S[i];
3274 if (!draw_at_arrival(i, job.cls)) {
3275 job.service = slot_snap(s.svc[job.cls].next_at(g_svc[i][job.cls], now),
"service time");
3276 job.remaining = job.service;
3282 if (!drew && s.preemptive && !s.resume && job.elapsed > 0.0) {
3283 job.service = slot_snap(s.svc[job.cls].next_at(g_svc[i][job.cls], now),
"service time");
3284 job.remaining = job.service;
3294 if (s.has_pools && !(s.preemptive && s.resume && job.elapsed > 0.0)) {
3295 const std::size_t ty = s.server_type[slot];
3296 if (s.type_has_svc[ty][job.cls]) {
3297 job.service = slot_snap(
3298 s.type_svc[ty][job.cls].next_at(g_hsvc[i][ty][job.cls], now),
3300 job.remaining = job.service;
3304 s.server[slot] = job;
3305 s.server_busy[slot] =
true;
3306 s.server_start[slot] = now;
3307 s.server_tag[slot] = ++ps_tag;
3313 double occupied = 1.0;
3314 if (s.has_parallelism && s.parallelism[job.cls] > 1) {
3315 const std::size_t need = s.parallelism[job.cls];
3316 std::vector<std::size_t> extras;
3317 for (std::size_t sl = 0; sl < s.nservers && extras.size() + 1 < need; ++sl) {
3318 if (sl == slot || s.server_busy[sl] || s.server_blocked[sl] || s.server_held[sl])
3320 s.server_busy[sl] =
true;
3321 extras.push_back(sl);
3323 if (extras.size() + 1 < need)
3324 throw InputError(
"SolverLDES (native engine): station '" + sn.stations[i].name +
3325 "' started a job holding " + std::to_string(need) +
3326 " servers with only " + std::to_string(extras.size() + 1) +
3328 s.held_extra[slot] = extras;
3329 occupied =
static_cast<double>(need);
3331 acc.update_busy(i, job.cls, now);
3334 acc.busy[i][job.cls] += occupied;
3341 const double scale0 = rate_scaling(i, job.cls);
3342 if (!(scale0 > 0.0))
return;
3344 e.t = now + job.remaining / scale0;
3349 e.tag = s.server_tag[slot];
3362 auto release_slots = [&](std::size_t i, std::size_t slot) ->
double {
3363 StationState& s = S[i];
3364 if (!s.has_parallelism || slot >= s.held_extra.size())
return 1.0;
3365 std::vector<std::size_t>& extras = s.held_extra[slot];
3366 if (extras.empty())
return 1.0;
3367 for (std::size_t j = 0; j < extras.size(); ++j) s.server_busy[extras[j]] =
false;
3368 const double n =
static_cast<double>(extras.size() + 1);
3381 auto serve_next_on_slot = [&](std::size_t i, std::size_t sl) {
3382 if (!buffer_has_for_slot(i, sl))
return;
3383 Job nextjob = buffer_pop_for_slot(i, sl);
3384 if (S[i].has_parallelism && free_slot_count(i) < S[i].parallelism[nextjob.cls]) {
3385 buffer_push(i, nextjob);
3388 start_service(i, sl, nextjob);
3398 auto preempt = [&](std::size_t i, std::size_t slot) {
3399 StationState& s = S[i];
3400 Job job = s.server[slot];
3401 const double served = now - s.server_start[slot];
3402 job.remaining = std::max(0.0, job.remaining - served);
3403 job.elapsed += served;
3404 s.server_busy[slot] =
false;
3405 s.server_tag[slot] = 0;
3406 const double freed = release_slots(i, slot);
3407 acc.update_busy(i, job.cls, now);
3408 acc.busy[i][job.cls] -= freed;
3411 buffer_push(i, job);
3422 std::function<void(std::size_t)> poll_serve;
3423 std::function<void(std::size_t)> poll_advance;
3435 auto poll_open_visit = [&](StationState& s) {
3436 std::size_t found = 0;
3437 for (
const Job& j : s.buffer)
3438 if (j.cls == s.poll_at) ++found;
3439 switch (s.poll_type) {
3443 s.poll_budget = std::numeric_limits<std::size_t>::max();
3444 s.poll_gate = found;
3447 default: s.poll_budget = std::numeric_limits<std::size_t>::max();
break;
3451 poll_serve = [&](std::size_t i) {
3452 StationState& s = S[i];
3453 if (s.server_busy[0] || s.server_blocked[0] || s.server_held[0])
return;
3455 std::size_t at = s.buffer.size();
3456 std::size_t waiting = 0;
3457 for (std::size_t j = 0; j < s.buffer.size(); ++j)
3458 if (s.buffer[j].cls == s.poll_at) {
3460 if (at == s.buffer.size() || s.buffer[j].t_arr < s.buffer[at].t_arr) at = j;
3464 const bool decremented =
3466 if (at == s.buffer.size() || s.poll_budget == 0 || decremented) {
3470 Job job = s.buffer[at];
3471 s.buffer.erase(s.buffer.begin() +
static_cast<std::ptrdiff_t
>(at));
3472 std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
3475 start_service(i, 0, job);
3478 poll_advance = [&](std::size_t i) {
3479 StationState& s = S[i];
3480 for (std::size_t step = 0; step < K; ++step) {
3481 const std::size_t from = s.poll_at;
3482 const std::size_t next = (from + 1) % K;
3485 if (s.has_switchover[from]) {
3486 s.poll_switching =
true;
3488 e.t = now + s.switchover[from].next(g_aux[i]);
3497 bool has_work =
false;
3498 for (
const Job& j : s.buffer)
3499 if (j.cls == next) {
3511 s.poll_parked =
true;
3526 std::function<void(std::size_t)> pas_reschedule = [&](std::size_t i) {
3527 StationState& s = S[i];
3528 if (s.pas_list.empty())
return;
3529 std::vector<std::size_t> seq;
3530 for (
const Job& j : s.pas_list) seq.push_back(j.cls);
3531 const double total = s.pas_rate(seq);
3532 if (!(total > 0.0))
return;
3533 s.pas_tag = ++ps_tag;
3535 e.t = now + (-std::log(uniform01(g_aux[i])) / total);
3540 e.slot = std::numeric_limits<std::size_t>::max();
3548 auto pas_complete = [&](std::size_t i) -> Job {
3549 StationState& s = S[i];
3550 std::vector<std::size_t> seq;
3551 for (
const Job& j : s.pas_list) seq.push_back(j.cls);
3552 std::vector<double> inc(seq.size(), 0.0);
3553 double prev = 0.0, total = 0.0;
3554 for (std::size_t p = 0; p < seq.size(); ++p) {
3555 std::vector<std::size_t> pre(seq.begin(),
3556 seq.begin() +
static_cast<std::ptrdiff_t
>(p + 1));
3557 const double cur = s.pas_rate(pre);
3558 inc[p] = std::max(0.0, cur - prev);
3562 std::size_t pos = 0;
3564 const double u = uniform01(g_routing) * total;
3566 for (std::size_t p = 0; p < inc.size(); ++p) {
3578 std::vector<std::size_t> chain(1, pos);
3579 std::size_t moving = seq[pos], cur = pos;
3581 std::size_t nxt = seq.size();
3582 for (std::size_t j = cur + 1; j < seq.size(); ++j) {
3583 const bool swappable =
3584 s.pas_swap.empty() ||
3585 (moving < s.pas_swap.size() && seq[j] < s.pas_swap[moving].size() &&
3586 s.pas_swap[moving][seq[j]]);
3592 if (nxt >= seq.size())
break;
3593 chain.push_back(nxt);
3597 Job departing = s.pas_list[chain.back()];
3598 std::vector<Job> arr = s.pas_list;
3599 for (std::size_t k = 0; k + 1 < chain.size(); ++k)
3600 arr[chain[k + 1]] = s.pas_list[chain[k]];
3601 std::vector<Job> survivors;
3602 for (std::size_t k = 0; k < arr.size(); ++k)
3603 if (k != chain.front()) survivors.push_back(arr[k]);
3604 s.pas_list = survivors;
3623 std::size_t hop_old_cls = K;
3624 std::function<void(std::size_t)> release_region;
3626 double pending_region_wait = 0.0;
3627 bool hop_internal =
false;
3630 bool pending_retry =
false;
3633 bool pending_orbit_exit =
false;
3636 bool pending_preset =
false;
3639 auto retrial_server_free = [&](std::size_t i, std::size_t cls) ->
bool {
3640 StationState& s = S[i];
3641 if (s.ps)
return s.ps_jobs.size() < s.nservers;
3642 return free_slot_for(i, cls) < s.nservers;
3644 auto orbit_tick = [&](std::size_t i) {
3645 StationState& s = S[i];
3646 const double dt = now - s.orbit_last;
3648 for (std::size_t r = 0; r < K; ++r) s.tot_orbit[r] += s.orbit_size[r] * dt;
3653 auto arm_retrial = [&](std::size_t i, std::size_t c, std::uint64_t id) {
3654 StationState& s = S[i];
3655 if (s.retrial_policy[c] == 2) {
3656 if (s.const_pending[c] || s.orbit[c].empty())
return;
3657 s.const_pending[c] =
true;
3661 e.t = now + s.retrial[c].next(g_svc[i][c]);
3668 auto orbit_join = [&](std::size_t i, Job job) {
3669 StationState& s = S[i];
3670 const std::size_t c = job.cls;
3672 s.orbit_size[c] += 1.0;
3676 job.deadline = now + classdeadline[c];
3678 acc.update_qlen(i, c, now);
3679 acc.qlen[i][c] += 1.0;
3680 acc.arrived[i][c] += 1.0;
3681 bp.track(i, c, +1, now);
3682 s.orbit[c].push_back(job);
3683 if (s.has_orbit_abandon[c]) {
3685 e.t = now + s.orbit_abandon[c].next(g_svc[i][c]);
3692 arm_retrial(i, c, job.id);
3696 auto orbit_leave = [&](std::size_t i, std::size_t c, std::size_t at,
bool leaves) {
3697 StationState& s = S[i];
3699 s.orbit_size[c] -= 1.0;
3700 s.orbit[c].erase(s.orbit[c].begin() +
static_cast<std::ptrdiff_t
>(at));
3703 acc.update_qlen(i, c, now);
3704 acc.qlen[i][c] -= 1.0;
3705 bp.track(i, c, -1, now);
3710 std::function<bool(std::size_t, Job, std::size_t)> admit =
3711 [&](std::size_t i, Job job, std::size_t from) ->
bool {
3712 StationState& s = S[i];
3713 const double rwait = pending_region_wait;
3714 pending_region_wait = 0.0;
3715 const bool is_retry = pending_retry;
3716 pending_retry =
false;
3717 const bool from_orbit = pending_orbit_exit;
3718 pending_orbit_exit =
false;
3719 const bool preset = pending_preset;
3720 pending_preset =
false;
3722 throw InputError(
"SolverLDES (native engine): a job of class '" +
3723 sn.classes[job.cls].name +
"' reached station '" +
3724 sn.stations[i].name +
"', which does not serve it");
3732 int rgi = region_of[i];
3733 if (rgi >= 0 && from < M && region_of[from] == rgi) {
3734 const std::size_t oc = hop_old_cls;
3739 hop_internal =
true;
3740 if (oc != job.cls) {
3741 const std::size_t rg =
static_cast<std::size_t
>(rgi);
3742 regions[rg].completed[oc] += 1.0;
3750 Region& R = regions[
static_cast<std::size_t
>(rgi)];
3751 if (R.would_exceed(job.cls)) {
3753 if (R.drops(job.cls)) {
3754 R.dropped[job.cls] += 1.0;
3755 dropped[i][job.cls] += 1.0;
3763 w.job.id = ++job_id;
3764 region_wait[
static_cast<std::size_t
>(rgi)].push_back(w);
3765 R.blocked[job.cls] += 1.0;
3774 if (s.role == Role::Queue && !from_orbit) {
3776 for (std::size_t r = 0; r < K; ++r) total += acc.qlen[i][r] - s.blocked_at[r];
3777 if (total + 1.0 > s.cap ||
3778 acc.qlen[i][job.cls] - s.blocked_at[job.cls] + 1.0 > s.classcap[job.cls]) {
3779 if (s.has_retrial[job.cls]) s.retrial_lost[job.cls] += 1.0;
3780 dropped[i][job.cls] += 1.0;
3784 acc.arrived[i][job.cls] += 1.0;
3789 if (s.has_retrial[job.cls] && !retrial_server_free(i, job.cls)) {
3802 if (!is_retry && !preset && !s.balk[job.cls].empty()) {
3804 for (std::size_t r = 0; r < K; ++r)
3805 total += acc.qlen[i][r] - (s.blocked_at.empty() ? 0.0 : s.blocked_at[r]);
3806 for (
const StationState::BalkRule& br : s.balk[job.cls])
3807 if (total >= br.min_jobs && (br.max_jobs < 0.0 || total <= br.max_jobs)) {
3808 if (uniform01(g_routing) < br.probability) {
3809 balked[i][job.cls] += 1.0;
3810 acc.arrived[i][job.cls] += 1.0;
3816 if (!from_orbit && !preset) job.t_arr = now;
3817 job.region_wait = rwait;
3818 job.priority = classprio[job.cls];
3824 if (draw_at_arrival(i, job.cls) && !job.preset) {
3825 job.service = slot_snap(S[i].svc[job.cls].next_at(g_svc[i][job.cls], now),
3827 job.remaining = job.service;
3830 job.rank = g_siro.aux.next_double();
3831 if (!from_orbit) job.deadline = (preset ? job.t_sys : now) + classdeadline[job.cls];
3834 Region& R = regions[
static_cast<std::size_t
>(rgi)];
3841 acc.update_qlen(i, job.cls, now);
3842 acc.qlen[i][job.cls] += 1.0;
3846 acc.arrived[i][job.cls] += 1.0;
3847 bp.track(i, job.cls, +1, now);
3850 if (s.role == Role::Delay) {
3851 job.due = now + job.service;
3852 if (track_delay_jobs) {
3854 delay_live[i].push_back(job);
3857 e.t = now + job.service;
3861 e.slot = std::numeric_limits<std::size_t>::max();
3869 if (has_gd) sd_reschedule(i);
3876 if (s.lps_limit > 0 && s.ps_jobs.size() >= s.lps_limit) {
3877 buffer_push(i, job);
3882 pj.priority = job.priority;
3883 pj.t_arr = job.t_arr;
3884 pj.region_wait = job.region_wait;
3885 pj.t_sys = job.t_sys;
3886 pj.total = job.preset ? job.elapsed + job.service : job.service;
3887 pj.remaining = job.service;
3888 pj.parent = job.parent;
3890 s.ps_jobs.push_back(pj);
3897 if (s.has_setup && !s.setup_on) {
3898 buffer_push(i, job);
3899 if (!s.setup_running) {
3900 s.setup_running =
true;
3902 e.t = now + s.setup_time.next(g_aux[i]);
3911 s.pas_list.push_back(job);
3912 acc.update_busy(i, job.cls, now);
3913 acc.busy[i][job.cls] += 1.0;
3918 buffer_push(i, job);
3919 if (s.poll_parked && !s.poll_switching) {
3920 s.poll_parked =
false;
3923 if (job.cls == s.poll_at) {
3929 }
else if (!s.server_busy[0] && !s.poll_switching) {
3938 buffer_push(i, job);
3939 if (!s.server_busy[0]) {
3940 s.server_busy[0] =
true;
3941 s.server_start[0] = now;
3942 acc.update_busy(i, job.cls, now);
3943 acc.busy[i][job.cls] += 1.0;
3949 const std::size_t sl = free_slot_for(i, job.cls);
3950 if (sl < s.nservers) {
3951 start_service(i, sl, job);
3966 (!s.has_parallelism || free_slot_count(i) + 1 >= s.parallelism[job.cls])) {
3967 std::vector<Job> held;
3968 std::vector<std::size_t> slot_of;
3969 std::vector<double> started;
3970 for (std::size_t sl = 0; sl < s.nservers; ++sl)
3971 if (s.server_busy[sl] && !s.server_blocked[sl]) {
3976 if (s.has_pools && !s.type_compat[s.server_type[sl]][job.cls])
continue;
3977 Job h = s.server[sl];
3981 h.remaining = std::max(0.0, h.remaining - (now - s.server_start[sl]));
3982 h.elapsed += now - s.server_start[sl];
3984 slot_of.push_back(sl);
3985 started.push_back(s.server_start[sl]);
3988 static_cast<double>(s.nservers), now);
3989 if (v < held.size()) {
3990 const std::size_t sl = slot_of[v];
3992 start_service(i, sl, job);
3997 buffer_push(i, job);
4002 if (s.has_patience[job.cls]) {
4004 e.t = now + s.patience[job.cls].next(g_svc[i][job.cls]);
4015 const std::map<std::size_t, double> kNoWeights;
4017 std::vector<std::size_t> tab_all;
4027 auto draw_prob = [&](
const std::vector<RouteDest>& tab) -> std::size_t {
4029 for (std::size_t k = 0; k < tab.size(); ++k) total += tab[k].mass;
4030 const double u = uniform01(g_routing) * total;
4032 std::size_t pick = 0;
4033 for (pick = 0; pick + 1 < tab.size(); ++pick) {
4034 cum += tab[pick].mass;
4035 if (u <= cum)
break;
4051 const bool has_sdr = !sn.sdr.empty() && !sn.sdr_nodes.empty();
4052 pfqn::SdrCoeff sdr_coeff;
4054 std::vector<double> sdr_pop, sdr_mass;
4066 auto draw_sdr = [&](
const std::vector<RouteDest>& tab) -> std::size_t {
4067 sdr_pop.assign(M, 0.0);
4068 for (std::size_t i = 0; i < M; ++i) {
4069 if (S[i].role == Role::Source)
continue;
4071 for (std::size_t q = 0; q < K; ++q) held += acc.qlen[i][q];
4077 sdr_mass.assign(tab.size(), 0.0);
4079 for (std::size_t i = 0; i < tab.size(); ++i) {
4080 const std::size_t dnode = tab[i].node - 1;
4082 for (std::size_t b = 1; b < sn.sdr_nodes.branch.size(); ++b)
4083 if (sn.sdr_nodes.entryOf[b] == dnode) p += Pb[b];
4084 if (sn.sdr_nodes.departure == dnode) p += Ped;
4087 if (p < 0.0) p = 0.0;
4096 if (!(total > 0.0))
return draw_prob(tab);
4097 const double u = uniform01(g_routing) * total;
4099 std::size_t pick = 0;
4100 for (pick = 0; pick + 1 < tab.size(); ++pick) {
4101 cum += sdr_mass[pick];
4102 if (u <= cum)
break;
4115 auto shortest_queue = [&](
const std::vector<RouteDest>& tab,
4116 const std::vector<std::size_t>& cand) -> std::size_t {
4117 double best = std::numeric_limits<double>::infinity();
4118 std::vector<std::size_t> tied;
4119 for (std::size_t c = 0; c < cand.size(); ++c) {
4120 const std::size_t k = cand[c];
4121 if (tab[k].station >= M || S[tab[k].station].role == Role::Source)
continue;
4123 for (std::size_t q = 0; q < K; ++q) held += acc.qlen[tab[k].station][q];
4128 }
else if (held == best) {
4132 if (tied.empty())
return cand.empty() ? 0 : cand[0];
4133 if (tied.size() == 1)
return tied[0];
4134 std::size_t at =
static_cast<std::size_t
>(uniform01(g_routing) *
4135 static_cast<double>(tied.size()));
4136 if (at >= tied.size()) at = tied.size() - 1;
4150 auto draw_node_route = [&](std::size_t inode, std::size_t r) -> RouteEntry {
4151 const std::vector<RouteDest>& tab = nroute[inode][r];
4153 throw InputError(
"SolverLDES (native engine): node '" + sn.nodes[inode - 1].name +
4154 "' has no routing for class '" + sn.classes[r].name +
"'");
4155 std::size_t pick = 0;
4156 if (tab.size() > 1) {
4157 if (tab_all.size() != tab.size()) {
4158 tab_all.resize(tab.size());
4159 for (std::size_t k = 0; k < tab.size(); ++k) tab_all[k] = k;
4161 switch (node_routing[inode][r]) {
4163 pick =
static_cast<std::size_t
>(uniform01(g_routing) *
4164 static_cast<double>(tab.size()));
4165 if (pick >= tab.size()) pick = tab.size() - 1;
4169 pick = rr_counter[inode][r] % tab.size();
4170 ++rr_counter[inode][r];
4179 const std::map<std::size_t, double>& w =
4180 (sn.nodes[inode - 1].routing_weights.size() > r)
4181 ? sn.nodes[inode - 1].routing_weights[r]
4183 std::vector<double> raw(tab.size(), 0.0);
4184 bool integral =
true;
4185 for (std::size_t k = 0; k < tab.size(); ++k) {
4186 const std::map<std::size_t, double>::const_iterator it = w.find(tab[k].node);
4187 raw[k] = (it == w.end()) ? 0.0 : it->second;
4188 if (raw[k] < 0.0 || raw[k] != std::floor(raw[k])) integral =
false;
4190 std::vector<std::size_t> quota(tab.size(), 0);
4191 std::size_t total = 0;
4192 for (std::size_t k = 0; k < tab.size(); ++k) {
4193 double v = integral ? raw[k] : raw[k] * 1000.0;
4194 std::size_t q = (v > 0.0) ?
static_cast<std::size_t
>(v) : 0;
4195 if (raw[k] > 0.0 && q < 1) q = 1;
4202 pick = rr_counter[inode][r] % tab.size();
4203 ++rr_counter[inode][r];
4206 std::size_t counter = rr_counter[inode][r] % total;
4207 ++rr_counter[inode][r];
4208 std::size_t cum = 0;
4209 for (pick = 0; pick + 1 < tab.size(); ++pick) {
4211 if (counter < cum)
break;
4223 pick = has_sdr ? draw_sdr(tab) : draw_prob(tab);
4227 pick = shortest_queue(tab, tab_all);
4236 std::size_t d = (sn.nodes[inode - 1].routing_param.size() > r &&
4237 sn.nodes[inode - 1].routing_param[r] > 0)
4238 ?
static_cast<std::size_t
>(
4239 sn.nodes[inode - 1].routing_param[r])
4241 std::vector<std::size_t> cand;
4242 for (std::size_t k = 0; k < d; ++k) {
4243 std::size_t at =
static_cast<std::size_t
>(
4244 uniform01(g_routing) *
static_cast<double>(tab.size()));
4245 if (at >= tab.size()) at = tab.size() - 1;
4248 pick = shortest_queue(tab, cand);
4252 pick = draw_prob(tab);
4257 const RouteDest& d = tab[pick];
4260 e.station = d.station;
4262 e.cls = d.cls[0].first;
4263 if (d.cls.size() > 1) {
4265 for (std::size_t k = 0; k < d.cls.size(); ++k) total += d.cls[k].second;
4266 const double u = uniform01(g_routing) * total;
4268 for (std::size_t k = 0; k < d.cls.size(); ++k) {
4269 cum += d.cls[k].second;
4270 if (u <= cum || k + 1 == d.cls.size()) {
4271 e.cls = d.cls[k].first;
4300 auto resolve_final_cls = [&](std::size_t node, std::size_t cls) -> std::size_t {
4301 std::size_t at_node = node, at_cls = cls;
4302 for (std::size_t hop = 0; hop <= nnodes && at_node != 0; ++hop) {
4303 const NodeType dt = sn.nodes[at_node - 1].nodetype;
4304 if (dt != NodeType::ClassSwitch && dt != NodeType::Router && dt != NodeType::Logger)
4306 const RouteEntry e = draw_node_route(at_node, at_cls);
4339 std::vector<std::pair<std::size_t, std::size_t>> fb_loggers;
4340 auto resolve_forced_dest = [&](
const RouteEntry&
re) -> RouteEntry {
4343 for (std::size_t hop = 0; hop <= nnodes; ++hop) {
4344 if (at.node == 0 || at.sink)
break;
4345 const NodeType dt = sn.nodes[at.node - 1].nodetype;
4346 if (dt != NodeType::ClassSwitch && dt != NodeType::Router && dt != NodeType::Logger)
4348 if (dt == NodeType::Logger) fb_loggers.push_back(std::make_pair(at.node, at.cls));
4349 const std::vector<RouteDest>& tab = nroute[at.node][at.cls];
4350 if (tab.size() != 1 || tab[0].cls.size() != 1)
break;
4352 nxt.node = tab[0].node;
4353 nxt.station = tab[0].station;
4354 nxt.sink = tab[0].sink;
4355 nxt.cls = tab[0].cls[0].first;
4366 if (o.busy_period_orders > 0) {
4367 std::vector<std::vector<std::size_t>> sets;
4368 std::vector<std::string> set_names;
4369 for (std::size_t i = 0; i < M; ++i) {
4370 if (S[i].role == Role::Source || S[i].role == Role::Synchronization)
continue;
4371 sets.push_back(std::vector<std::size_t>(1, i));
4372 set_names.push_back(sn.stations[i].name);
4374 for (
const std::vector<std::size_t>& sub : o.busy_period_subnets) {
4376 for (std::size_t s2 : sub) {
4377 if (s2 >= M || S[s2].role == Role::Source)
4378 throw InputError(
"SolverLDES (native engine): station " +
4379 std::to_string(s2) +
4380 " of a busy period subnetwork is not a service station");
4381 nm += (nm.empty() ?
"" :
"+") + sn.stations[s2].name;
4383 sets.push_back(sub);
4384 set_names.push_back(nm);
4386 std::vector<BpTarget> targets;
4387 for (std::size_t si = 0; si < sets.size(); ++si) {
4389 t.stations = sets[si];
4391 t.name = set_names[si];
4392 targets.push_back(t);
4393 for (std::size_t r = 0; r < K; ++r) {
4395 tc.stations = sets[si];
4396 tc.job_class =
static_cast<int>(r);
4397 tc.name = set_names[si] +
":" + sn.classes[r].name;
4398 targets.push_back(tc);
4401 bp.init(targets, o.busy_period_orders, M);
4414 std::function<void(std::size_t)> release_blocked = [&](std::size_t j) {
4415 bool progress =
true;
4418 std::size_t best_i = M, best_sl = 0;
4419 double oldest = std::numeric_limits<double>::infinity();
4420 for (std::size_t a = 0; a < M; ++a) {
4424 if (S[a].role != Role::Queue || S[a].ps || S[a].server_blocked.empty())
continue;
4425 for (std::size_t sl = 0; sl < S[a].nservers; ++sl)
4426 if (S[a].server_blocked[sl] && S[a].blocked_dest[sl] == j &&
4427 S[a].blocked_job[sl].t_arr < oldest) {
4428 oldest = S[a].blocked_job[sl].t_arr;
4433 if (best_i >= M)
break;
4434 const std::size_t dcls = S[best_i].blocked_dest_cls[best_sl];
4435 if (!dest_has_room(j, dcls))
break;
4437 StationState& src = S[best_i];
4438 Job moved = src.blocked_job[best_sl];
4439 src.server_blocked[best_sl] =
false;
4440 src.server_busy[best_sl] =
false;
4441 src.server_tag[best_sl] = 0;
4442 acc.update_qlen(j, dcls, now);
4443 acc.qlen[j][dcls] -= 1.0;
4444 S[j].blocked_at[dcls] -= 1.0;
4447 acc.update_busy(best_i, moved.cls, now);
4448 if (acc.busy[best_i][moved.cls] > 0.0) acc.busy[best_i][moved.cls] -= 1.0;
4452 fresh.t_sys = moved.t_sys;
4453 fresh.parent = moved.parent;
4454 admit(j, fresh, best_i);
4458 serve_next_on_slot(best_i, best_sl);
4477 release_region = [&](std::size_t rg) {
4478 std::vector<RegionWaiter>& q = region_wait[rg];
4479 while (!q.empty()) {
4481 const std::size_t w = 0;
4483 if (R.would_exceed(q[w].job.cls))
break;
4484 RegionWaiter taken = q[w];
4486 R.blocked[taken.job.cls] -= 1.0;
4487 q.erase(q.begin() +
static_cast<std::ptrdiff_t
>(w));
4489 fresh.cls = taken.job.cls;
4490 fresh.t_sys = taken.job.t_sys;
4493 fresh.parent = taken.job.parent;
4494 fresh.call = taken.job.call;
4498 pending_region_wait = now - taken.job.t_arr;
4499 admit(taken.station, fresh, M);
4507 auto signal_held = [&](std::size_t i) -> std::size_t {
4508 const StationState& s = S[i];
4509 if (s.role == Role::Delay)
return delay_live[i].size();
4510 if (s.pas)
return s.pas_list.size();
4511 std::size_t n = s.buffer.size();
4512 if (s.ps)
return n + s.ps_jobs.size();
4513 for (std::size_t sl = 0; sl < s.nservers; ++sl)
4514 if (s.server_busy[sl]) ++n;
4535 StationState& s = S[i];
4538 if (s.role == Role::Delay) {
4539 std::vector<Job>& live = delay_live[i];
4540 if (live.empty())
return K;
4543 for (std::size_t j = 1; j < live.size(); ++j)
4544 if (live[j].t_arr < live[at].t_arr) at = j;
4546 for (std::size_t j = 1; j < live.size(); ++j)
4547 if (live[j].t_arr > live[at].t_arr) at = j;
4549 at =
static_cast<std::size_t
>(uniform01(g_routing) *
4550 static_cast<double>(live.size()));
4551 if (at >= live.size()) at = live.size() - 1;
4553 const std::size_t rc = live[at].cls;
4554 live.erase(live.begin() +
static_cast<std::ptrdiff_t
>(at));
4562 if (s.pas_list.empty())
return K;
4565 for (std::size_t j = 1; j < s.pas_list.size(); ++j)
4566 if (s.pas_list[j].t_arr < s.pas_list[at].t_arr) at = j;
4568 for (std::size_t j = 1; j < s.pas_list.size(); ++j)
4569 if (s.pas_list[j].t_arr > s.pas_list[at].t_arr) at = j;
4571 at =
static_cast<std::size_t
>(uniform01(g_routing) *
4572 static_cast<double>(s.pas_list.size()));
4573 if (at >= s.pas_list.size()) at = s.pas_list.size() - 1;
4575 const std::size_t rc = s.pas_list[at].cls;
4576 s.pas_list.erase(s.pas_list.begin() +
static_cast<std::ptrdiff_t
>(at));
4577 acc.update_busy(i, rc, now);
4578 if (acc.busy[i][rc] > 0.0) acc.busy[i][rc] -= 1.0;
4583 const std::size_t waiting = s.buffer.size();
4588 const std::size_t inservice = s.ps_jobs.size();
4589 if (waiting + inservice == 0)
return K;
4592 const std::size_t pick =
static_cast<std::size_t
>(
4593 uniform01(g_routing) *
static_cast<double>(waiting + inservice));
4594 from_wait = pick < waiting;
4596 from_wait = waiting > 0;
4601 for (std::size_t j = 1; j < s.buffer.size(); ++j)
4602 if (s.buffer[j].t_arr < s.buffer[at].t_arr) at = j;
4604 for (std::size_t j = 1; j < s.buffer.size(); ++j)
4605 if (s.buffer[j].t_arr > s.buffer[at].t_arr) at = j;
4607 at = std::min(waiting - 1,
4608 static_cast<std::size_t
>(uniform01(g_routing) *
4609 static_cast<double>(waiting)));
4611 const std::size_t rc = s.buffer[at].cls;
4612 s.buffer.erase(s.buffer.begin() +
static_cast<std::ptrdiff_t
>(at));
4613 if (s.sched != SchedStrategy::FSP)
4614 std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
4620 for (std::size_t j = 1; j < s.ps_jobs.size(); ++j)
4621 if (s.ps_jobs[j].t_arr < s.ps_jobs[at].t_arr) at = j;
4623 for (std::size_t j = 1; j < s.ps_jobs.size(); ++j)
4624 if (s.ps_jobs[j].t_arr > s.ps_jobs[at].t_arr) at = j;
4626 at = std::min(inservice - 1,
4627 static_cast<std::size_t
>(uniform01(g_routing) *
4628 static_cast<double>(inservice)));
4630 const std::size_t rc = s.ps_jobs[at].cls;
4631 s.ps_jobs.erase(s.ps_jobs.begin() +
static_cast<std::ptrdiff_t
>(at));
4633 if (s.lps_limit > 0 && !s.buffer.empty() && s.ps_jobs.size() < s.lps_limit) {
4634 Job nextjob = buffer_pop(i);
4636 pj.cls = nextjob.cls;
4637 pj.priority = nextjob.priority;
4638 pj.t_arr = nextjob.t_arr;
4639 pj.region_wait = nextjob.region_wait;
4640 pj.t_sys = nextjob.t_sys;
4641 pj.total = nextjob.preset ? nextjob.elapsed + nextjob.service : nextjob.service;
4642 pj.remaining = nextjob.service;
4643 pj.parent = nextjob.parent;
4644 pj.ext = nextjob.ext;
4645 s.ps_jobs.push_back(pj);
4654 std::size_t inservice = 0;
4655 for (std::size_t sl = 0; sl < s.nservers; ++sl)
4656 if (s.server_busy[sl]) ++inservice;
4657 if (waiting + inservice == 0)
return K;
4660 const std::size_t pick =
static_cast<std::size_t
>(
4661 uniform01(g_routing) *
static_cast<double>(waiting + inservice));
4662 from_wait = pick < waiting;
4664 from_wait = waiting > 0;
4670 for (std::size_t j = 1; j < s.buffer.size(); ++j)
4671 if (s.buffer[j].t_arr < s.buffer[at].t_arr) at = j;
4673 for (std::size_t j = 1; j < s.buffer.size(); ++j)
4674 if (s.buffer[j].t_arr > s.buffer[at].t_arr) at = j;
4676 at = std::min(waiting - 1,
4677 static_cast<std::size_t
>(uniform01(g_routing) *
4678 static_cast<double>(waiting)));
4680 const std::size_t rc = s.buffer[at].cls;
4681 s.buffer.erase(s.buffer.begin() +
static_cast<std::ptrdiff_t
>(at));
4682 if (s.sched != SchedStrategy::FSP)
4683 std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
4687 std::size_t slot = s.nservers;
4688 for (std::size_t sl = 0; sl < s.nservers; ++sl) {
4689 if (!s.server_busy[sl])
continue;
4690 if (slot == s.nservers) {
4695 if (s.server[sl].t_arr < s.server[slot].t_arr) slot = sl;
4697 if (s.server[sl].t_arr > s.server[slot].t_arr) slot = sl;
4701 std::size_t pick = std::min(inservice - 1,
4702 static_cast<std::size_t
>(uniform01(g_routing) *
4703 static_cast<double>(inservice)));
4704 for (std::size_t sl = 0; sl < s.nservers; ++sl)
4705 if (s.server_busy[sl]) {
4713 if (slot == s.nservers)
return K;
4715 const std::size_t rc = s.server[slot].cls;
4718 s.server_busy[slot] =
false;
4719 s.server_tag[slot] = 0;
4720 const double freed = release_slots(i, slot);
4721 acc.update_busy(i, rc, now);
4722 acc.busy[i][rc] -= freed;
4727 }
else if (!s.server_blocked[slot] && !s.server_held[slot]) {
4728 serve_next_on_slot(i, slot);
4750 auto apply_signal = [&](std::size_t i,
const Job& sig) {
4751 StationState& s = S[i];
4752 const std::size_t held = signal_held(i);
4754 std::size_t to_remove = 1;
4755 if (sig.cls < sn.signaltype.size() &&
4758 }
else if (sig.cls < sn.signalremdist.size() && !sn.signalremdist[sig.cls].empty()) {
4759 const std::vector<T>& pmf = sn.signalremdist[sig.cls];
4760 const double u = uniform01(g_routing);
4763 for (; b + 1 < pmf.size(); ++b) {
4764 cum +=
static_cast<double>(pmf[b]);
4767 to_remove = std::min(b, held);
4769 to_remove = std::min<std::size_t>(1, held);
4772 if (sig.cls < sn.signalrempolicy.size()) policy = sn.signalrempolicy[sig.cls];
4773 for (std::size_t rep = 0; rep < to_remove; ++rep) {
4774 if (signal_held(i) == 0)
break;
4775 const std::size_t rc = signal_remove_one(i, policy);
4777 last_lost[i][rc] = now;
4779 acc.update_qlen(i, rc, now);
4780 acc.qlen[i][rc] -= 1.0;
4781 bp.track(i, rc, -1, now);
4785 if (region_of[i] >= 0) {
4786 const std::size_t rg =
static_cast<std::size_t
>(region_of[i]);
4787 regions[rg].update(now);
4788 regions[rg].leave(rc);
4797 sys_resp_sum[sig.cls] += now - sig.t_sys;
4798 sys_resp_cnt[sig.cls] += 1.0;
4799 sys_completed[sig.cls] += 1.0;
4816 std::uint64_t total_completions = 0;
4828 std::function<void(CacheState&, std::size_t, std::size_t)> release_delayed_hits;
4844 std::function<void()> spn_settle;
4845 std::function<void(std::size_t, Job, std::size_t)> deliver =
4846 [&](std::size_t node, Job job, std::size_t from) {
4847 const NodeType nt = sn.nodes[node - 1].nodetype;
4856 if (nt == NodeType::Place) {
4857 place_deposit(place_of[node], job.cls, 1.0);
4861 if (nt == NodeType::Transition) {
4866 if (nt == NodeType::Sink) {
4867 sys_resp_sum[job.cls] += now - job.t_sys;
4868 sys_resp_cnt[job.cls] += 1.0;
4869 sys_completed[job.cls] += 1.0;
4873 if (nt == NodeType::Cache) {
4874 auto ci = caches.find(node);
4875 if (ci == caches.end())
4876 throw InputError(
"SolverLDES (native engine): Cache node '" +
4877 sn.nodes[node - 1].name +
"' carries no parameters");
4878 CacheState& cs = ci->second;
4886 if (cs.has_retrieval) {
4887 const int fetched = cs.fetches_item(job.cls);
4889 const std::size_t fi =
static_cast<std::size_t
>(fetched);
4891 cache_miss(cs, fi, uniform01(g_routing), uniform01(g_routing));
4892 cs.in_flight[fi] = 0;
4893 cs.total_fetch_time += now - cs.fetch_start[fi];
4894 cs.completed_fetches += 1.0;
4895 release_delayed_hits(cs, fi, node);
4896 if (cs.miss_class[job.cls] >= 0)
4897 out.cls =
static_cast<std::size_t
>(cs.miss_class[job.cls]);
4898 const RouteEntry
re = draw_node_route(node, out.cls);
4900 if (
re.node != 0 && sn.nodes[
re.node - 1].nodetype == NodeType::Sink)
4901 ++total_completions;
4902 deliver(
re.node, out, M);
4907 const std::vector<double>& pop = cs.popularity[job.cls];
4909 throw InputError(
"SolverLDES (native engine): class '" +
4910 sn.classes[job.cls].name +
"' reads cache '" +
4911 sn.nodes[node - 1].name +
"' with no item popularity");
4912 const double u = uniform01(g_routing);
4913 std::size_t item = 0;
4914 while (item + 1 < pop.size() && u >= pop[item]) ++item;
4916 const int at = cs.find(item);
4918 cs.hits[job.cls] += 1.0;
4919 cs.hits_per_list[job.cls][
static_cast<std::size_t
>(at)] += 1.0;
4921 cache_hit(cs, item,
static_cast<std::size_t
>(at), uniform01(g_routing));
4922 if (cs.hit_class[job.cls] >= 0)
4923 out.cls =
static_cast<std::size_t
>(cs.hit_class[job.cls]);
4924 }
else if (cs.has_retrieval && cs.retrieval_of(item, job.cls) > 0) {
4925 if (cs.in_flight[item]) {
4929 CacheState::HeldRequest hr;
4932 hr.t_sys = job.t_sys;
4933 cs.held[item].push_back(hr);
4939 cs.misses[job.cls] += 1.0;
4940 cs.in_flight[item] = 1;
4941 cs.fetch_start[item] = now;
4942 out.cls = cs.retrieval_of(item, job.cls) - 1;
4944 cs.misses[job.cls] += 1.0;
4946 cache_miss(cs, item, uniform01(g_routing), uniform01(g_routing));
4947 if (cs.miss_class[job.cls] >= 0)
4948 out.cls =
static_cast<std::size_t
>(cs.miss_class[job.cls]);
4952 const RouteEntry e = draw_node_route(node, out.cls);
4956 if (e.node != 0 && sn.nodes[e.node - 1].nodetype == NodeType::Sink)
4957 ++total_completions;
4958 deliver(e.node, out, M);
4962 if (nt == NodeType::ClassSwitch || nt == NodeType::Router || nt == NodeType::Logger) {
4978 std::size_t at = node;
4979 const std::uint64_t walk = logger_traces.empty() ? 0 : log_walk_id++;
4980 if (nt == NodeType::Logger) log_passage(node, out.cls, walk);
4981 for (std::size_t hop = 0; hop <= nnodes; ++hop) {
4982 const RouteEntry e = draw_node_route(at, out.cls);
4984 if (e.node == 0)
return;
4985 const NodeType dt = sn.nodes[e.node - 1].nodetype;
4986 if (dt == NodeType::ClassSwitch || dt == NodeType::Router ||
4987 dt == NodeType::Logger) {
4988 if (dt == NodeType::Logger) log_passage(e.node, out.cls, walk);
4992 deliver(e.node, out, from);
4995 throw InputError(
"SolverLDES (native engine): the routing out of node '" +
4996 sn.nodes[node - 1].name +
4997 "' cycles through pass-through nodes without reaching a station");
5000 if (nt == NodeType::Fork) {
5004 const std::vector<RouteDest>& tab = nroute[node][job.cls];
5005 if (tab.empty())
return;
5006 std::vector<RouteEntry> links;
5007 for (std::size_t di = 0; di < tab.size(); ++di)
5008 for (std::size_t ci = 0; ci < tab[di].cls.size(); ++ci) {
5010 e.node = tab[di].node;
5011 e.station = tab[di].station;
5012 e.sink = tab[di].sink;
5013 e.cls = tab[di].cls[ci].first;
5016 const qn::NodeDef& fk = sn.nodes[node - 1];
5017 const qn::ForkParam<double>* fp = sn.fork_param_of(node);
5018 const int per_link = std::max(1,
static_cast<int>(fk.tasks_per_link + 0.5));
5026 std::vector<int> per_branch(links.size(), per_link);
5028 const bool variable = fork_is_variable(node);
5029 for (std::size_t li = 0; li < links.size(); ++li) {
5031 const std::size_t k0 = links[li].node - 1, r0 = job.cls;
5032 const double p = fp->fan_out_prob(k0, r0);
5033 if (p < 1.0 && g_fork.aux.next_double() >= p) {
5035 }
else if (!fp->fan_out_dist[k0][r0].disabled) {
5036 per_branch[li] = sample_fork_degree(fp->fan_out_dist[k0][r0]);
5038 per_branch[li] =
static_cast<int>(fp->fan_out_link(k0, r0) + 0.5);
5041 total += per_branch[li];
5044 throw InputError(
"SolverLDES (native engine): Fork '" + fk.name +
5045 "' emitted no sibling; at least one branch must be certain "
5046 "to emit at least one task");
5050 fs.t_sys = job.t_sys;
5051 fs.outer = job.parent;
5053 fs.pair = fj_pair_of_fork[node];
5055 fs.required = fs.total;
5056 const int jn = join_of_fork[node];
5058 auto jd = sn.joindecl.find(
static_cast<std::size_t
>(jn));
5059 if (jd != sn.joindecl.end() &&
5065 fs.required = std::min(fs.total,
5066 std::max(1,
static_cast<int>(jd->second.quorum + 0.5)));
5069 const std::uint64_t pid = ++next_parent;
5070 fork_sync[pid] = fs;
5071 for (std::size_t li = 0; li < links.size(); ++li)
5072 for (
int t = 0; t < per_branch[li]; ++t) {
5074 sib.cls = links[li].cls;
5075 sib.t_sys = job.t_sys;
5077 deliver(links[li].node, sib, M);
5082 if (nt == NodeType::Join) {
5083 auto it = fork_sync.find(job.parent);
5084 if (it == fork_sync.end()) {
5096 const std::size_t jd = node_to_station[node];
5098 acc.arrived[jd][job.cls] += 1.0;
5099 acc.join_dropped[jd][job.cls] += 1.0;
5103 ForkSync& fs = it->second;
5112 const std::size_t js = node_to_station[node];
5114 acc.update_qlen(js, job.cls, now);
5115 acc.qlen[js][job.cls] += 1.0;
5116 acc.arrived[js][job.cls] += 1.0;
5118 fs.siblings.push_back(std::make_pair(job.cls, now));
5119 if (
static_cast<int>(fs.siblings.size()) < fs.required)
return;
5123 for (std::size_t si = 0; si < fs.siblings.size(); ++si) {
5124 const std::size_t sc = fs.siblings[si].first;
5125 acc.update_qlen(js, sc, now);
5126 acc.qlen[js][sc] -= 1.0;
5129 acc.tally_resp(js, sc, now - fs.siblings[si].second);
5131 acc.completed[js][fs.cls] += 1.0;
5134 merged.cls = fs.cls;
5135 merged.t_sys = fs.t_sys;
5136 merged.parent = fs.outer;
5137 if (fs.pair >= 0 && fs.cls < K) {
5138 const std::size_t p =
static_cast<std::size_t
>(fs.pair);
5139 fj_resp_sum[p][fs.cls] += now - fs.t_fork;
5140 fj_resp_cnt[p][fs.cls] += 1.0;
5141 fj_release[p][fs.cls].push_back(now);
5143 fork_sync.erase(it);
5144 const RouteEntry e = draw_node_route(node, merged.cls);
5146 deliver(e.node, merged, M);
5151 const std::size_t j = node_to_station[node];
5155 if (has_removal_signal && j < M && is_removal_signal[job.cls]) {
5156 apply_signal(j, job);
5169 if (has_sync_call && j < M && is_reply_signal[job.cls]) {
5170 acc.completed[j][job.cls] += 1.0;
5171 typename std::map<std::uint64_t, PendingCall>::iterator pc =
5172 pending_reply.find(job.call);
5173 if (pc != pending_reply.end()) {
5174 const std::size_t bi = pc->second.station, bs = pc->second.slot;
5175 const std::size_t bc = pc->second.cls;
5176 pending_reply.erase(pc);
5177 StationState& bst = S[bi];
5178 bst.server_held[bs] =
false;
5179 acc.update_busy(bi, bc, now);
5180 if (acc.busy[bi][bc] > 0.0) acc.busy[bi][bc] -= 1.0;
5181 acc.update_qlen(bi, bc, now);
5182 acc.held[bi][bc] -= 1.0;
5183 if (!bst.server_blocked[bs]) {
5184 serve_next_on_slot(bi, bs);
5187 release_blocked(bi);
5189 const RouteEntry re2 = draw_node_route(node, job.cls);
5191 onward.cls = re2.cls;
5192 onward.t_sys = job.t_sys;
5193 onward.parent = job.parent;
5194 deliver(re2.node, onward, M);
5202 j + 1 == sn.classes[job.cls].refstat) {
5203 sys_resp_sum[job.cls] += now - job.t_sys;
5204 sys_resp_cnt[job.cls] += 1.0;
5205 sys_completed[job.cls] += 1.0;
5209 throw InputError(
"SolverLDES (native engine): a closed job was dropped at "
5210 "station '" + sn.stations[j].name +
5211 "'; a closed class cannot lose population");
5214 release_delayed_hits = [&](CacheState& cs, std::size_t item, std::size_t node) {
5215 if (item >= cs.held.size() || cs.held[item].empty())
return;
5216 std::vector<CacheState::HeldRequest> freed;
5217 freed.swap(cs.held[item]);
5218 for (std::size_t i = 0; i < freed.size(); ++i) {
5219 const CacheState::HeldRequest& hr = freed[i];
5220 cs.delayed[hr.cls] += 1.0;
5221 cs.delayed_wait += now - hr.hold_time;
5223 out.cls = (cs.hit_class[hr.cls] >= 0) ?
static_cast<std::size_t
>(cs.hit_class[hr.cls])
5225 out.t_sys = hr.t_sys;
5226 const RouteEntry e = draw_node_route(node, out.cls);
5228 if (e.node != 0 && sn.nodes[e.node - 1].nodetype == NodeType::Sink)
5229 ++total_completions;
5230 deliver(e.node, out, M);
5242 auto flat_marking = [&]() {
5243 std::vector<double> tok(place_nodes.size() * K, 0.0);
5244 for (std::size_t p = 0; p < place_nodes.size(); ++p)
5245 for (std::size_t r = 0; r < K; ++r) tok[p * K + r] = marking[p][r];
5259 auto flat_avail = [&]() {
5260 std::vector<double> tok(place_nodes.size() * K, 0.0);
5261 for (std::size_t p = 0; p < place_nodes.size(); ++p)
5262 for (std::size_t r = 0; r < K; ++r) tok[p * K + r] = avail[p][r];
5275 auto spn_apply = [&](
const SpnMode& m) {
5276 for (std::size_t p = 0; p < place_nodes.size(); ++p)
5277 for (std::size_t r = 0; r < K; ++r) {
5278 const double e = m.enabling[p * K + r];
5279 if (e > 0.0) place_take(p, r, e);
5281 for (std::size_t p = 0; p < place_nodes.size(); ++p)
5282 for (std::size_t r = 0; r < K; ++r) {
5283 const double f = m.firing[p * K + r];
5284 if (f > 0.0) place_deposit(p, r, f);
5288 for (
const auto& ex : m.exits) {
5289 const std::size_t nd = std::get<0>(ex), r = std::get<1>(ex);
5290 const double n = std::get<2>(ex);
5291 if (sn.nodes[nd - 1].nodetype == NodeType::Sink) {
5292 sys_completed[r] += n;
5295 for (
double t = 0.0; t + 0.5 < n; t += 1.0) {
5299 deliver(nd, job, M);
5313 std::uint64_t fire_tag = 0;
5315 std::vector<std::vector<std::vector<std::uint64_t>>> spn_clocks(transitions.size());
5316 for (std::size_t t = 0; t < transitions.size(); ++t)
5317 spn_clocks[t].resize(transitions[t].modes.size());
5318 spn_settle = [&]() {
5319 if (transitions.empty())
return;
5320 for (
int guard = 0; guard < 1000000; ++guard) {
5321 bool fired_any =
false;
5322 std::vector<double> tok = flat_avail();
5323 for (std::size_t t = 0; t < transitions.size(); ++t) {
5325 if (k < 0)
continue;
5326 const SpnMode& m = transitions[t].modes[
static_cast<std::size_t
>(k)];
5330 transitions[t].fired[
static_cast<std::size_t
>(k)] += 1.0;
5331 ++total_completions;
5335 if (!fired_any)
break;
5342 const std::vector<double> tok = flat_avail();
5343 const std::vector<double> mk = flat_marking();
5344 for (std::size_t t = 0; t < transitions.size(); ++t)
5345 for (std::size_t k = 0; k < transitions[t].modes.size(); ++k) {
5346 const SpnMode& m = transitions[t].modes[k];
5347 if (m.immediate)
continue;
5348 std::vector<std::uint64_t>& live = spn_clocks[t][k];
5349 const bool on =
spn_enabled(m, tok) && !spn_law[t][k].disabled();
5355 const std::size_t target =
5357 while (live.size() > target) {
5358 const std::size_t pick = std::min(
5359 live.size() - 1,
static_cast<std::size_t
>(uniform01(g_spn) * live.size()));
5360 live.erase(live.begin() +
static_cast<std::ptrdiff_t
>(pick));
5362 while (live.size() < target) {
5364 e.t = now + spn_law[t][k].next(g_spn) / g;
5369 live.push_back(e.tag);
5400 auto job_state_check = [&](std::size_t i) {
5401 const StationState& s = S[i];
5403 if (!place_nodes.empty() || !regions.empty() || !sn.fj.empty() || has_spawn)
5404 why =
"Petri-net places, finite capacity regions, fork-join or spawning";
5405 else if (has_removal_signal || has_sync_call)
5406 why =
"G-network removal signals or synchronous calls";
5407 else if (has_gd || slotted)
5408 why =
"global state dependence or slotted time";
5409 else if (!s.lld.empty() || s.has_cd || s.state_dependent)
5410 why =
"load-, class- or state-dependent service";
5411 else if (s.pas || s.polling || s.bmsp)
5412 why =
"pass-and-swap, polling or batch service";
5413 else if (s.has_setup || s.has_breakdown || s.has_pools || s.has_parallelism)
5414 why =
"setup, breakdowns, server pools or server parallelism";
5415 else if (s.role == Role::Queue && !s.ps && s.sched != SchedStrategy::FCFS)
5417 " (supported: FCFS, the PS family and Delay)";
5419 for (std::size_t r = 0; r < K && why.empty(); ++r) {
5420 if (s.off[r])
continue;
5421 if ((r < s.has_retrial.size() && s.has_retrial[r]) ||
5422 (r < s.has_patience.size() && s.has_patience[r]) || s.svc[r].time_varying())
5423 why =
"retrial, reneging or a time-varying service law";
5427 "SolverLDES (native engine): init_jobs/track_jobs do not support station '" +
5428 sn.stations[i].name +
"': " + why);
5447 const Matrix<double>& J = o.init_jobs;
5449 throw InputError(
"SolverLDES (native engine): init_jobs needs the columns [station, class, "
5450 "remaining, elapsed, qage, sage, order(, id)], got " +
5451 std::to_string(J.cols()) +
" columns");
5452 std::vector<std::size_t> rows(J.rows());
5453 std::iota(rows.begin(), rows.end(), std::size_t(0));
5454 std::stable_sort(rows.begin(), rows.end(), [&](std::size_t a, std::size_t b) {
5455 if (J(a, 0) != J(b, 0)) return J(a, 0) < J(b, 0);
5456 return J(a, 6) < J(b, 6);
5458 std::size_t last_station = M;
5459 double last_arrival = -std::numeric_limits<double>::infinity();
5460 for (std::size_t row : rows) {
5461 const double sv = J(row, 0), kv = J(row, 1);
5462 if (!(sv >= 0.0 && sv <
static_cast<double>(M) && kv >= 0.0 && kv <
static_cast<double>(K)))
5463 throw InputError(
"SolverLDES (native engine): init_jobs row " + std::to_string(row) +
5464 " names station " + std::to_string(sv) +
", class " + std::to_string(kv) +
5465 "; the model has " + std::to_string(M) +
" stations and " +
5466 std::to_string(K) +
" classes");
5467 const std::size_t i =
static_cast<std::size_t
>(sv), k =
static_cast<std::size_t
>(kv);
5468 if (S[i].role != Role::Queue && S[i].role != Role::Delay)
5469 throw InputError(
"SolverLDES (native engine): init_jobs row " + std::to_string(row) +
5470 " places a job at station '" + sn.stations[i].name +
5471 "', which is not a Queue or Delay");
5473 const double remaining = J(row, 2), el = J(row, 3), qage = J(row, 4), sage = J(row, 5);
5474 if (!(el >= 0.0 && qage >= 0.0 && sage >= qage))
5475 throw InputError(
"SolverLDES (native engine): init_jobs row " + std::to_string(row) +
5476 " needs elapsed >= 0 and sage >= qage >= 0");
5479 job.t_sys = now - sage;
5480 if (i != last_station) {
5482 last_arrival = -std::numeric_limits<double>::infinity();
5484 job.t_arr = std::max(now - qage,
5485 std::nextafter(last_arrival, std::numeric_limits<double>::infinity()));
5486 last_arrival = job.t_arr;
5487 job.ext = (J.cols() >= 8) ?
static_cast<std::int64_t
>(J(row, 7)) : -1;
5488 if (remaining > 0.0) {
5489 job.service = remaining;
5490 }
else if (remaining == -2.0) {
5491 const int max_tries = 1000000;
5493 for (; t < max_tries; ++t) {
5494 const double v = S[i].svc[k].next_at(g_svc[i][k], now);
5496 job.service = v - el;
5501 throw InputError(
"SolverLDES (native engine): init_jobs could not draw a residual for "
5502 "class '" + sn.classes[k].name +
"' with elapsed = " +
5503 std::to_string(el) +
" (P(S > elapsed) is below 1e-6); give it explicitly");
5504 }
else if (remaining != -1.0) {
5505 throw InputError(
"SolverLDES (native engine): init_jobs row " + std::to_string(row) +
5506 " has remaining = " + std::to_string(remaining) +
5507 "; use a positive residual, -1 (draw a full service) or -2 (draw given elapsed)");
5509 if (job.service > 0.0) {
5511 job.remaining = job.service;
5514 pending_preset =
true;
5515 if (!admit(i, job, M))
5516 throw InputError(
"SolverLDES (native engine): init_jobs row " + std::to_string(row) +
5517 " does not fit station '" + sn.stations[i].name +
"'");
5519 for (std::size_t r = 0; r < K; ++r) {
5522 for (std::size_t i = 0; i < M; ++i) held += acc.qlen[i][r];
5523 const double want = sn.classes[r].population;
5524 if (std::fabs(held - want) > 0.5)
5525 throw InputError(
"SolverLDES (native engine): init_jobs holds " +
5526 std::to_string(
static_cast<long>(held + 0.5)) +
" jobs of class '" +
5527 sn.classes[r].name +
"' against a population of " +
5528 std::to_string(
static_cast<long>(want + 0.5)));
5532 const std::vector<double> placement =
5533 init_jobs_on ? std::vector<double>()
5535 const bool warm_start = !init_jobs_on && !placement.empty() && place_nodes.empty();
5537 for (std::size_t i = 0; i < M; ++i) {
5538 if (S[i].role != Role::Queue && S[i].role != Role::Delay)
continue;
5539 for (std::size_t k = K; k-- > 0;) {
5540 const std::size_t idx = i * K + k;
5541 if (idx >= placement.size())
continue;
5542 const double v = placement[idx];
5543 if (!(v > 0.0))
continue;
5544 const std::size_t count =
static_cast<std::size_t
>(v);
5545 for (std::size_t j = 0; j < count; ++j) {
5549 if (!admit(i, job, M))
5550 throw InputError(
"SolverLDES (native engine): the warm-start placement "
5551 "of class '" + sn.classes[k].name +
"' does not fit "
5552 "station '" + sn.stations[i].name +
"'");
5558 for (std::size_t r = 0; r < K; ++r) {
5561 for (std::size_t i = 0; i < M; ++i) held += acc.qlen[i][r];
5562 const double want = sn.classes[r].population;
5563 if (std::fabs(held - want) > 0.5)
5564 throw InputError(
"SolverLDES (native engine): the warm-start placement holds " +
5565 std::to_string(
static_cast<long>(held + 0.5)) +
" jobs of class '" +
5566 sn.classes[r].name +
"' against a population of " +
5567 std::to_string(
static_cast<long>(want + 0.5)));
5573 for (std::size_t i = 0; i < M; ++i) {
5574 if (!S[i].has_breakdown)
continue;
5576 e.t = S[i].failure_time.next_at(g_aux[i], 0.0);
5582 std::vector<Sampler> arrival(K);
5583 std::vector<double> lambda(K, 0.0);
5592 std::vector<Sampler> arr_batch(K);
5593 std::vector<bool> has_arr_batch(K,
false);
5594 std::vector<Rng> g_arr_batch;
5595 g_arr_batch.reserve(K);
5596 for (std::size_t k = 0; k < K; ++k)
5597 g_arr_batch.push_back(Rng(seed_ll, 991LL + 97LL * (
static_cast<long long>(k) + 1LL)));
5598 if (source_st < M) {
5599 const auto& sb = sn.stations[source_st].arrival_batch;
5600 for (std::size_t r = 0; r < K && r < sb.size(); ++r) {
5601 if (sb[r].disabled)
continue;
5602 arr_batch[r] = Sampler(sb[r],
"the arrival batch of class '" + sn.classes[r].name +
5603 "' at station '" + sn.stations[source_st].name +
"'");
5604 has_arr_batch[r] =
true;
5621 std::vector<std::size_t> mark_class;
5622 std::vector<double> mark_rate;
5623 std::size_t mark_carrier = K;
5624 if (source_st < M) {
5625 for (std::size_t c : sn.stations[source_st].marked_classes)
5626 if (c >= 1 && c <= K) mark_class.push_back(c - 1);
5635 auto nominal_of = [](
const std::vector<T>& bp,
const std::vector<Matrix<T>>& segs) {
5636 const T zero = num_traits<T>::from_int(0);
5637 Matrix<T> acc(segs[0].rows(), segs[0].cols(), zero);
5639 for (std::size_t k = 0; k + 1 < bp.size(); ++k) total += T(bp[k + 1] - bp[k]);
5640 for (std::size_t k = 0; k < segs.size(); ++k) {
5641 const T w = T(T(bp[k + 1] - bp[k]) / total);
5642 for (std::size_t a = 0; a < acc.rows(); ++a)
5643 for (std::size_t b = 0; b < acc.cols(); ++b) acc(a, b) += w * segs[k](a, b);
5658 auto mean_batch_of = [&](
const lang::Distrib<T>& cd) ->
double {
5659 std::vector<Matrix<T>> blocks;
5660 std::vector<std::size_t> sizes;
5662 for (std::size_t b = 0; b < cd.Dmark.size(); ++b) {
5663 blocks.push_back(cd.Dmark[b]);
5664 sizes.push_back(b + 1);
5666 }
else if (cd.has_batch_schedule()) {
5667 for (std::size_t c = 0; c < cd.sched_Dbatch.size(); ++c)
5668 for (std::size_t b = 0; b < cd.sched_Dbatch[c].size(); ++b) {
5669 blocks.push_back(nominal_of(cd.sched_bp, cd.sched_Dbatch[c][b]));
5670 sizes.push_back(b + 1);
5673 if (blocks.empty())
return 1.0;
5679 double epochs = 0.0, jobs = 0.0;
5680 for (std::size_t b = 0; b < blocks.size(); ++b) {
5681 const double v = num_traits<T>::to_double(lam[b]);
5683 jobs +=
static_cast<double>(sizes[b]) * v;
5685 return (epochs > 0.0) ? (jobs / epochs) : 1.0;
5687 if (!mark_class.empty()) {
5688 mark_carrier = mark_class[0];
5689 const lang::Distrib<T>& cd = sn.service[source_st][mark_carrier];
5690 if (cd.Dmark.size() < mark_class.size())
5691 throw InputError(
"SolverLDES (native engine): the Source binds " +
5692 std::to_string(mark_class.size()) +
5693 " marks but its arrival process carries only " +
5694 std::to_string(cd.Dmark.size()) +
" marked blocks");
5695 if (cd.has_batch_schedule()) {
5700 const std::size_t B = cd.sched_Dbatch[0].size();
5701 std::vector<Matrix<T>> blocks;
5702 for (std::size_t c = 0; c < cd.sched_Dbatch.size(); ++c)
5703 for (std::size_t b = 0; b < B; ++b)
5704 blocks.push_back(nominal_of(cd.sched_bp, cd.sched_Dbatch[c][b]));
5710 for (std::size_t k = 0; k < mark_class.size(); ++k) {
5712 for (std::size_t b = 0; b < B; ++b)
5713 jobs +=
static_cast<double>(b + 1) *
5714 num_traits<T>::to_double(lam[k * B + b]);
5715 mark_rate.push_back(jobs);
5723 for (std::size_t k = 0; k < mark_class.size(); ++k)
5724 mark_rate.push_back(num_traits<T>::to_double(lam[k]));
5728 const std::vector<std::size_t>& mkc = mark_class;
5729 auto marked_non_carrier = [&mkc, &mark_carrier](std::size_t r) {
5730 if (mkc.empty() || r == mark_carrier)
return false;
5731 for (std::size_t c : mkc)
5732 if (c == r)
return true;
5736 for (std::size_t r = 0; r < K; ++r) {
5739 throw InputError(
"SolverLDES (native engine): an open class with no Source");
5740 if (S[source_st].off[r])
continue;
5741 arrival[r] = S[source_st].svc[r];
5745 lambda[r] = mean_batch_of(sn.service[source_st][r]) / arrival[r].mean();
5746 if (has_arr_batch[r]) lambda[r] = arr_batch[r].mean() / arrival[r].mean();
5747 for (std::size_t k = 0; k < mark_class.size(); ++k)
5748 if (mark_class[k] == r) lambda[r] = mark_rate[k];
5750 if (marked_non_carrier(r))
continue;
5752 e.t = slot_snap(arrival[r].next_at(g_arr[r], 0.0),
"interarrival time");
5754 e.station = source_st;
5758 if (warm_start || init_jobs_on)
continue;
5759 const double n = sn.classes[r].population;
5760 if (!(n >= 0.0) || !std::isfinite(n))
5761 throw InputError(
"SolverLDES (native engine): class '" + sn.classes[r].name +
5762 "' has a population that is not a finite count");
5763 const std::size_t refst = sn.classes[r].refstat;
5764 if (refst == 0 || refst > M)
5765 throw InputError(
"SolverLDES (native engine): class '" + sn.classes[r].name +
5766 "' has no reference station");
5776 if (sn.stations[refst - 1].nodetype == NodeType::Place)
continue;
5777 const std::size_t count =
static_cast<std::size_t
>(n + 0.5);
5778 for (std::size_t j = 0; j < count; ++j) {
5782 if (!admit(refst - 1, job, M))
5783 throw InputError(
"SolverLDES (native engine): the initial population of "
5784 "class '" + sn.classes[r].name +
5785 "' does not fit its reference station's capacity");
5792 cnvg.init(M, K, o, max_events);
5793 std::vector<std::size_t> servers_of(M, 1);
5794 std::vector<std::vector<bool>> off_of(M, std::vector<bool>(K,
true));
5795 for (std::size_t i = 0; i < M; ++i) {
5796 servers_of[i] = S[i].nservers;
5797 for (std::size_t r = 0; r < K; ++r)
5798 off_of[i][r] = S[i].off[r] || S[i].role == Role::Source ||
5799 S[i].role == Role::Synchronization;
5801 std::uint64_t last_cnvg_events = 0;
5802 bool converged =
false;
5807 for (std::size_t p = 0; p < place_nodes.size(); ++p)
5808 if (qplace[p]) place_try_start(p);
5813 Observations obs(M, K, o.tranfilter ==
"mser5", o.mserbatch > 0 ? o.mserbatch : 5);
5814 for (std::size_t i = 0; i < M; ++i)
5815 if (S[i].role == Role::Source || S[i].role == Role::Synchronization) obs.in_mser[i] = 0;
5816 std::uint64_t mser_interval = max_events / 1000ULL;
5817 if (mser_interval < 1) mser_interval = 1;
5818 std::uint64_t last_mser_events = 0;
5831 const bool transient_run = o.has_timespan && std::isfinite(o.t1) && o.t1 > o.t0;
5832 const double horizon = transient_run ? o.t1 : std::numeric_limits<double>::infinity();
5835 const double tran_points = (o.tranobs > 0) ?
static_cast<double>(o.tranobs) : 1000.0;
5836 const double tran_interval = transient_run ? (o.t1 - o.t0) / tran_points : 0.0;
5837 double next_tran_sample = transient_run ? o.t0 + tran_interval : 0.0;
5838 double last_tran_time = transient_run ? o.t0 : 0.0;
5839 std::vector<double> tran_times;
5840 std::vector<std::vector<std::vector<double>>> tran_q(M, std::vector<std::vector<double>>(K));
5841 std::vector<std::vector<std::vector<double>>> tran_u(M, std::vector<std::vector<double>>(K));
5842 std::vector<std::vector<std::vector<double>>> tran_t(M, std::vector<std::vector<double>>(K));
5843 std::vector<std::vector<double>> last_tran_q(M, std::vector<double>(K, 0.0));
5844 std::vector<std::vector<double>> last_tran_b(M, std::vector<double>(K, 0.0));
5845 std::vector<std::vector<double>> last_tran_c(M, std::vector<double>(K, 0.0));
5856 std::map<std::vector<int>,
double> histogram;
5857 std::vector<std::pair<double, std::vector<int>>> trajectory;
5858 double hist_last = 0.0;
5861 const bool want_traj = transient_run || o.export_trajectory;
5862 const bool want_hist = transient_run || o.export_histogram || want_traj;
5864 auto joint_state = [&]() {
5865 std::vector<int> row(M * K, 0);
5866 for (std::size_t i = 0; i < M; ++i)
5867 for (std::size_t r = 0; r < K; ++r)
5868 row[i * K + r] =
static_cast<int>(acc.qlen[i][r] + acc.held[i][r] + 0.5);
5871 auto hist_accumulate = [&]() {
5872 if (!want_hist)
return;
5873 const double dt = now - hist_last;
5875 const std::vector<int> row = joint_state();
5876 histogram[row] += dt;
5877 if (want_traj) trajectory.push_back(std::make_pair(hist_last, row));
5882 auto tran_sample = [&]() {
5883 const double dt = now - last_tran_time;
5884 tran_times.push_back(now);
5885 for (std::size_t i = 0; i < M; ++i)
5886 for (std::size_t r = 0; r < K; ++r) {
5888 tran_q[i][r].push_back((acc.tot_qlen[i][r] - last_tran_q[i][r]) / dt);
5889 const double c = S[i].util_peak;
5890 tran_u[i][r].push_back(
5891 (S[i].role == Role::Delay)
5892 ? (acc.tot_qlen[i][r] - last_tran_q[i][r]) / dt
5893 : (acc.tot_busy[i][r] - last_tran_b[i][r]) / (dt * c));
5894 tran_t[i][r].push_back((acc.completed[i][r] - last_tran_c[i][r]) / dt);
5896 tran_q[i][r].push_back(acc.qlen[i][r] + acc.held[i][r]);
5897 tran_u[i][r].push_back(0.0);
5898 tran_t[i][r].push_back(0.0);
5900 last_tran_q[i][r] = acc.tot_qlen[i][r];
5901 last_tran_b[i][r] = acc.tot_busy[i][r];
5902 last_tran_c[i][r] = acc.completed[i][r];
5904 last_tran_time = now;
5908 while (!evq.empty() && (transient_run || total_completions < max_events)) {
5909 const Event ev = evq.top();
5910 if (transient_run && ev.t > horizon)
break;
5917 if (want_hist) hist_accumulate();
5918 if (transient_run) {
5919 while (now >= next_tran_sample && next_tran_sample <= horizon) {
5920 const double save = now;
5921 now = next_tran_sample;
5922 for (std::size_t a = 0; a < M; ++a) {
5923 if (S[a].ps) ps_advance(a);
5924 for (std::size_t r = 0; r < K; ++r) {
5925 acc.update_qlen(a, r, now);
5926 acc.update_busy(a, r, now);
5930 next_tran_sample += tran_interval;
5935 if (ev.kind == EV_BREAKDOWN) {
5940 StationState& s = S[ev.station];
5942 ps_advance(ev.station);
5944 sd_advance(ev.station);
5947 ps_reschedule(ev.station);
5949 sd_reschedule(ev.station);
5959 nxt.t = now + (s.up ? s.failure_time.next_at(g_aux[ev.station], now)
5960 : s.repair_time.next_at(g_aux[ev.station], now));
5962 nxt.station = ev.station;
5967 if (ev.kind == EV_FIRING) {
5969 if (ev.station >= transitions.size())
continue;
5970 std::vector<std::uint64_t>& live = spn_clocks[ev.station][ev.cls];
5971 const auto it = std::find(live.begin(), live.end(), ev.tag);
5972 if (it == live.end())
continue;
5974 SpnTransition& tr = transitions[ev.station];
5975 const SpnMode& m = tr.modes[ev.cls];
5978 tr.fired[ev.cls] += 1.0;
5979 ++total_completions;
5994 if (ev.kind == EV_PLACE_SVC) {
5995 const std::size_t p = ev.station;
5996 const std::size_t r = ev.cls;
5997 const std::size_t ist = place_station[p];
5999 acc.update_busy(ist, r, now);
6000 place_insvc[p][r] -= 1.0;
6001 acc.busy[ist][r] = place_insvc[p][r];
6010 if (ev.kind == EV_SWITCHOVER) {
6011 StationState& sp = S[ev.station];
6012 sp.poll_switching =
false;
6013 sp.poll_at = ev.cls;
6014 poll_open_visit(sp);
6016 for (
const Job& j : sp.buffer)
6017 if (j.cls == sp.poll_at) {
6022 poll_serve(ev.station);
6024 poll_advance(ev.station);
6028 if (ev.kind == EV_RETRIAL) {
6029 StationState& sr = S[ev.station];
6030 const std::size_t c = ev.cls;
6031 if (ev.tag == 0) sr.const_pending[c] =
false;
6032 std::deque<Job>& ob = sr.orbit[c];
6033 std::size_t at = ob.size();
6035 if (!ob.empty()) at = 0;
6037 for (std::size_t j = 0; j < ob.size(); ++j)
6038 if (ob[j].
id == ev.tag) {
6043 if (at == ob.size())
continue;
6044 if (retrial_server_free(ev.station, c)) {
6045 const Job job = ob[at];
6046 orbit_leave(ev.station, c, at,
false);
6047 sr.retried[c] += 1.0;
6050 pending_retry =
true;
6051 pending_orbit_exit =
true;
6052 admit(ev.station, job, M);
6054 ob[at].attempts += 1;
6055 if (sr.max_attempts[c] > 0 && ob[at].attempts >= sr.max_attempts[c]) {
6056 orbit_leave(ev.station, c, at,
true);
6057 sr.retrial_lost[c] += 1.0;
6058 dropped[ev.station][c] += 1.0;
6059 }
else if (sr.retrial_policy[c] == 2) {
6061 const Job jb = ob[at];
6062 ob.erase(ob.begin() +
static_cast<std::ptrdiff_t
>(at));
6065 arm_retrial(ev.station, c, ob[at].id);
6068 if (sr.retrial_policy[c] == 2) arm_retrial(ev.station, c, 0);
6072 if (ev.kind == EV_ORBIT_ABANDON) {
6073 StationState& sr = S[ev.station];
6074 const std::size_t c = ev.cls;
6075 std::deque<Job>& ob = sr.orbit[c];
6076 std::size_t at = ob.size();
6077 for (std::size_t j = 0; j < ob.size(); ++j)
6078 if (ob[j].
id == ev.tag) {
6082 if (at == ob.size())
continue;
6083 orbit_leave(ev.station, c, at,
true);
6084 dropped[ev.station][c] += 1.0;
6085 if (sr.retrial_policy[c] == 2) arm_retrial(ev.station, c, 0);
6089 if (ev.kind == EV_SETUP) {
6090 StationState& ss = S[ev.station];
6095 for (std::size_t r = 0; r < K; ++r) held += acc.qlen[ev.station][r];
6096 if (held > 0.0 || now + 1e-12 < ss.delayoff_at)
continue;
6097 ss.setup_on =
false;
6098 ss.delayoff_at = std::numeric_limits<double>::infinity();
6101 if (ss.setup_on)
continue;
6102 ss.setup_running =
false;
6104 ss.delayoff_at = std::numeric_limits<double>::infinity();
6105 for (std::size_t sl = 0; sl < ss.nservers && !ss.buffer.empty(); ++sl)
6106 if (!ss.server_busy[sl] && !ss.server_blocked[sl] && !ss.server_held[sl] &&
6107 buffer_has_for_slot(ev.station, sl)) {
6108 Job nextjob = buffer_pop_for_slot(ev.station, sl);
6109 start_service(ev.station, sl, nextjob);
6111 sd_reschedule(ev.station);
6115 if (ev.kind == EV_RENEGE) {
6119 StationState& sr = S[ev.station];
6120 std::size_t at = sr.buffer.size();
6121 for (std::size_t j = 0; j < sr.buffer.size(); ++j)
6122 if (sr.buffer[j].id == ev.tag) {
6126 if (at == sr.buffer.size())
continue;
6127 const std::size_t rc = sr.buffer[at].cls;
6128 sd_advance(ev.station);
6129 acc.update_qlen(ev.station, rc, now);
6130 acc.qlen[ev.station][rc] -= 1.0;
6131 bp.track(ev.station, rc, -1, now);
6132 sr.buffer.erase(sr.buffer.begin() +
static_cast<std::ptrdiff_t
>(at));
6133 if (sr.sched != SchedStrategy::FSP)
6134 std::make_heap(sr.buffer.begin(), sr.buffer.end(), sr.cmp);
6135 reneged[ev.station][rc] += 1.0;
6136 last_lost[ev.station][rc] = now;
6137 sd_reschedule(ev.station);
6141 if (ev.kind == EV_ARRIVAL) {
6146 const int mark = arrival[ev.cls].last_mark();
6151 std::size_t batch = 1;
6152 if (arrival[ev.cls].batched()) {
6153 const int b = arrival[ev.cls].last_batch();
6154 if (b > 0) batch =
static_cast<std::size_t
>(b);
6155 }
else if (has_arr_batch[ev.cls]) {
6156 const double sampled = arr_batch[ev.cls].next(g_arr_batch[ev.cls]);
6157 const double rounded = std::floor(sampled + 0.5);
6158 if (!(rounded >= 1.0) || std::fabs(sampled - rounded) > 1e-9)
6160 "SolverLDES (native engine): the arrival batch of class '" +
6161 sn.classes[ev.cls].name +
"' sampled " + std::to_string(sampled) +
6162 ", which is not a positive integer; a batch-size law must be "
6163 "supported on {1,2,...}");
6164 batch =
static_cast<std::size_t
>(rounded);
6167 const double gap = slot_snap(arrival[ev.cls].next_at(g_arr[ev.cls], now),
6168 "interarrival time");
6171 if (!(gap > 0.0) && arrival[ev.cls].time_varying())
continue;
6174 nxt.station = ev.station;
6180 std::size_t arv_cls = ev.cls;
6181 if (mark > 0 &&
static_cast<std::size_t
>(mark) <= mark_class.size())
6182 arv_cls = mark_class[
static_cast<std::size_t
>(mark) - 1];
6189 for (std::size_t b = 0; b < batch; ++b) {
6190 const RouteEntry e = draw_node_route(sn.station_to_node[ev.station], arv_cls);
6194 deliver(e.node, job, M);
6200 const std::size_t i = ev.station;
6201 StationState& s = S[i];
6215 if (s.bmsp_tag != ev.tag)
continue;
6216 const std::size_t r = s.bmsp_cls;
6217 const std::size_t want = s.bmsp_pending;
6218 std::vector<Job> served;
6219 for (std::size_t d = 0; d < want && !s.buffer.empty(); ++d)
6220 served.push_back(buffer_pop(i));
6221 for (std::size_t d = 0; d < served.size(); ++d) {
6222 const Job& sj = served[d];
6223 acc.update_qlen(i, sj.cls, now);
6224 acc.qlen[i][sj.cls] -= 1.0;
6225 bp.track(i, sj.cls, -1, now);
6226 acc.completed[i][sj.cls] += 1.0;
6227 acc.tally_resp(i, sj.cls, now - sj.t_arr + sj.region_wait);
6228 if (want_respt) resp_samples[i][sj.cls].push_back(now - sj.t_arr + sj.region_wait);
6229 ++total_completions;
6234 if (!s.buffer.empty()) {
6237 acc.update_busy(i, r, now);
6238 if (acc.busy[i][r] > 0.0) acc.busy[i][r] -= 1.0;
6239 s.server_busy[0] =
false;
6244 for (std::size_t d = 0; d < served.size(); ++d) {
6245 const Job& sj = served[d];
6246 const RouteEntry re2 = draw_node_route(sn.station_to_node[i], sj.cls);
6247 if (region_of[i] >= 0) {
6248 const std::size_t rg =
static_cast<std::size_t
>(region_of[i]);
6249 const bool inside = !re2.sink && re2.station < M &&
6250 region_of[re2.station] == region_of[i];
6252 regions[rg].update(now);
6253 regions[rg].leave(sj.cls);
6254 regions[rg].completed[sj.cls] += 1.0;
6259 moved.cls = re2.cls;
6260 moved.t_sys = sj.t_sys;
6261 moved.parent = sj.parent;
6262 moved.call = sj.call;
6263 deliver(re2.node, moved, i);
6265 if (total_completions >= max_events)
break;
6269 if (s.role == Role::Delay) {
6271 if (track_delay_jobs) {
6277 std::size_t at = delay_live[i].size();
6278 for (std::size_t j = 0; j < delay_live[i].size(); ++j)
6279 if (delay_live[i][j].
id == job.id) {
6283 if (at == delay_live[i].size())
continue;
6284 if (delay_live[i][at].tag != ev.tag)
continue;
6285 delay_live[i].erase(delay_live[i].begin() +
static_cast<std::ptrdiff_t
>(at));
6289 std::size_t idx = s.ps_jobs.size();
6290 for (std::size_t j = 0; j < s.ps_jobs.size(); ++j)
6291 if (s.ps_jobs[j].tag == ev.tag) {
6298 if (idx == s.ps_jobs.size())
continue;
6299 const PsJob& pj = s.ps_jobs[idx];
6301 job.t_arr = pj.t_arr;
6302 job.region_wait = pj.region_wait;
6303 job.t_sys = pj.t_sys;
6304 job.parent = pj.parent;
6306 s.ps_jobs.erase(s.ps_jobs.begin() +
static_cast<std::ptrdiff_t
>(idx));
6310 if (s.pas_tag != ev.tag || s.pas_list.empty())
continue;
6311 job = pas_complete(i);
6312 acc.update_busy(i, job.cls, now);
6313 if (acc.busy[i][job.cls] > 0.0) acc.busy[i][job.cls] -= 1.0;
6318 if (!s.server_busy[ev.slot] || s.server_tag[ev.slot] != ev.tag)
continue;
6320 job = s.server[ev.slot];
6321 s.server_busy[ev.slot] =
false;
6322 s.server_tag[ev.slot] = 0;
6323 const double freed = release_slots(i, ev.slot);
6324 acc.update_busy(i, job.cls, now);
6325 acc.busy[i][job.cls] -= freed;
6329 acc.update_qlen(i, job.cls, now);
6330 acc.qlen[i][job.cls] -= 1.0;
6331 bp.track(i, job.cls, -1, now);
6332 acc.completed[i][job.cls] += 1.0;
6333 acc.tally_resp(i, job.cls, now - job.t_arr + job.region_wait);
6334 if (want_respt) resp_samples[i][job.cls].push_back(now - job.t_arr + job.region_wait);
6335 ++total_completions;
6336 if (track_jobs && job.ext >= 0)
6337 job_departure_log.push_back({
static_cast<double>(job.ext),
static_cast<double>(i),
6338 static_cast<double>(job.cls), now});
6341 if (s.lps_limit > 0 && !s.buffer.empty() && s.ps_jobs.size() < s.lps_limit) {
6342 Job nextjob = buffer_pop(i);
6344 pj.cls = nextjob.cls;
6345 pj.priority = nextjob.priority;
6346 pj.t_arr = nextjob.t_arr;
6347 pj.region_wait = nextjob.region_wait;
6348 pj.t_sys = nextjob.t_sys;
6349 pj.total = nextjob.preset ? nextjob.elapsed + nextjob.service : nextjob.service;
6350 pj.remaining = nextjob.service;
6351 pj.parent = nextjob.parent;
6352 pj.ext = nextjob.ext;
6353 s.ps_jobs.push_back(pj);
6358 const RouteEntry
re = draw_node_route(sn.station_to_node[i], job.cls);
6359 const bool to_sink =
re.sink;
6360 const std::size_t dst =
re.station, dcls =
re.cls;
6374 const RouteEntry fb = has_immfeed ? resolve_forced_dest(re) :
re;
6375 const std::size_t fcls = fb.cls;
6376 if (has_immfeed && !fb.sink && fb.station == i && s.role == Role::Queue && !s.ps &&
6377 !s.pas && ev.slot < s.nservers && fcls < sn.immfeed[i].size() && sn.immfeed[i][fcls]) {
6379 if (!fb_loggers.empty()) {
6380 const std::uint64_t walk = log_walk_id++;
6381 for (std::size_t q = 0; q < fb_loggers.size(); ++q)
6382 log_passage(fb_loggers[q].first, fb_loggers[q].second, walk);
6386 fed.t_sys = job.t_sys;
6387 fed.parent = job.parent;
6389 fed.priority = classprio[fcls];
6390 if (draw_at_arrival(i, fcls)) {
6391 fed.service = slot_snap(s.svc[fcls].next_at(g_svc[i][fcls], now),
"service time");
6392 fed.remaining = fed.service;
6395 fed.rank = g_siro.aux.next_double();
6396 fed.deadline = now + classdeadline[fcls];
6399 acc.update_qlen(i, fcls, now);
6400 acc.qlen[i][fcls] += 1.0;
6406 acc.arrived[i][fcls] += 1.0;
6407 bp.track(i, fcls, +1, now);
6408 start_service(i, ev.slot, fed);
6410 if (total_completions >= max_events)
break;
6429 if (has_sync_call && sync_reply[job.cls] < K &&
6430 !is_reply_signal[resolve_final_cls(
re.node, dcls)] &&
6431 s.role == Role::Queue && !s.ps && !s.pas && ev.slot < s.nservers) {
6432 const std::uint64_t call = ++next_call;
6438 pending_reply[call] = pc;
6439 s.server_held[ev.slot] =
true;
6440 s.held_cls[ev.slot] = job.cls;
6441 acc.update_busy(i, job.cls, now);
6442 acc.busy[i][job.cls] += 1.0;
6445 acc.update_qlen(i, job.cls, now);
6446 acc.held[i][job.cls] += 1.0;
6449 moved.t_sys = job.t_sys;
6450 moved.parent = job.parent;
6452 deliver(
re.node, moved, i);
6453 if (total_completions >= max_events)
break;
6462 if (!to_sink && dst < M && s.role == Role::Queue && !s.ps && ev.slot < s.nservers &&
6463 !dest_has_room(dst, dcls)) {
6468 s.server_blocked[ev.slot] =
true;
6469 s.blocked_job[ev.slot] = held;
6470 s.blocked_dest[ev.slot] = dst;
6471 s.blocked_dest_cls[ev.slot] = dcls;
6478 acc.update_qlen(dst, dcls, now);
6479 acc.qlen[dst][dcls] += 1.0;
6480 S[dst].blocked_at[dcls] += 1.0;
6486 acc.update_busy(i, job.cls, now);
6487 acc.busy[i][job.cls] += 1.0;
6488 blocked_count[i][job.cls] += 1.0;
6500 }
else if (s.role == Role::Queue && !s.ps && !s.pas && !s.server_blocked.empty() &&
6501 ev.slot < s.nservers && !s.server_blocked[ev.slot] &&
6502 !s.server_held[ev.slot]) {
6503 serve_next_on_slot(i, ev.slot);
6510 std::uint64_t spawn_parent = 0;
6511 if (has_spawn && spawn_of[job.cls] < K && spawn_joins_at[spawn_of[job.cls]]) {
6512 spawn_parent = job.parent;
6519 const int dep_rg = region_of[i];
6521 regions[
static_cast<std::size_t
>(dep_rg)].update(now);
6522 regions[
static_cast<std::size_t
>(dep_rg)].leave(job.cls);
6528 if (has_spawn && spawn_of[job.cls] < K) {
6529 const std::size_t scls = spawn_of[job.cls];
6532 spawned.t_sys = now;
6537 spawned.parent = spawn_parent;
6538 admit(i, spawned, M);
6543 moved.t_sys = job.t_sys;
6544 moved.parent = job.parent;
6548 moved.call = job.call;
6549 hop_old_cls = job.cls;
6550 hop_internal =
false;
6551 deliver(
re.node, moved, i);
6557 const bool stays_inside = hop_internal;
6558 hop_internal =
false;
6559 if (!stays_inside) {
6560 regions[
static_cast<std::size_t
>(dep_rg)].completed[job.cls] += 1.0;
6561 release_region(
static_cast<std::size_t
>(dep_rg));
6566 if (s.has_setup && s.setup_on) {
6568 for (std::size_t r = 0; r < K; ++r) held += acc.qlen[i][r];
6569 if (!(held > 0.0)) {
6571 s.delayoff_time.disabled() ? 0.0 : s.delayoff_time.next(g_aux[i]);
6572 s.delayoff_at = now + idle;
6574 e.t = s.delayoff_at;
6584 if (total_completions >= max_events)
break;
6585 if (cnvg.enabled() && (total_completions - last_cnvg_events) >= cnvg.interval()) {
6586 for (std::size_t a = 0; a < M; ++a) {
6587 if (S[a].ps) ps_advance(a);
6588 for (std::size_t r = 0; r < K; ++r) {
6589 acc.update_qlen(a, r, now);
6590 acc.update_busy(a, r, now);
6593 cnvg.finalize_batch(acc, servers_of, now);
6594 last_cnvg_events = total_completions;
6595 if (cnvg.converged(off_of)) {
6600 if ((total_completions - last_mser_events) >= mser_interval) {
6601 for (std::size_t a = 0; a < M; ++a) {
6602 if (S[a].ps) ps_advance(a);
6603 for (std::size_t r = 0; r < K; ++r) {
6604 acc.update_qlen(a, r, now);
6605 acc.update_busy(a, r, now);
6608 obs.collect(acc, now);
6609 last_mser_events = total_completions;
6614 for (std::size_t i = 0; i < M; ++i) {
6615 if (S[i].ps) ps_advance(i);
6616 for (std::size_t r = 0; r < K; ++r) {
6617 acc.update_qlen(i, r, now);
6618 acc.update_busy(i, r, now);
6623 const Truncation tr = obs.truncate();
6624 const double sim_time = now - tr.warmup_end;
6625 const double elapsed = tr.applied ? (now - obs.time[tr.index]) : sim_time;
6627 const double t_meas = tr.applied ? obs.time[tr.index] : 0.0;
6628 std::vector<std::vector<bool>> lost_unserved(M, std::vector<bool>(K,
false));
6629 for (std::size_t i = 0; i < M; ++i)
6630 for (std::size_t r = 0; r < K; ++r) lost_unserved[i][r] = last_lost[i][r] >= t_meas;
6636 res.nchains = sn.nchains;
6637 for (std::size_t i = 0; i < M; ++i) res.station_names.push_back(sn.stations[i].name);
6638 for (std::size_t r = 0; r < K; ++r) res.class_names.push_back(sn.classes[r].name);
6639 res.QN = Matrix<double>(M, K, 0.0);
6640 res.UN = Matrix<double>(M, K, 0.0);
6641 res.RN = Matrix<double>(M, K, 0.0);
6642 res.TN = Matrix<double>(M, K, 0.0);
6643 res.CN = Matrix<double>(1, K, 0.0);
6644 res.XN = Matrix<double>(1, K, 0.0);
6645 res.AN = Matrix<double>(M, K, 0.0);
6646 res.WN = Matrix<double>(M, K, 0.0);
6649 Matrix<double> u_law(M, K, 0.0);
6651 for (std::size_t i = 0; i < M; ++i) {
6652 for (std::size_t r = 0; r < K; ++r) {
6653 if (S[i].role == Role::Source) {
6654 res.TN(i, r) = lambda[r];
6657 if (S[i].off[r])
continue;
6658 if (elapsed > 0.0) {
6659 const double q0 = tr.applied ? obs.qt[i][r][tr.index] : 0.0;
6660 const double b0 = tr.applied ? obs.bt[i][r][tr.index] : 0.0;
6661 const double c0 = tr.applied ? obs.cmp[i][r][tr.index] : 0.0;
6662 res.QN(i, r) = (acc.tot_qlen[i][r] - q0) / elapsed;
6663 res.TN(i, r) = (acc.completed[i][r] - c0) / elapsed;
6664 if (S[i].role == Role::Delay) {
6673 res.UN(i, r) = res.TN(i, r) * S[i].class_mean[r] / S[i].gd_peak[r];
6674 u_law(i, r) = S[i].class_mean[r] / S[i].gd_peak[r];
6675 }
else if (sn.stations[i].nodetype == NodeType::Place && !sn.is_queueing_place(i)) {
6677 res.UN(i, r) = res.QN(i, r);
6678 if (res.TN(i, r) > 0.0) u_law(i, r) = res.QN(i, r) / res.TN(i, r);
6679 }
else if (S[i].has_cd || has_gd) {
6699 (S[i].has_cd && r < S[i].cd_peak_class.size()) ? S[i].cd_peak_class[r]
6700 : (S[i].has_cd || !S[i].lld.empty()) ? S[i].util_peak : 1.0;
6701 const double peak = base * S[i].gd_peak[r];
6703 (peak > 0.0) ? res.TN(i, r) * S[i].class_mean[r] / peak : 0.0;
6704 if (peak > 0.0) u_law(i, r) = S[i].class_mean[r] / peak;
6706 res.UN(i, r) = (acc.tot_busy[i][r] - b0) / (elapsed * S[i].util_peak);
6707 if (S[i].util_peak > 0.0) u_law(i, r) = S[i].class_mean[r] / S[i].util_peak;
6710 if (lost_unserved[i][r]) {
6714 res.RN(i, r) = (res.TN(i, r) > 0.0) ? res.QN(i, r) / res.TN(i, r) : 0.0;
6715 }
else if (acc.resp_cnt[i][r] > 0.0) {
6716 res.RN(i, r) = acc.resp_sum[i][r] / acc.resp_cnt[i][r];
6717 }
else if (sn.stations[i].nodetype == NodeType::Place && res.TN(i, r) > 0.0) {
6725 res.RN(i, r) = res.QN(i, r) / res.TN(i, r);
6732 const double a0 = tr.applied ? obs.arv[i][r][tr.index] : 0.0;
6733 if (elapsed > 0.0) res.AN(i, r) = (acc.arrived[i][r] - a0) / elapsed;
6734 res.WN(i, r) = res.RN(i, r);
6742 if (!sn.fj.empty()) {
6743 res.DropRateJoin = Matrix<double>(M, K, 0.0);
6744 for (std::size_t i = 0; i < M; ++i)
6745 for (std::size_t r = 0; r < K; ++r) {
6746 const double d0 = tr.applied ? obs.drp[i][r][tr.index] : 0.0;
6748 res.DropRateJoin(i, r) = (acc.join_dropped[i][r] - d0) / elapsed;
6753 if (!fj_pairs.empty()) {
6754 const std::size_t P = fj_pairs.size();
6755 res.FJQLen = Matrix<double>(P, K, 0.0);
6756 res.FJRespT = Matrix<double>(P, K, 0.0);
6757 res.FJTput = Matrix<double>(P, K, 0.0);
6758 const double t_trunc = tr.applied ? obs.time[tr.index] : -1.0;
6759 for (std::size_t p = 0; p < P; ++p) {
6760 res.fj_forks.push_back(sn.nodes[fj_pairs[p].first - 1].name);
6761 res.fj_joins.push_back(sn.nodes[fj_pairs[p].second - 1].name);
6762 for (std::size_t r = 0; r < K; ++r) {
6763 if (fj_resp_cnt[p][r] > 0.0) res.FJRespT(p, r) = fj_resp_sum[p][r] / fj_resp_cnt[p][r];
6765 for (double tr_k : fj_release[p][r])
6766 if (tr_k > t_trunc) nrel += 1.0;
6767 if (elapsed > 0.0) res.FJTput(p, r) = nrel / elapsed;
6768 res.FJQLen(p, r) = res.FJTput(p, r) * res.FJRespT(p, r);
6772 for (std::size_t r = 0; r < K; ++r) {
6773 if (sim_time > 0.0) res.XN(0, r) = sys_completed[r] / sim_time;
6774 if (sys_resp_cnt[r] > 0.0) res.CN(0, r) = sys_resp_sum[r] / sys_resp_cnt[r];
6784 const std::size_t rs = sn.classes[r].refstat;
6785 if (res.XN(0, r) == 0.0 && rs != 0 && rs <= M &&
6786 sn.stations[rs - 1].nodetype == NodeType::Place)
6787 res.XN(0, r) = res.TN(rs - 1, r);
6791 for (std::size_t ti = 0; ti < bp.targets().size(); ++ti) {
6792 LdesResult::BusyPeriodTarget out;
6793 out.name = bp.targets()[ti].name;
6794 out.stations = bp.targets()[ti].stations;
6795 out.job_class = bp.targets()[ti].job_class;
6796 for (std::size_t oi = 0; oi < static_cast<std::size_t>(bp.orders()); ++oi) {
6797 out.mean.push_back(bp.mean(ti, oi));
6798 out.count.push_back(bp.count(ti, oi));
6800 res.busy_periods.push_back(out);
6804 if (transient_run) {
6805 now = std::min(now, horizon);
6808 res.QNt.assign(M, std::vector<Matrix<double>>(K));
6809 res.UNt.assign(M, std::vector<Matrix<double>>(K));
6810 res.TNt.assign(M, std::vector<Matrix<double>>(K));
6811 for (std::size_t i = 0; i < M; ++i)
6812 for (std::size_t r = 0; r < K; ++r) {
6813 const std::size_t n = tran_times.size();
6816 Matrix<double> q(n, 2, 0.0), u(n, 2, 0.0), t(n, 2, 0.0);
6817 for (std::size_t k = 0; k < n; ++k) {
6818 q(k, 0) = tran_q[i][r][k];
6819 q(k, 1) = tran_times[k];
6820 u(k, 0) = tran_u[i][r][k];
6821 u(k, 1) = tran_times[k];
6822 t(k, 0) = tran_t[i][r][k];
6823 t(k, 1) = tran_times[k];
6829 res.stopping_reason =
"max_time";
6833 res.respTimeSamples.assign(M, std::vector<std::vector<double>>(K));
6834 for (std::size_t i = 0; i < M; ++i)
6835 for (std::size_t r = 0; r < K; ++r) res.respTimeSamples[i][r] = resp_samples[i][r];
6843 if (want_hist && !transient_run) hist_accumulate();
6844 if (!histogram.empty()) {
6845 res.histogram_space = Matrix<double>(histogram.size(), M * K, 0.0);
6846 res.histogram_time = Matrix<double>(histogram.size(), 1, 0.0);
6847 std::size_t row = 0;
6848 for (const auto& kv : histogram) {
6849 for (std::size_t c = 0; c < kv.first.size(); ++c)
6850 res.histogram_space(row, c) = kv.first[c];
6851 res.histogram_time(row, 0) = kv.second;
6855 if (!trajectory.empty()) {
6856 res.traj_space = Matrix<double>(trajectory.size(), M * K, 0.0);
6857 res.traj_time = Matrix<double>(trajectory.size(), 1, 0.0);
6858 for (std::size_t k = 0; k < trajectory.size(); ++k) {
6859 for (std::size_t c = 0; c < trajectory[k].second.size(); ++c)
6860 res.traj_space(k, c) = trajectory[k].second[c];
6861 res.traj_time(k, 0) = trajectory[k].first;
6866 for (
auto& kv : caches) {
6867 CacheState& cs = kv.second;
6868 LdesCacheMetrics cm;
6869 cm.hit = Matrix<double>(1, K, 0.0);
6870 cm.miss = Matrix<double>(1, K, 0.0);
6875 if (cs.has_retrieval) cm.delayed = Matrix<double>(1, K, 0.0);
6876 for (std::size_t r = 0; r < K; ++r) {
6877 const double dl = cs.has_retrieval ? cs.delayed[r] : 0.0;
6878 const double tot = cs.hits[r] + cs.misses[r] + dl;
6880 cm.hit(0, r) = cs.hits[r] /
tot;
6881 cm.miss(0, r) = cs.misses[r] /
tot;
6882 if (cs.has_retrieval) cm.delayed(0, r) = dl /
tot;
6902 if (cs.has_retrieval)
6903 for (std::size_t r = 0; r < K; ++r) reads += cs.hits[r] + cs.misses[r] + cs.delayed[r];
6905 cm.latency = Matrix<double>(1, K, 0.0);
6906 const double mean_wait = (cs.total_fetch_time + cs.delayed_wait) / reads;
6907 for (std::size_t r = 0; r < K; ++r)
6908 cm.latency(0, r) = (cs.hits[r] + cs.misses[r] + cs.delayed[r] > 0.0)
6910 : std::numeric_limits<double>::quiet_NaN();
6916 const std::size_t H = cs.lists.size();
6917 const double nan = std::numeric_limits<double>::quiet_NaN();
6918 cm.hitList = Matrix<double>(K, H, nan);
6919 for (std::size_t r = 0; r < K; ++r) {
6920 const double tot = cs.hits[r] + cs.misses[r] + (cs.has_retrieval ? cs.delayed[r] : 0.0);
6922 for (std::size_t l = 0; l < H; ++l) cm.hitList(r, l) = cs.hits_per_list[r][l] /
tot;
6925 cm.itemProb = Matrix<double>(cs.nitems, H + 1, nan);
6927 for (std::size_t it = 0; it < cs.nitems; ++it) {
6928 double in_cache = 0.0;
6929 for (std::size_t l = 0; l < H; ++l) {
6930 const double p = cs.item_level_time[it][l] / now;
6931 cm.itemProb(it, l + 1) = p;
6934 cm.itemProb(it, 0) = std::max(0.0, 1.0 - in_cache);
6936 if (!cs.item_size.empty()) {
6937 cm.listCost = Matrix<double>(1, H, 0.0);
6938 for (std::size_t l = 0; l < H; ++l)
6939 for (std::size_t it = 0; it < cs.nitems; ++it)
6940 cm.listCost(0, l) += cs.size_of(it) * cm.itemProb(it, l + 1);
6943 res.cache_metrics[sn.nodes[cs.node - 1].name] = cm;
6947 if (!regions.empty()) {
6948 const std::size_t NR = regions.size();
6950 res.QNfcr = Matrix<double>(NR, K, 0.0);
6951 res.TNfcr = Matrix<double>(NR, K, 0.0);
6952 res.WeightNfcr = Matrix<double>(NR, K, 0.0);
6953 res.MemOccNfcr = Matrix<double>(NR, K, 0.0);
6954 res.DropRateNfcr = Matrix<double>(NR, K, 0.0);
6955 for (std::size_t g = 0; g < NR; ++g) {
6956 regions[g].update(now);
6957 for (std::size_t r = 0; r < K; ++r) {
6958 if (sim_time > 0.0) {
6962 res.QNfcr(g, r) = regions[g].tot_jobs[r] / sim_time;
6963 res.WeightNfcr(g, r) = regions[g].tot_weight[r] / sim_time;
6964 res.MemOccNfcr(g, r) = regions[g].tot_mem[r] / sim_time;
6965 res.TNfcr(g, r) = regions[g].completed[r] / sim_time;
6966 res.DropRateNfcr(g, r) = regions[g].dropped[r] / sim_time;
6977 res.total_simulated_events =
static_cast<long long>(total_completions);
6982 double any_balk = 0.0, any_renege = 0.0;
6983 for (std::size_t i = 0; i < M; ++i)
6984 for (std::size_t r = 0; r < K; ++r) {
6985 any_balk += balked[i][r];
6986 any_renege += reneged[i][r];
6990 bool any_orbit =
false;
6991 for (std::size_t i = 0; i < M; ++i)
6992 for (std::size_t r = 0; r < K; ++r) any_orbit = any_orbit || S[i].has_retrial[r];
6994 res.retriedCustomers = Matrix<double>(M, K, 0.0);
6995 res.retrialDropped = Matrix<double>(M, K, 0.0);
6996 res.avgOrbitSize = Matrix<double>(M, K, 0.0);
6997 for (std::size_t i = 0; i < M; ++i) {
6998 const double dt = now - S[i].orbit_last;
6999 for (std::size_t r = 0; r < K; ++r) {
7000 double tot = S[i].tot_orbit[r];
7001 if (dt > 0.0)
tot += S[i].orbit_size[r] * dt;
7002 res.retriedCustomers(i, r) = S[i].retried[r];
7003 res.retrialDropped(i, r) = S[i].retrial_lost[r];
7004 if (sim_time > 0.0) res.avgOrbitSize(i, r) =
tot / sim_time;
7009 if (any_balk > 0.0 || any_renege > 0.0) {
7010 res.balkedCustomers = Matrix<double>(M, K, 0.0);
7011 res.renegedCustomers = Matrix<double>(M, K, 0.0);
7012 res.renegingRate = Matrix<double>(M, K, 0.0);
7013 for (std::size_t i = 0; i < M; ++i)
7014 for (std::size_t r = 0; r < K; ++r) {
7015 res.balkedCustomers(i, r) = balked[i][r];
7016 res.renegedCustomers(i, r) = reneged[i][r];
7017 if (sim_time > 0.0) res.renegingRate(i, r) = reneged[i][r] / sim_time;
7029 if (transient_run) now = horizon;
7030 std::int64_t next_id = 0;
7031 for (std::size_t r = 0; r < o.init_jobs.rows() && o.init_jobs.cols() >= 8; ++r)
7032 next_id = std::max(next_id,
static_cast<std::int64_t
>(o.init_jobs(r, 7)) + 1);
7033 std::vector<std::array<double, 8>> out;
7034 auto emit = [&](std::size_t i,
const Job& j,
double remaining,
double el,
double& order) {
7035 out.push_back({
static_cast<double>(i),
static_cast<double>(j.cls),
7036 std::max(remaining, std::numeric_limits<double>::min()), el, now - j.t_arr,
7037 now - j.t_sys, order,
static_cast<double>(j.ext)});
7040 for (std::size_t i = 0; i < M; ++i) {
7041 StationState& s = S[i];
7042 if (s.role != Role::Queue && s.role != Role::Delay)
continue;
7043 const std::size_t before = out.size();
7045 if (s.role == Role::Delay) {
7046 for (
const Job& d : delay_live[i]) emit(i, d, d.due - now, d.elapsed + d.service - (d.due - now), order);
7049 for (
const PsJob& pj : s.ps_jobs) {
7055 emit(i, j, pj.remaining, pj.total - pj.remaining, order);
7058 std::vector<std::size_t> busy;
7059 for (std::size_t sl = 0; sl < s.nservers; ++sl) {
7060 if (s.server_blocked[sl] || s.server_held[sl])
7061 throw UnsupportedError(ErrorKind{
"InitJobs"},
7062 "SolverLDES (native engine): track_jobs cannot export station '" +
7063 sn.stations[i].name +
"', whose server is blocked or held");
7064 if (s.server_busy[sl]) busy.push_back(sl);
7066 std::sort(busy.begin(), busy.end(), [&](std::size_t a, std::size_t b) {
7067 return s.server[a].t_arr < s.server[b].t_arr;
7069 for (std::size_t sl : busy) {
7070 const Job& j = s.server[sl];
7071 const double done = now - s.server_start[sl];
7072 emit(i, j, j.remaining - done, j.elapsed + done, order);
7076 std::vector<Job> waiting = s.buffer;
7077 for (std::size_t n = waiting.size(); n > 0; --n) {
7078 std::pop_heap(waiting.begin(), waiting.begin() +
static_cast<std::ptrdiff_t
>(n), s.cmp);
7079 const Job& j = waiting[n - 1];
7080 emit(i, j, j.remaining, j.elapsed, order);
7082 if (out.size() > before) job_state_check(i);
7084 for (
const auto& row : out) next_id = std::max(next_id,
static_cast<std::int64_t
>(row[7]) + 1);
7085 res.finalJobs = Matrix<double>(out.size(), 8, 0.0);
7086 for (std::size_t r = 0; r < out.size(); ++r) {
7087 if (out[r][7] < 0.0) out[r][7] =
static_cast<double>(next_id++);
7088 for (std::size_t c = 0; c < 8; ++c) res.finalJobs(r, c) = out[r][c];
7090 res.jobDepartures = Matrix<double>(job_departure_log.size(), 4, 0.0);
7091 for (std::size_t r = 0; r < job_departure_log.size(); ++r)
7092 for (std::size_t c = 0; c < 4; ++c) res.jobDepartures(r, c) = job_departure_log[r][c];
7094 res.converged = converged;
7095 res.convergence_batches = cnvg.batches();
7096 if (!transient_run) res.stopping_reason = converged ?
"convergence" :
"max_events";
7097 res.engine =
"native";
7098 res.method = o.method;