populations.h

Go to the documentation of this file.

/*
* Copyright (C) 2020-2026 MEmilio
*
* Authors: Jan Kleinert, 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_EPI_POPULATIONS_H
#define MIO_EPI_POPULATIONS_H
 
#include "memilio/config.h" // IWYU pragma: keep
#include "memilio/math/eigen.h" // IWYU pragma: keep
#include "memilio/math/math_utils.h"
#include "memilio/utils/custom_index_array.h"
#include "memilio/utils/uncertain_value.h"
 
#include <concepts>
#include <numeric>
 
namespace mio
{
 
template <typename FP, class... Categories>
class Populations : public CustomIndexArray<UncertainValue<FP>, Categories...>
{
public:
    using Base  = CustomIndexArray<UncertainValue<FP>, Categories...>;
    using Index = typename Base::Index;
 
    template <class... Ts>
        requires std::is_constructible_v<UncertainValue<FP>, Ts...>
    explicit Populations(Index const& sizes, Ts... args)
        : Base(sizes, args...)
    {
    }
 
    explicit Populations(Base&& array)
        : Base(std::move(array))
    {
    }
 
    template <class OtherType>
        requires std::convertible_to<typename Base::Type::Type, OtherType>
    Populations<OtherType, Categories...> convert() const
    {
        return Populations<OtherType, Categories...>(
            Base::template convert<typename Base::Type::Type>().template convert<UncertainValue<OtherType>>());
    }
 
    size_t get_num_compartments() const
    {
        return this->numel();
    }
 
    inline Eigen::VectorX<FP> get_compartments() const
    {
        return this->array().template cast<FP>();
    }
 
    template <class Arr>
    decltype(auto) get_from(Arr&& y, Index const& cats) const
    {
        static_assert(std::is_lvalue_reference_v<Arr>, "get_from is disabled for temporary arrays.");
        return y[this->get_flat_index(cats)];
    }
 
    template <class T>
    FP get_group_total(mio::Index<T> group_idx) const
    {
        auto const s = this->template slice<T>({(size_t)group_idx, 1});
        return std::accumulate(s.begin(), s.end(), FP(0.0), [](const FP& a, const UncertainValue<FP>& b) {
            return evaluate_intermediate<FP>(a + b);
        });
    }
 
    template <class T>
    void set_group_total(mio::Index<T> group_idx, FP value)
    {
        using std::fabs;
        FP current_population = get_group_total(group_idx);
        auto s                = this->template slice<T>({(size_t)group_idx, 1});
 
        if (fabs(current_population) < 1e-12) {
            for (auto& v : s) {
                v = value / s.numel();
            }
        }
        else {
            for (auto& v : s) {
                v *= value / current_population;
            }
        }
    }
 
    FP get_total() const
    {
        return this->array().template cast<FP>().sum();
    }
 
    template <class T>
    void set_difference_from_group_total(Index const& midx, FP total_group_population)
 
    {
        auto group_idx        = mio::get<T>(midx);
        FP current_population = get_group_total(group_idx);
        size_t idx            = this->get_flat_index(midx);
        current_population -= this->array()[idx];
 
        assert(current_population <= total_group_population + 1e-10);
 
        this->array()[idx] = total_group_population - current_population;
    }
 
    void set_total(FP value)
    {
        using std::fabs;
        FP current_population = get_total();
        if (fabs(current_population) < 1e-12) {
            FP ysize = static_cast<FP>(this->array().size());
            for (size_t i = 0; i < this->get_num_compartments(); i++) {
                this->array()[(Eigen::Index)i] = value / ysize;
            }
        }
        else {
            for (size_t i = 0; i < this->get_num_compartments(); i++) {
                this->array()[(Eigen::Index)i] *= value / current_population;
            }
        }
    }
 
    void set_difference_from_total(Index midx, FP total_population)
    {
        FP current_population = get_total();
        size_t idx            = this->get_flat_index(midx);
        current_population -= this->array()[idx];
 
        assert(current_population <= total_population);
 
        this->array()[idx] = total_population - current_population;
    }
 
    bool apply_constraints()
    {
        bool corrected = false;
        for (int i = 0; i < this->array().size(); i++) {
            if (this->array()[i] < 0.0) {
                if (this->array()[i] > -1e-10) {
                    log_warning("Constraint check: Compartment number {} changed from {} to {}", i, this->array()[i],
                                0);
                }
                else {
                    log_error("Constraint check: Compartment number {} changed from {} to {}", i, this->array()[i], 0);
                }
                this->array()[i] = 0.0;
                corrected        = true;
            }
        }
        return corrected;
    }
 
    bool check_constraints() const
    {
        for (int i = 0; i < this->array().size(); i++) {
            FP value = this->array()[i];
            if (value < 0.0) {
                log_error("Constraint check: Compartment number {} is {} and smaller {}", i, value, 0);
                return true;
            }
        }
        return false;
    }
 
    template <class IOContext>
    static IOResult<Populations> deserialize(IOContext& io)
    {
        return Base::deserialize(io, Tag<Populations>{});
    }
};
 
template <typename FP, class... Categories>
class Populations<FP, Index<Categories...>> : public Populations<FP, Categories...>
{
    using Populations<FP, Categories...>::Populations;
};
 
} // namespace mio
 
#endif // MIO_EPI_POPULATIONS_H