Line data    Source code

       1             : /* 
       2             : * Copyright (C) 2020-2024 MEmilio
       3             : *
       4             : * Authors: Daniel Abele, Elisabeth Kluth, Khoa Nguyen, David Kerkmann, Rene Schmieding
       5             : *
       6             : * Contact: Martin J. Kuehn <Martin.Kuehn@DLR.de>
       7             : *
       8             : * Licensed under the Apache License, Version 2.0 (the "License");
       9             : * you may not use this file except in compliance with the License.
      10             : * You may obtain a copy of the License at
      11             : *
      12             : *     http://www.apache.org/licenses/LICENSE-2.0
      13             : *
      14             : * Unless required by applicable law or agreed to in writing, software
      15             : * distributed under the License is distributed on an "AS IS" BASIS,
      16             : * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
      17             : * See the License for the specific language governing permissions and
      18             : * limitations under the License.
      19             : */
      20             : 
      21             : #include "abm/model_functions.h"
      22             : #include "abm/location.h"
      23             : #include "abm/person.h"
      24             : #include "abm/random_events.h"
      25             : #include "abm/infection.h"
      26             : #include "abm/virus_variant.h"
      27             : #include "memilio/epidemiology/age_group.h"
      28             : #include "memilio/utils/logging.h"
      29             : 
      30             : namespace mio
      31             : {
      32             : namespace abm
      33             : {
      34             : 
      35        1980 : ScalarType daily_transmissions_by_contacts(const ContactExposureRates& rates, const CellIndex cell_index,
      36             :                                            const VirusVariant virus, const AgeGroup age_receiver,
      37             :                                            const LocalInfectionParameters& params)
      38             : {
      39        1980 :     assert(age_receiver < rates.size<AgeGroup>());
      40        1980 :     ScalarType prob = 0;
      41       13860 :     for (AgeGroup age_transmitter(0); age_transmitter < rates.size<AgeGroup>(); ++age_transmitter) {
      42       11880 :         prob +=
      43       23760 :             rates[{cell_index, virus, age_transmitter}] * params.get<ContactRates>()[{age_receiver, age_transmitter}];
      44             :     }
      45        1980 :     return prob;
      46             : }
      47             : 
      48        1980 : ScalarType daily_transmissions_by_air(const AirExposureRates& rates, const CellIndex cell_index,
      49             :                                       const VirusVariant virus, const Parameters& global_params)
      50             : {
      51        3960 :     return rates[{cell_index, virus}] * global_params.get<AerosolTransmissionRates>()[{virus}];
      52             : }
      53             : 
      54        2817 : void interact(PersonalRandomNumberGenerator& personal_rng, Person& person, const Location& location,
      55             :               const AirExposureRates& local_air_exposure, const ContactExposureRates& local_contact_exposure,
      56             :               const TimePoint t, const TimeSpan dt, const Parameters& global_parameters)
      57             : {
      58             :     // make sure all dimensions are set correctly and all indices are valid
      59        2817 :     assert(location.get_cells().size() == local_air_exposure.size<CellIndex>().get());
      60        2817 :     assert(location.get_cells().size() == local_contact_exposure.size<CellIndex>().get());
      61        2817 :     assert(local_contact_exposure.size<VirusVariant>() == local_air_exposure.size<VirusVariant>());
      62        2817 :     assert(local_contact_exposure.size<VirusVariant>() == VirusVariant::Count);
      63        2817 :     assert(local_contact_exposure.size<AgeGroup>().get() == global_parameters.get_num_groups());
      64        2817 :     assert(person.get_age() < local_contact_exposure.size<AgeGroup>());
      65        5634 :     assert(std::all_of(person.get_cells().begin(), person.get_cells().end(), [&](const auto& cell) {
      66             :         return cell < location.get_cells().size();
      67             :     }));
      68             : 
      69        2817 :     if (person.get_infection_state(t) == InfectionState::Susceptible) {
      70        1980 :         auto& local_parameters = location.get_infection_parameters();
      71             :         // TODO: we need to define what a cell is used for, as the loop may lead to incorrect results for multiple cells
      72        1980 :         auto age_receiver          = person.get_age();
      73        1980 :         ScalarType mask_protection = person.get_mask_protective_factor(global_parameters);
      74        1980 :         assert(person.get_cells().size() && "Person is in multiple cells. Interact logic is incorrect at the moment.");
      75        3960 :         for (auto cell_index : person.get_cells()) {
      76        1980 :             std::pair<VirusVariant, ScalarType> local_indiv_trans_prob[static_cast<uint32_t>(VirusVariant::Count)];
      77        3960 :             for (uint32_t v = 0; v != static_cast<uint32_t>(VirusVariant::Count); ++v) {
      78        1980 :                 VirusVariant virus = static_cast<VirusVariant>(v);
      79        1980 :                 ScalarType local_indiv_trans_prob_v =
      80        1980 :                     (std::min(local_parameters.get<MaximumContacts>(),
      81        3960 :                               daily_transmissions_by_contacts(local_contact_exposure, cell_index, virus, age_receiver,
      82        1980 :                                                               local_parameters)) +
      83        3960 :                      daily_transmissions_by_air(local_air_exposure, cell_index, virus, global_parameters)) *
      84        3960 :                     dt.days() * (1 - mask_protection) * (1 - person.get_protection_factor(t, virus, global_parameters));
      85             : 
      86        1980 :                 local_indiv_trans_prob[v] = std::make_pair(virus, local_indiv_trans_prob_v);
      87             :             }
      88             :             VirusVariant virus =
      89        1980 :                 random_transition(personal_rng, VirusVariant::Count, dt,
      90        1980 :                                   local_indiv_trans_prob); // use VirusVariant::Count for no virus submission
      91        1980 :             if (virus != VirusVariant::Count) {
      92          54 :                 person.add_new_infection(Infection(personal_rng, virus, age_receiver, global_parameters, t + dt / 2,
      93          27 :                                                    mio::abm::InfectionState::Exposed, person.get_latest_protection(),
      94             :                                                    false)); // Starting time in first approximation
      95             :             }
      96             :         }
      97             :     }
      98        2817 :     person.add_time_at_location(dt);
      99        2817 : }
     100             : 
     101        2889 : void add_exposure_contribution(AirExposureRates& local_air_exposure, ContactExposureRates& local_contact_exposure,
     102             :                                const Person& person, const Location& location, const TimePoint t, const TimeSpan dt)
     103             : {
     104        2889 :     if (person.get_location() != location.get_id()) {
     105           0 :         mio::log_debug("In add_exposure_contribution: Person {} is not at Location {}", person.get_id().get(),
     106           0 :                        location.get_id().get());
     107             :     }
     108             : 
     109        2889 :     if (person.is_infected(t)) {
     110         846 :         auto& infection = person.get_infection();
     111         846 :         auto virus      = infection.get_virus_variant();
     112         846 :         auto age        = person.get_age();
     113             :         // average infectivity over the time step to second order accuracy using midpoint rule
     114        1692 :         for (CellIndex cell : person.get_cells()) {
     115         846 :             if (location.get_infection_parameters().get<UseLocationCapacityForTransmissions>()) {
     116           0 :                 local_air_exposure[{cell, virus}] +=
     117           0 :                     infection.get_infectivity(t + dt / 2) *
     118           0 :                     location.get_cells()[cell.get()].compute_space_per_person_relative();
     119             :             }
     120             :             else {
     121        1692 :                 local_air_exposure[{cell, virus}] += infection.get_infectivity(t + dt / 2);
     122             :             }
     123        1692 :             local_contact_exposure[{cell, virus, age}] += infection.get_infectivity(t + dt / 2);
     124             :         }
     125             :     }
     126        2889 : }
     127             : 
     128         360 : bool change_location(Person& person, const Location& destination, const TransportMode mode,
     129             :                      const std::vector<uint32_t>& cells)
     130             : {
     131         738 :     assert(std::all_of(cells.begin(), cells.end(), [&](const auto& cell) {
     132             :         return cell < destination.get_cells().size();
     133             :     })); // make sure cell indices are valid
     134             : 
     135         360 :     if (person.get_location() != destination.get_id()) {
     136         351 :         person.set_location(destination.get_type(), destination.get_id());
     137         351 :         person.get_cells() = cells;
     138         351 :         person.set_last_transport_mode(mode);
     139             : 
     140         351 :         return true;
     141             :     }
     142             :     else {
     143          18 :         mio::log_debug("In change_location: Person {} already is at Location {}", person.get_id().get(),
     144          18 :                        destination.get_id().get());
     145           9 :         return false;
     146             :     }
     147             : }
     148             : 
     149             : } // namespace abm
     150             : } // namespace mio