Basic R count

src/default_vector.jl

@@ -0,0 +1,31 @@
+#=
+deafult_vector:
+- Julia version: 
+- Author: jacek
+- Date: 2022-03-18
+=#
+
+mutable struct DefaultVector{T}
+    defaultValue::T
+    core::Vector{T}
+
+    DefaultVector(defaultValue_) = new{typeof(defaultValue_)}(defaultValue_, Vector{typeof(defaultValue_)}(undef, 0))
+    DefaultVector(defaultValue_, K) = new{K}(defaultValue_, Vector{K}(undef, 0))
+end
+
+Base.setindex!(self::DefaultVector{T}, value::T, index::Int64) where T = _ensure_index_in_range!(self, index) do
+    self.core[index] = value
+end
+
+Base.getindex(self::DefaultVector{T}, index::Int64) where T = _ensure_index_in_range!(self, index) do
+    self.core[index]
+end
+
+Base.convert(::Type{Vector}, self::DefaultVector{T}) where T = copy(self.core)
+
+function _ensure_index_in_range!(next_func::Function, self::DefaultVector{T}, index::Int64) where T
+    while index > length(self.core)
+        push!(self.core, self.defaultValue)
+    end
+    next_func()
+end

src/event.jl

@@ -57,6 +57,7 @@ struct Event
   Event(::Val{OutsideInfectionEvent}, time::Real, subject::Integer, strain::StrainKind) = new(time, subject, 0, OutsideInfectionEvent, UInt8(OutsideContact), strain)
   Event(::Val{TransmissionEvent}, ::Real, ::Integer) = error("source and contact kind needed for transmission event")
   Event(::Val{TransmissionEvent}, time::Real, subject::Integer, source::Integer, contact_kind::ContactKind, strain::StrainKind) = new(time, subject, source, TransmissionEvent, UInt8(contact_kind), strain)
+  Event(::Val{TransmissionEvent}; time::Real, subject::Integer, source::Integer, contact_kind::ContactKind, strain::StrainKind) = new(time, subject, source, TransmissionEvent, UInt8(contact_kind), strain)
   Event(::Val{QuarantinedEvent}, time::Real, subject::Integer, extension::Bool) = new(time, subject, 0, QuarantinedEvent, UInt8(extension), NullStrain)
   Event(::Val{TracedEvent}, ::Real, ::Integer) = error("source and tracing kind must be given for TracedEvent")
   Event(::Val{TracedEvent}, time::Real, subject::Integer, source::Integer, tracing_kind::TracingKind) = new(time, subject, source, TracedEvent, UInt8(tracing_kind), NullStrain)

src/r_count.jl

@@ -0,0 +1,96 @@
+#=
+r_count:
+- Julia version: 
+- Author: jacek
+- Date: 2022-03-07
+=#
+using DataStructures: Deque
+using Debugger
+include("../src/default_vector.jl")
+
+function count_infected_per_time(forest::RobinForest, max_days::Int, time_unit=1.0)
+    infected_per_time = Array{Int}(undef, max_days)
+    fill!(infected_per_time, 0)
+    N = size(forest.inedges)
+    for e in forest.inedges
+        if kind(e) == TransmissionEvent
+            day_number = (time(e) / time_unit) |> floor |> Int
+            infected_per_time[day_number] += 1
+            end
+        end
+        infected_per_time
+    end
+
+function count_r_per_time(forest::RobinForest, max_days::Int, time_unit=1.0)
+    r_per_time = DefaultVector(-1.0, Float64)
+    infected_per_time = DefaultVector(0, Int)
+    edge_springsoff_per_time = DefaultVector(0, Int)
+
+    for e in forest.inedges
+        println(e)
+        end
+        println()
+        println("outdegrees")
+    for e in forest.outdegrees
+        println(e)
+        end
+        println()
+        println("outedgedict")
+    for e in forest.outedgedict
+        println(e)
+        end
+
+    for p in forwardinfections(forest, 2)
+        println("2 -> $(subject(p)) $p")
+        end
+
+
+    init_node = 1
+    d = Deque{Int64}()
+    pushfirst!(d, 1)
+    pushfirst!(d, 2)
+
+    println(pop!(d))
+    println(pop!(d))
+
+    person_infecting = zeros(Int, length(forest.inedges))
+
+    person_infected_at = Vector{Int}(undef, length(forest.inedges))
+    person_infected_at[1] = -1
+    person_infecting[1] = 1
+    pushfirst!(d, 1)
+    @bp
+
+    while !isempty(d)
+        n = pop!(d)
+        for e::Event in forwardinfections(forest, n)
+            n2 = subject(e)
+            pushfirst!(d, n2)
+            t = (time(e) / time_unit) |> floor |> Int
+            println("n = $n, n2 = $n2, t = $t")
+            infected_per_time[t] += 1
+            person_infected_at[n2] = t
+            person_infecting[n] = t
+#             if person_infected_at[n] + 1 == person_infected_at[n2]
+#                 edge_springsoff_per_time[t] += 1
+                edge_springsoff_per_time[t] += 1
+#             end
+            println("person_infected_at[$t] = $(person_infected_at[t])")
+        end
+    end
+
+    person_infecting_at_time = zeros(Int, length(forest.inedges))
+    for i in eachindex(person_infecting)
+        if person_infecting[i] > 0
+                person_infecting_at_time[person_infecting[i]] += 1
+            end
+    end
+
+
+#     return convert(Vector, r_per_time)
+    println("edge_springsoff_per_time = $(convert(Vector, edge_springsoff_per_time))")
+    println("infected_per_time = $(convert(Vector, infected_per_time))")
+    println("person_infecting = $(convert(Vector, person_infecting))")
+    println("person_infecting_at_time = $(convert(Vector, person_infecting_at_time))")
+    return [(springs_off / per_day) for (springs_off, per_day) in zip(convert(Vector, infected_per_time), convert(Vector, person_infecting_at_time))]
+    end

test/test_default_vector.jl

@@ -0,0 +1,30 @@
+#=
+test_default_vector:
+- Julia version: 
+- Author: jacek
+- Date: 2022-03-18
+=#
+using Test
+
+include("../src/default_vector.jl")
+
+
+function test_default_vector()
+    for v in [DefaultVector(1), DefaultVector(1, Int)]
+        v[1] = 5
+        v[3] = 7
+        @test v[2] == 1
+        @test convert(Vector, v) == [5, 1, 7]
+
+        v[7] += 1
+
+        @test convert(Vector, v) == [5, 1, 7, 1, 1, 1, 2]
+    end
+
+    for v in [DefaultVector('c'), DefaultVector('c', typeof('c'))]
+        v[3] = 'a'
+        @test convert(Vector, v) == ['c', 'c', 'a']
+    end
+end
+
+test_default_vector()

test/test_forest.jl

@@ -0,0 +1,86 @@
+#=
+test_forest:
+- Julia version: 
+- Author: jacek
+- Date: 2022-02-28
+=#
+
+using Test
+using DataStructures
+using FixedPointNumbers
+const PersonIdx=UInt32
+const TimePoint = Fixed{Int32, 16}
+include("../src/enums.jl")
+include("../src/event.jl")
+include("../src/robin_forest.jl")
+
+function trans(; src, subject, t=0.0)
+    return Event(Val(TransmissionEvent), time=t, subject=subject, source=src, contact_kind=NoContact, strain=NullStrain)
+end
+
+@testset "ForestTesting" begin
+    @testset "R count 1" begin
+        include("../src/r_count.jl")
+        f = RobinForest(10)
+        push!(f, trans(src=1, subject=2, t=1.1))
+        push!(f, trans(src=2, subject=3, t=2.1))
+        push!(f, trans(src=2, subject=4, t=2.2))
+        push!(f, trans(src=2, subject=5, t=2.2))
+        push!(f, trans(src=5, subject=6, t=3.2))
+        push!(f, trans(src=6, subject=7, t=3.3))
+        @test 2 + 2 == 4
+
+        println("backwardinfection $(backwardinfection(f, 2) )")
+        println("backwardinfection subject $(  subject(backwardinfection(f, 2)) )")
+        println("backwardinfection source $(  source(backwardinfection(f, 2)) )")
+
+        for x in forwardinfections(f, 2)
+            println("forwardinfections -> $x")
+        end
+
+        for x in f.inedges
+            println("inedges -> $x")
+        end
+
+        infected_per_day = count_infected_per_time(f, 3)
+        println("infected_per_day = $infected_per_day")
+
+
+        r_per_time = count_r_per_time(f, 3)
+        println("infected_per_day = $r_per_time")
+
+    end
+
+    @testset "Small Tree" begin
+        f = RobinForest(10)
+        push!(f, trans(src=1, subject=2))
+        push!(f, trans(src=2, subject=3))
+        push!(f, trans(src=2, subject=4))
+        @test 2 + 2 == 4
+
+        println("backwardinfection $(backwardinfection(f, 2) )")
+        println("backwardinfection subject $(  subject(backwardinfection(f, 2)) )")
+        println("backwardinfection source $(  source(backwardinfection(f, 2)) )")
+
+        for x in forwardinfections(f, 2)
+            println("forwardinfections -> $x")
+        end
+
+        for x in f.inedges
+            println("inedges -> $x")
+        end
+
+    end
+
+    @testset "ForestForward" begin
+
+        f = RobinForest(10)
+
+        e = Event(Val(TransmissionEvent), time=0.0, subject=2, source=1, contact_kind=NoContact, strain=NullStrain)
+
+        push!(f, e)
+        push!(f, trans(src=1, subject=2))
+        @test 2 + 2 == 4
+
+    end
+end

test/test_r_count.jl

@@ -0,0 +1,72 @@
+#=
+test_r_count:
+- Julia version: 
+- Author: jacek
+- Date: 2022-03-07
+=#
+
+using Test
+using DataStructures
+using FixedPointNumbers
+const PersonIdx=UInt32
+const TimePoint = Fixed{Int32, 16}
+include("../src/enums.jl")
+include("../src/event.jl")
+include("../src/robin_forest.jl")
+include("../src/r_count.jl")
+
+function trans(; src, subject, t=0.0)
+    return Event(Val(TransmissionEvent), time=t, subject=subject, source=src, contact_kind=NoContact, strain=NullStrain)
+end
+
+function nan2zero(x)
+    if isnan(x)
+        return 0.0
+        end
+        return abs(x)
+    end
+
+    function max_abs_diff(arr1, arr2)
+        maximum(map(nan2zero, arr1 .- arr2))
+        end
+
+        EPS = 0.00001
+
+
+function test1()
+    f = RobinForest(10)
+    push!(f, trans(src=1, subject=2, t=1.1))  # day 1
+
+    push!(f, trans(src=2, subject=3, t=2.1))  # day 2
+    push!(f, trans(src=2, subject=4, t=2.2))  # day 2
+    push!(f, trans(src=2, subject=5, t=2.2))  # day 2
+
+    push!(f, trans(src=5, subject=6, t=3.2))  # day 3
+    push!(f, trans(src=6, subject=7, t=3.3))  # day 3
+
+    max_day = maximum(f.inedges) do e
+        time(e) |> floor |> Int
+    end
+
+#         infected_per_day = count_infected_per_time(f, max_day)
+#         println("Number of infected people per each day ($max_day days total) = $infected_per_day")
+#         @test infected_per_day == [1, 3, 2]
+
+
+    r_per_time = count_r_per_time(f, max_day)
+    println("R value per each day ($max_day days total) = $r_per_time")
+#     @test max_abs_diff(r_per_time, [1, 3, 0.33, 0.5] ) < EPS
+end
+
+test1()
+
+#=
+On Jacek machine it runs like:
+
+Number of infected people per each day (3 days total) = [1, 3, 2]
+R value per each day (3 days total) = [NaN, 0.75, 0.3333333333333333]
+Test Summary: | Pass  Total
+RCountTesting |    2      2
+
+Process finished with exit code 0
+=#