Added data used by the model and scripts to generate them

model/data/buffer.csv

+tile
+542
+543
+544
+545
+546
+547
+582
+587
+622
+627
+662
+667
+702
+707
+742
+747
+782
+787
+822
+823
+824
+825
+826
+827

model/data/create_tiles.sh

+#!/bin/bash
+for tile in `cat list_of_tiles.txt` 
+do
+	python extract_area.py $tile
+done

model/data/extract_area.py

+## Usage: python extract_area.py TILE [RESOLUTION THRESHOLD ALPHA_BOAR_BOAR ALPHA_BOAR_STATIC ALPHA_OTHER_PATHWAYS N1 N2...]
+
+## assuming that RESOLUTION >= max relevant distances
+## i.e. >= Ni (neighbourhood distances)
+## and >= alpha * sqrt(ln(1/epsilon)) for all alphas (cut-off distance)
+## with this assumption, works fine with 8 Moore neighbours
+## otherwise revise how adjacent tiles are retrieved but careful with computational cost
+
+import itertools as it
+import sys
+import time
+import numpy as np
+import pandas as pd
+import scipy.sparse as sp
+from scipy.spatial import distance_matrix
+
+TARGET_TILE = int(sys.argv[1])
+if len(sys.argv) > 2:
+    C = int(sys.argv[2])
+    epsilon = float(sys.argv[3])
+    alpha_boar_boar = int(sys.argv[4])
+    alpha_boar_static = int(sys.argv[5])
+    alpha_other_pathways = int(sys.argv[6])
+    neighbourhoods = [ int(sys.argv[i]) for i in range(7, len(sys.argv))]
+else:
+    C = 15000
+    alpha_boar_boar = 2000
+    alpha_boar_static = 1000
+    alpha_other_pathways = 1000
+    epsilon = 1e-9
+    neighbourhoods = [3000, 10000, 15000]
+
+
+print("Extracting tile", TARGET_TILE, "with resolution =", C, "m and pruning threshold ", epsilon)
+print("Alpha values:", alpha_boar_boar, alpha_boar_static, alpha_boar_static)
+print("Neighbour distances (m):", *neighbourhoods)
+
+start = time.perf_counter()
+print('Loading pig herd info')
+# load pig herds info
+units = pd.read_csv('herds.csv', sep=',', header=0).drop(columns=['biosecurity','is_commercial','is_outdoor','multisite','production','size'])
+nb_pigherds = units.shape[0]
+
+print('Loading wildboar info')
+# load wildboars info
+units = pd.concat([units, pd.read_csv('boars.csv', sep=',', header=0).drop(columns=['Zone'])],
+                  ignore_index = True)
+units = units.astype('int32')
+nb_wildboars = units.shape[0] - nb_pigherds
+
+print('Filtering points')
+# prepare the list of points
+points = units[['X', 'Y']].values
+
+Xmin, Xmax = units['X'].min(), units['X'].max()
+Ymin, Ymax = units['Y'].min(), units['Y'].max()
+
+W = int(np.ceil((Xmax - Xmin) / C))
+H = int(np.ceil((Ymax - Ymin) / C))
+
+units['row'] = (np.floor((units['Y'] - Ymin) / C)).astype(int)
+units['col'] = (np.floor((units['X'] - Xmin) / C)).astype(int)
+units['tile'] = units['col'] + W * units['row']
+
+print('Loading target tiles')
+COL = TARGET_TILE % W
+ROW = TARGET_TILE // W
+
+ADJACENT_TILES = [(COL + i) + (ROW + j) * W
+                  for j in [-1, 0, 1]
+                  for i in [-1, 0, 1]
+                  if (i != 0 or j != 0)]
+
+tile_idx,  = (units['tile'] == TARGET_TILE).to_numpy().nonzero()
+ph_idx = tile_idx[tile_idx < nb_pigherds]
+wb_idx = tile_idx[tile_idx >= nb_pigherds]
+nb_ph_tile, = ph_idx.shape
+
+tile = points[tile_idx, :]
+tile_ph = points[ph_idx, :]
+tile_wb = points[wb_idx, :]
+
+adjacent_idx = np.concatenate([tile_idx] +\
+                              [(units['tile'] == other).to_numpy().nonzero()[0]
+                               for other in ADJACENT_TILES])
+nph_idx = adjacent_idx[adjacent_idx < nb_pigherds]
+nwb_idx = adjacent_idx[adjacent_idx >= nb_pigherds]
+nb_ph_adj, = nph_idx.shape
+
+neighbour_tiles = points[adjacent_idx, :]
+n_tiles_ph = points[nph_idx, :]
+n_tiles_wb = points[nwb_idx, :]
+
+### WARNING: slow operation !!! (about 3.10^6 coordinates)
+# coords = tuple(zip(*it.product(tile_idx, adjacent_idx)))
+irow, icol = np.meshgrid(tile_idx, adjacent_idx)
+irow, icol = irow.T, icol.T
+
+irowph, icolph = np.meshgrid(ph_idx, nph_idx)
+irowph, icolph = irowph.T, icolph.T
+
+irowwb, icolwb = np.meshgrid(wb_idx-nb_pigherds, nwb_idx-nb_pigherds)
+irowwb, icolwb = irowwb.T, icolwb.T
+
+nrowph, ncolph = np.meshgrid(ph_idx, nph_idx)
+nrowph, ncolph = nrowph.T, ncolph.T
+
+print('Computing kernel other_pathways (PH/PH)')
+# Kernel for herd-herd indirect transmission
+d_tile = distance_matrix(tile_ph, n_tiles_ph)
+kernel = np.exp(- (d_tile / alpha_other_pathways)**2)
+kernel[kernel < epsilon] = 0
+np.fill_diagonal(kernel, 0)
+kernel_other_whole = sp.lil_matrix((nb_pigherds,nb_pigherds))
+kernel_other_whole[irowph, icolph] = kernel
+
+print('Computing kernel boar-to-boar (WB/WB)')
+# Kernel for direct transmission between I and S alive boars
+d_tile = distance_matrix(tile_wb, n_tiles_wb)
+kernel = np.exp(- (d_tile / alpha_boar_boar)**2)
+kernel[kernel < epsilon] = 0
+np.fill_diagonal(kernel, 0)
+kernel_boars_whole = sp.lil_matrix((nb_wildboars,nb_wildboars))
+kernel_boars_whole[irowwb, icolwb] = kernel
+
+print('Computing kernel boar-static(PH+WB/PH+WB)')
+# Kernel boar-to-static epid unit (I/S wild boar <--> boar carcass or pig herd)
+d_tile = distance_matrix(tile, neighbour_tiles)
+kernel = np.exp(- (d_tile / alpha_boar_static)**2)
+kernel[kernel < epsilon] = 0
+np.fill_diagonal(kernel, 0)
+kernel_static_whole = sp.lil_matrix((nb_pigherds+nb_wildboars,nb_pigherds+nb_wildboars))
+kernel_static_whole[irow, icol] = kernel
+
+
+print('Computing neighbours')
+# Distance matrix between pig herds on tile and pig herds in neighbourhood
+d_tile = distance_matrix(tile_ph, n_tiles_ph)
+neighbours = np.zeros(d_tile.shape, dtype='int32')
+for dist in sorted(neighbourhoods, reverse=True):
+    neighbours[d_tile <= dist] = dist
+np.fill_diagonal(neighbours, 0)
+neighbours_whole = sp.lil_matrix((nb_pigherds,nb_pigherds), dtype='int32')
+neighbours_whole[nrowph, ncolph] = neighbours
+
+
+print('Saving')
+sp.save_npz('tiles_{}/kernel_boars_{}_{}'.format(C, alpha_boar_boar, TARGET_TILE), kernel_boars_whole.tocsr())
+sp.save_npz('tiles_{}/kernel_static_{}_{}'.format(C, alpha_boar_static, TARGET_TILE), kernel_static_whole.tocsr())
+sp.save_npz('tiles_{}/kernel_other_{}_{}'.format(C, alpha_other_pathways, TARGET_TILE), kernel_other_whole.tocsr())
+
+sp.save_npz('tiles_{}/neighbours_{}'.format(C, TARGET_TILE), neighbours_whole.tocsr())
+
+stop = time.perf_counter()
+print('Computation finished in {:.2f} s'.format(stop-start))

model/data/fences.csv

+tile
+583
+584
+585
+586
+623
+624
+625
+626
+663
+664
+665
+666
+703
+704
+705
+706
+743
+744
+745
+746
+783
+784
+785
+786

model/data/list_of_tiles.txt

+11
+12
+13
+14
+15
+16
+17
+18
+19
+51
+52
+53
+54
+55
+56
+57
+58
+81
+83
+84
+85
+87
+88
+89
+90
+91
+92
+93
+94
+95
+96
+97
+98
+120
+121
+122
+123
+124
+125
+126
+127
+128
+129
+130
+131
+132
+133
+134
+135
+136
+137
+138
+160
+161
+162
+163
+164
+165
+166
+167
+168
+169
+170
+171
+172
+173
+174
+175
+176
+177
+178
+200
+201
+202
+203
+204
+205
+206
+207
+208
+209
+210
+211
+212
+213
+214
+215
+216
+217
+218
+240
+241
+242
+243
+244
+245
+246
+247
+248
+249
+250
+251
+252
+253
+254
+255
+256
+257
+258
+259
+260
+281
+282
+283
+284
+285
+286
+287
+288
+289
+290
+291
+292
+293
+294
+295
+296
+297
+298
+299
+300
+301
+302
+321
+322
+323
+324
+325
+326
+327
+328
+329
+330
+331
+332
+333
+334
+335
+336
+337
+338
+339
+340
+341
+342
+343
+344
+345
+360
+361
+362
+363
+364
+365
+366
+367
+368
+369
+370
+371
+372
+373
+374
+375
+376
+377
+378
+379
+380
+381
+382
+383
+384
+385
+386
+400
+401
+402
+403
+404
+405
+406
+407
+408
+409
+410
+411
+412
+413
+414
+415
+416
+417
+418
+419
+420
+421
+422
+423
+424
+425
+426
+440
+441
+442
+443
+444
+445
+446
+447
+448
+449
+450
+451
+452
+453
+454
+455
+456
+457
+458
+459
+460
+461
+462
+463
+464
+465
+466
+467
+480
+481
+482
+483
+484
+485
+486
+487
+488
+489
+490
+491
+492
+493
+494
+495
+496
+497
+498
+499
+500
+501
+502
+503
+504
+505
+506
+507
+525
+526
+527
+528
+529
+530
+531
+532
+533
+534
+535
+536
+537
+538
+539
+540
+541
+542
+543
+544
+545
+546
+547
+550
+551
+552
+566
+567
+568
+569
+570
+571
+572
+573
+574
+575
+576
+577
+578
+579
+580
+581
+582
+583
+584
+585
+586
+587
+588
+589
+590
+591
+592
+606
+607
+608
+609
+610
+611
+612
+613
+614
+615
+616
+617
+618
+619
+620
+621
+622
+623
+624
+625
+626
+627
+628
+629
+630
+631
+647
+648
+649
+650
+651
+652
+653
+654
+655
+656
+657
+658
+659
+660
+661
+662
+663
+664
+665
+666
+667
+668
+669
+670
+671
+672
+688
+689
+690
+691
+692
+693
+694
+695
+696
+697
+698
+699
+700
+701
+702
+703
+704
+705
+706
+707
+708
+709
+710
+711
+712
+729
+730
+731
+732
+733
+734
+735
+736
+737
+738
+739
+740
+741
+742
+743
+744
+745
+746
+747
+748
+749
+750
+751
+752
+753