ssa_mlmc.h

Go to the documentation of this file.

#ifndef LINE_SOLVERS_SSA_MLMC_H
#define LINE_SOLVERS_SSA_MLMC_H
 
#include <algorithm>
#include <cmath>
#include <functional>
#include <numeric>
#include <limits>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
#include <boost/math/distributions/normal.hpp>
 
namespace line { namespace ssa {
 
/** Generic fixed-increment jump process. Rates/rewards include signed auxiliary states. */
struct MlmcProblem {
    std::vector<double> x0;
    std::vector<std::vector<double>> nu; // dimension by channels
    std::function<std::vector<double>(const std::vector<double>&)> rates, observe;
};
struct MlmcOptions {
    std::vector<double> absTol{0.1};
    double confidenceLevel=0.95;
    int levels=3, baseSteps=4, refinementFactor=2, pilotPaths=64, minPaths=64;
    int maxPaths=1000000, maxEvents=1000000;
    unsigned long seed=23000;
};
struct MlmcResult {
    std::vector<double> estimate, standardError, steps, costsPerSample;
    std::vector<std::vector<double>> confidenceInterval, levelMeans, levelVariances, pilotVariances;
    std::vector<int> pathsPerLevel;
    std::vector<bool> samplingToleranceMet;
    bool exactCorrectionIncluded=false;
    double time=0, confidenceLevel=0.95;
    std::string biasStatus, method;
    std::vector<std::string> names;
};
namespace mlmc_detail {
class Engine {
    using Vec=std::vector<double>;
    const MlmcProblem& p;
    MlmcOptions o;
    double time,z;
    std::size_t d,k,q,nl;
    std::vector<int> nsteps;
    Vec steps,tol,initial;
    std::mt19937_64 rng;
    double work=0;
    [[noreturn]] static void fail(const char* s) { throw std::invalid_argument(s); }
    Vec rates(const Vec& x) {
        Vec a=p.rates(x); double sum=0;
        if(a.size()!=k) fail("MLMC rate dimension changed");
        for(double v:a) { if(!std::isfinite(v)||v<0) fail("MLMC rates must be finite and nonnegative on auxiliary states"); sum+=v; }
        if(!std::isfinite(sum)) fail("MLMC total rate overflow");
        return a;
    }
    Vec observe(const Vec& x) {
        Vec f=p.observe(x);
        if(f.empty()) fail("Define at least one MLMC reward");
        for(double v:f) if(!std::isfinite(v)) fail("MLMC rewards must be finite on auxiliary states");
        return f;
    }
    double uniform() { return std::generate_canonical<double,53>(rng); }
    double poisson(double mean) {
        if(mean==0) return 0;
        if(!std::isfinite(mean)||mean>std::numeric_limits<int>::max()/4.0)
            fail("MLMC Poisson mean exceeds supported range; refine the grid");
        return static_cast<double>(std::poisson_distribution<int>(mean)(rng));
    }
    void jump(Vec& x,std::size_t channel,double count) {
        for(std::size_t i=0;i<d;++i) x[i]+=p.nu[i][channel]*count;
    }
    Vec difference(const Vec& x,const Vec& c) {
        Vec a=observe(x),b=observe(c);
        if(a.size()!=q||b.size()!=q) fail("MLMC reward dimension changed");
        for(std::size_t i=0;i<q;++i) a[i]-=b[i];
        return a;
    }
    Vec sample(int level) {
        Vec x=p.x0,c=p.x0; work=0;
        if(level==0) {
            for(int j=0;j<nsteps[0];++j) {
                Vec a=rates(x);
                for(std::size_t r=0;r<k;++r) jump(x,r,poisson(a[r]*steps[0]));
                work+=std::max<std::size_t>(k,1);
            }
            return observe(x);
        }
        if(level<=o.levels) {
            Vec ac; double h=steps[level];
            for(int j=0;j<nsteps[level];++j) {
                if(j%o.refinementFactor==0) ac=rates(c);
                Vec af=rates(x);
                for(std::size_t r=0;r<k;++r) {
                    double shared=std::min(af[r],ac[r]),common=poisson(shared*h);
                    jump(x,r,common+poisson((af[r]-shared)*h));
                    jump(c,r,common+poisson((ac[r]-shared)*h));
                }
                work+=3*std::max<std::size_t>(k,1);
            }
            return difference(x,c);
        }
        double t=0,h=steps.back(); int events=0;
        for(int j=1;j<=nsteps.back();++j) {
            Vec ac=rates(c); double boundary=j*h;
            while(t<boundary) {
                Vec ax=rates(x),a(3*k); double total=0;
                for(std::size_t r=0;r<k;++r) {
                    a[r]=std::min(ax[r],ac[r]); a[k+r]=ax[r]-a[r]; a[2*k+r]=ac[r]-a[r];
                    total+=ax[r]+ac[r]-a[r];
                }
                work+=3*std::max<std::size_t>(k,1);
                if(!std::isfinite(total)) fail("MLMC split rate overflow");
                if(total==0) { t=boundary; break; }
                double dt=-std::log(std::max(std::numeric_limits<double>::min(),uniform()))/total;
                if(t+dt>=boundary) { t=boundary; break; }
                if(t+dt==t) fail("MLMC event time cannot advance");
                double u=uniform()*total,cs=0; std::size_t channel=0;
                for(;channel+1<a.size();++channel) { cs+=a[channel]; if(u<cs) break; }
                if(channel<2*k) jump(x,channel%k,1);
                if(channel<k||channel>=2*k) jump(c,channel%k,1);
                for(double v:x) if(v<0) fail("MLMC exact process left the nonnegative domain");
                t+=dt;
                if(++events>o.maxEvents) fail("MLMC exact correction exceeds maxEvents");
            }
        }
        return difference(x,c);
    }
    struct Batch { Vec mean,var; double cost; };
    Batch batch(int level,int n) {
        Batch b{Vec(q),Vec(q),0};
        for(int i=1;i<=n;++i) {
            Vec y=sample(level); b.cost+=work;
            for(double v:y) if(!std::isfinite(v)) fail("MLMC reward difference overflow");
            if(y.size()!=q) fail("MLMC reward dimension changed");
            for(std::size_t r=0;r<q;++r) {
                double delta=y[r]-b.mean[r]; b.mean[r]+=delta/i; b.var[r]+=delta*(y[r]-b.mean[r]);
            }
        }
        for(double v:b.mean) if(!std::isfinite(v)) fail("MLMC moment overflow");
        for(double& v:b.var) {
            if(!std::isfinite(v)) fail("MLMC moment overflow");
            v=std::max(0.0,v/(n-1));
        }
        b.cost=std::max(1.0,b.cost/n); return b;
    }
public:
    Engine(const MlmcProblem& p_,double time_,const MlmcOptions& o_):p(p_),o(o_),time(time_),rng(o.seed) {
        if(!std::isfinite(time)||time<0||o.levels<0||o.levels>20||o.baseSteps<1||o.refinementFactor<2
           ||o.pilotPaths<2||o.minPaths<2||o.maxPaths<2||o.maxEvents<1||!(o.confidenceLevel>0&&o.confidenceLevel<1))
            fail("Invalid MLMC time, hierarchy or sampling options");
        d=p.x0.size();
        if(d==0||p.nu.size()!=d) fail("Invalid MLMC state dimension");
        k=p.nu[0].size();
        for(std::size_t i=0;i<d;++i) {
            if(!std::isfinite(p.x0[i])||p.x0[i]<0||p.x0[i]!=std::floor(p.x0[i])||p.nu[i].size()!=k)
                fail("MLMC needs nonnegative integer initial counts");
            for(double v:p.nu[i]) if(!std::isfinite(v)||v!=std::floor(v)) fail("MLMC needs integer increments");
        }
        initial=observe(p.x0); q=initial.size(); nl=o.levels+2;
        if(o.absTol.size()!=1&&o.absTol.size()!=q) fail("absTol must be scalar or one per reward");
        for(std::size_t r=0;r<q;++r) {
            double v=o.absTol[o.absTol.size()==1?0:r];
            if(!std::isfinite(v)||v<=0) fail("absTol must be positive");
            tol.push_back(v);
        }
        long long n=o.baseSteps;
        for(int l=0;l<=o.levels;++l) {
            if(n>o.maxEvents) fail("MLMC hierarchy exceeds maxEvents");
            nsteps.push_back(static_cast<int>(n)); steps.push_back(time/n);
            if(time>0 && steps.back()==0) fail("MLMC time steps underflow; use a coarser grid"); n*=o.refinementFactor;
        }
        z=-boost::math::quantile(boost::math::normal(),(1-o.confidenceLevel)/2);
    }
    MlmcResult run() {
        MlmcResult e; e.time=time; e.confidenceLevel=o.confidenceLevel; e.steps=steps;
        e.exactCorrectionIncluded=true;
        e.levelMeans.assign(q,Vec(nl)); e.levelVariances=e.levelMeans; e.pilotVariances=e.levelMeans;
        e.pathsPerLevel.assign(nl,0); e.costsPerSample.assign(nl,0);
        e.estimate.assign(q,0); e.standardError.assign(q,0); e.confidenceInterval.assign(q,Vec(2));
        e.samplingToleranceMet.assign(q,false);
        if(time>0) {
            if(static_cast<long long>(nl)*(static_cast<long long>(o.pilotPaths)+o.minPaths)>o.maxPaths)
                fail("MLMC minimum sampling allocation exceeds maxPaths");
            for(std::size_t l=0;l<nl;++l) { Batch b=batch(static_cast<int>(l),o.pilotPaths); e.costsPerSample[l]=b.cost;
                for(std::size_t r=0;r<q;++r) e.pilotVariances[r][l]=b.var[r]; }
            std::fill(e.pathsPerLevel.begin(),e.pathsPerLevel.end(),o.minPaths);
            for(std::size_t r=0;r<q;++r) {
                double sum=0;
                for(std::size_t l=0;l<nl;++l) sum+=std::sqrt(e.pilotVariances[r][l]*e.costsPerSample[l]);
                for(std::size_t l=0;l<nl;++l) {
                    double n=std::ceil(sum/std::pow(tol[r]/z,2)*std::sqrt(e.pilotVariances[r][l]/e.costsPerSample[l]));
                    if(!std::isfinite(n)||n>o.maxPaths) fail("MLMC pilot allocation exceeds maxPaths");
                    e.pathsPerLevel[l]=std::max(e.pathsPerLevel[l],static_cast<int>(n));
                }
            }
            long long total=nl*o.pilotPaths; for(int n:e.pathsPerLevel) total+=n;
            if(total>o.maxPaths) fail("MLMC pilot allocation exceeds maxPaths");
            for(std::size_t l=0;l<nl;++l) { Batch b=batch(static_cast<int>(l),e.pathsPerLevel[l]);
                for(std::size_t r=0;r<q;++r) {
                    e.levelMeans[r][l]=b.mean[r]; e.levelVariances[r][l]=b.var[r];
                    e.estimate[r]+=b.mean[r]; e.standardError[r]+=b.var[r]/e.pathsPerLevel[l];
                }
            }
        } else e.estimate=initial;
        for(std::size_t r=0;r<q;++r) {
            e.standardError[r]=std::sqrt(e.standardError[r]);
            e.confidenceInterval[r][0]=e.estimate[r]-z*e.standardError[r];
            e.confidenceInterval[r][1]=e.estimate[r]+z*e.standardError[r];
            e.samplingToleranceMet[r]=z*e.standardError[r]<=tol[r];
        }
        e.biasStatus=time==0?"zero":"zero-discretization-bias";
        return e;
    }
};
} // mlmc_detail
 
/** Independent pilots, fixed production counts, pointwise normal intervals. */
inline MlmcResult ssa_mlmc(const MlmcProblem& p,double time,const MlmcOptions& o=MlmcOptions()) {
    return mlmc_detail::Engine(p,time,o).run();
}
}} // line::ssa
#endif