In [1]:
## Contact networks and simple viral spreading.
In [2]:
##
using Random
Random.seed!(0) # ensure repeatable results...
Out[2]:
MersenneTwister(UInt32[0x00000000], Random.DSFMT.DSFMT_state(Int32[748398797, 1073523691, -1738140313, 1073664641, -1492392947, 1073490074, -1625281839, 1073254801, 1875112882, 1073717145  …  943540191, 1073626624, 1091647724, 1073372234, -1273625233, -823628301, 835224507, 991807863, 382, 0]), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], UInt128[0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000  …  0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000, 0x00000000000000000000000000000000], 1002, 0)
In [3]:
## Generate a simple spatial graph model to look at.\
using NearestNeighbors, Distributions, SparseArrays
function spatial_graph_edges(n::Integer,d::Integer;degreedist=LogNormal(log(4),1))
  xy = rand(d,n)
  T = BallTree(xy)
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    deg = min(ceil(Int,rand(degreedist)),n-1)
    idxs, dists = knn(T, xy[:,i], deg+1)
    for j in idxs
      if i != j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  return xy, ei, ej
end
function spatial_network(n::Integer, d::Integer; degreedist=LogNormal(log(3),1))
  xy, ei, ej = spatial_graph_edges(n, d;degreedist=degreedist)
  A = sparse(ei,ej,1,n,n)
  return max.(A,A'), xy
end
A,xy = spatial_network(1000, 2)
Out[3]:
(
  [16 ,    1]  =  1
  [192,    1]  =  1
  [382,    1]  =  1
  [436,    1]  =  1
  [511,    1]  =  1
  [576,    1]  =  1
  [598,    1]  =  1
  [838,    1]  =  1
  [145,    2]  =  1
  [227,    2]  =  1
  [246,    2]  =  1
  [497,    2]  =  1
  â‹®
  [500,  999]  =  1
  [522,  999]  =  1
  [624,  999]  =  1
  [713,  999]  =  1
  [79 , 1000]  =  1
  [129, 1000]  =  1
  [194, 1000]  =  1
  [279, 1000]  =  1
  [289, 1000]  =  1
  [417, 1000]  =  1
  [486, 1000]  =  1
  [652, 1000]  =  1
  [933, 1000]  =  1, [0.8236475079774124 0.16456579813368521 … 0.22375320653445185 0.41147647589610803; 0.9103565379264364 0.17732884646626457 … 0.3371041633752143 0.6063082992397912])
In [4]:
##
using Plots
function plotgraph(A::SparseMatrixCSC,xy::AbstractArray{T,2};kwargs...) where T
  px,py = zeros(T,0),zeros(T,0)
  P = [px,py]
  rows = rowvals(A)
  skip = NaN.*xy[:,begin] # first row
  for j=1:size(A,2) # for each column
    for nzi in nzrange(A, j)
      i = rows[nzi]
      if i > j
        push!.(P, @view xy[:,i])
        push!.(P, @view xy[:,j])
        push!.(P, skip)
      end
    end
  end
  plot(px,py;framestyle=:none,legend=false,kwargs...)
end
plotgraph(A,xy,alpha=0.25); scatter!(xy[1,:],xy[2,:],
  markersize=2, markerstrokewidth=0, color=1)
Out[4]:
In [5]:
## Run a contact simulation
p = 0.2 # probability of infection

# The result we use is that your probability of infection is
# 1 - probability of not being infected
# and your probability of not being infected is that you
# didn't get infected from any neighbor.
# https://www.khanacademy.org/math/ap-statistics/probability-ap/probability-multiplication-rule/a/probabilities-involving-at-least-one-success
function evolve(x::Vector,  p::Real, A::AbstractMatrix)
  log_not_infected = log.(1 .- p.*x)
  y = 1 .- exp.(A*log_not_infected)
  y = max.(y, x)
end

x = zeros(size(A,1))
x[2] = 1.0
anim = @animate for i=1:100
  global x = evolve(x, p, A)
  plotgraph(A,xy,alpha=0.25, size=(600,600))
  scatter!(xy[1,:],xy[2,:],
    markersize=5, markerstrokewidth=0,
    color=1, marker_z = x, clim=(0,1))
  title!("$i")
end
gif(anim, "viral-wrong.gif", fps=10)
┌ Info: Saved animation to 
│   fn = /Users/dgleich/Dropbox/courses/cs515-2020/web/input/julia/viral-wrong.gif
â”” @ Plots /Users/dgleich/.julia/packages/Plots/ViMfq/src/animation.jl:104
Out[5]:
In [6]:
## Okay, fix one bug, you need to scale by the probability
# you weren't infected at the previous step.

# 1 - probability of not being infected
# and your probability of not being infected is that you
# didn't get infected from any neighbor.
# https://www.khanacademy.org/math/ap-statistics/probability-ap/probability-multiplication-rule/a/probabilities-involving-at-least-one-success
function evolve_self(x::Vector,  p::Real, A::AbstractMatrix)
  log_not_infected = log.(1 .- p.*x)
  y = (1 .- exp.(A*log_not_infected).*(1 .- x))
  y = max.(y, x)
end
p = 0.2
x = zeros(size(A,1))
x[2] = 1.0
anim = @animate for i=1:100
  global x = evolve_self(x, p, A)
  plotgraph(A,xy,alpha=0.25, size=(600,600))
  scatter!(xy[1,:],xy[2,:],
    markersize=5, markerstrokewidth=0,
    color=1, marker_z = x, clim=(0,1))
  title!("$i")
end
gif(anim, "viral-self.gif", fps=10)
┌ Info: Saved animation to 
│   fn = /Users/dgleich/Dropbox/courses/cs515-2020/web/input/julia/viral-self.gif
â”” @ Plots /Users/dgleich/.julia/packages/Plots/ViMfq/src/animation.jl:104
Out[6]:
In [7]:
##
function evolve_self(x::Vector,  p::Real, A::AbstractMatrix)
  log_not_infected = log.(1 .- p.*x)
  y = (1 .- exp.(A*log_not_infected).*(1 .- x))
  y = max.(y, x)
end
p = 0.01
x = zeros(size(A,1))
x[2] = 1.0
anim = @animate for i=1:500
  global x = evolve_self(x, p, A)
  plotgraph(A,xy,alpha=0.25, size=(600,600))
  scatter!(xy[1,:],xy[2,:],
    markersize=5, markerstrokewidth=0,
    color=1, marker_z = x, clim=(0,1))
  title!("$i")
end
gif(anim, "viral-self-slow-long.gif", fps=10)
┌ Info: Saved animation to 
│   fn = /Users/dgleich/Dropbox/courses/cs515-2020/web/input/julia/viral-self-slow-long.gif
â”” @ Plots /Users/dgleich/.julia/packages/Plots/ViMfq/src/animation.jl:104
Out[7]: