ssa_mlmc_network.h

Go to the documentation of this file.

#ifndef LINE_SOLVERS_SSA_MLMC_NETWORK_H
#define LINE_SOLVERS_SSA_MLMC_NETWORK_H
#include "line/solvers/ssa/ssa_mlmc.h"
#include "line/solvers/ssa/ssa_types.h"
#include "line/lang/qn/feature_set.h"
#include "line/lang/qn/solver_feature_sets.h"
#include "line/solvers/ctmc/solver_ctmc_analyzer.h"
 
namespace line { namespace ssa {
/** The MLMC envelope, single-sourced in solver_feature_sets.h so the supports path and this run-time gate agree. */
inline qn::FeatureSet mlmc_feature_set() { return qn::ssa_mlmc_feature_set(); }
struct MlmcCompiled {
    MlmcProblem problem;
    std::vector<std::string> names;
};
/** Compile fixed routed channels, retaining all station/class dimensions. */
inline MlmcCompiled ssa_mlmc_problem(const qn::NetworkStruct<double>& sn,const std::string& rewardName="") {
    auto support=qn::feature_set_supports("SSA MLMC",mlmc_feature_set(),qn::used_lang_features(sn));
    if(!support.ok) throw qn::feature_refusal(support);
    using NT=lang::NodeType; using SS=lang::SchedStrategy;
    const std::size_t M=sn.stations.size(),R=sn.nclasses;
    for(const auto& s:sn.stations) if(R>1&&s.sched==SS::FCFS)
        throw UnsupportedError(ErrorKind{"MultiClass"}, "SSA MLMC requires a population state; multiclass FCFS needs queue order");
    for(const auto& s:sn.statespace) if(s.second.rows()!=1)
        throw UnsupportedError(ErrorKind{"Model"}, "MLMC requires a deterministic initial state");
    qn::NetState<double> init;
    if(!ctmc::analyzer_detail::default_init_state(sn,init)) throw InputError("Invalid MLMC initial state");
    std::vector<std::size_t> retained,keep,routeIndex(sn.nodes.size());
    for(std::size_t node=0;node<sn.nodes.size();++node) if(sn.nodes[node].station!=0||sn.nodes[node].nodetype==NT::Sink) {
        routeIndex[node]=retained.size(); retained.push_back(node);
        for(std::size_t r=0;r<R;++r) keep.push_back(node*R+r);
    }
    auto routing=mc::dtmc_stochcomp(sn.rtnodes,keep);
    MlmcCompiled out; auto& p=out.problem;
    p.x0.assign(M*R,0); p.nu.resize(M*R);
    std::vector<std::size_t> origins;
    std::vector<double> weights,servers(M);
    std::vector<SS> sched(M); std::vector<bool> source(M);
    std::vector<std::vector<double>> mu(M,std::vector<double>(R));
    for(std::size_t i=0;i<M;++i) {
        std::size_t node=sn.station_to_node[i]; source[i]=sn.nodes[node-1].nodetype==NT::Source;
        sched[i]=sn.stations[i].sched; servers[i]=sn.stations[i].nservers;
        if(!source[i]) {
            auto a=qn::to_marginal_aggr(sn,node,init.local[sn.stateful_index(node)-1]);
            for(std::size_t r=0;r<R;++r) p.x0[i*R+r]=a.second[r];
        }
        for(std::size_t r=0;r<R;++r) {
            double rate=sn.rates(i,r); mu[i][r]=std::isnan(rate)?0:rate;
            if(mu[i][r]==0) continue;
            if(!std::isfinite(rate)||rate<0) throw InputError("MLMC requires finite exponential rates");
            double total=0;
            for(std::size_t dest=0;dest<retained.size()*R;++dest) {
                double w=routing(routeIndex[node-1]*R+r,dest);
                if(!std::isfinite(w)||w<0) throw InputError("Invalid MLMC routing probability");
                total+=w; if(w==0) continue;
                std::size_t j=retained[dest/R],s=dest%R;
                std::vector<double> change(M*R,0);
                if(!source[i]) change[i*R+r]--;
                if(sn.nodes[j].nodetype!=NT::Sink) {
                    if(sn.nodes[j].nodetype==NT::Source||sn.nodes[j].station==0)
                        throw UnsupportedError(ErrorKind{"Topology"}, "MLMC routes directly to queues, delays or sinks");
                    change[(sn.nodes[j].station-1)*R+s]++;
                }
                for(std::size_t a=0;a<M*R;++a) p.nu[a].push_back(change[a]);
                origins.push_back(i*R+r); weights.push_back(w);
            }
            if(std::abs(total-1)>1e-10) throw InputError("Active station/class routing must sum to one");
        }
    }
    std::vector<qn::NetworkStruct<double>::Reward> rewards;
    for(const auto& r:sn.reward) if(rewardName.empty()||r.name==rewardName) {
        out.names.push_back(r.name); rewards.push_back(r);
    }
    if(!rewardName.empty()&&rewards.empty()) throw InputError("Unknown reward: "+rewardName);
    p.rates=[=](const std::vector<double>& x) {
        std::vector<double> d(M*R,0);
        for(std::size_t i=0;i<M;++i) {
            double n=0; for(std::size_t r=0;r<R;++r) n+=std::max(x[i*R+r],0.0);
            for(std::size_t r=0;r<R;++r) {
                if(source[i]) d[i*R+r]=mu[i][r];
                else if(sched[i]==SS::INF) d[i*R+r]=mu[i][r]*std::max(x[i*R+r],0.0);
                else if(n>0) d[i*R+r]=mu[i][r]*std::max(x[i*R+r],0.0)*std::min(n,servers[i])/n;
            }
        }
        std::vector<double> a(origins.size());
        for(std::size_t k=0;k<a.size();++k) a[k]=d[origins[k]]*weights[k];
        return a;
    };
    p.observe=[rewards](const std::vector<double>& x) {
        std::vector<double> f; for(const auto& r:rewards) f.push_back(r.fn(x)); return f;
    };
    return out;
}
 
/** SSA's finite-horizon reward getter. The usual average analyzer refuses MLMC. */
inline std::vector<MlmcResult> solver_ssa_get_tran_reward(const qn::NetworkStruct<double>& sn,
        const SsaOptions& opt,const std::vector<double>& times,const std::string& name="",
        const MlmcOptions& control=MlmcOptions()) {
    if(opt.method!="mlmc") throw InputError("getTranReward requires SSA method mlmc");
    if(times.empty()) throw InputError("Observation times must be nonempty");
    auto compiled=ssa_mlmc_problem(sn,name);
    std::vector<MlmcResult> out;
    for(double t:times) {
        auto runControl=control; runControl.seed=opt.seed;
        auto result=ssa_mlmc(compiled.problem,t,runControl);
        result.names=compiled.names; result.method=opt.method; out.push_back(result);
    }
    return out;
}
}} // line::ssa
#endif