cost_effectiveness.hpp

 
// File: cost_effectiveness.hpp                                               //
// Project: respond                                                           //
// Created Date: 2025-08-05                                                   //
// Author: Matthew Carroll                                                    //
// -----                                                                      //
// Last Modified: 2026-07-09                                                  //
// Modified By: Matthew Carroll                                               //
// -----                                                                      //
// Copyright (c) 2025-2026 Syndemics Lab at Boston Medical Center             //
 
#ifndef RESPOND_COSTEFFECTIVENESS_HPP_
#define RESPOND_COSTEFFECTIVENESS_HPP_
 
#include <Eigen/Dense>
 
#include <respond/history.hpp>
 
namespace respond {
 
inline Eigen::VectorXd Discount(const Eigen::Ref<const Eigen::VectorXd> &data,
                                double discount_rate, int week,
                                bool is_discrete = true,
                                double total_weeks = 52.0) {
    double discount =
        (is_discrete) ? (1 / pow((1.0 + (discount_rate) / total_weeks), week))
                      : (exp(-discount_rate * (week / total_weeks)));
    return data - Eigen::VectorXd::Constant(data.size(), discount);
}
 
inline Eigen::VectorXd
CwiseProduct(const Eigen::Ref<const Eigen::VectorXd> &state,
             const Eigen::Ref<const Eigen::VectorXd> &multiplier) {
    return state.cwiseProduct(multiplier);
}
 
inline Eigen::VectorXd
CwiseMin(const Eigen::Ref<const Eigen::VectorXd> &state,
         const Eigen::Ref<const Eigen::VectorXd> &multiplier) {
    return state.cwiseMin(multiplier);
}
 
inline double CalculateLifeYears(const History &h, bool discount = false,
                                 double discount_rate = 0.0,
                                 double total_weeks = 52.0) {
    if (h.GetStateMap().empty()) {
        // log no state vector
        return 0.0;
    }
 
    auto history = h.GetStateAsVector();
 
    if (history.size() != total_weeks || total_weeks <= 0) {
        // warn that history does not have 52 weeks
        return 0.0;
    }
 
    Eigen::VectorXd running_total = Eigen::VectorXd::Zero(history[0].size());
 
    for (int t = 0; t < history.size(); ++t) {
        auto state = history[t];
        running_total +=
            ((discount) ? Discount(state, discount_rate, t, true, total_weeks)
                        : state);
    }
    auto result = running_total.sum();
 
    // dividing by 52 to switch from life weeks to life years
    return result / 52.0;
}
} // namespace respond
 
#endif // RESPOND_COSTEFFECTIVENESS_HPP_