euler_maruyama.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Martin J. Kuehn, Daniel Abele
*
* Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
*     http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#ifndef MIO_MATH_EULER_MARUYAMA_H
#define MIO_MATH_EULER_MARUYAMA_H
 
#include "memilio/math/eigen.h" // IWYU pragma: keep
#include "memilio/math/integrator.h"
#include "memilio/math/math_utils.h"
#include "memilio/utils/logging.h"
 
namespace mio
{
 
template <typename FP>
class EulerMaruyamaIntegratorCore : public SdeIntegratorCore<FP>
{
public:
    EulerMaruyamaIntegratorCore()
        : SdeIntegratorCore<FP>(FP{}, FP{})
    {
    }
 
    std::unique_ptr<SdeIntegratorCore<FP>> clone() const override
    {
        return std::make_unique<EulerMaruyamaIntegratorCore>(*this);
    }
 
    bool step(const DerivFunction<FP>& f, const DerivFunction<FP>& noise_f, Eigen::Ref<const Eigen::VectorX<FP>> yt,
              FP& t, FP& dt, Eigen::Ref<Eigen::VectorX<FP>> ytp1) const override
    {
        using std::sqrt;
 
        f(yt, t, ytp1);
 
        Eigen::VectorX<FP> noise = Eigen::VectorX<FP>::Zero(yt.size());
        noise_f(yt, t, noise);
 
        ytp1 = yt + dt * ytp1 + sqrt(dt) * noise;
 
        if (!map_to_nonnegative(ytp1)) [[unlikely]] {
            mio::log_error("Failed to rescale values in Euler Maruyama step, total population is negative. "
                           "Redistributing absolute values evenly.");
        }
 
        t += dt;
        return true;
    }
};
 
} // namespace mio
 
#endif // MIO_MATH_EULER_MARUYAMA_H