module classroom_transmission use :: iso_c_binding use :: fgsl use :: utilities use :: prob use :: sort use :: graphs, only : ContactGraph use :: classroom_params use :: classroom_member use :: classroom_records implicit none private integer, parameter :: dp = c_double integer, parameter :: i4 = c_int integer, parameter :: i8 = c_int64_t public :: simulate type, public :: ClassRoom type(Parameters) :: Params type(Pupil), dimension(:), allocatable :: pupils type(Teacher), dimension(:), allocatable :: teachers integer(kind=c_int) :: time_step integer(kind=c_int64_t) :: day type(Records) :: result type(ContactGraph), dimension(:), allocatable :: contact_graph logical :: new_unwell logical :: class_quarantine contains private procedure, pass(self) :: addPupil generic, public :: add_pupil => addPupil procedure, pass(self) :: addTeacher generic, public :: add_teacher => addTeacher procedure, pass(self) :: addTempTeacher generic, public :: add_temp_teacher => addTempTeacher procedure, pass(self) :: replaceTempTeacher generic, public :: replace_temp_teacher => replaceTempTeacher procedure, pass(self) :: removeTempTeacher generic, public :: remove_temp_teacher => removeTempTeacher procedure, pass(self) :: getCounts generic, public :: counts => getCounts procedure, pass(self) :: getInClassCount generic, public :: inclass_count => getInClassCount procedure, pass(self) :: getSusceptibleCount generic, public :: susceptible_count => getSusceptibleCount procedure, pass(self) :: getExposedCount generic, public :: exposed_count => getExposedCount procedure, pass(self) :: getUnwellCount generic, public :: unwell_count => getUnwellCount procedure, pass(self) :: getQuarantinedCount generic, public :: quarantined_count => getQuarantinedCount procedure, pass(self) :: getRecoveredCount generic, public :: recovered_count => getRecoveredCount procedure, pass(self) :: getInClass generic, public :: get_inclass => getInClass procedure, pass(self) :: getInfected generic, public :: get_infected => getInfected procedure, pass(self) :: getNotInfected generic, public :: get_not_infected => getNotInfected procedure, public, pass(self) :: get_num_infected_pupils procedure, pass(self) :: new_infection_i4, new_infection_i8 generic, public :: new_infection => new_infection_i4, new_infection_i8 procedure, pass(self) :: getSusceptible generic, public :: get_susceptible => getSusceptible procedure, pass(self) :: getExposed generic, public :: get_exposed => getExposed procedure, pass(self) :: getUnwell generic, public :: get_unwell => getUnwell procedure, pass(self) :: getQuarantined generic, public :: get_quarantined => getQuarantined procedure, pass(self) :: getRecovered generic, public :: get_recovered => getRecovered procedure, pass(self) :: updateDisease generic, public :: update_disease => updateDisease procedure, pass(self) :: updateQuarantine_pupils generic, public :: update_pupils_quarantine => updateQuarantine_pupils procedure, pass(self) :: updateQuarantine_teachers generic, public :: update_teachers_quarantine => updateQuarantine_teachers procedure, pass(self) :: QuarantineClass generic, public :: quarantine_class => QuarantineClass procedure, pass(self) :: updateQuarantineClass generic, public :: update_quarantine_class => updateQuarantineClass procedure, pass(self) :: updateInClass generic, public :: update_in_class => updateInClass procedure, pass(self) :: updateCommunity generic, public :: update_community => updateCommunity procedure, pass(self) :: updateRapidTest generic, public :: update_rapid_test => updateRapidTest procedure, pass(self) :: updatePCRTest generic, public :: update_pcr_test => updatePCRTest procedure, pass(self) :: updateClassroom generic, public :: update => updateClassroom final :: ClassRoom_destruct end type ClassRoom interface ClassRoom module procedure init_ClassRoom_Params end interface ClassRoom contains function init_ClassRoom_Params(Params) result(self) ! Constructor for ClassRoom object ! Initialize ClassRoom object of pupils and teacher using parameters Params. ! Inputs: ! Params - Parameters object containing parameter settings implicit none type(ClassRoom) :: self type(Parameters), intent(in) :: Params ! Indexing integer(kind=c_int) :: j integer(kind=c_int64_t) :: jj ! Counters integer(kind=c_int) :: num_susceptible ! Number in Susceptible compartment integer(kind=c_int) :: num_exposed ! Number in Exposed compartment integer(kind=c_int) :: num_unwell ! Number in Unwell compartment integer(kind=c_int) :: num_recovered ! Number is Recovered compartment integer(kind=c_int) :: num_quarantined ! Number in Quarantine compartment integer(kind=c_int) :: num_in_class ! Number in class integer(kind=c_int) :: num_infected ! Number currently infected integer(kind=c_int) :: num_infected_pupils ! Number of pupils currently infected integer(kind=c_int) :: num_classroom_transmission ! Cumulative number of classroom transmissions integer(kind=c_int) :: num_community_transmission ! Cumulative number of community transmissions logical :: infected_teacher ! Is the teacher infected? integer(kind=c_int64_t) :: N ! Class size N = Params%class_size ! Set class size from Params self%Params = Params ! Copy Params into ClassRoom object ! Initialize counters -- will go to initialize self%results num_susceptible = 0 num_exposed = 0 num_unwell = 0 num_quarantined = 0 num_recovered = 0 num_infected = 0 num_in_class = 0 num_infected_pupils = 0 num_classroom_transmission = 0 num_community_transmission = 0 ! Initialize time step self%time_step = 1 self%day = Params%start_day ! Get starting day of the week from Params ! Add pupils to the ClassRoom object do j=1,Params%class_size call self%add_pupil() end do ! Add teacher to the ClassRoom object call self%add_teacher() ! If the teacher is currently in quarantine, need to add a temporary teacher if (self%teachers(1)%quarantined) then call self%add_temp_teacher() end if infected_teacher = self%teachers(1)%infected self%new_unwell = .False. ! On initialization, the teacher cannot by unwell self%class_quarantine = .False. ! On initialization, the class is not in Quarantine ! Update counters ! First loop over pupils... do j=1,size(self%pupils) if (self%pupils(j)%susceptible) num_susceptible = num_susceptible + 1 if (self%pupils(j)%exposed) num_exposed = num_exposed + 1 if (self%pupils(j)%unwell) num_unwell = num_unwell + 1 if (self%pupils(j)%recovered) num_recovered = num_recovered + 1 if (self%pupils(j)%quarantined) num_quarantined = num_quarantined + 1 if (self%pupils(j)%infected) then num_infected = num_infected + 1 num_infected_pupils = num_infected_pupils + 1 num_community_transmission = num_community_transmission + 1 ! On initialization, any infections must have come from community transmission end if if (self%pupils(j)%in_class) num_in_class = num_in_class + 1 end do ! ...and now over the teachers do j=1,size(self%teachers) if (self%teachers(j)%susceptible) num_susceptible = num_susceptible + 1 if (self%teachers(j)%exposed) num_exposed = num_exposed + 1 if (self%teachers(j)%unwell) num_unwell = num_unwell + 1 if (self%teachers(j)%recovered) num_recovered = num_recovered + 1 if (self%teachers(j)%quarantined) num_quarantined = num_quarantined + 1 if (self%teachers(j)%infected) then num_infected = num_infected + 1 num_community_transmission = num_community_transmission + 1 end if end do ! Initialize contact_graph allocate(self%contact_graph(self%Params%days+1)) ! Allocate memory for contact graphs ! Initialize records self%result = Records(self%Params%days+1, & self%Params%class_size, & num_in_class, & num_infected, & num_infected_pupils, & infected_teacher, & num_susceptible, & num_exposed, & num_unwell, & num_quarantined, & num_recovered, & num_community_transmission) ! Add pupils to the appropriate SEUQR compartments do jj=1,size(self%pupils) if (self%pupils(jj)%susceptible) call push(self%result%cat_lists(1)%susceptible, jj) if (self%pupils(jj)%exposed) call push(self%result%cat_lists(1)%exposed, jj) if (self%pupils(jj)%unwell) call push(self%result%cat_lists(1)%unwell, jj) if (self%pupils(jj)%quarantined) call push(self%result%cat_lists(1)%quarantined, jj) if (self%pupils(jj)%recovered) call push(self%result%cat_lists(1)%recovered, jj) end do ! Add teachers to the appropriate SEUQR compartments do jj=1,size(self%teachers) if (self%teachers(jj)%susceptible) call push(self%result%cat_lists(1)%susceptible, jj+N) if (self%teachers(jj)%exposed) call push(self%result%cat_lists(1)%exposed, jj+N) if (self%teachers(jj)%unwell) call push(self%result%cat_lists(1)%unwell, jj+N) if (self%teachers(jj)%quarantined) call push(self%result%cat_lists(1)%quarantined, jj+N) if (self%teachers(jj)%recovered) call push(self%result%cat_lists(1)%recovered, jj+N) end do return end function init_ClassRoom_Params subroutine addPupil(self) ! Add a pupil to the ClassRoom object implicit none class(ClassRoom), intent(inout) :: self type(Pupil) :: new_pupil ! Pupil object to contain the new pupil type(Pupil), dimension(:), allocatable :: pupils ! Temporary array of Pupil objects to allow new pupil to be added logical :: infected ! Is the new pupil infected? logical :: vaccinated ! Is the new pupil vaccinated? logical :: symptomatic ! Will the new pupil be symptomatic if infected? integer :: N ! number of pupils real(kind=c_double) :: prevalence ! community infection prevalence for children real(kind=c_double) :: p_symptomatic ! probability of symptomatic infection for children ! Is the pupil vaccinated? vaccinated = pdists%bernoulli_rvs(self%Params%child_P_vaccinated) ! Perform a Bernoulli trial using proportion of child vaccination ! Is the pupil infected? prevalence = self%Params%timeseries%child_prevalence(self%time_step) ! Child prevalence at the time pupil is added if (vaccinated) prevalence = prevalence*self%Params%Vaccine_Susceptible_prob_factor ! Reduce probability of infection for vaccinated person infected = pdists%bernoulli_rvs(prevalence) ! Perform Bernoulli trial using prevalence data ! Will this pupil be symptomatic? p_symptomatic = self%Params%child_P_symptomatic ! Probability of symptomatic infection for children if (vaccinated) p_symptomatic = p_symptomatic*self%Params%Vaccine_Symptomatic_prob_factor ! Reduce probability of symptoms for vaccinated person symptomatic = pdists%bernoulli_rvs(p_symptomatic) ! Perform Bernoulli trial using probability of symptoms for children ! Initialize this Pupil new_pupil = Pupil(self%Params, infected, symptomatic, vaccinated) ! Add new pupil to pupils array in ClassRoom object if (allocated(self%pupils)) then ! This is not the first pupil added, so... N = size(self%pupils) allocate(pupils(N+1)) ! allocate temporary array pupils(1:N) = self%pupils ! Copy existing pupils into temporary array pupils(N+1) = new_pupil ! Add new pupil call move_alloc(pupils, self%pupils) ! Copy (and destroy) temporary array to ClassRoom object else ! This is the first pupil added allocate(self%pupils(1)) ! Allocate Pupils array in ClassRoom object self%pupils(1) = new_pupil ! Add the first pupil end if return end subroutine addPupil subroutine addTeacher(self) ! Add a teacher to the ClassRoom object implicit none class(ClassRoom), intent(inout) :: self type(Teacher) :: new_teacher ! Teacher object to contain the new teacher type(Teacher), dimension(:), allocatable :: teachers ! Temporary array of Teacher objects to allow new teacher to be added logical :: infected ! Is the new teacher infected? logical :: vaccinated ! Is the new teacher vaccinated? logical :: symptomatic ! Will the new teacher be symptomatic if infected? integer :: N ! number of teachers (could be temporary) real(kind=c_double) :: prevalence ! community infection prevalence for adults real(kind=c_double) :: p_symptomatic ! probability of symptomatic infection for adults ! Is the teacher vaccinated? vaccinated = pdists%bernoulli_rvs(self%Params%adult_P_vaccinated) ! Perform a Bernoulli trial using proportion of adult vaccination ! Is the teacher infected? prevalence = self%Params%timeseries%adult_prevalence(self%time_step) ! Adult prevalence at the time teacher is added if (vaccinated) prevalence = prevalence*self%Params%Vaccine_Susceptible_prob_factor ! Reduce probability of infection for vaccinated person infected = pdists%bernoulli_rvs(prevalence) ! Perform Bernoulli trial using prevalence data ! Will this teacher be symptomatic? p_symptomatic = self%Params%adult_P_symptomatic ! Probability of symptomatic infection for adults if (vaccinated) p_symptomatic = p_symptomatic*self%Params%Vaccine_Symptomatic_prob_factor ! Reduce probability of symptoms for vaccinated person symptomatic = pdists%bernoulli_rvs(p_symptomatic) ! Perform Bernoulli trial using probability of symptoms for adults ! Initialize this Teacher new_teacher = Teacher(self%Params, infected, symptomatic, vaccinated, temporary=.False.) ! Add new teacher to teachers array in ClassRoom object if (allocated(self%teachers)) then ! This is not the first teacher added, so... N = size(self%teachers) allocate(teachers(N+1)) ! allocate temporary array teachers(1:N) = self%teachers ! Copy existing teachers into temporary array teachers(N+1) = new_teacher ! Add new teacher call move_alloc(teachers, self%teachers) ! Copy (and destroy) temporary array to ClassRoom object else ! This is the first teacher added allocate(self%teachers(1)) ! Allocate Teachers array in ClassRoom object self%teachers(1) = new_teacher ! Add the first teacher end if return end subroutine addTeacher subroutine addTempTeacher(self) ! Add a temporary teacher to the ClassRoom object implicit none class(ClassRoom), intent(inout) :: self type(Teacher) :: new_teacher ! Teacher object to contain the new teacher type(Teacher), dimension(:), allocatable :: teachers ! Temporary array of Teacher objects to allow new teacher to be added ! We cannot add a temporary teacher who is known to be infected (either in Unwell compartment, or through testing) logical :: bad_new_teacher ! Is this an unacceptable temporary teacher? logical :: vaccinated ! Is the new teacher infected? logical :: infected ! Is the new teacher vaccinated? logical :: symptomatic ! Will the new teacher be symptomatic if infected? real(kind=c_double) :: prevalence ! community infection prevalence for adults real(kind=c_double) :: p_symptomatic ! probability of symptomatic infection for adults integer(kind=c_int64_t) :: N N = self%Params%class_size ! Is the teacher vaccinated? vaccinated = pdists%bernoulli_rvs(self%Params%adult_P_vaccinated) ! Perform a Bernoulli trial using proportion of adult vaccination ! Is the teacher infected? prevalence = self%Params%timeseries%adult_prevalence(self%time_step) ! Adult prevalence at the time teacher is added if (vaccinated) prevalence = prevalence*self%Params%Vaccine_Susceptible_prob_factor ! Reduce probability of infection for vaccinated person infected = pdists%bernoulli_rvs(prevalence) ! Perform Bernoulli trial using prevalence data ! Will this teacher be symptomatic? p_symptomatic = self%Params%adult_P_symptomatic ! Probability of symptomatic infection for adults if (vaccinated) p_symptomatic = p_symptomatic*self%Params%Vaccine_Symptomatic_prob_factor ! Reduce probability of symptoms for vaccinated person symptomatic = pdists%bernoulli_rvs(p_symptomatic) ! Perform Bernoulli trial using probability of symptoms for adults ! Initialize this Teacher, making sure they are not in quarantine bad_new_teacher = .True. do while (bad_new_teacher) ! Keep initiating potential temporary teachers until an acceptable temporary teacher is found new_teacher = Teacher(self%Params, infected, symptomatic, vaccinated, temporary=.True.) ! Instantiate a Teacher bad_new_teacher = new_teacher%quarantined ! Unacceptable teacher if in Quarantine when instantiated end do ! If we need a temporary teacher, either: ! - there is currently a temporary teacher who should be replaced ! or - the permanent teacher has been Quarantined, so a temporary is needed allocate(teachers(2)) ! Allocate temporary teachers array teachers(1) = self%teachers(1) ! Copy permanent teacher teachers(2) = new_teacher ! Add temporary teacher call move_alloc(teachers, self%teachers) ! Copy (and destroy) temporary array to ClassRoom object ! The new teacher temporary could be infected when added. ! If so, need to update infection records if (new_teacher%infected) then call self%new_infection(N+2, .False.) end if return end subroutine addTempTeacher subroutine removeTempTeacher(self) ! Remove temporary teacher from the ClassRoom object implicit none class(ClassRoom), intent(inout) :: self type(Teacher), dimension(:), allocatable :: teachers ! Temporary array of Teacher objects ! Check if trying to remove temporary teacher that doesn't exist. ! This is a fatal error. if (size(self%teachers)/=2) then call FatalErrorMessage("Trying to remove temp teacher but none are added") end if allocate(teachers(1)) ! Allocate temporary array for copying teachers(1) = self%teachers(1) ! Copy existing permanent teacher call move_alloc(teachers, self%teachers) ! Copy (and destroy) temporary array to ClassRoom object return end subroutine removeTempTeacher subroutine replaceTempTeacher(self) ! Replace temporary teacher in the ClassRoom object implicit none class(ClassRoom), intent(inout) :: self ! Replace the temporary teach by first removing... call self%remove_temp_teacher() ! ... then adding a new temporary teacher. call self%add_temp_teacher() return end subroutine replaceTempTeacher subroutine ClassRoom_destruct(self) ! Destructor for ClassRoom object ! Not much is done here -- allocated array are simply deallocated implicit none type(ClassRoom) :: self if (allocated(self%pupils)) deallocate(self%pupils) if (allocated(self%teachers)) deallocate(self%teachers) if (allocated(self%contact_graph)) deallocate(self%contact_graph) end subroutine ClassRoom_destruct function get_num_infected_pupils(self) result(num_infected_pupils) ! Get the number of infected pupils in the ClassRoom object ! Output: ! num_infected_pupils [c_int] -- the number of infected pupils implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: num_infected_pupils integer(kind=c_int64_t), dimension(:), allocatable :: infected ! Temporary array holding indices of infected individuals (including teachers) call self%get_infected(infected) ! Get indices of infected pupils. if (allocated(infected)) then ! If there are any infected pupils, then infected will be allocated num_infected_pupils = size(infected) ! Get size of infected array if (in_array(infected, self%Params%class_size+1)) then ! Take away teacher, if they are infected num_infected_pupils = num_infected_pupils - 1 end if if (in_array(infected, self%Params%class_size+2)) then ! Take away temporary teacher, if they are infected num_infected_pupils = num_infected_pupils - 1 end if else ! If infected is not allocated, then there are no infected pupils num_infected_pupils = 0 end if deallocate(infected) ! Deallocate array of indices of infected individuals return end function get_num_infected_pupils subroutine getCounts(self) ! Update counters and indices lists in the ClassRoom object at current time step implicit none class(ClassRoom), intent(inout) :: self integer(kind=c_int) :: num_infected_pupils ! Number of pupils currently infected integer(kind=c_int) :: num_in_class ! Number in class integer(kind=c_int) :: infected_teacher ! Is the teacher infected? integer(kind=c_int) :: susceptible ! Number in Susceptible compartment integer(kind=c_int) :: exposed ! Number in Exposed compartment integer(kind=c_int) :: unwell ! Number in Unwell compartment integer(kind=c_int) :: quarantined ! Number in Quarantine compartment integer(kind=c_int) :: recovered ! Number is Recovered compartment integer :: j ! Initialize counters num_infected_pupils = 0 num_in_class = 0 infected_teacher = 0 susceptible = 0 exposed = 0 unwell = 0 quarantined = 0 recovered = 0 ! Update ClassRoom objects' lists of indices of agents in SEUQR compartment at current time step call self%get_susceptible(self%result%cat_lists(self%time_step+1)%susceptible) call self%get_exposed(self%result%cat_lists(self%time_step+1)%exposed) call self%get_unwell(self%result%cat_lists(self%time_step+1)%unwell) call self%get_quarantined(self%result%cat_lists(self%time_step+1)%quarantined) call self%get_recovered(self%result%cat_lists(self%time_step+1)%recovered) ! Loop over pupils and count number of infected pupils and number in class do j=1, size(self%pupils) if (self%pupils(j)%infected) num_infected_pupils = num_infected_pupils + 1 if (self%pupils(j)%in_class) num_in_class = num_in_class + 1 end do ! Add teachers do j=1, size(self%teachers) if (self%teachers(j)%infected) infected_teacher = infected_teacher + 1 end do ! Get number in Susceptible compartment if (allocated(self%result%cat_lists(self%time_step+1)%susceptible)) then susceptible = size(self%result%cat_lists(self%time_step+1)%susceptible) else susceptible = 0 end if call push(self%result%susceptible, susceptible) ! Update the ClassRoom object susceptible list for current time step ! Get number in Exposed compartment if (allocated(self%result%cat_lists(self%time_step+1)%exposed)) then exposed = size(self%result%cat_lists(self%time_step+1)%exposed) else exposed = 0 end if call push(self%result%exposed, exposed) ! Update the ClassRoom object exposed list for current time step if (allocated(self%result%cat_lists(self%time_step+1)%unwell)) then unwell = size(self%result%cat_lists(self%time_step+1)%unwell) else unwell = 0 end if call push(self%result%unwell, unwell) if (allocated(self%result%cat_lists(self%time_step+1)%quarantined)) then quarantined = size(self%result%cat_lists(self%time_step+1)%quarantined) else quarantined = 0 end if call push(self%result%quarantined, quarantined) if (allocated(self%result%cat_lists(self%time_step+1)%recovered)) then recovered = size(self%result%cat_lists(self%time_step+1)%recovered) else recovered = 0 end if call push(self%result%recovered, recovered) call push(self%result%num_infected_pupils, num_infected_pupils) call push(self%result%num_in_class, num_in_class) call push(self%result%infected_teacher, infected_teacher) return end subroutine getCounts pure function getInClassCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%in_class) count = count + 1 end do return end function getInClassCount pure function getSusceptibleCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%susceptible) count = count + 1 end do do j=1,size(self%teachers) if (self%teachers(j)%susceptible) count = count + 1 end do return end function getSusceptibleCount pure function getExposedCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%exposed) count = count + 1 end do do j=1,size(self%teachers) if (self%teachers(j)%exposed) count = count + 1 end do return end function getExposedCount pure function getUnwellCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%unwell) count = count + 1 end do do j=1,size(self%teachers) if (self%teachers(j)%unwell) count = count + 1 end do return end function getUnwellCount pure function getQuarantinedCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%quarantined) count = count + 1 end do do j=1,size(self%teachers) if (self%teachers(j)%quarantined) count = count + 1 end do return end function getQuarantinedCount pure function getRecoveredCount(self) result(count) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int) :: count integer :: j count=0 do j=1,size(self%pupils) if (self%pupils(j)%recovered) count = count + 1 end do do j=1,size(self%teachers) if (self%teachers(j)%recovered) count = count + 1 end do end function getRecoveredCount function getInClass(self) result(inclass) implicit none class(Classroom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable :: inclass integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%in_class) call push(inclass, j) end do do j=1,size(self%teachers) if (self%teachers(j)%in_class) call push(inclass, j+N) end do return end function getInClass subroutine getSusceptible(self, susceptible) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: susceptible integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%susceptible) then call push(susceptible, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%susceptible) then call push(susceptible, j+N) end if end do return end subroutine getSusceptible subroutine getExposed(self, exposed) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: exposed integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%exposed) then call push(exposed, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%exposed) then call push(exposed, j+N) end if end do return end subroutine getExposed subroutine getUnwell(self, unwell) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: unwell integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%unwell) then call push(unwell, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%unwell) then call push(unwell, j+N) end if end do return end subroutine getUnwell subroutine getQuarantined(self, quarantined) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: quarantined integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%quarantined) then call push(quarantined, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%quarantined) then call push(quarantined, j+N) end if end do return end subroutine getQuarantined subroutine getRecovered(self, recovered) implicit none class(ClassRoom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: recovered integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%recovered) then call push(recovered, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%recovered) then call push(recovered, j+N) end if end do return end subroutine getRecovered subroutine getInfected(self, infected) implicit none class(Classroom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable, intent(out) :: infected integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (self%pupils(j)%infected) then call push(infected, j) end if end do do j=1,size(self%teachers) if (self%teachers(j)%infected) then call push(infected, j+N) end if end do return end subroutine getInfected function getNotInfected(self) result(not_infected) implicit none class(Classroom), intent(in) :: self integer(kind=c_int64_t), dimension(:), allocatable :: not_infected integer(kind=c_int64_t) :: j, N N = size(self%pupils) do j=1,N if (.not. self%pupils(j)%infected) then call push(not_infected, j) end if end do do j=1,size(self%teachers) if (.not. self%teachers(j)%infected) then call push(not_infected, j+N) end if end do return end function getNotInfected subroutine new_infection_i4(self, indx, class_transmission) ! Update member of ClassRoom for new infection implicit none class(ClassRoom), intent(inout) :: self integer(kind=c_int), intent(in) :: indx ! Index of member with new infection logical, intent(in) :: class_transmission ! Is this a classroom transmission? integer(kind=c_int) :: N N = size(self%pupils) if (indx<=N) then ! indx corresponds to pupil ! Mark as infected self%pupils(indx)%infected = .True. ! Move from Susceptible to Exposed self%pupils(indx)%susceptible = .False. self%pupils(indx)%exposed = .True. ! update total infected pupils self%result%total_infected_pupils = self%result%total_infected_pupils + 1 ! start secondary infections counter self%result%secondary_infections(indx) = 0_i8 else ! indx corresponds to teacher self%teachers(indx-N)%infected = .True. ! Move from Susceptible to Exposed self%teachers(indx-N)%susceptible = .False. self%teachers(indx-N)%exposed = .True. ! start secondary infections counter self%result%secondary_infections(indx) = 0_i8 end if ! update total infected self%result%total_infected = self%result%total_infected + 1 if (class_transmission) then self%result%num_classroom_transmission = self%result%num_classroom_transmission + 1 else self%result%num_community_transmission = self%result%num_community_transmission + 1 end if return end subroutine new_infection_i4 subroutine new_infection_i8(self, indx, class_transmission) ! Update member of ClassRoom for new infection implicit none class(ClassRoom), intent(inout) :: self integer(kind=c_int64_t), intent(in) :: indx ! Index of member with new infection logical, intent(in) :: class_transmission ! Is this a classroom transmission? integer(kind=c_int64_t) :: N N = size(self%pupils) if (indx<=N) then ! indx corresponds to pupil ! Mark as infected self%pupils(indx)%infected = .True. ! Move from Susceptible to Exposed self%pupils(indx)%susceptible = .False. self%pupils(indx)%exposed = .True. ! update total infected pupils self%result%total_infected_pupils = self%result%total_infected_pupils + 1 ! start secondary infections counter self%result%secondary_infections(indx) = 0_i8 else ! indx corresponds to teacher self%teachers(indx-N)%infected = .True. ! Move from Susceptible to Exposed self%teachers(indx-N)%susceptible = .False. self%teachers(indx-N)%exposed = .True. ! start secondary infections counter self%result%secondary_infections(indx) = 0_i8 end if ! update total infected self%result%total_infected = self%result%total_infected + 1 if (class_transmission) then self%result%num_classroom_transmission = self%result%num_classroom_transmission + 1 else self%result%num_community_transmission = self%result%num_community_transmission + 1 end if return end subroutine new_infection_i8 subroutine updateDisease(self) implicit none class(Classroom), intent(inout) :: self integer :: j type(Pupil) :: pupil type(Teacher) :: teacher self%new_unwell = .False. ! Loop over pupils do j=1,size(self%pupils) pupil = self%pupils(j) ! Is member infected? if (pupil%infected) then ! Yes -- advance their days_infected pupil%days_infected = pupil%days_infected+1 end if ! Go through compartments and advance as appropriate ! If in exposed... if (pupil%exposed) then ! and symptomatic, then might need to move to unwell if (pupil%symptomatic) then ! if they have been infected for the incubation time if (pupil%days_infected > pupil%incubation_time) then pupil%exposed = .False. pupil%unwell = .True. self%new_unwell = .True. end if else ! not symptomatic, so might need to move to recovered ! if the days_infected is greater than infection clearing time if (pupil%days_infected > pupil%infection_clearing_time) then pupil%exposed = .False. pupil%unwell = .False. pupil%recovered = .True. pupil%infected = .False. pupil%immune = .True. end if end if ! If in unwell... else if (pupil%unwell) then ! must be in symptomatic compartment, so don't need to check pupil%days_unwell = pupil%days_unwell + 1 ! If infected for longer than infection clearing time if (pupil%days_infected > pupil%infection_clearing_time) then ! move to the recovered compartment pupil%exposed = .False. pupil%unwell = .False. pupil%recovered = .True. pupil%infected = .False. pupil%immune = .True. end if end if self%pupils(j) = pupil end do ! Loop over teachers do j=1,size(self%teachers) teacher = self%teachers(j) ! Is member infected? if (teacher%infected) then ! Yes -- advance their days_infected teacher%days_infected = teacher%days_infected+1 end if ! Go through compartments and advance as appropriate ! If in exposed... if (teacher%exposed) then ! and symptomatic, then might need to move to unwell if (teacher%symptomatic) then ! if they have been infected for the incubation time if (teacher%days_infected > teacher%incubation_time) then teacher%exposed = .False. teacher%unwell = .True. end if else ! not symptomatic, so might need to move to recovered ! if the days_infected is greater than infection clearing time if (teacher%days_infected > teacher%infection_clearing_time) then teacher%exposed = .False. teacher%unwell = .False. teacher%recovered = .True. teacher%infected = .False. teacher%immune = .True. end if end if ! If in unwell... else if (teacher%unwell) then ! must be in symptomatic compartment, so don't need to check teacher%days_unwell = teacher%days_unwell + 1 ! If infected for longer than infection clearing time if (teacher%days_infected > teacher%infection_clearing_time) then ! move to the recovered compartment teacher%exposed = .False. teacher%unwell = .False. teacher%recovered = .True. teacher%infected = .False. teacher%immune = .True. end if end if self%teachers(j) = teacher end do return end subroutine updateDisease subroutine updateQuarantine_pupils(self) implicit none class(Classroom), intent(inout) :: self integer :: j type(Pupil) :: p ! Loop over pupils do j=1,size(self%pupils) p = self%pupils(j) ! Is member in quarantine? if (p%quarantined) then p%days_quarantined = p%days_quarantined + 1 if ((self%day /= 6) .and. (self%day /= 0)) then self%result%pupil_absent_days(j) = self%result%pupil_absent_days(j) + 1 end if ! If they have been quarantined for long enough, release them if ((.not. p%unwell) .and. (p%days_quarantined > self%Params%quarantine_days)) then p%quarantined = .False. p%pcr_test_positive = .False. end if else ! If they are not quarantined, they could be if they have detectable symptoms if (p%symptomatic .and. p%unwell) then p%quarantined = pdists%bernoulli_rvs(self%Params%Detection_prob_sympt) p%days_quarantined = 0 end if end if p%in_class = .not. p%quarantined self%pupils(j) = p end do end subroutine updateQuarantine_pupils subroutine updateQuarantine_teachers(self) implicit none class(ClassRoom), intent(inout) :: self type(Teacher) :: perm, temp ! Loop over teachers -- there is always one and potentially two: one permanent, the other temporary perm = self%teachers(1) ! this is the permanent teacher ! Is the permanent teacher in quarantine? if (perm%quarantined) then ! Yes, so there is a temporary teacher too temp = self%teachers(2) perm%days_quarantined = perm%days_quarantined + 1 ! If they have been quarantined for long enough, release them if ((.not. perm%unwell) .and. (perm%days_quarantined > self%Params%quarantine_days)) then perm%quarantined = .False. call self%remove_temp_teacher() else ! Permanent teacher still in quarantine, so advance temp teacher. ! The temp cannot be currently in quarantine ! but may be entering quarantine, if they are symptomatic if (temp%symptomatic .and. temp%unwell) then call self%replace_temp_teacher() end if end if else ! If they are not quarantined, they could be if they have detectable symptoms if (perm%unwell .and. perm%in_class) then perm%quarantined = .True. perm%days_quarantined = 0 call self%add_temp_teacher() end if end if perm%in_class = .not. perm%quarantined self%teachers(1) = perm return end subroutine updateQuarantine_teachers subroutine QuarantineClass(self) implicit none class(Classroom), intent(inout) :: self integer :: j type(Pupil) :: p ! Loop over pupils do j=1,size(self%pupils) p = self%pupils(j) p%quarantined = .True. p%days_quarantined = 0 p%in_class = .False. self%pupils(j) = p end do self%teachers(1)%quarantined = .True. self%teachers(1)%days_quarantined = 0 self%teachers(1)%in_class = .False. self%class_quarantine = .True. if (size(self%teachers)==2) then call self%remove_temp_teacher() end if return end subroutine QuarantineClass subroutine updateQuarantineClass(self) implicit none class(ClassRoom), intent(inout) :: self integer :: j type(Pupil) :: p type(Teacher) :: t do j=1,size(self%pupils) p = self%pupils(j) p%days_quarantined = p%days_quarantined + 1 if ((self%day /= 6) .and. (self%day /= 0)) then self%result%pupil_absent_days(j) = self%result%pupil_absent_days(j) + 1 end if if (p%days_quarantined < self%Params%quarantine_days+1) then p%in_class = .False. p%quarantined = .True. else if (.not. p%unwell) then p%in_class = .True. p%quarantined = .False. end if end if self%pupils(j) = p end do t = self%teachers(1) t%days_quarantined = t%days_quarantined + 1 if (t%days_quarantined < self%Params%quarantine_days+1) then t%quarantined = .True. t%in_class = .False. else self%class_quarantine = .False. if (.not. t%unwell) then t%in_class = .True. t%quarantined = .False. else t%in_class = .False. t%quarantined = .True. call self%add_temp_teacher() end if end if self%teachers(1) = t return end subroutine updateQuarantineClass subroutine updateInClass(self) implicit none class(Classroom), intent(inout) :: self integer(kind=c_int64_t) :: indxA, indxB class(ClassMember), allocatable :: memberA, memberB integer(kind=c_int64_t), dimension(:), allocatable :: inclass integer(kind=c_int64_t), dimension(:), allocatable :: infected integer(kind=c_int64_t), dimension(:), allocatable :: exposed integer(kind=c_int64_t), dimension(:), allocatable :: inclass_degrees real(kind=c_double) :: transmission_factor real(kind=c_double) :: transmission_prob integer(kind=c_int64_t) :: j logical :: transmission logical :: any_infected type(ContactGraph) :: ContactNet integer(kind=c_int64_t) :: N N = size(self%pupils) call self%get_infected(infected) call self%get_exposed(exposed) ! Who is still in the class? ! These will be the nodes in the contact graph. inclass = self%get_inclass() ! Is anyone exposed? if (size(exposed)>0) then any_infected=.True. else any_infected=.False. end if ! Get the degrees (i.e. number of contacts) of the nodes allocate(inclass_degrees(size(inclass))) do j=1,size(inclass) if (j<=N) then ! inclass_degrees(j) = self%pupils(j)%reduced_daily_contacts inclass_degrees(j) = self%pupils(j)%daily_contacts else ! inclass_degrees(j) = self%teachers(j-N)%reduced_daily_contacts inclass_degrees(j) = self%teachers(j-N)%daily_contacts end if end do if ((any_infected) .or. (self%Params%contact_net)) then ! Build the contact network ContactNet = ContactGraph(inclass, inclass_degrees) ! Loop through the edges of the contact network do j=1, ContactNet%size transmission = .False. ! Default to no transmission if (allocated(memberA)) deallocate(memberA) if (allocated(memberB)) deallocate(memberB) ! Get indices of the contacts indxA = ContactNet%edges(j,1) indxB = ContactNet%edges(j,2) if (indxA<=N) then allocate(memberA, source=self%pupils(indxA)) else allocate(memberA, source=self%teachers(indxA-N)) end if if (indxB<=N) then allocate(memberB, source=self%pupils(indxB)) else allocate(memberB, source=self%teachers(indxB-N)) end if ! If it is a self-loop, then skip if (indxA==indxB) cycle ! If neither of the agents is infected, then skip if ((.not. memberA%infected) .and. (.not. memberB%infected)) then if (self%Params%contact_net) then call self%result%add_contact_record(self%time_step, indxA, indxB) end if cycle end if ! If both of the pupils are infected, then skip if ((memberA%infected) .and. (memberB%infected)) then if (self%Params%contact_net) then call self%result%add_contact_record(self%time_step, indxA, indxB) end if cycle end if ! If we reach here, only one of indxA and indxB are infected, so transmission is possible ! Swap indices so that indxA is the infected individual, indxB is not if (memberB%infected) call swap(indxA, indxB) if (allocated(memberA)) deallocate(memberA) if (allocated(memberB)) deallocate(memberB) if (indxA<=N) then allocate(memberA, source=self%pupils(indxA)) else allocate(memberA, source=self%teachers(indxA-N)) end if if (indxB<=N) then allocate(memberB, source=self%pupils(indxB)) else allocate(memberB, source=self%teachers(indxB-N)) end if ! Now individual memberA is infected. ! For transmission, memberB must be susceptible ! If indxB is not susceptible, then skip if (.not. memberB%susceptible) then if (self%Params%contact_net) then call self%result%add_contact_record(self%time_step, indxA, indxB) end if cycle else ! infected -- susceptible contact. Possibility of transmission. ! Determine the transmission probability = transmission_prob*transmission_factor ! where transmission_factor depends on the length of time memberA has been infected ! relative to the tost_time transmission_factor = memberA%rel_trans%factor(memberA%days_infected-memberA%incubation_time-memberA%tost_time) transmission_prob = memberA%mitigation_factor * memberA%transmission_prob * transmission_factor ! If susceptible indxB is vaccinated, reduce transmission prob if (memberB%vaccinated) transmission_prob = transmission_prob * self%Params%Vaccine_Susceptible_prob_factor ! Perform Bernoulli trial for transmission in contact transmission = pdists%bernoulli_rvs(transmission_prob) if (transmission) then ! This contact results in transmission so update contacted individual call self%new_infection(indxB, .True.) self%result%secondary_infections(indxA) = self%result%secondary_infections(indxA) + 1 end if end if ! Update the records if we are storing the contact network if (self%Params%contact_net) then call self%result%add_contact_record(self%time_step, & indxA, & indxB, & transmission=transmission) end if end do end if return end subroutine updateInClass subroutine updateCommunity(self) implicit none class(ClassRoom), intent(inout) :: self integer :: j, N real(kind=c_double) :: transmission_prob logical :: infected N = size(self%pupils) ! Loop through pupils members do j=1,N ! If pupil is currently susceptible, they could become infected through community transmission if (self%pupils(j)%susceptible) then ! Get probability of infection from community incidence at time step transmission_prob = self%Params%timeseries%community_incidence(self%time_step) ! If vaccinated, reduce the transmission_prob if (self%pupils(j)%vaccinated) transmission_prob = transmission_prob * self%Params%Vaccine_Susceptible_prob_factor ! Draw from Bernoulli using community incidence infected = pdists%bernoulli_rvs(transmission_prob) if (infected) then ! New infection call self%new_infection(j, .False.) end if end if end do ! Loop through teachers members do j=1,size(self%teachers) ! If teacher is currently susceptible, they could become infected through community transmission if (self%teachers(j)%susceptible) then ! Get probability of infection from community incidence at time step transmission_prob = self%Params%timeseries%community_incidence(self%time_step) ! If vaccinated, reduce the transmission_prob if (self%teachers(j)%vaccinated) transmission_prob = transmission_prob * self%Params%Vaccine_Susceptible_prob_factor ! Draw from Bernoulli using community incidence infected = pdists%bernoulli_rvs(transmission_prob) if (infected) then ! New infection call self%new_infection(j+N, .False.) end if end if end do return end subroutine updateCommunity subroutine updateRapidTest(self) implicit none class(ClassRoom), intent(inout) :: self integer :: j type(Pupil) :: p type(Teacher) :: t logical :: rapid_test_positive ! Loop through pupils members do j=1,size(self%pupils) p = self%pupils(j) if ((.not. p%unwell) .and. (.not. p%pcr_test_positive)) then if (p%infected) then rapid_test_positive = pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_sensitivity) else rapid_test_positive = .not. pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_specificity) end if if (rapid_test_positive) then p%rapid_test_positive = .True. p%rapid_positive_day = self%time_step p%quarantined = .True. p%days_quarantined = 0 p%in_class = .False. self%pupils(j) = p if ((.not. self%Params%mitigations%pcr_test) .and. (self%Params%mitigations%class_quarantine)) call self%quarantine_class() end if end if end do ! Is there a temporary teacher if (size(self%teachers)==2) then t = self%teachers(2) ! The temporary teacher cannot be unwell so skip this check if (t%infected) then rapid_test_positive = pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_sensitivity) else rapid_test_positive = .not. pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_specificity) end if if (rapid_test_positive) then call self%replace_temp_teacher() else t%rapid_test_positive = .False. self%teachers(2) = t end if end if ! Consider now the permanent teacher t = self%teachers(1) if ((.not. t%unwell) .and. (.not. t%pcr_test_positive)) then if (t%infected) then rapid_test_positive = pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_sensitivity) else rapid_test_positive = .not. pdists%bernoulli_rvs(self%Params%mitigations%rapid_test_specificity) end if if (rapid_test_positive) then if (t%in_class) then call self%add_temp_teacher() end if t%rapid_test_positive = .True. t%rapid_positive_day = self%time_step t%quarantined = .True. t%days_quarantined = 0 t%in_class = .False. self%teachers(1) = t end if end if return end subroutine updateRapidTest subroutine updatePCRTest(self) implicit none class(ClassRoom), intent(inout) :: self integer :: j type(Pupil) :: p type(Teacher) :: t logical :: pcr_test_positive do j=1,size(self%pupils) p = self%pupils(j) if ((p%rapid_test_positive) .and. (self%time_step - p%rapid_positive_day==2)) then if (p%infected) then pcr_test_positive = pdists%bernoulli_rvs(self%Params%mitigations%pcr_test_sensitivity) else pcr_test_positive = .not. pdists%bernoulli_rvs(self%Params%mitigations%pcr_test_specificity) end if p%pcr_test_positive = pcr_test_positive if (pcr_test_positive) then if (self%Params%mitigations%class_quarantine) call self%quarantine_class() else p%rapid_test_positive = .False. p%quarantined = .False. p%days_quarantined = 0 p%in_class = .True. end if self%pupils(j) = p end if end do ! consider permanent teacher t = self%teachers(1) if ((t%rapid_test_positive) .and. (self%time_step - t%rapid_positive_day==2)) then if (t%infected) then pcr_test_positive = pdists%bernoulli_rvs(self%Params%mitigations%pcr_test_sensitivity) else pcr_test_positive = .not. pdists%bernoulli_rvs(self%Params%mitigations%pcr_test_specificity) end if t%pcr_test_positive = pcr_test_positive if (pcr_test_positive) then if (self%Params%mitigations%class_quarantine) call self%quarantine_class() else t%rapid_test_positive = .False. t%quarantined = .False. t%days_quarantined = 0 t%in_class = .True. if (size(self%teachers)==2) then call self%remove_temp_teacher() end if end if self%teachers(1) = t end if ! TODO: what to do with temporary teacher? end subroutine updatePCRTest subroutine updateClassroom(self) implicit none class(ClassRoom), intent(inout) :: self ! update from community transmissions call self%update_community() if (self%Params%mitigations%rapid_test) then if (any(self%Params%mitigations%rapid_test_days == self%day)) then call self%update_rapid_test() end if end if if (self%Params%mitigations%pcr_test) then call self%update_pcr_test() end if if (.not. self%class_quarantine) then ! If it is a school day, update classroom transmissions if ((self%day /= 6) .and. (self%day /= 0)) then call self%update_in_class() end if end if ! Update disease model call self%update_disease() if (.not. self%class_quarantine) then if ((self%new_unwell) .and. (self%Params%mitigations%class_quarantine)) then call self%quarantine_class() else call self%update_pupils_quarantine() call self%update_teachers_quarantine() end if else call self%update_quarantine_class() end if ! update records call self%counts() ! Step in time self%time_step = self%time_step + 1 self%day = mod(self%day+1, int(7, kind=c_int)) return end subroutine updateClassroom subroutine simulate(Params,results) implicit none type(Parameters), intent(in) :: Params type(Records), intent(out) :: results type(ClassRoom) :: cs integer :: j cs = ClassRoom(Params) do j=1,Params%days call cs%update() end do results = cs%result%copy() return end subroutine simulate end module classroom_transmission