state_age_function.h

Go to the documentation of this file.

/* 
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Anna Wendler, Lena Ploetzke
*
* 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 STATEAGEFUNCTION_H
#define STATEAGEFUNCTION_H
 
#include "memilio/config.h"
#include "memilio/utils/compiler_diagnostics.h"
#include "memilio/math/smoother.h"
#include "memilio/math/floating_point.h"
 
#include "boost/math/distributions/gamma.hpp"
#include "boost/math/distributions/lognormal.hpp"
#include "memilio/utils/logging.h"
 
namespace mio
{
/**************************
* Define StateAgeFunction *
***************************/
 
template <typename FP>
struct StateAgeFunction {
 
    StateAgeFunction(FP init_distribution_parameter, FP init_location = 0, FP init_scale = 1)
        : m_distribution_parameter{init_distribution_parameter}
        , m_location{init_location}
        , m_scale{init_scale}
        , m_mean{-1.0} // Initialize mean as not set.
        , m_support_max{-1.0} // Initialize support maximum as not set.
    {
        if (m_scale <= 0) {
            log_error("The scale parameter of a StateAgeFunction has to be positive. Set scale to 1.");
            m_scale = 1.0;
        }
    }
 
    virtual ~StateAgeFunction() = default;
 
    StateAgeFunction(const StateAgeFunction<FP>& other) = default;
 
    StateAgeFunction(StateAgeFunction<FP>&& other) = default;
 
    StateAgeFunction<FP>& operator=(const StateAgeFunction<FP>& other) = default;
 
    StateAgeFunction<FP>& operator=(StateAgeFunction<FP>&& other) = default;
 
    bool operator==(const StateAgeFunction<FP>& other) const
    {
        return (typeid(*this).name() == typeid(other).name() &&
                m_distribution_parameter == other.get_distribution_parameter() && m_location == other.get_location() &&
                m_scale == other.get_scale());
    }
 
    virtual FP eval(FP state_age) = 0;
 
    FP get_distribution_parameter() const
    {
        return m_distribution_parameter;
    }
 
    void set_distribution_parameter(FP new_distribution_parameter)
    {
        m_distribution_parameter = new_distribution_parameter;
        m_support_max            = -1.;
        m_mean                   = -1.;
    }
 
    FP get_location() const
    {
        return m_location;
    }
 
    void set_location(FP new_location)
    {
        m_location    = new_location;
        m_support_max = -1.;
        m_mean        = -1.;
    }
 
    FP get_scale() const
    {
        return m_scale;
    }
 
    void set_scale(FP new_scale)
    {
        if (new_scale <= 0) {
            log_error("The scale parameter of a StateAgeFunction has to be positive. Set scale to 1.");
            new_scale = 1;
        }
        m_scale       = new_scale;
        m_support_max = -1.;
        m_mean        = -1.;
    }
 
    virtual FP get_support_max(FP dt, FP tol = 1e-10)
    {
        FP support_max = 0.0;
 
        if (!floating_point_equal<FP>(m_support_tol, tol, 1e-14) ||
            floating_point_equal<FP>(m_support_max, -1.0, 1e-14)) {
            while (eval(support_max) >= tol) {
                support_max += dt;
            }
 
            m_support_max = support_max;
            m_support_tol = tol;
        }
 
        return m_support_max;
    }
 
    virtual FP get_mean(FP dt = 1.0, FP tol = 1e-10)
    {
        using std::ceil;
        if (!floating_point_equal<FP>(m_mean_tol, tol, 1e-14) || floating_point_equal<FP>(m_mean, -1., 1e-14)) {
            // Integration using trapezoidal rule.
            // Compute value for i=0.
            FP mean         = 0.5 * dt * eval(FP(0 * dt));
            FP supp_max_idx = ceil(get_support_max(dt, tol) / dt);
            // We start with i=1 since the value for i=0 was already considered above when defining the variable mean.
            // Note that we do not consider indices i>=supp_max_idx since it holds eval(i*dt)<tol for all i>=supp_max_idx
            // by definition of the support_max.
            for (int i = 1; i < supp_max_idx; i++) {
                mean += dt * eval(FP(i * dt));
            }
 
            m_mean     = mean;
            m_mean_tol = tol;
        }
        return m_mean;
    }
 
    std::string get_state_age_function_type() const
    {
        return typeid(*this).name();
    }
 
    std::unique_ptr<StateAgeFunction<FP>> clone() const
    {
        return std::unique_ptr<StateAgeFunction<FP>>(clone_impl());
    }
 
protected:
    virtual StateAgeFunction<FP>* clone_impl() const = 0;
 
    FP m_distribution_parameter; 
    FP m_location; 
    FP m_scale; 
    FP m_mean; 
    FP m_mean_tol{-1.0}; 
    FP m_support_max; 
    FP m_support_tol{-1.0}; 
};
 
/**************************************
* Derived classes of StateAgeFunction *
***************************************/
 
template <typename FP>
struct ExponentialSurvivalFunction : public StateAgeFunction<FP> {
 
    ExponentialSurvivalFunction(FP init_distribution_parameter, FP init_location = 0)
        : StateAgeFunction<FP>(init_distribution_parameter, init_location)
    {
    }
 
    FP eval(FP state_age) override
    {
        using std::exp;
        if (state_age <= this->m_location) {
            return 1;
        }
        return exp(-this->m_distribution_parameter * (state_age - this->m_location));
    }
 
    FP get_mean(FP dt = 1.0, FP tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        return 1.0 / this->m_distribution_parameter + this->m_location;
    }
 
protected:
    StateAgeFunction<FP>* clone_impl() const override
    {
        return new ExponentialSurvivalFunction<FP>(*this);
    }
};
 
template <typename FP>
struct SmootherCosine : public StateAgeFunction<FP> {
 
    SmootherCosine(FP init_distribution_parameter, FP init_location = 0)
        : StateAgeFunction<FP>(init_distribution_parameter, init_location)
    {
    }
 
    FP eval(FP state_age) override
    {
        if (state_age <= this->m_location) {
            return 1.0;
        }
        return smoother_cosine<FP>(state_age - this->m_location, 0.0, this->m_distribution_parameter, 1.0, 0.0);
    }
 
    FP get_support_max(FP dt, FP tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        this->m_support_max = this->m_distribution_parameter + this->m_location;
        return this->m_support_max;
    }
 
    FP get_mean(FP dt = 1.0, FP tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        return 0.5 * this->m_distribution_parameter + this->m_location;
    }
 
protected:
    StateAgeFunction<FP>* clone_impl() const override
    {
        return new SmootherCosine<FP>(*this);
    }
};
 
struct GammaSurvivalFunction : public StateAgeFunction<ScalarType> {
 
    GammaSurvivalFunction(ScalarType init_shape = 1, ScalarType init_location = 0, ScalarType init_scale = 1)
        : StateAgeFunction<ScalarType>(init_shape, init_location, init_scale)
    {
    }
 
    ScalarType eval(ScalarType state_age) override
    {
        if (state_age <= this->m_location) {
            return 1;
        }
        boost::math::gamma_distribution<ScalarType, boost::math::policies::policy<>> gamma(
            this->m_distribution_parameter, this->m_scale);
        return boost::math::cdf(boost::math::complement(gamma, state_age - this->m_location));
    }
 
    ScalarType get_mean(ScalarType dt = 1.0, ScalarType tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        return this->m_distribution_parameter * this->m_scale + this->m_location;
    }
 
protected:
    StateAgeFunction<ScalarType>* clone_impl() const override
    {
        return new GammaSurvivalFunction(*this);
    }
};
 
struct LognormSurvivalFunction : public StateAgeFunction<ScalarType> {
 
    LognormSurvivalFunction(ScalarType init_distribution_parameter, ScalarType init_location = 0,
                            ScalarType init_scale = 1)
        : StateAgeFunction<ScalarType>(init_distribution_parameter, init_location, init_scale)
    {
    }
 
    ScalarType eval(ScalarType state_age) override
    {
        if (state_age < this->m_location) {
            return 1;
        }
        boost::math::lognormal_distribution<ScalarType, boost::math::policies::policy<>> logn(
            0.0, this->m_distribution_parameter);
        return boost::math::cdf(boost::math::complement(logn, (state_age - this->m_location) / this->m_scale));
    }
 
    // Closed form for the mean value is kind of complex, use default implementation.
    // For testing purposes, a class must exist anyway that uses the default implementation.
 
protected:
    StateAgeFunction<ScalarType>* clone_impl() const override
    {
        return new LognormSurvivalFunction(*this);
    }
};
 
template <typename FP>
struct ConstantFunction : public StateAgeFunction<FP> {
    ConstantFunction(FP init_distribution_parameter)
        : StateAgeFunction<FP>(init_distribution_parameter)
    {
    }
 
    FP eval(FP state_age) override
    {
        unused(state_age);
        return this->m_distribution_parameter;
    }
 
    FP get_support_max(FP dt, FP tol = 1e-10) override
    {
        // In case of a ConstantFunction we would have support_max = infinity
        // This type of function is not suited to be a TransitionDistribution
        // Log error and return -2.
 
        unused(dt);
        unused(tol);
        this->m_support_max = -2.0;
 
        log_error("This function is not suited to be a TransitionDistribution. Do not call in case of "
                  "StateAgeFunctions of type b); see documentation of StateAgeFunction Base class.");
 
        return this->m_support_max;
    }
 
    FP get_mean(FP dt = 1.0, FP tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        log_warning("Attention: This function is not suited to be a TransitionDistribution. Do not call in case of "
                    "StateAgeFunctions of type b); see documentation of StateAgeFunction Base class.");
        return this->m_distribution_parameter;
    }
 
protected:
    StateAgeFunction<FP>* clone_impl() const override
    {
        return new ConstantFunction<FP>(*this);
    }
};
 
struct ErlangDensity : public StateAgeFunction<ScalarType> {
 
    ErlangDensity(unsigned int init_shape, ScalarType init_scale)
        : StateAgeFunction<ScalarType>(init_shape, 0.0, init_scale)
    {
    }
 
    ScalarType eval(ScalarType state_age) override
    {
        using std::exp;
        using std::pow;
        if (state_age < 0) {
            return 0;
        }
        int shape = static_cast<int>(this->m_distribution_parameter);
        return pow(state_age / this->m_scale, shape - 1) /
               (this->m_scale * boost::math::factorial<ScalarType>(shape - 1)) * exp(-state_age / this->m_scale);
    }
 
    ScalarType get_support_max(ScalarType dt, ScalarType tol = 1e-10) override
    {
        // We are looking for the smallest time value t where function(tau)=0 for all tau>t. Thus support max is bigger
        // than the mean.
        // We use, that the density is monotonically decreasing for tau>mean here.
        ScalarType mean        = this->m_distribution_parameter * this->m_scale;
        ScalarType support_max = dt * static_cast<int>(mean / dt);
 
        if (!floating_point_equal<ScalarType>(this->m_support_tol, tol, 1e-14) ||
            floating_point_equal<ScalarType>(this->m_support_max, -1.0, 1e-14)) {
            while (eval(support_max) >= tol) {
                support_max += dt;
            }
 
            this->m_support_max = support_max;
            this->m_support_tol = tol;
        }
 
        return this->m_support_max;
    }
 
    ScalarType get_mean(ScalarType dt = 1.0, ScalarType tol = 1e-10) override
    {
        unused(dt);
        unused(tol);
        log_warning("Attention: This function is not suited to be a TransitionDistribution. Do not call in case of "
                    "StateAgeFunctions of type b); see documentation of StateAgeFunction Base class.");
        return this->m_distribution_parameter * this->m_scale;
    }
 
protected:
    StateAgeFunction<ScalarType>* clone_impl() const override
    {
        return new ErlangDensity(*this);
    }
};
 
/*********************************
* Define StateAgeFunctionWrapper *
*********************************/
 
template <typename FP>
struct StateAgeFunctionWrapper {
 
    StateAgeFunctionWrapper()
        : m_function(mio::SmootherCosine<FP>(1.0).clone())
    {
    }
    StateAgeFunctionWrapper(StateAgeFunction<FP>& init_function)
        : m_function(init_function.clone())
    {
    }
 
    StateAgeFunctionWrapper(StateAgeFunctionWrapper<FP> const& other)
        : m_function(other.m_function->clone())
    {
    }
 
    StateAgeFunctionWrapper(StateAgeFunctionWrapper<FP>&& other) = default;
 
    StateAgeFunctionWrapper<FP>& operator=(StateAgeFunctionWrapper<FP> const& other)
    {
        m_function = other.m_function->clone();
        return *this;
    }
 
    StateAgeFunctionWrapper<FP>& operator=(StateAgeFunctionWrapper<FP>&& other) = default;
 
    ~StateAgeFunctionWrapper() = default;
 
    bool operator==(const StateAgeFunctionWrapper<FP>& other) const
    {
        return (m_function->get_state_age_function_type() == other.get_state_age_function_type() &&
                m_function->get_distribution_parameter() == other.get_distribution_parameter() &&
                m_function->get_location() == other.get_location() && m_function->get_scale() == other.get_scale());
    }
 
    void set_state_age_function(StateAgeFunction<FP>& new_function)
    {
        m_function = new_function.clone();
    }
 
    std::string get_state_age_function_type() const
    {
        return m_function->get_state_age_function_type();
    }
 
    FP eval(FP state_age) const
    {
        return m_function->eval(state_age);
    }
 
    FP get_distribution_parameter() const
    {
        return m_function->get_distribution_parameter();
    }
 
    void set_distribution_parameter(FP new_distribution_parameter)
    {
        m_function->set_distribution_parameter(new_distribution_parameter);
    }
 
    FP get_location() const
    {
        return m_function->get_location();
    }
 
    void set_location(FP new_location)
    {
        m_function->set_location(new_location);
    }
    FP get_scale() const
    {
        return m_function->get_scale();
    }
 
    void set_scale(FP new_scale)
    {
        m_function->set_scale(new_scale);
    }
 
    FP get_support_max(FP dt, FP tol = 1e-10) const
    {
        return m_function->get_support_max(dt, tol);
    }
 
    FP get_mean(FP dt = 1.0, FP tol = 1e-10) const
    {
        return m_function->get_mean(dt, tol);
    }
 
private:
    std::unique_ptr<StateAgeFunction<FP>> m_function; 
};
 
} // namespace mio
 
#endif // STATEAGEFUNCTION_H