Custom implementation of weighted sampling

[Tortar]

This helps the simulation time, shaving another 15% from the timings of https://github.com/bancaditalia/BeforeIT.jl/blob/main/examples/benchmark_w_matlab.jl#L24. Unfortunately it seems to slightly change the results, this should be only because the implementation is different from StatsBase.wsample

[Tortar]

Fixed the issue with the tests, now everything works the same as before

[Devetak]

Great!

[AldoGl]

Before merging I would actually make sure we find a non negligible speedup.
I did run the "benchmark_with_matlabl" script in the examples folder and I could not measure it

[Tortar]

The problem is that the GC runs unreliably between runs, so to measure realibly I think we should use @benchmark or @btime, will come back with the results using them

[Tortar]

With

using BeforeIT

# We will run the model without any output to avoid the overhead of printing the results.
function run(parameters, initial_conditions, T; multi_threading = false)
    model = BeforeIT.init_model(parameters, initial_conditions, T)
    data = BeforeIT.init_data(model);
    
    for _ in 1:T
        BeforeIT.one_epoch!(model; multi_threading = multi_threading)
        BeforeIT.update_data!(data, model)
    end
    return model, data
end

parameters = BeforeIT.AUSTRIA2010Q1.parameters
initial_conditions = BeforeIT.AUSTRIA2010Q1.initial_conditions
T = 1

@benchmark run($parameters, $initial_conditions, $T) 

This is what I see

julia> @benchmark run(parameters, initial_conditions, T) # this PR
BenchmarkTools.Trial: 135 samples with 1 evaluation.
 Range (min … max):  35.383 ms …  41.835 ms  ┊ GC (min … max): 4.44% … 3.82%
 Time  (median):     36.941 ms               ┊ GC (median):    4.67%
 Time  (mean ± σ):   37.141 ms ± 944.112 μs  ┊ GC (mean ± σ):  5.08% ± 1.09%

           ▂  ▂▁▅▇ ▅▄ ▁▅▁ █▁                                    
  ▃▁▁▁▅▅▅▃▃█▆▆███████████▆██▁█▃▁▃▃▃▅▅▆▃▁▁▆▃▃▁▁▁▃▃▁▁▃▁▁▁▁▃▁▁▁▁▃ ▃
  35.4 ms         Histogram: frequency by time         40.4 ms <

 Memory estimate: 122.47 MiB, allocs estimate: 130124.

julia> @benchmark run(parameters, initial_conditions, T) # master
BenchmarkTools.Trial: 117 samples with 1 evaluation.
 Range (min … max):  40.476 ms … 47.697 ms  ┊ GC (min … max): 6.15% … 14.14%
 Time  (median):     42.503 ms              ┊ GC (median):    7.65%
 Time  (mean ± σ):   42.956 ms ±  1.570 ms  ┊ GC (mean ± σ):  8.75% ±  2.40%

              █      ▂  ▃                                      
  ▃▁▁▁▁▆▅▅▆▇▆███▇▇▅█▆█▇▅█▇▅▅▁▃▃▃▅▃▁▁▁▁▆▅▁▃▁▁▁▇▅▆▁▁▅▃▃▁▃▃▁▃▅▁▅ ▃
  40.5 ms         Histogram: frequency by time        46.7 ms <

 Memory estimate: 139.23 MiB, allocs estimate: 678914.

[Devetak]

Confirm @Tortar. With the example from main.jl without printing the speed up for me, with a single thread is 50%.
One thing before merging. In general the function you provide should work best for small-medium sized weights. @AldoGl do we have a larger test case from Poledna's paper?

[AldoGl]

Unfortunately there is no medium-size default parametrisation yet. The standard one has around 8000 agents. We probably need to add one with around 100 000, which btw could be a target size for all experimentations. To generate one would be sufficient to use the calibration script with a custom "scale" parameter

[Tortar]

I guess this should always be the best algorithm, because I don't think there could exist a case where this function is slower if the simulation is dynamic and we can't reuse the sum of weights between different runs (which is the case here as far as I can tell)

[Devetak]

You are absolutely correct! I suggest we merge since it makes the code faster.
I benchmarked the custom weighted sampling with the one from StatsBase and they take the same amount. So I am not sure how the BeforeIT version is so much faster.

function test(N)
       a = rand(N)
       b = rand(N)
       @btime wsample_single(a, b)
       @btime wsample(a, b)
       end

[Tortar]

I think because this version is non-allocating so the GC should be triggered less often, it is probably possible to improve the version of StatsBase, but I'm a bit unsure the time it would take, because maintenance is rather slow there (I have some PRs floating around there from months)

[Tortar]

I have actually pushed some more improvements exploting the unique properties which are satisfied in this case, now it should be even a bit faster

[Tortar]

The last commit broke the tests but it shouldn't have any logical difference between the previous method, it just select differently but with the same probability structure

[Tortar]

@Devetak what do you think? It seems to me that we merge this one after adjusting tests

[Devetak]

@Devetak what do you think? It seems to me that we merge this one after adjusting tests

Few sparse comments:
(a) see also comment on Issue #25, I think some of the later changes make the code less readable. So I will keep deleteat! or put a comment
(b) the test that fail are the ones that compare with matlab. hence the most important ones! We should not change those. Those are the most important ones.

Btw the bug is the change of deleteat! I think. Since

deleteat!([a,b,c,d], 2) = [a,c,d] but

F_g[e], F_g[end] = F_g[end], F_g[e] # [a,d,c,e]
pop!(F_g) # [a,d,c]

In practice they are equivalent, but making it impossible to compare to Matlab is in my opinion not worth it. That is really the most important check! How much speed up there is from that?

[Tortar]

It has a big impact (something like 30 seconds in the big models with scale=0.01) but let's keep this simple and discuss possible improvements later, so I will just remove that improvement. I will just add a function for that for the record. I think that nonetheless we could reproduce the results "in expectation"

[Tortar]

Now it should be okay

[Devetak]

An idea could be to change swap pop to removeat for the deterministic run. This would then mirror the Matlab code.

[Tortar]

you are totally right, we will have to change the some search_and_matching tests anyway but they are not as critical I believe, right?

[Devetak]

Which ones? Also make a test for swap_pop.

[Tortar]

you can see them here: https://github.com/bancaditalia/BeforeIT.jl/actions/runs/11643028585/job/32423294518?pr=20

Will add a test

[Devetak]

I would ask @AldoGl how where the tests constructed and see from there. Definitely important to have comprehensive tests if we are going to work on performance.

[Tortar]

Closing in favour of #30