write tests for contact vectors and degree distributions finally

CovidAlertVaccinationModel/prem_mixing_approx.jl

-using CovidAlertVaccinationModel:abm,ModelSolution, get_parameters, neighbors, GraphEdge, get_weight
+using CovidAlertVaccinationModel:abm,ModelSolution, get_parameters, neighbors, GraphEdge, get_weight, contact_time_distributions
 using LightGraphs
 function approximate_mixing_matricies_1(sol)
     mean_mixing_degree = zeros(3,3)
@@ -62,5 +62,8 @@ sol = abm(p,nothing)
 approximate_mixing_matricies_1(sol);
 approximate_mixing_matricies_2(sol);
 println("dasdlkjasdas")
+
+
+
 #rerun with weights
 #
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/contact_vectors.jl

@@ -17,6 +17,12 @@ We repeat this function until `csum == 0`.
     j_to_i_contacts[j_index] = sample_list_j[k]    
     return csum
 end
+@inline function reindex!(k,csum,index_list_i,i_to_i_contacts,sample_list_i)
+    i_index = index_list_i[k]  
+    csum +=  i_to_i_contacts[i_index] -  sample_list_i[k]
+    i_to_i_contacts[i_index] = sample_list_i[k]   
+    return csum
+end
 
 
 """
@@ -53,4 +59,25 @@ function generate_contact_vectors!(ij_dist,ji_dist,i_to_j_contacts::Vector{T}, j
         end
     end
     return nothing
+end
+function generate_contact_vectors!(ii_dist,i_to_i_contacts::Vector{T}) where T
+    
+    rand!(Random.default_rng(Threads.threadid()),ii_dist,i_to_i_contacts)
+    l_i = length(i_to_i_contacts)
+
+    csum = sum(i_to_i_contacts)
+    inner_iter = 10
+    index_list_i = @MVector zeros(Int,inner_iter)
+    sample_list_i = @MVector zeros(T,inner_iter)
+    
+    while csum % 2 != 0
+        sample!(Random.default_rng(Threads.threadid()),1:l_i,index_list_i)
+        rand!(Random.default_rng(Threads.threadid()),ii_dist,sample_list_i)
+        @inbounds for i = 1:inner_iter
+            if csum != 0
+                csum = reindex!(i,csum,index_list_i,i_to_i_contacts,sample_list_i)
+            end
+        end
+    end
+    return nothing
 end
\ No newline at end of file

CovidAlertVaccinationModel/src/ABM/mixing_graphs.jl

@@ -87,20 +87,34 @@ function create_mixing_edges(demographic_index_vectors,mixing_matrix,weights_dis
     contact_array =[Vector{GraphEdge}() for i in 1:length(demographic_index_vectors),j in 1:length(demographic_index_vectors)]
     tot = 0
     for i in 1:3, j in 1:i  #diagonal
-        num_degrees_ij = similar(demographic_index_vectors[i])
-        num_degrees_ji = similar(demographic_index_vectors[j])
-        generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji)
-        m = sum(num_degrees_ij)
-        stubs_i = Vector{Int}(undef,m)
-        stubs_j = Vector{Int}(undef,m)
-        if m>0
-            sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[i],Weights(num_degrees_ij./m),stubs_i)
-            sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[j],Weights(num_degrees_ji./m),stubs_j)
-            tot += m
-        end    
-
-        contact_array[j,i] = GraphEdge.(stubs_i,stubs_j)
-
+        if i != j 
+            num_degrees_ij = similar(demographic_index_vectors[i])
+            num_degrees_ji = similar(demographic_index_vectors[j])
+            generate_contact_vectors!(mixing_matrix[i,j],mixing_matrix[j,i],num_degrees_ij,num_degrees_ji)
+            num_edges = sum(num_degrees_ij)
+            stubs_i = Vector{Int}(undef,num_edges)
+            stubs_j = Vector{Int}(undef,num_edges)
+            if num_edges>0
+                sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[i],Weights(num_degrees_ij./num_edges),stubs_i)
+                sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[j],Weights(num_degrees_ji./num_edges),stubs_j)
+                tot += num_edges
+            end    
+
+            contact_array[j,i] = GraphEdge.(stubs_i,stubs_j)
+        else #from one group to itself we need another algorithm
+            num_degrees_ii = similar(demographic_index_vectors[i])
+            generate_contact_vectors!(mixing_matrix[i,j],num_degrees_ii)
+            m = sum(num_degrees_ii)
+            num_edges= div(m,2)
+            stubs_i = Vector{Int}(undef,num_edges)
+            stubs_j = Vector{Int}(undef,num_edges)
+            if m>0
+                sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[i],Weights(num_degrees_ii./m),stubs_i)
+                sample!(Random.default_rng(Threads.threadid()),demographic_index_vectors[i],Weights(num_degrees_ii./m),stubs_j)
+                tot += num_edges
+            end 
+            contact_array[j,i] = GraphEdge.(stubs_i,stubs_j)
+        end
     end
     return MixingEdges(tot,contact_array,weights_distribution_matrix)
 end

CovidAlertVaccinationModel/test/ABM/mixing_test.jl

 
-using CovidAlertVaccinationModel
-using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,AgentStatus,get_u_0,get_parameters,solve!, DebugRecorder, WeightedGraph
+using CovidAlertVaccinationModel:generate_contact_vectors!
+using CovidAlertVaccinationModel:ModelSolution,AgentDemographic,mean,get_weight
+using CovidAlertVaccinationModel:AgentStatus,get_u_0,get_parameters,solve!, DebugRecorder, WeightedGraph, contact_time_distributions,GraphEdge
+using LightGraphs
 using Random
 using Plots
+using UnPack
 using ThreadsX
 using OnlineStats
+using StatsBase
 const model_sizes = [1000,5000]
  const dem_cat = AgentDemographic.size -1
- const samples = 10
+ const samples = 1
 
 #WRITE MORE TESTS
 
-@testset "mixing matrices, size: $sz" for sz in model_sizes
-    for rep = 1:samples
-        m = ModelSolution(100,get_parameters(),sz)
-        ws_dist = m.ws_matrix_tuple
-        r_dist = m.rest_matrix_tuple 
-        index_vec =m.index_vectors
-        @testset "workschool" for i in dem_cat, j in dem_cat
-            for t in 1:length(ws_dist)
-                @test mean(ws_dist[t][i,j])*length(index_vec[i]) == mean(ws_dist[t][j,i])*length(index_vec[j])
-            end
-        end
-        @testset "rest" for i in dem_cat, j in dem_cat
-            for t in 1:length(ws_dist)
-                @test mean(r_dist[t][i,j])*length(index_vec[i]) == mean(r_dist[t][j,i])*length(index_vec[j])
-            end
-        end
-    end
-end
+# @testset "mixing matrices, size: $sz" for sz in model_sizes
+#     for rep = 1:samples
+#         m = ModelSolution(100,get_parameters(),sz)
+#         ws_dist = m.ws_matrix_tuple
+#         r_dist = m.rest_matrix_tuple 
+#         index_vec =m.index_vectors
+#         @testset "workschool" for i in dem_cat, j in dem_cat
+#             for t in 1:length(ws_dist)
+#                 @test mean(ws_dist[t][i,j])*length(index_vec[i]) == mean(ws_dist[t][j,i])*length(index_vec[j])
+#             end
+#         end
+#         @testset "rest" for i in dem_cat, j in dem_cat
+#             for t in 1:length(ws_dist)
+#                 @test mean(r_dist[t][i,j])*length(index_vec[i]) == mean(r_dist[t][j,i])*length(index_vec[j])
+#             end
+#         end
+#     end
+# end
 
 function approximate_mixing_matricies_weights_dict(p,sol)
     mean_mixing_degree = zeros(3,3)
@@ -69,18 +73,69 @@ function approximate_mixing_matricies_graph(p,sol)
         mean_mixing_weighted_degree ./= (p.sim_length .* length.(sol.index_vectors) * 2)
     return mean_mixing_degree, mean_mixing_weighted_degree
 end
-@testset "weights dict and graph match" for _ in samples
-    p = get_parameters()
-    sol = abm(p,nothing)
-    graph_deg_matrix,graph_wdeg_matrix = approximate_mixing_matricies_graph(p,sol)
-    dict_deg_matrix,dict_wdeg_matrix = approximate_mixing_matricies_weights_dict(p,sol)
+# @testset "weights dict and graph match" for _ in samples
+#     p = get_parameters()
+#     sol = abm(p,nothing)
+#     graph_deg_matrix,graph_wdeg_matrix = approximate_mixing_matricies_graph(p,sol)
+#     dict_deg_matrix,dict_wdeg_matrix = approximate_mixing_matricies_weights_dict(p,sol)
 
-    @test all(graph_deg_matrix .≈ dict_deg_matrix)
-    @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
-end
+#     @test all(graph_deg_matrix .≈ dict_deg_matrix)
+#     @test all(graph_wdeg_matrix .≈ dict_wdeg_matrix)
+# end
+
+using OnlineStats
 
 
+# @testset "contact vector generation" for _ in samples
+#     model = ModelSolution(1, get_parameters(), 1_000_000) #big
+#     @unpack demographics, index_vectors, rest_matrix_tuple,ws_matrix_tuple = model
+#     ThreadsX.map((ws_matrix_tuple...,rest_matrix_tuple...)) do dist
+#         for i in 1:3, j in 1:i  #diagonal
+#             if i != j 
+#                 num_degrees_ij = similar(index_vectors[i])
+#                 num_degrees_ji = similar(index_vectors[j])
+#                 generate_contact_vectors!(dist[i,j],dist[j,i],num_degrees_ij,num_degrees_ji) 
+#                 @test sum(num_degrees_ji) == sum(num_degrees_ij)
+#                 @test mean(num_degrees_ij) ≈ mean(dist[i,j]) atol = 1e-2
+#                 @test mean(num_degrees_ji) ≈ mean(dist[j,i]) atol = 1e-2
+#                 @test var(num_degrees_ij) ≈ var(dist[i,j]) atol = 1e-1
+#                 @test var(num_degrees_ji) ≈ var(dist[j,i]) atol = 1e-1
+#                 @test skewness(num_degrees_ij) ≈ skewness(dist[i,j]) atol = 5e-1
+#                 @test skewness(num_degrees_ji) ≈ skewness(dist[j,i]) atol = 5e-1
+#             else
+#                 num_degrees_ii = similar(index_vectors[i])
+#                 generate_contact_vectors!(dist[i,j],num_degrees_ii) 
+#                 @test iseven(sum(num_degrees_ii))
+#                 @test mean(num_degrees_ii) ≈ mean(dist[i,j]) atol = 1e-2
+#                 @test var(num_degrees_ii) ≈ var(dist[i,j]) atol = 1e-1
+#                 @test skewness(num_degrees_ii) ≈ skewness(dist[i,j]) atol = 5e-1
+#             end
+#         end
+#     end
+# end
 
-@testset "degree distribution generation" for _ in samples
-    
+@testset "degree distribution generation"  for _ in samples
+    model = ModelSolution(1, get_parameters(), 100_000)
+    @unpack demographics, index_vectors, rest_matrix_tuple, ws_matrix_tuple = model
+    for dist in (ws_matrix_tuple...,rest_matrix_tuple...)
+        g = WeightedGraph(demographics,index_vectors,dist,contact_time_distributions.ws)
+        mixing_dist = [Variance() for _ in 1:3, _ in 1:3]
+        for v in vertices(g.g)  
+            demo_v = demographics[v]
+            degs = zeros(3)
+            for w in LightGraphs.neighbors(g.g,v)
+                demo_w = demographics[w]
+                degs[Int(demo_w)] += 1
+            end
+            for (j,d) in enumerate(degs)
+                fit!(mixing_dist[Int(demo_v), j],d)
+            end
+        end
+        for i in eachindex(mixing_dist)
+            @test mean(mixing_dist[i]) ≈ mean(dist[i])  atol = 0.2
+        end
+        display(var.(mixing_dist))
+        display(var.(dist))
+        display("====")
+    end
 end
\ No newline at end of file

CovidAlertVaccinationModel/test/runtests.jl

+using OnlineStats: neighbors
 
 using RandomNumbers.Xorshifts
 using Test