Commit e91ab4c1 authored by Sébastien Picault's avatar Sébastien Picault
Browse files

Added model files and updated README

parent ee98c925
# ASF Challenge 2020
ASF Challenge 2020
========================
Public repository for the model and documents used in the ASF Challenge (Aug 2020-Jan 2021), coordinated by Pauline Ezanno, Sébastien Picault and Timothée Vergne (INRAE).
\ No newline at end of file
This is the public repository for the model and documents used in the
[ASF Challenge](https://www6.inrae.fr/asfchallenge/) (Aug 2020-Jan
2021), coordinated by Pauline Ezanno, Sébastien Picault and Timothée
Vergne (INRAE).
<!-- markdown-toc start - Don't edit this section. Run M-x markdown-toc-refresh-toc -->
**Table of Contents**
- [ASF Challenge 2020](#asf-challenge-2020)
- [Content](#content)
- [Running model M0](#running-model-m0)
- [Data provided to ASF Challenge participants ("players")](#data-provided-to-asf-challenge-participants-players)
<!-- markdown-toc end -->
Content
-------
Running model M0
----------------
Data provided to ASF Challenge participants ("players")
-------------------------------------------------------
"""This file is aimed at providing specific code add-on for the PPA
model (ppa.yaml). It is based on the code skeletton generated
automatically by EMULSION.
emulsion run ppa.yaml --output-dir outputs/test --no-count -r 1 -t 10 -p nb_infected_wild_boars=1 -p transmission_boar_pig=0.03
"""
import sys
from pathlib import Path
import csv
import numpy as np
import scipy.sparse as sp
from scipy.spatial import distance_matrix
import pandas as pd
from emulsion.agent.managers import CompartProcessManager
from emulsion.agent.managers import MetapopProcessManager
from emulsion.agent.views import AutoStructuredView
from emulsion.tools.functions import random_bool
#===============================================================
# CLASS Metapop (LEVEL 'metapop')
#===============================================================
class Metapop(MetapopProcessManager):
"""
Level of the spatial population.
"""
#----------------------------------------------------------------
# Level initialization
#----------------------------------------------------------------
def initialize_level(self, **others):
"""Initialize a metapop composed of pig herds and wild boars.
BEWARE: data files for herds and boars must have an index
column starting at 0 and incremented sequentially
"""
# set of tiles surrounded by fences
self.surrounded_tiles = set()
# set of tiles outside fences but just adjacent to those within fences
self.surrounding_tiles = set()
# array of indexes of wildboars within fences
self.wildboars_within_fences = np.array([], dtype='int32')
# array of indexes of wildboars within fences
self.wildboars_within_buffer = np.array([], dtype='int32')
## dictionary {step -> list of epidunits} for reactivating agents at specific time
self.activation_schedule = {}
## dictionary {step -> list of epidunits} for preventive culling
self.culling_schedule = {}
self.init_epid_unit_data()
self.create_initial_epidunits()
self.init_trade()
## init fences location from data
self.setup_fences()
def init_epid_unit_data(self):
print('Loading pig herd info')
# load pig herds info
herd_file = self.model.input_dir.joinpath('herds.csv')
self.herds = pd.read_csv(herd_file, sep=',', header=0)
self.nb_herds = self.herds.shape[0]
print('Nb pig herds:', self.nb_herds)
print('Loading wildboar info')
# load wildboars info
boar_file = self.model.input_dir.joinpath('boars.csv')
self.wildboars = pd.read_csv(boar_file, sep=',', header=0)
self.nb_wildboars = self.wildboars.shape[0]
print('Nb wildboars:', self.nb_wildboars)
# index -> healthy carcass there ?
self.boars_natural_carcass = np.zeros(self.nb_wildboars)
# index -> no more wildboar here (hunt or carcass removal)
self.removed_wildboars = np.zeros(self.nb_wildboars)
## use pig_herd sizes to know how many susceptible animals are
## present in each herd, and other associated pre-computed
## information used to calculate infection probabilities
print("Structuring herd biosecurity/outdoor information")
self.susceptibles = self.herds['size'].values.reshape((1, self.nb_herds))
self.biorisk = sp.csr_matrix(1 - self.herds['biosecurity'].values)
self.biorisk_outdoor = sp.csr_matrix(self.herds['is_outdoor'].values).multiply(self.biorisk)
print("Computing tiles")
self.points = np.concatenate([self.herds[['X', 'Y']].values, self.wildboars[['X', 'Y']].values])
self.Xmin, self.Ymin = self.points.min(axis=0)
self.Xmax, self.Ymax = self.points.max(axis=0)
self.C = int(self.get_model_value('resolution'))
self.W = int(np.ceil((self.Xmax - self.Xmin) / self.C))
self.H = int(np.ceil((self.Ymax - self.Ymin) / self.C))
self.herds['tile'] = self.W * (np.floor((self.herds['Y'] - self.Ymin) / self.C)).astype(int) +\
(np.floor((self.herds['X'] - self.Xmin) / self.C)).astype(int)
self.wildboars['tile'] = self.W * (np.floor((self.wildboars['Y'] - self.Ymin) / self.C)).astype(int) +\
(np.floor((self.wildboars['X'] - self.Xmin) / self.C)).astype(int)
tiles = set(self.herds['tile'])
tiles |= set(self.wildboars['tile'])
print(len(tiles), "distinct tiles in simulation environment")
self.tiles_info = pd.DataFrame({'tile': self.wildboars['tile'].unique()})
self.tiles_info['active'] = 0
self.tiles_info.set_index(['tile'], inplace=True)
# initial kernels (0)
self.kernel_other_pathways = sp.csr_matrix((self.nb_herds, self.nb_herds))
self.kernel_boar_static = sp.csr_matrix((self.nb_herds+self.nb_wildboars, self.nb_herds+self.nb_wildboars))
self.kernel_boar_to_boar = sp.csr_matrix((self.nb_wildboars, self.nb_wildboars))
# initial matrix to find neighbours within a given radius
self.neighbours = sp.lil_matrix((self.nb_herds, self.nb_herds), dtype='int32')
# date where suspicion due to trade contacts ends (population_index -> date to end suspicion)
self.traced_end = sp.lil_matrix((1, self.nb_herds))
# date where protection area ends (population_index -> date to end suspicion)
self.protection_end = sp.lil_matrix((1, self.nb_herds))
# date where surveillance area ends (population_index -> date to end suspicion)
self.surveillance_end = sp.lil_matrix((1, self.nb_herds))
# ## reset hunting bag if first run
# self.hunting_bag_name= Path(self.log_path()).with_name('hunting_bag.txt')
# if self.hunting_bag_name.exists() and self.statevars.simu_id==0:
# self.hunting_bag_name.unlink()
## reset log of removed trade movements if first run
self.trade_log_name= Path(self.log_path()).with_name('trade_log.txt')
if self.trade_log_name.exists() and self.statevars.simu_id==0:
self.trade_log_name.unlink()
## reset log of herd status (P protection, S surveillance, T traced) if first run
self.status_log_name= Path(self.log_path()).with_name('herd_status.txt')
if self.status_log_name.exists() and self.statevars.simu_id==0:
self.status_log_name.unlink()
## reset detailed WB log if first run
self.wblog_name= Path(self.log_path()).with_name('wildboars.txt')
if self.wblog_name.exists() and self.statevars.simu_id==0:
self.wblog_name.unlink()
## initialize schedules for active wildboar carcass search
self.search_schedule={} # step -> list of population_ID of centers of research areas
self.wb_points = self.points[self.nb_herds:, :]
def tile(self, col, row):
"""Return the number of the tile at column *col* and row *row*."""
return col + self.W * row
def coords(self, tile):
"""Return the coordinates (column, row) of the specified tile
number.
"""
return (tile % self.W, tile // self.W)
def is_adjacent(self, tile1, tile2):
"""Test if *tile1* and *tile2* are adjacent."""
col1, row1 = self.coords(tile1)
col2, row2 = self.coords(tile2)
return max(abs(col1 - col2), abs(row1 - row2)) == 1
def customize_pigherd_prototype(self, pop_id, nb_exposed):
proto = self.model.get_prototype('epid_unit', 'pig_herd')
proto['initial_exposed'] = nb_exposed
proto['is_infected'] = 1
proto['is_traced'] = self.traced_pig_herds[0, pop_id-1]
proto['is_in_protection_area'] = self.protection_pig_herds[0, pop_id-1]
proto['is_in_surveillance_area'] = self.surveillance_pig_herds[0, pop_id-1]
proto['date_traced_end'] = self.traced_end[0, pop_id-1]
proto['date_protection_end'] = self.protection_end[0, pop_id-1]
proto['date_surveillance_end'] = self.surveillance_end[0, pop_id-1]
proto['changed'] = 1
herd = self.herds[self.herds['population_id'] == pop_id]
for field in ['population_id', 'X', 'Y', 'size', 'biosecurity', 'is_outdoor', 'is_commercial', 'tile']:
proto[field] = float(herd[field])
return proto
def customize_wildboar_prototype(self, pop_id):
proto = self.model.get_prototype('epid_unit', 'wild_boar')
proto['initial_exposed'] = 1
proto['is_infected'] = 1
proto['size'] = 1
proto['changed'] = 1
boars = self.wildboars[self.wildboars['population_id'] == pop_id]
for field in ['population_id', 'X', 'Y', 'tile']:
proto[field] = float(boars[field])
return proto
def create_initial_epidunits(self):
# list of active tiles (loaded)
self.tiles = set()
# list of tiles for which the pig herd dist matrix has been loaded
self.neigh_tiles = set()
## INIT VECTORS AND OTHER VARIABLES USED TO REPRESENT RELEVANT HEALTH STATES/PREVALENCES
self.statevars.initial_nb_boars_alive = self.get_model_value('nb_wildboars')
self.reset_status()
self.contribution_boar_to_pig = sp.csr_matrix((1, self.nb_herds))
self.contribution_other_pathways = sp.csr_matrix((1, self.nb_herds))
self.contribution_boar_to_boar = sp.csr_matrix((1, self.nb_wildboars))
self.contribution_pig_to_boar = sp.csr_matrix((1, self.nb_wildboars))
self.detected_pig_herds = sp.lil_matrix((1, self.nb_herds), dtype='int32')
self.protection_pig_herds = sp.lil_matrix((1, self.nb_herds), dtype='int32')
self.surveillance_pig_herds = sp.lil_matrix((1, self.nb_herds), dtype='int32')
self.traced_pig_herds = sp.lil_matrix((1, self.nb_herds), dtype='int32')
# CREATE INITIALLY INFECTED EPID UNITS
nb_herds = int(self.get_model_value('nb_infected_pig_herds'))
new_herds = [] if nb_herds == 0 else\
[(int(self.get_model_value('infected_pig_herd_ID')), 1)]
nb_boars = int(self.get_model_value('nb_infected_wild_boars'))
new_boars = [] if nb_boars == 0 else\
[int(self.get_model_value('infected_boar_ID'))]
self.add_new_epidunits(new_herds, new_boars)
# print('Calculating initial contributions')
# self.compute_contributions_from_kernels()
def add_new_epidunits(self, herd_list, boar_list):
new_epid_units = [self.new_atom(custom_prototype=self.customize_pigherd_prototype(herd_id, nb_E), execute_actions=False)
for (herd_id, nb_E) in herd_list]
new_epid_units += [self.new_atom(custom_prototype=self.customize_wildboar_prototype(boar_id), execute_actions=False)
for boar_id in boar_list]
self.add_atoms(new_epid_units)
new_tiles = set()
for agent in new_epid_units:
agent.reapply_prototype(execute_actions=True)
# agent.log_custom_vars()
new_tiles.add(agent.statevars.tile)
tiles_to_add = new_tiles - self.tiles
for tile in tiles_to_add:
self.load_tile(tile)
def load_tile(self, new_tile):
# load kernels and integrate them to existing data
# print('Loading tile', new_tile)
# OTHER PATHWAYS (between herds through indirect contacts)
print('Loading new tile', new_tile)
alpha = self.get_model_value('alpha_other_pathways')
kop_file = self.model.input_dir.joinpath('tiles_{}'.format(self.C),
'kernel_other_{}_{}.npz'.format(int(alpha), int(new_tile)))
# print('Processing', kop_file)
self.kernel_other_pathways += sp.load_npz(kop_file)
# new_kernel = sp.load_npz(kop_file).tolil()
# idx = new_kernel.nonzero()
# self.kernel_other_pathways[idx] = new_kernel[idx]
# BOAR-PIG HERD RELATIONSHIPS
alpha = self.get_model_value('alpha_boar_static')
kbs_file = self.model.input_dir.joinpath('tiles_{}'.format(self.C),
'kernel_static_{}_{}.npz'.format(int(alpha), int(new_tile)))
# print('Processing', kbs_file)
self.kernel_boar_static += sp.load_npz(kbs_file)
# new_kernel = sp.load_npz(kbs_file).tolil()
# idx = new_kernel.nonzero()
# self.kernel_boar_static[idx] = new_kernel[idx]
# BOAR-BOAR RELATIONSHIPS
alpha = self.get_model_value('alpha_boar_boar')
kbb_file = self.model.input_dir.joinpath('tiles_{}'.format(self.C),
'kernel_boars_{}_{}.npz'.format(int(alpha), int(new_tile)))
# print('Processing', kbb_file)
self.kernel_boar_to_boar += sp.load_npz(kbb_file)
# new_kernel = sp.load_npz(kbb_file).tolil()
# idx = new_kernel.nonzero()
# self.kernel_boar_to_boar[idx] = new_kernel[idx]
self.tiles.add(new_tile)
self.tiles_info['active'][new_tile] = 1
# if fences have been installed, if new tile is the neighbour of a sourrounded tile, modify the corresponding kernel
if self.get_model_value('use_fences') and self.statevars.date_first_detection > 0 and self.time >= self.get_model_value('delay_use_fences') + self.statevars.date_first_detection:
if new_tile in self.surrounding_tiles:
print('Adapting kernel for new tile in buffer area', new_tile)
self.kernel_update_fences([new_tile], self.surrounded_tiles)
if new_tile in self.surrounded_tiles:
print('Adapting kernel for new tile within fences', new_tile)
self.kernel_update_fences([new_tile], self.surrounding_tiles)
def load_neighbours(self, new_tile):
print('Loading neighbourhood for tile', new_tile)
if new_tile not in self.neigh_tiles:
#### LOAD neighbours only when needed
# NEIGHBOURHOOD RELATIONSHIPS
neighb_file = self.model.input_dir.joinpath('tiles_{}'.format(self.C),
'neighbours_{}.npz'.format(int(new_tile)))
print('Processing', neighb_file)
#### remember that matrix is triangular in the part concerning
#### each tile (we only need to find who is at distance D of
#### another epidunit => look at one row/col in matrix without
#### making the matrix symmetric)
self.neighbours += sp.load_npz(neighb_file)
self.neigh_tiles.add(new_tile)
def init_trade(self):
## TRADE
# read a CSV data file for moves:
# date of movement, source herd, destination herd, quantity
# and restructure it according to origin_date and delta_t:
# {step: {source_id: [(dest_id, qty), ...],
# ...},
# ...}
origin = self.model.origin_date
step_duration = self.model.step_duration
moves = {}
move_file = Path(self.model.input_dir, 'moves.csv')
print("Processing", move_file)
with open(move_file) as csvfile:
# read the CSV file
csvreader = csv.DictReader(csvfile, delimiter=',')
for row in csvreader:
# read date and incorporate the offset between user
# time (0=beginning of hunting period) and simulation
# time (0 = a few days before)
day = int(row['date'])
# if day < 0:
# # ignore dates before simulation
# continue
# convert dates into simulation steps
# step = (day - origin) // step_duration
step = day // step_duration.days
# group information by step and source herd
if step not in moves:
moves[step] = {}
src, dest, qty = int(row['source']), int(row['dest']), int(row['qty'])
if src != dest:
if src not in moves[step]:
moves[step][src] = []
moves[step][src].append([dest, qty])
self.moves = moves
def log_active_search(self, cases, healthy):
"""Store custom information for WB carcasses found in active search.
*cases* and *healthy* represent indexes of respectively infected and healthy carcasses
"""
## logged variables:
## simu_id, step, population_id, is_infected, is_detected, date_detected, is_traced, date_traced_end, is_in_protection_area, date_protection_end, is_in_surveillance_area, date_surveillance_end, is_hunted, is_tested, S, E, I, C, D
dur_search = int(self.get_model_value('duration_carcass_active_search'))
message = ''
values = [
[ self.statevars.simu_id,
self.statevars.step,
pop_idx + self.nb_herds + 1,
1, ## is_infected,
0, ## is_detected,
1, ## is_searched,
0, ## is_culled
self.statevars.step - np.random.random_integers(0, dur_search), ## date_detected,
0, ## is_traced,
0, ## date_traced_end,
0, ## is_in_protection_area,
0, ## date_protection_end,
0, ## is_in_surveillance_area,
0, ## date_surveillance_end,
0, ## is_hunted,
0, ## is_tested,
0, ## total_S,
0, ## total_E,
0, ## total_I,
1, ## total_C,
0, ## total_D,
] for pop_idx in cases
]
values += [
[ self.statevars.simu_id,
self.statevars.step,
pop_idx + self.nb_herds + 1,
0, ## is_infected,
0, ## is_detected,
1, ## is_searched,
0, ## is_culled
self.statevars.step - np.random.random_integers(0, dur_search), ## date_detected,
0, ## is_traced,
0, ## date_traced_end,
0, ## is_in_protection_area,
0, ## date_protection_end,
0, ## is_in_surveillance_area,
0, ## date_surveillance_end,
0, ## is_hunted,
0, ## is_tested,
0, ## total_S,
0, ## total_E,
0, ## total_I,
0, ## total_C,
1, ## total_D,
] for pop_idx in healthy
]
if any(values):
nb_fields = len(values[0])
message = '\n'.join((','.join(['{}'] * nb_fields)).format(*row) for row in values)
with open(self.log_path(), 'a') as logfile:
logfile.write('{}\n'.format(message))
def log_preventive_culling(self, healthy_idx):
"""Store custom information for healthy pig herds (list of indices)
culled preventively (infected pig herds culled preventively
even when not detected are in charge of logging their own
variables since they are running in the simulation).
"""
## logged variables:
## simu_id, step, population_id, is_infected, is_detected, date_detected, is_traced, date_traced_end, is_in_protection_area, date_protection_end, is_in_surveillance_area, date_surveillance_end, is_hunted, is_tested, S, E, I, C, D
message = ''
values = [
[ self.statevars.simu_id,
self.statevars.step,
idx+1, # population_id
0, ## is_infected,
0, ## is_detected,
0, ## is_searched,
1, ## is_culled,
self.statevars.step, ## date culled,
self.traced_pig_herds[0, idx], ## is_traced,
self.traced_end[0, idx], ## date_traced_end,
self.protection_pig_herds[0, idx], ## is_in_protection_area,
self.protection_end[0, idx], ## date_protection_end,
self.surveillance_pig_herds[0, idx], ## is_in_surveillance_area,
self.surveillance_end[0, idx], ## date_surveillance_end,
0, ## is_hunted,
0, ## is_tested,
0, ## total_S,
0, ## total_E,
0, ## total_I,
0, ## total_C,
0, ## total_D,
] for idx in healthy_idx
]
if any(values):
nb_fields = len(values[0])
message = '\n'.join((','.join(['{}'] * nb_fields)).format(*row) for row in values)
with open(self.log_path(), 'a') as logfile:
logfile.write('{}\n'.format(message))
def log_hunted_wildboars(self, healthy_idx):
"""Store custom information for tested hunted wildboars (list of
indices).
"""
## logged variables:
## simu_id, step, population_id, is_infected, is_detected, date_detected, is_traced, date_traced_end, is_in_protection_area, date_protection_end, is_in_surveillance_area, date_surveillance_end, is_hunted, is_tested, S, E, I, C, D
nb_PH = self.nb_herds
message = ''
values = [
[ self.statevars.simu_id,
self.statevars.step,
idx + nb_PH + 1, # population_id
0, ## is_infected,
0, ## is_detected,
0, ## is_searched,
0, ## is_culled,
self.statevars.step, ## date hunted,
0, ## is_traced,
0, ## date_traced_end,
0, ## is_in_protection_area,
0, ## date_protection_end,
0, ## is_in_surveillance_area,
0, ## date_surveillance_end,
1, ## is_hunted,
1, ## is_tested,
0, ## total_S,
0, ## total_E,
0, ## total_I,
0, ## total_C,
0, ## total_D,
] for idx in healthy_idx
]
if any(values):
nb_fields = len(values[0])
message = '\n'.join((','.join(['{}'] * nb_fields)).format(*row) for row in values)
with open(self.log_path(), 'a') as logfile:
logfile.write('{}\n'.format(message))
#----------------------------------------------------------------
# Statevars
#----------------------------------------------------------------
@property
def nb_boars_alive(self):
return self.statevars.initial_nb_boars_alive - self.statevars.nb_deceased_wildboars
#----------------------------------------------------------------
# Processes
#----------------------------------------------------------------
def change_seed(self):
""" Use a new seed for the rest of the simulation
"""
if self.get_model_value('use_second_seed') and self.statevars.date_first_detection > 0 and self.time == self.get_model_value('delay_change_seed') + self.statevars.date_first_detection:
np.random.seed(int(self.get_model_value('second_seed')))
def setup_fences(self):
"""Identify location of future physical barriers between tiles specified in data (fences.csv and buffer.csv).
"""
if self.get_model_value('use_fences'):
# read explicit list of tiles to enclose
fences_file = self.model.input_dir.joinpath('fences.csv')
self.surrounded_tiles = set(pd.read_csv(fences_file, sep=',', header=0)['tile'])
buffer_file = self.model.input_dir.joinpath('buffer.csv')
self.surrounding_tiles = set(pd.read_csv(buffer_file, sep=',', header=0)['tile'])
# identify indexes of wildboars within tiles
self.wildboars_within_fences = self.wildboars[self.wildboars['tile'].isin(self.surrounded_tiles)].index
self.wildboars_within_buffer = self.wildboars[self.wildboars['tile'].isin(self.surrounding_tiles)].index
def install_fences(self):
"""Install physical barriers between predefined surrounded/surrounding tiles. Kernel values can only be modified for surrounded/surrounding tiles already active. When new tiles are loaded due to infection spread, kernels must be modified consequently.
"""
### TODO: prevent installing fences while not yet detected (temporarily allowed for test purpose)
if self.get_model_value('use_fences') and self.statevars.date_first_detection > 0 and self.time == self.get_model_value('delay_use_fences') + self.statevars.date_first_detection: # and self.statevars.date_first_detection > 0
nb_active_tiles = len(self.surrounded_tiles)
area = nb_active_tiles * (self.get_model_value('resolution') / 1000)**2
print('INSTALLING FENCES at t =', self.statevars.step, 'for', nb_active_tiles, 'tiles =', area, 'km2')
print('SURROUNDED TILES:', self.surrounded_tiles)
print('SURROUNDING TILES:', self.surrounding_tiles)
self.kernel_update_fences(self.surrounded_tiles, self.surrounding_tiles)
self.kernel_update_fences(self.surrounding_tiles, self.surrounded_tiles)
self.statevars.changed = True
def check_if_finished(self):
"""Stop simulation if time >= date detection + game duration (enough data
for players) if a game duration > 0 is set.
"""
game_duration = self.get_model_value('game_duration')
if game_duration > 0 and self.statevars.date_first_detection > 0 and self.time >= self.statevars.date_first_detection + game_duration:
print('HALTING SIMULATION AT', self.time)
self.deactivate()
def kernel_update_fences(self, from_tiles, to_tiles):
"""Modify boar_static and boar_to_boar kernels between tiles contained
in *from_tiles* and tiles contained in *to_tiles*, to account
for newly installed or existing fences. The operation is not
symmetric. Two first calls (symmetrical) of this method must be done when
installing fences, then when a new tile is loaded, if it is in the list
of surrounded (resp. surrounding) tiles, the method must be
called with from_tiles = new tile and to_tiles = surrounding (resp.
surrounded) tiles.
"""
# retrieve parameter value
efficacy = self.get_model_value('fences_efficacy')
## iteration on source tiles to limit memory footprint !
for source_tile in from_tiles:
## skip non-active source tiles
if not self.tiles_info['active'][source_tile]:
continue
## skip non-adjacent tiles since kernels between
## non-adjacent tiles is already 0
adjacent_tiles = set(dest_tile for dest_tile in to_tiles