授業の課題で、アーラン分布における待ち行列のシミュレーションを書くプログラム。途中で力尽きた。

erlang_distribution_simulation.cpp

#include <algorithm>//copy
#include <cassert>//assert
#include <iostream>//cout
#include <iterator>//iterator
#include <limits>//numeric
#include <memory>//shared_ptr
#include <numeric>//accumelate
#include <random>//mt19937, distribution
#include <sstream>//stringstream;
#include <type_traits>//is_arithmetic
#include <vector>//vector

//実体はなくてよい
extern void* enabler;

template <
    typename T = double,//乱数の型
    typename std::enable_if<std::is_arithmetic<T>::value>::type*& = enabler//Tが算術型かどうか判断する
>
class erlang_distribution
{
    //アーラン分布
    //このクラスのインターフェイスは標準の分布クラスに従っている

public:
    typedef T result_type;

    struct param_type
    {
        typedef erlang_distribution<T> distribution_type;

        param_type(T lambda, std::size_t dim = 5U)
            : lambda(lambda)
        {

            //kは1以上の正整数値でなければならない
            assert(1U <= dim && "K must be greater than or equal to 1");
            this->dim = dim;
        }

        bool operator==(const param_type& right) const
        {
            return dim == right.dim && lambda == right.lambda;
        }

        bool operator!=(const param_type& right) const
        {
            return !(*this == right);
        }

        std::size_t dim;
        T lambda;
    };

    explicit erlang_distribution(T lambda = static_cast<T>(1), std::size_t dim = 9U)
        : par(lambda, dim)
    {
        //アーラン分布のパラメータkとλを設定する
        //課題で指定されたデフォルト引数9を入れておく
    }

    explicit erlang_distribution(param_type par)
        : par(par)
    {
    }

    std::size_t dimension() const
    {
        return par.dim;
    }

    T lambda() const
    {
        return par.lambda;
    }

    param_type param() const
    {
        return par;
    }

    void param(const param_type& par)
    {
        this->par = par;
    }

    result_type(min)() const
    {
        return static_cast<result_type>(0);
    }

    result_type(max)() const
    {
        return (std::numeric_limits<result_type>::max)();
    }

    template <typename Engine>
    result_type operator()(Engine& engine) const
    {
        return execute(engine, par);
    }

    template <typename Engine>
    result_type operator()(Engine& engine, const param_type& param) const
    {
        return execute(engine, param);
    }

private:
    template <typename Engine>
    result_type execute(Engine& engine, const param_type& param) const
    {
        typedef std::exponential_distribution<result_type> exp_type;

        random<exp_type> rand(param.lambda);
        auto it = rand.begin();
        auto end = rand.end(param.dim);

        return std::accumulate(it, end, static_cast<result_type>(0));
    }

    param_type par;
};

template <typename Result>
class random_iterator
    : public std::iterator<std::input_iterator_tag, Result>
{
    //乱数をイテレータで扱えるようにする

public:
    typedef Result result_type;
    typedef std::iterator<std::input_iterator_tag, result_type> base_iterator;
    typedef typename base_iterator::value_type value_type;
    typedef typename base_iterator::iterator_category iterator_category;
    typedef typename base_iterator::pointer pointer;
    typedef typename base_iterator::reference reference;
    typedef typename base_iterator::difference_type difference_type;

    //rand_holder_baseとrand_holder<Rand>は、
    //型消去(Type Erasure)というテクニックを使用
    class rand_holder_base
    {
    public:
        rand_holder_base()
        {
        }

        virtual ~rand_holder_base()
        {
        }

        virtual result_type operator()() = 0;
    };

    template <typename Rand>
    class rand_holder : public rand_holder_base
    {
    public:
        rand_holder(const Rand& rand)
            : rand(rand)
        {
        }

        virtual result_type operator()() override
        {
            return rand();
        }

    private:
        Rand rand;
    };

    template <typename Random>
    random_iterator(const Random& rnd, std::size_t max = 0U)
        : rand(std::make_shared < rand_holder < Random >> (rnd)),
        max(max),
        counter(0U)
    {
        value = (*rand)();
    }

    bool operator!=(const random_iterator& right) const
    {
        //2つイテレータが不等かどうか返す
        return counter < right.max;
    }

    bool operator==(const random_iterator& right) const
    {
        //2つのイテレータが等しいかどうか返す
        return !(*this != right);
    }

    value_type operator*() const
    {
        //イテレータの指す要素を返す
        return value;
    }

    value_type operator++()
    {
        //イテレータを１つ進める
        value = (*rand)();
        ++counter;
        return value;
    }

private:
    std::shared_ptr<rand_holder_base> rand;
    value_type value;
    std::size_t max;
    std::size_t counter;
};

template <
    typename Dist,//乱数の分布
    typename Engine = std::mt19937//乱数生成アルゴリズム(デフォルトはメルセンヌ・ツイスター)
>
class random
{
    //指定した乱数生成アルゴリズムと分布と型に基づいた乱数を生成するファンクタ

public:
    typedef typename Dist::result_type result_type;
    typedef typename random_iterator<result_type> iterator;

    template <typename... Args>
    random(Args... args)
        : dist(std::forward<Args>(args)...),
        engine(std::random_device()())
    {
    }

    iterator begin() const
    {
        //開始イテレータを返す
        return iterator(*this);
    }

    iterator end(std::size_t max = 10U) const
    {
        //末尾イテレータを返す
        return iterator(*this, max);
    }

    template <typename... Args>
    result_type operator()(Args... args)
    {
        if (0U < sizeof...(Args))
        {
            //パラメータが設定されている時は更新する
            typename Dist::param_type p(std::forward<Args>(args)...);
            dist.param(p);
        }
        return dist(engine);
    }

private:
    Dist dist;//乱数分布
    Engine engine;//乱数生成アルゴリズム
};

struct fcfs
{
    //サービス規律の先着順ファンクタ

    typedef std::size_t size_type;

    fcfs() : s(0U)
    {
        //コンストラクタ
    }

    void push()
    {
        //キューに１つ追加
        ++s;
    }

    void pop()
    {
        //キューから１つ削除
        assert(!empty() && "!empty()");
        --s;
    }

    size_type size() const
    {
        //現在の長さを返す
        return s;
    }

    bool empty() const
    {
        //キューが空かどうか返す
        return s == 0U;
    }

private:
    size_type s;//待ち行列の客数
};

template <typename T>
struct typeof
{
    typedef T type;
};

#define TYPEOF(value) typeof<decltype(value)>::type

template <
    typename A,  //到着間隔分布
    typename B,	//サービス時間分布
    std::size_t	X,//窓口数
    std::size_t Y = UINT_MAX,//系内客数上限(デフォルトは∞、std::size_tの上限)
    typename Z = fcfs,//サービス規律(デフォルトは先着順)
    typename std::enable_if<1 <= X>::type*& = enabler//窓口数は１以上でなければならない
>
class simulation
{
    //待ち行列のシミュレーションをするファンクタ

public:
    typedef A arrival_dist_type;
    typedef typename arrival_dist_type::result_type arrival_type;
    typedef B service_dist_type;
    typedef typename service_dist_type::result_type service_type;
    typedef Z rule_type;

    void operator()(std::ostream& os, arrival_type lambda)
    {
        arrival_dist_type arrival_rand;
        arrival_type arrival_time = static_cast<arrival_type>(0);
        service_dist_type serivce_rand;
        service_type service_time = static_cast<service_type>(0);
        rule_type queue;
        decltype(arrival_time + service_time) total_time = static_cast<decltype(arrival_time + service_time)>(0);
        typedef std::vector<decltype(total_time)> customer_type;
        customer_type number_of_customers(DEFAULT_VECTOR_RESERVE_SIZE);
        const double MU = 9.0;

        for (int i = 0; i < MAX_SIMULATION_TIME; i++)
        {
            //到着までの時間(無記憶性により再計算して良い)
            arrival_time = arrival_rand(lambda);

            if (X == 1U)
            {
                //退去までの時間(窓口数1)
                service_time = serivce_rand(MU);
            }
            else
            {
                //退去までの時間(窓口数2以上)
                if (queue.size() <= 1)
                {
                    service_time = serivce_rand(MU);
                }
                else
                {
                    service_time = serivce_rand(2.0 * MU);
                }
            }

            if (queue.empty() || arrival_time < service_time)
            {
                //系内客数が0か、到着までの時間が短い場合は次に到着が発生

                total_time += arrival_time;

                if (queue.size() >= DEFAULT_VECTOR_RESERVE_SIZE)
                {
                    number_of_customers.resize(queue.size() + DEFAULT_VECTOR_RESERVE_SIZE);
                    number_of_customers[queue.size()] += arrival_time;
                }
                else
                {
                    number_of_customers[queue.size()] += arrival_time;
                }

                //系内客数は1増える
                queue.push();
            }
            else
            {
                //系内客数が1以上で、退去までの時間が短い場合は次に退去が発生

                total_time += service_time;

                if (queue.size() >= DEFAULT_VECTOR_RESERVE_SIZE)
                {
                    number_of_customers.resize(queue.size() + DEFAULT_VECTOR_RESERVE_SIZE);
                    number_of_customers[queue.size()] += arrival_time;
                }
                else
                {
                    number_of_customers[queue.size()] += arrival_time;
                }

                //系内客数は1減る
                queue.pop();
            }
        }
        for (auto& it : number_of_customers)
        {
            it /= total_time;
        }

        typename customer_type::value_type average = static_cast<typename customer_type::value_type>(0);

        for (std::size_t i = 0U; i < number_of_customers.size(); i++)
        {
            average += i * number_of_customers[i];
        }
        os << average << std::endl;
    }

private:
    static const std::size_t MAX_SIMULATION_TIME = 1000000U;//1度のシュミレーションの試行回数
    static const std::size_t DEFAULT_VECTOR_RESERVE_SIZE = 1 << 8U;
};

//指数分布
typedef random<std::exponential_distribution<double>> M;
//アーラン分布
typedef random<erlang_distribution<double>> E;

int main(int argc, char** argv)
{
    if (argc != 2)
    {
        std::cout << "argc != 2" << std::endl;
        return 1;
    }

    double lambda;
    std::stringstream ss;

    ss << argv[1];
    ss >> lambda;

    simulation<M, E, 1> s;
    s(std::cout, lambda);
}