compile time repressilator simulation depending on sprout

constexpr-repressilator.cpp

#include <sprout/config.hpp>
#include <sprout/container/traits.hpp>
#include <sprout/container/functions.hpp>
#include <sprout/container/indexes.hpp>
#include <sprout/algorithm/fixed/unfold.hpp>
#include <sprout/utility/pair.hpp>
#include <sprout/array.hpp>

constexpr static double dt = 0.1;

template<typename realT>
struct transcript
{
    constexpr static realT a  = 200.0;
    constexpr static realT a0 = 0.2;

    constexpr static realT
    invoke(const realT rna, const realT pro)
    {
        return a0 - rna + a / (1 + pro * pro);
    }
};

template<typename realT>
struct translate
{
    constexpr static realT b = 3.0;

    constexpr static realT
    invoke(const realT pro, const realT rna)
    {
        return b * (rna - pro);
    }
};

template<typename realT, typename reactT>
inline constexpr realT
euler(realT self, realT reactor)
{
    return self + reactT::invoke(self, reactor) * dt;
}

template<typename stateT, typename indexT>
inline constexpr typename sprout::container_traits<stateT>::value_type
next_state(const stateT st, indexT i)
{
    return (i == 0) ?
        euler<typename sprout::container_traits<stateT>::value_type,
              translate<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[3]) :
           (i == 1) ?
        euler<typename sprout::container_traits<stateT>::value_type,
              translate<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[4]) :
           (i == 2) ?
        euler<typename sprout::container_traits<stateT>::value_type,
              translate<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[5]) :
           (i == 3) ?
        euler<typename sprout::container_traits<stateT>::value_type,
              transcript<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[1]) :
           (i == 4) ?
        euler<typename sprout::container_traits<stateT>::value_type,
              transcript<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[2]) :
          // i == 5
        euler<typename sprout::container_traits<stateT>::value_type,
              transcript<typename sprout::container_traits<stateT>::value_type>
          >(st[i], st[0]);
}

template<typename stateT, sprout::index_t ... Is>
inline constexpr stateT
next_states(const stateT& st, sprout::index_tuple<Is...>)
{
    return sprout::remake<stateT>(st, sprout::size(st), next_state(st, Is)...);
}

template<typename stateT>
inline constexpr stateT next_states(const stateT& st)
{
    return next_states(st, sprout::container_indexes<stateT>::make());
}

struct gen_next_state
{
    template<typename stateT>
    constexpr sprout::pair<stateT, stateT>
    operator()(const stateT& st) const
    {
        return sprout::pair<stateT, stateT>(st, next_states(st));
    }
};

template<typename stateT, std::size_t N>
inline constexpr sprout::array<stateT, N>
time_evolution(const stateT& init)
{
    return sprout::fixed::unfold<sprout::array<stateT, N>>(gen_next_state(), init);
}

int main()
{
    typedef double real_type;
    typedef sprout::array<real_type, 6> state_type;
    constexpr std::size_t total_step = 1000;

    constexpr state_type init{{10., 10., 10., 100., 80., 50.}};
    constexpr auto states = time_evolution<state_type, total_step>(init);

    for(auto&& st : states)
    {
        for(auto&& conc : st)
            std::cout << conc << " ";
        std::cout << std::endl;
    }
    return 0;
}