solve.hpp

Go to the documentation of this file.

/// @file solve/solve.hpp
/// @brief Unified solve() dispatcher -- the single entry point for all solvers.
#pragma once
 
#include "ode/implicit.hpp"
#include "ode/ode.hpp"
#include "solve/algorithms.hpp"
#include "solve/problems.hpp"
#include "stochastic/mcmc.hpp"
#include <concepts>
#include <random>
 
namespace num {
 
template<typename P>
concept IsODEProblem = requires(const P& p, double t, const Vector& y, Vector& dy) {
    p.f(t, y, dy);
    { p.u0 } -> std::convertible_to<const Vector&>;
    { p.t0 } -> std::convertible_to<double>;
    { p.tf } -> std::convertible_to<double>;
};
 
template<typename A>
concept IsExplicitODEAlg = std::same_as<std::remove_cvref_t<A>, Euler>
                           || std::same_as<std::remove_cvref_t<A>, RK4>
                           || std::same_as<std::remove_cvref_t<A>, RK45>;
 
template<typename A>
concept IsImplicitODEAlg = requires(const A& a) {
    { a.solver } -> std::convertible_to<LinearSolver>;
    { a.dt } -> std::convertible_to<double>;
    { a.nstep } -> std::convertible_to<int>;
};
 
template<typename A>
concept IsMCMCAlg = std::same_as<std::remove_cvref_t<A>, Metropolis>;
 
template<IsODEProblem P>
ODEResult solve(const P& prob, const RK45& alg, ObserverFn obs = nullptr) {
    ODEParams p{.t0 = prob.t0,
                .tf = prob.tf,
                .h = alg.h,
                .rtol = alg.rtol,
                .atol = alg.atol,
                .max_steps = alg.max_steps};
    return ode_rk45(prob.f, prob.u0, p, obs);
}
 
template<IsODEProblem P>
ODEResult solve(const P& prob, const RK4& alg, ObserverFn obs = nullptr) {
    return ode_rk4(prob.f, prob.u0, {.t0 = prob.t0, .tf = prob.tf, .h = alg.h}, obs);
}
 
template<IsODEProblem P>
ODEResult solve(const P& prob, const Euler& alg, ObserverFn obs = nullptr) {
    return ode_euler(prob.f, prob.u0, {.t0 = prob.t0, .tf = prob.tf, .h = alg.h}, obs);
}
 
template<VecField F>
void solve(F& u, const BackwardEuler& alg) {
    ode::advance(u, alg.solver, {alg.nstep, alg.dt});
}
 
template<VecField F, typename Observer>
void solve(F& u, const BackwardEuler& alg, Observer&& obs) {
    ode::advance(u, alg.solver, {alg.nstep, alg.dt}, std::forward<Observer>(obs));
}
 
/// @brief Run equilibration and measurement sweeps, then return the mean.
template<IsMCMCAlg A, typename MeasureFn, typename RNG>
double solve(const MCMCProblem& prob, const A& alg, MeasureFn&& measure, RNG& rng) {
    for (int s = 0; s < alg.equilibration; ++s)
        markov::metropolis_sweep_prob(prob.n_sites, prob.accept_prob, prob.propose, rng);
    double acc = 0.0;
    for (int s = 0; s < alg.measurements; ++s) {
        markov::metropolis_sweep_prob(prob.n_sites, prob.accept_prob, prob.propose, rng);
        acc += measure();
    }
    return acc / alg.measurements;
}
 
} // namespace num