Plot kills Julia

Author: gdalle

docs/plots_csbm.jl

@@ -8,35 +8,38 @@ using ProgressMeter
 
 BLAS.set_num_threads(1)
 
-function compute_fig1(; α, d, μ, ρ, N_values, λ_values, trials)
+function compute_fig1_csbm(; α, d, μ, ρ, N_values, λ_values, trials)
     rng = default_rng()
     I, J, K = length(N_values), length(λ_values), trials
     qu_values = zeros(I, J, K)
     qv_values = zeros(I, J, K)
     converged_values = zeros(Bool, I, J, K)
     prog = Progress(I * J * K; desc="CSBM - Fig 1")
-    @threads for i in 1:I, j in 1:J, k in 1:K
-        N, λ = N_values[i], λ_values[j]
-        P = ceil(Int, N / α)
-        csbm = CSBM(; N, P, d, μ, λ, ρ)
-        (qu, qv, converged) = evaluate_amp(rng, csbm)
-        qu_values[i, j, k] = qu
-        qv_values[i, j, k] = qv
-        converged_values[i, j, k] = converged
-        next!(prog)
+    for i in 1:I
+        for j in 1:J
+            @threads for k in 1:K
+                N, λ = N_values[i], λ_values[j]
+                P = ceil(Int, N / α)
+                csbm = CSBM(; N, P, d, μ, λ, ρ)
+                (qu, qv, converged) = evaluate_amp(rng, csbm)
+                qu_values[i, j, k] = qu
+                qv_values[i, j, k] = qv
+                converged_values[i, j, k] = converged
+                next!(prog)
+            end
+        end
     end
     return (; qu_values, qv_values)
 end
 
-function plot_fig1(res; N_values, λ_values)
+function plot_fig1_csbm(res; N_values, λ_values)
     qu_values, qv_values = res
     f = Figure()
     for (i, N) in enumerate(N_values)
-        ax = Axis(f[i, 1]; title=" N = $N", xlabel="λ", ylabel="qᵤ", limits=(0, 2, 0, 1))
-        qu_means = dropdims(mean(qu_values[i, :, :]; dims=1); dims=1)
-        qu_stds = dropdims(std(qu_values[i, :, :]; dims=1); dims=1)
+        ax = Axis(f[i, 1]; title="N = $N", xlabel="λ", ylabel="q_u", limits=(0, 2, 0, 1))
+        qu_means = dropdims(mean(qu_values[i, :, :]; dims=2); dims=2)
+        qu_stds = dropdims(std(qu_values[i, :, :]; dims=2); dims=2)
         lines!(ax, λ_values, qu_means)
-        scatter!(ax, λ_values, qu_means)
         errorbars!(ax, λ_values, qu_means, qu_stds)
     end
     return f
@@ -45,10 +48,10 @@ end
 α = 10
 d = 5
 μ = 2
-ρ = 0.1
+ρ = 0.0
 N_values = reverse([3 * 10^3, 10^4, 3 * 10^4, 10^5])
 N_values = reverse([3 * 10^3, 10^4])
 λ_values = 0:0.1:2
 trials = 10
-res1 = compute_fig1(; α, d, μ, ρ, N_values, λ_values, trials)
-plot_fig1(res1; N_values, λ_values)
+res1 = compute_fig1_csbm(; α, d, μ, ρ, N_values, λ_values, trials)
+plot_fig1_csbm(res1; N_values, λ_values)

docs/plots_glmsbm.jl

@@ -0,0 +1,63 @@
+using Base.Threads
+using CairoMakie
+using LinearAlgebra
+using Random: default_rng
+using Statistics
+using StochasticBlockModelVariants
+using ProgressMeter
+
+BLAS.set_num_threads(1)
+
+function compute_fig1_glmsbm(; N, c, ρ, Pʷ, α_values, λ_values, trials)
+    rng = default_rng()
+    I, J, K = length(α_values), length(λ_values), trials
+    qs_values = zeros(I, J, K)
+    qw_values = zeros(I, J, K)
+    converged_values = zeros(Bool, I, J, K)
+    prog = Progress(I * J * K; desc="CSBM - Fig 1")
+    for i in 1:I
+        for j in 1:J
+            @threads for k in 1:K
+                α, λ = α_values[i], λ_values[j]
+                M = ceil(Int, N / α)
+                glmsbm = GLMSBM(; N, M, c, λ, ρ, Pʷ)
+                (qs, qw, converged) = evaluate_amp(rng, glmsbm)
+                qs_values[i, j, k] = qs
+                qw_values[i, j, k] = qw
+                converged_values[i, j, k] = converged
+                next!(prog)
+            end
+        end
+    end
+    return (; qs_values, qw_values)
+end
+
+function plot_fig1_glmsbm(res; α_values, λ_values)
+    qs_values, qw_values = res
+    f = Figure()
+    ax1 = Axis(f[1, 1]; title="Communities", xlabel="λ", ylabel="q_s", limits=(0, 2, 0, 1))
+    ax2 = Axis(f[1, 2]; title="Weights", xlabel="λ", ylabel="q_w", limits=(0, 2, 0, 1))
+    for (i, α) in enumerate(α_values)
+        qs_means = dropdims(mean(qs_values[i, :, :]; dims=2); dims=2)
+        qs_stds = dropdims(std(qs_values[i, :, :]; dims=2); dims=2)
+        lines!(ax1, λ_values, qs_means; label="α=$α")
+        errorbars!(ax1, λ_values, qs_means, qs_stds; label=nothing)
+
+        qw_means = dropdims(mean(qs_values[i, :, :]; dims=2); dims=2)
+        qw_stds = dropdims(std(qs_values[i, :, :]; dims=2); dims=2)
+        lines!(ax2, λ_values, qw_means; label="α=$α")
+        errorbars!(ax2, λ_values, qw_means, qw_stds; label=nothing)
+    end
+    return f
+end
+
+N = 10^4
+c = 5
+ρ = 0.0
+Pʷ = GaussianWeightPrior()
+α_values = reverse([0.3, 1, 3, 10, 30])
+α_values = reverse([0.3, 1])
+λ_values = 0:0.1:2
+trials = 10
+res1 = compute_fig1_glmsbm(; N, c, ρ, Pʷ, α_values, λ_values, trials)  # kills Julia
+plot_fig1_glmsbm(res1; α_values, λ_values)

src/csbm_inference.jl

@@ -129,7 +129,7 @@ function run_amp(
     init_std=1e-3,
     max_iterations=200,
     convergence_threshold=1e-3,
-    recent_past=10,
+    recent_past=5,
     damping=0.0,
     show_progress=false,
 )

src/glmsbm_inference.jl

@@ -170,7 +170,7 @@ function run_amp(
     init_std=1e-3,
     max_iterations=100,
     convergence_threshold=1e-3,
-    recent_past=10,
+    recent_past=5,
     damping=0.5,
     show_progress=false,
 )

src/precompile.jl

@@ -5,6 +5,6 @@
         GLMSBM(; N=10^2, M=10^2, c=5, λ=2, ρ=0.0, Pʷ=GaussianWeightPrior()),
         GLMSBM(; N=10^2, M=10^2, c=5, λ=2, ρ=0.1, Pʷ=RademacherWeightPrior()),
     ]
-        evaluate_amp(rng, sbm)
+        evaluate_amp(rng, sbm; max_iterations=2)
     end
 end