==============================================================================
THE SAME DISEASE ON FOUR DIFFERENT NETWORKS
Science Journaling Club, Volume 1 Issue 2, Winter 2025
==============================================================================
Master seed                 : 20251215
Nodes per network N         : 10000
Target mean degree <k>      : 6
Recovery probability gamma  : 0.200 per step (mean infectious period 5.0 steps)
Large-outbreak cutoff       : 100 nodes (1.0% of N)
Runs per configuration      : 600 (main sweep), 600 (threshold refinement)
Mode                        : full
numpy                       : 2.4.2

==============================================================================
PART 1.  BUILDING THE FOUR NETWORKS
==============================================================================
Every graph below is a single fixed realisation drawn from the master seed.
The same four graphs are reused for every transmission rate, so differences
between networks are not confounded with graph-to-graph noise.

network      edges      <k>     <k^2>   k_min   k_max  transit. mean_path    giant
------------------------------------------------------------------------------
lattice      30000   6.0000    36.000       6       6    0.4000    38.891   1.0000
ER           30000   6.0000    42.025       0      18    0.0004     5.342   0.9977
WS           30000   6.0000    36.292       3       9    0.5088    11.215   1.0000
BA           29994   5.9988   114.291       3     287    0.0029     4.267   1.0000

Degree histograms (count of nodes at each degree; BA truncated at 20):
  lattice  6:10000
  ER       0:23 1:136 2:459 3:896 4:1353 5:1645 6:1583 7:1327 8:991 9:723 10:442 11:216 12:128 13:44 14:23 15:7 16:2 18:2
  WS       3:2 4:51 5:1202 6:7526 7:1132 8:83 9:4
  BA       3:4032 4:2011 5:1105 6:708 7:460 8:344 9:238 10:182 11:145 12:100 13:94 14:82 15:63 16:53 17:37 18:32 19:20 20:29  >20:265

BA degree tail, nodes with degree >= 40:
  count = 78, top ten degrees = [287, 269, 250, 201, 182, 179, 162, 132, 130, 126]

-- Graph-to-graph variability of <k^2> (4 extra seeds) -----------------------
The error bars later in this study come from run-to-run variation on ONE
fixed graph. This table shows how much the graph itself would move if we
redrew it, which is a source of uncertainty we do NOT propagate.

network       seed A     seed B     seed C     seed D
------------------------------------------------------------------------------
ER            41.993     42.044     42.209     42.022
WS            36.299     36.291     36.299     36.279
BA           110.519    112.402    107.931    110.950
lattice   deterministic, <k^2> = 36.000 by construction

Parallel workers            : 8
Run seeds are assigned by task index, so the numbers below do not depend
on how many cores happen to be available.

==============================================================================
PART 2.  VALIDATION A - THE MASS-ACTION CASE AGAINST THE ANALYTIC ANSWER
==============================================================================
No network at all. Each infective makes Poisson(lambda) contacts per step
with uniformly chosen members of the population; R0 = lambda / gamma.
In the N -> infinity limit the final attack rate z among those ever infected
solves the Kermack-McKendrick final size equation

        z = 1 - exp(-R0 * z)

and the threshold sits exactly at R0 = 1. This holds for ANY infectious
period distribution, because only the total infectious person-time enters.
Our infectious period is geometric, so this is a real test of the code and
not a tautology.

R0       n_large  measured   analytic       diff         SE   P(large)  branching
                    attack          z  meas-pred              measured   1-q pred
------------------------------------------------------------------------------
0.50           0         -     0.0000          -          -     0.0000     0.0000
0.70           0         -     0.0000          -          -     0.0000     0.0000
0.90           8    0.0185     0.0000     0.0185     0.0037     0.0133     0.0000
1.00          44    0.0361     0.0000     0.0361     0.0049     0.0733     0.0000
1.10          62    0.1415     0.1761    -0.0347     0.0095     0.1033     0.1010
1.25         142    0.3640     0.3714    -0.0074     0.0056     0.2367     0.2220
1.50         221    0.5832     0.5828     0.0003     0.0011     0.3683     0.3696
1.75         288    0.7116     0.7127    -0.0011     0.0007     0.4800     0.4747
2.00         319    0.7977     0.7968     0.0009     0.0005     0.5317     0.5533
2.50         387    0.8928     0.8926     0.0001     0.0003     0.6450     0.6626
3.00         432    0.9403     0.9405    -0.0002     0.0002     0.7200     0.7348
4.00         495    0.9801     0.9802    -0.0001     0.0001     0.8250     0.8234

-- Verdict on validation A ---------------------------------------------------
Largest absolute discrepancy between measured attack rate and the analytic
final size equation, over all R0 >= 1.25 : 0.00741
Tolerance we set before running          : 0.01500
PASS

Threshold check. Below and at R0 = 1 the analytic final size is exactly zero
and the epidemic probability is zero; a finite population still throws up the
occasional run above the 1% cutoff. Measured P(large) at R0 <= 1:
   R0 = 0.50   P(large) = 0.0000 +/- 0.0000   mean final size over ALL runs = 2.2 nodes
   R0 = 0.70   P(large) = 0.0000 +/- 0.0000   mean final size over ALL runs = 2.9 nodes
   R0 = 0.90   P(large) = 0.0133 +/- 0.0047   mean final size over ALL runs = 8.5 nodes
   R0 = 1.00   P(large) = 0.0733 +/- 0.0106   mean final size over ALL runs = 31.7 nodes

And just above: 
   R0 = 1.10   P(large) = 0.1033 +/- 0.0124   analytic 1-q = 0.1010
   R0 = 1.25   P(large) = 0.2367 +/- 0.0174   analytic 1-q = 0.2220
   R0 = 1.50   P(large) = 0.3683 +/- 0.0197   analytic 1-q = 0.3696

Note that P(large) and the attack rate are DIFFERENT numbers here. With a
Poisson offspring distribution they would coincide. Our offspring count is
Poisson(lambda*L) with L geometric, which is over-dispersed, so extinction is
more likely than the naive z would suggest. The measured P(large) tracks the
branching-process prediction, not z. That agreement is the second validation.

-- Mean prevalence curve, mass action at R0 = 2.0 (319 large outbreaks) ------
step : I(t)/N
    0 : 0.00010
    4 : 0.00033
    8 : 0.00072
   12 : 0.00148
   16 : 0.00304
   20 : 0.00623
   24 : 0.01251
   28 : 0.02410
   32 : 0.04307
   36 : 0.06909
   40 : 0.09758
   44 : 0.11991
   48 : 0.12862
   52 : 0.12247
   56 : 0.10487
   60 : 0.08227
   64 : 0.06020
   68 : 0.04217
   72 : 0.02849
   76 : 0.01877
   80 : 0.01210
   84 : 0.00769
   88 : 0.00478
   92 : 0.00297
   96 : 0.00187
  100 : 0.00120
  104 : 0.00079
  108 : 0.00056
  112 : 0.00044
  116 : 0.00035
  120 : 0.00029
  124 : 0.00025
  128 : 0.00021
  132 : 0.00017
  136 : 0.00016
  140 : 0.00015
  144 : 0.00017
  148 : 0.00020
  152 : 0.00020
  156 : 0.00015

==============================================================================
PART 3.  MAIN SWEEP - THE SAME DISEASE ON FOUR NETWORKS
==============================================================================
Transmissibility T is the probability that a single edge ever carries the
infection. It is the honest common axis: two networks at the same T have
identical per-contact disease parameters. beta = T*gamma / (1 - T*(1-gamma)).

Well-mixed reference: a homogeneous population with the same mean degree has
R0 = T * <k> = 6T, so its threshold sits at T = 1/6 = 0.166667.

  lattice  done   (42 s elapsed)
  ER       done   (56 s elapsed)
  WS       done   (72 s elapsed)
  BA       done   (88 s elapsed)

-- lattice  -  Triangular lattice, 100x100 torus, degree exactly 6 -----------
     T     beta P(large)      +/-SE  n_large  finalfrac    +/-SE  peakfrac    +/-SE    t_peak
------------------------------------------------------------------------------
 0.020   0.0041   0.0000     0.0000        0          -        -         -        -         -
 0.040   0.0083   0.0000     0.0000        0          -        -         -        -         -
 0.060   0.0126   0.0000     0.0000        0          -        -         -        -         -
 0.080   0.0171   0.0000     0.0000        0          -        -         -        -         -
 0.100   0.0217   0.0000     0.0000        0          -        -         -        -         -
 0.140   0.0315   0.0000     0.0000        0          -        -         -        -         -
 0.180   0.0421   0.0000     0.0000        0          -        -         -        -         -
 0.220   0.0534   0.0000     0.0000        0          -        -         -        -         -
 0.280   0.0722   0.0500     0.0089       30     0.0134   0.0006   0.00162  0.00007     47.73
 0.350   0.0972   0.4283     0.0202      257     0.0643   0.0042   0.00342  0.00009     92.86
 0.420   0.1265   0.7300     0.0181      438     0.8703   0.0040   0.03136  0.00027    220.46
 0.500   0.1667   0.8450     0.0148      507     0.9722   0.0001   0.06021  0.00016    136.64
 0.650   0.2708   0.9650     0.0075      579     0.9978   0.0000   0.09599  0.00014     88.13
 0.800   0.4444   0.9917     0.0037      595     0.9999   0.0000   0.12685  0.00012     65.98

-- ER  -  Erdos-Renyi G(N,M), M = 30000 --------------------------------------
     T     beta P(large)      +/-SE  n_large  finalfrac    +/-SE  peakfrac    +/-SE    t_peak
------------------------------------------------------------------------------
 0.020   0.0041   0.0000     0.0000        0          -        -         -        -         -
 0.040   0.0083   0.0000     0.0000        0          -        -         -        -         -
 0.060   0.0126   0.0000     0.0000        0          -        -         -        -         -
 0.080   0.0171   0.0000     0.0000        0          -        -         -        -         -
 0.100   0.0217   0.0000     0.0000        0          -        -         -        -         -
 0.140   0.0315   0.0033     0.0024        2     0.0106   0.0005   0.00100  0.00010     39.50
 0.180   0.0421   0.1383     0.0141       83     0.1040   0.0081   0.00608  0.00042    108.60
 0.220   0.0534   0.3083     0.0189      185     0.4367   0.0036   0.04048  0.00052     86.75
 0.280   0.0722   0.5050     0.0204      303     0.6811   0.0006   0.12090  0.00040     51.40
 0.350   0.0972   0.6933     0.0188      416     0.8223   0.0003   0.22041  0.00034     34.87
 0.420   0.1265   0.7733     0.0171      464     0.8953   0.0002   0.31234  0.00033     26.04
 0.500   0.1667   0.8367     0.0151      502     0.9410   0.0001   0.40697  0.00029     20.18
 0.650   0.2708   0.9450     0.0093      567     0.9785   0.0001   0.55612  0.00028     13.28
 0.800   0.4444   0.9800     0.0057      588     0.9918   0.0000   0.67717  0.00035      9.46

-- WS  -  Watts-Strogatz ring, k=6, rewiring p = 0.05 ------------------------
     T     beta P(large)      +/-SE  n_large  finalfrac    +/-SE  peakfrac    +/-SE    t_peak
------------------------------------------------------------------------------
 0.020   0.0041   0.0000     0.0000        0          -        -         -        -         -
 0.040   0.0083   0.0000     0.0000        0          -        -         -        -         -
 0.060   0.0126   0.0000     0.0000        0          -        -         -        -         -
 0.080   0.0171   0.0000     0.0000        0          -        -         -        -         -
 0.100   0.0217   0.0000     0.0000        0          -        -         -        -         -
 0.140   0.0315   0.0000     0.0000        0          -        -         -        -         -
 0.180   0.0421   0.0000     0.0000        0          -        -         -        -         -
 0.220   0.0534   0.0000     0.0000        0          -        -         -        -         -
 0.280   0.0722   0.0333     0.0073       20     0.0162   0.0011   0.00195  0.00008     53.10
 0.350   0.0972   0.2983     0.0187      179     0.1808   0.0097   0.00837  0.00042    130.49
 0.420   0.1265   0.6217     0.0198      373     0.7241   0.0023   0.06072  0.00040    101.63
 0.500   0.1667   0.8200     0.0157      492     0.9228   0.0005   0.14528  0.00039     61.30
 0.650   0.2708   0.9667     0.0073      580     0.9954   0.0001   0.30087  0.00036     33.30
 0.800   0.4444   0.9967     0.0024      598     0.9999   0.0000   0.43601  0.00032     21.49

-- BA  -  Barabasi-Albert preferential attachment, m = 3 ---------------------
     T     beta P(large)      +/-SE  n_large  finalfrac    +/-SE  peakfrac    +/-SE    t_peak
------------------------------------------------------------------------------
 0.020   0.0041   0.0000     0.0000        0          -        -         -        -         -
 0.040   0.0083   0.0033     0.0024        2     0.0122   0.0018   0.00245  0.00025     22.00
 0.060   0.0126   0.0117     0.0044        7     0.0128   0.0006   0.00326  0.00043     26.29
 0.080   0.0171   0.0250     0.0064       15     0.0370   0.0044   0.00561  0.00056     45.20
 0.100   0.0217   0.0750     0.0108       45     0.0821   0.0047   0.01225  0.00066     41.56
 0.140   0.0315   0.1683     0.0153      101     0.2249   0.0023   0.03665  0.00073     40.40
 0.180   0.0421   0.3000     0.0187      180     0.3695   0.0013   0.07321  0.00060     34.01
 0.220   0.0534   0.4133     0.0201      248     0.4967   0.0010   0.11683  0.00055     28.24
 0.280   0.0722   0.5183     0.0204      311     0.6496   0.0007   0.18974  0.00051     22.96
 0.350   0.0972   0.6400     0.0196      384     0.7744   0.0005   0.27304  0.00046     19.02
 0.420   0.1265   0.7333     0.0181      440     0.8576   0.0003   0.35383  0.00040     15.37
 0.500   0.1667   0.8283     0.0154      497     0.9198   0.0002   0.44060  0.00035     12.78
 0.650   0.2708   0.9417     0.0096      565     0.9772   0.0001   0.58459  0.00031      9.46
 0.800   0.4444   0.9733     0.0066      584     0.9963   0.0000   0.71033  0.00038      7.11

-- Head-to-head at four transmissibilities (final attack fraction, large outbreaks) 
       T      lattice           ER           WS           BA   well-mixed
------------------------------------------------------------------------------
   0.100   (no epid.)   (no epid.)   (no epid.)       0.0821       0.0000
   0.180   (no epid.)       0.1040   (no epid.)       0.3695       0.1444
   0.280       0.0134       0.6811       0.0162       0.6496       0.6820
   0.500       0.9722       0.9410       0.9228       0.9198       0.9405
   0.800       0.9999       0.9918       0.9999       0.9963       0.9914

-- Head-to-head: fraction of runs that fizzle out (final size < 100 nodes) ---
       T      lattice           ER           WS           BA
------------------------------------------------------------------------------
   0.020       1.0000       1.0000       1.0000       1.0000
   0.040       1.0000       1.0000       1.0000       0.9967
   0.060       1.0000       1.0000       1.0000       0.9883
   0.080       1.0000       1.0000       1.0000       0.9750
   0.100       1.0000       1.0000       1.0000       0.9250
   0.140       1.0000       0.9967       1.0000       0.8317
   0.180       1.0000       0.8617       1.0000       0.7000
   0.220       1.0000       0.6917       1.0000       0.5867
   0.280       0.9500       0.4950       0.9667       0.4817
   0.350       0.5717       0.3067       0.7017       0.3600
   0.420       0.2700       0.2267       0.3783       0.2667
   0.500       0.1550       0.1633       0.1800       0.1717
   0.650       0.0350       0.0550       0.0333       0.0583
   0.800       0.0083       0.0200       0.0033       0.0267

-- Head-to-head: peak prevalence and time to peak at T = 0.50 ----------------
network        peak I/N           SE       t_peak           SE
------------------------------------------------------------------------------
lattice         0.06021      0.00016       136.64         0.42
ER              0.40697      0.00029        20.18         0.11
WS              0.14528      0.00039        61.30         0.40
BA              0.44060      0.00035        12.78         0.10

-- Mean prevalence curves at T = 0.50 (I(t)/N, averaged over large outbreaks) 
(curves averaged over all large outbreaks, zero-padded after an outbreak ends)
 step      lattice           ER           WS           BA
------------------------------------------------------------------------------
    0      0.00010      0.00010      0.00010      0.00010
    4      0.00064      0.00097      0.00059      0.00589
    8      0.00143      0.00735      0.00117      0.13605
   12      0.00236      0.05029      0.00199      0.37622
   16      0.00333      0.21409      0.00320      0.31681
   20      0.00450      0.34782      0.00524      0.17234
   24      0.00587      0.26942      0.00842      0.08067
   28      0.00717      0.15259      0.01340      0.03557
   32      0.00861      0.07444      0.02093      0.01521
   36      0.01010      0.03345      0.03167      0.00641
   40      0.01170      0.01437      0.04635      0.00269
   44      0.01345      0.00609      0.06433      0.00112
   48      0.01512      0.00257      0.08438      0.00047
   52      0.01678      0.00108      0.10273      0.00020
   56      0.01867      0.00043      0.11528      0.00008
   60      0.02046      0.00017      0.11903      0.00003
   64      0.02241      0.00007      0.11361      0.00002
   68      0.02415      0.00003      0.10087      0.00000
   72      0.02597      0.00001      0.08419      0.00000
   76      0.02764      0.00001      0.06647      0.00000
   80      0.02969      0.00000      0.05044      0.00000
   84      0.03158      0.00000      0.03679      0.00000
   88      0.03366      0.00000      0.02629      0.00000
   92      0.03570      0.00000      0.01834      0.00000
   96      0.03747      0.00000      0.01267      0.00000
  100      0.03922      0.00000      0.00879      0.00000
  104      0.04132      0.00000      0.00595      0.00000
  108      0.04318      0.00000      0.00400      0.00000
  112      0.04511      0.00000      0.00263      0.00000
  116      0.04722      0.00000      0.00169      0.00000
  120      0.04948      0.00000      0.00108      0.00000
  124      0.05149      0.00000      0.00066      0.00000
  128      0.05329      0.00000      0.00038      0.00000
  132      0.05460      0.00000      0.00024      0.00000
  136      0.05534      0.00000      0.00016      0.00000
  140      0.05475      0.00000      0.00011      0.00000
  144      0.05254      0.00000      0.00008      0.00000
  148      0.04869      0.00000      0.00006      0.00000
  152      0.04294      0.00000      0.00004      0.00000
  156      0.03579      0.00000      0.00002      0.00000
  160      0.02842      0.00000      0.00002      0.00000
  164      0.02140      0.00000      0.00001      0.00000
  168      0.01523      0.00000      0.00001      0.00000
  172      0.01022      0.00000      0.00000      0.00000
  176      0.00658      0.00000      0.00000      0.00000
  180      0.00411      0.00000      0.00000      0.00000
  184      0.00243      0.00000      0.00000      0.00000
  188      0.00140      0.00000      0.00000      0.00000
  192      0.00079      0.00000      0.00000      0.00000
  196      0.00042      0.00000      0.00000      0.00000
  200      0.00022      0.00000      0.00000      0.00000
  204      0.00011      0.00000      0.00000      0.00000
  208      0.00005      0.00000      0.00000      0.00000
  212      0.00003      0.00000      0.00000      0.00000
  216      0.00001      0.00000      0.00000      0.00000
  220      0.00001      0.00000      0.00000      0.00000
  224      0.00000      0.00000      0.00000      0.00000
  228      0.00000      0.00000      0.00000      0.00000
  232      0.00000      0.00000      0.00000      0.00000

==============================================================================
PART 4.  VALIDATION B - EPIDEMIC THRESHOLDS, MEASURED AGAINST PREDICTED
==============================================================================
Two standard predictions, both in transmissibility:

  Molloy-Reed / Newman (locally tree-like bond percolation):
        T_c = <k> / (<k^2> - <k>)
  Heterogeneous mean field (Pastor-Satorras and Vespignani):
        T_c = <k> / <k^2>

The well-mixed prediction is T_c = 1/<k> = 0.166667.

Estimator 1 (primary): the susceptibility chi = <s^2>/<s> over the final
sizes of runs that fizzle. In percolation this peaks at the critical point.
We refine the peak with a parabola through the three surrounding grid points
and bootstrap the runs for a 95% interval.
Estimator 2: extrapolate the conditional attack fraction linearly to zero
over the window where it lies between 0.03 and 0.25.

  lattice  refinement done   (118 s elapsed)
  ER       refinement done   (124 s elapsed)
  WS       refinement done   (139 s elapsed)
  BA       refinement done   (140 s elapsed)


-- lattice refinement sweep --------------------------------------------------
       T   P(large)      +/-SE    n_large  cond.attack chi=<s2>/<s>
------------------------------------------------------------------------------
  0.2800     0.0617     0.0098         37       0.0154        45.64
  0.3000     0.0967     0.0121         58       0.0164        52.65
  0.3200     0.2517     0.0177        151       0.0249        54.16
  0.3400     0.3583     0.0196        215       0.0435        55.58
  0.3600     0.5083     0.0204        305       0.1137        55.70
  0.3800     0.6100     0.0199        366       0.4355        53.66
  0.4000     0.6750     0.0191        405       0.7577        44.01
  0.4200     0.7400     0.0179        444       0.8685        37.60
  0.4400     0.7717     0.0171        463       0.9178        19.50
  0.4700     0.8267     0.0155        496       0.9539         4.99

-- ER refinement sweep -------------------------------------------------------
       T   P(large)      +/-SE    n_large  cond.attack chi=<s2>/<s>
------------------------------------------------------------------------------
  0.1100     0.0000     0.0000          0            -        11.21
  0.1300     0.0017     0.0017          1            -        22.29
  0.1500     0.0183     0.0055         11       0.0189        31.81
  0.1600     0.0367     0.0077         22       0.0232        35.69
  0.1700     0.0700     0.0104         42       0.0533        41.37
  0.1800     0.1233     0.0134         74       0.1107        27.51
  0.1900     0.1650     0.0152         99       0.2065        29.47
  0.2100     0.2650     0.0180        159       0.3818        24.49
  0.2300     0.3700     0.0197        222       0.4949        10.66
  0.2600     0.4317     0.0202        259       0.6201         8.89

-- WS refinement sweep -------------------------------------------------------
       T   P(large)      +/-SE    n_large  cond.attack chi=<s2>/<s>
------------------------------------------------------------------------------
  0.2400     0.0017     0.0017          1            -        26.76
  0.2700     0.0317     0.0071         19       0.0136        36.34
  0.2900     0.0550     0.0093         33       0.0188        39.89
  0.3100     0.0883     0.0116         53       0.0250        43.87
  0.3300     0.1817     0.0157        109       0.0677        43.76
  0.3500     0.2700     0.0181        162       0.1741        37.50
  0.3800     0.4417     0.0203        265       0.4921        44.10
  0.4100     0.5950     0.0200        357       0.6782        28.92
  0.4500     0.7050     0.0186        423       0.8277        27.33
  0.5000     0.8383     0.0150        503       0.9215        19.83

-- BA refinement sweep -------------------------------------------------------
       T   P(large)      +/-SE    n_large  cond.attack chi=<s2>/<s>
------------------------------------------------------------------------------
  0.0200     0.0000     0.0000          0            -         1.41
  0.0300     0.0000     0.0000          0            -         3.68
  0.0400     0.0000     0.0000          0            -         8.93
  0.0500     0.0000     0.0000          0            -        18.59
  0.0600     0.0067     0.0033          4       0.0171        13.64
  0.0700     0.0250     0.0064         15       0.0253        30.47
  0.0850     0.0400     0.0080         24       0.0552        15.30
  0.1000     0.0667     0.0102         40       0.0842        19.88
  0.1200     0.0967     0.0121         58       0.1501        17.81
  0.1500     0.1800     0.0157        108       0.2619         8.80

-- Threshold summary ---------------------------------------------------------
network      T_c meas        95% CI  T_c extrap  MolloyReed HMF <k>/<k2>   wellmixed   peak?
------------------------------------------------------------------------------
lattice        0.3511  [0.3078,0.3751]      0.3276      0.2000      0.1667      0.1667 interior
ER             0.1679  [0.1610,0.1705]      0.1639      0.1665      0.1428      0.1667 interior
WS             0.3741  [0.3013,0.3778]      0.3173      0.1981      0.1653      0.1667 interior
BA             0.0728  [0.0507,0.1198]      0.0665      0.0554      0.0525      0.1667 interior

"peak?" says whether the susceptibility maximum was interior to the grid.
An "EDGE!" there means the grid failed to bracket the peak and the number
to its left is a grid boundary, not a measurement. We print this because a
silent boundary fallback is exactly how Validation C first went wrong.

Differences, measured minus predicted:
network    vs MolloyReed         vs HMF  vs well-mixed
------------------------------------------------------------------------------
lattice           0.1511         0.1844         0.1844
ER                0.0014         0.0251         0.0012
WS                0.1760         0.2088         0.2074
BA                0.0174         0.0203        -0.0939

Ratios, measured divided by predicted (1.00 is perfect):
network    vs MolloyReed         vs HMF  vs well-mixed
------------------------------------------------------------------------------
lattice            1.755          2.106          2.106
ER                 1.008          1.176          1.007
WS                 1.889          2.263          2.245
BA                 1.314          1.387          0.437

-- The scale-free case in full, because it is the one the theory is loudest about 
  <k>                                  = 5.9988
  <k^2>                                = 114.2906
  <k^2>/<k>                            = 19.0522
  HMF prediction   T_c = <k>/<k^2>     = 0.05249
  Molloy-Reed      T_c = <k>/(<k^2>-<k>) = 0.05539
  Club measured    T_c                 = 0.07281   95% CI [0.05068, 0.11975]
  measured - HMF                       = +0.02032
  measured - Molloy-Reed               = +0.01741
  measured / well-mixed (1/6)          = 0.437

  Reading: the scale-free threshold sits at roughly one 2.3 th of the well-mixed
  value. A disease far too weak to spread in a well-mixed population of the
  same average contact rate still takes off on this network.

==============================================================================
PART 5.  VALIDATION C - THE LATTICE AGAINST AN EXACTLY KNOWN NUMBER
==============================================================================
With a CONSTANT infectious period of 5 steps, the transmission events on the
edges leaving one node become independent, and the SIR process is exactly bond
percolation with bond probability T = 1 - (1-beta)^5.
The bond percolation threshold of the triangular lattice is known in closed
form: p_c = 2 sin(pi/18) = 0.3472963553  (Sykes and Essam 1964).
This is the single hardest number in the study to fudge, so we run it.

A note on the cutoff, because it decides the answer. The susceptibility is
taken over runs that did NOT take off, so we must say what taking off means.
At p_c on a 2-D lattice the incipient cluster is a fractal of mass ~L^(91/48),
which for L = 100 is about 6190 nodes, or 62% of N. A 1%-of-N cutoff therefore
counts ordinary critical clusters as epidemics and censors the peak downwards.
For this validation we use 1000 nodes (10% of N) and print both.

       T      beta   P(large)    n_large  cond.attack      chi @1%     chi @10%
------------------------------------------------------------------------------
  0.2600   0.05844     0.0417         25       0.0142        45.25        71.79
  0.2800   0.06359     0.1567         94       0.0164        52.52       134.18
  0.3000   0.06885     0.3350        201       0.0290        55.76       329.03
  0.3150   0.07288     0.4650        279       0.0499        55.44       466.16
  0.3300   0.07697     0.5700        342       0.1522        47.68       498.24
  0.3450   0.08114     0.7233        434       0.4665        47.95       471.52
  0.3600   0.08539     0.8000        480       0.7376        38.67       486.05
  0.3750   0.08972     0.8533        512       0.8372        26.09        80.75
  0.4000   0.09712     0.8983        539       0.9015        19.97        19.97
  0.4300   0.10633     0.9317        559       0.9402         3.29         3.29

  With the 1%-of-N cutoff      p_c = 0.30546  (interior)
  With the 10%-of-N cutoff     p_c = 0.33068  (interior)

  Club measured p_c (constant period)  = 0.33068   95% CI [0.32181, 0.35612]
  Exact triangular bond percolation    = 0.34730
  measured - exact                     = -0.01661
  measured / exact                     = 0.9522
  bootstrap SE on the measurement      = 0.00875
  discrepancy in standard errors       = 1.90 SE
  Tolerance set before running         = 0.02000
  PASS

-- Is the residual gap a finite-lattice effect? Finite-size scaling ----------
Percolation theory says the pseudo-critical point of an L x L system sits at
p_c(L) = p_c(inf) + a * L^(-1/nu) with nu = 4/3 in two dimensions, so the gap
should shrink like L^(-3/4). If our simulator is right, p_c(L) must march
toward 0.34730 as L grows. If it is wrong, it will march somewhere else.
Same constant-period process, 300 runs per point, cutoff at 10% of each N.

     L        N      p_c(L) minus exact     peak?
------------------------------------------------------------------------------
    40     1600     0.31787    -0.02942  interior
    70     4900     0.32222    -0.02508  interior
   100    10000     0.32769    -0.01960  interior
   140    19600     0.33145    -0.01585  interior

  Fit p_c(L) = p_inf + a * L^(-3/4):  a = -0.34650, p_inf = 0.33870
  Exact value                                        = 0.34730
  extrapolated p_inf minus exact                     = -0.00859
  The gap at each L shrinks with L and the slope a is negative, which is the
  sign percolation predicts. The simulator is measuring the right critical
  point; a 100 x 100 lattice is simply too small to sit on it.

  For comparison, the GEOMETRIC-period lattice threshold measured in Part 4
  was 0.35107. The two differ by +0.02039. That gap is real and it is the price of
  a variable infectious period: transmissions out of one node are correlated
  through the shared period, so the process is no longer exactly percolation.
  We report it rather than quietly picking whichever number matched better.

  Tree-like theory would have put the lattice threshold at <k>/(<k^2>-<k>) =
  0.20000. It is out by a factor of 1.65. Clustering is why: on this lattice
  most of a node's neighbours are neighbours of each other, so a chain of
  infection keeps running into people it has already met.

==============================================================================
PART 6.  MONTE CARLO CONVERGENCE AND UNCERTAINTY
==============================================================================
Everything above is an average over independent runs, so every number carries
a Monte Carlo error that shrinks like 1/sqrt(n). We check that directly: one
configuration, 1500 runs, running mean with a 95% interval as trials accumulate.

Configuration: ER network, T = 0.28 (R0 well-mixed equivalent = 1.68).

Running mean of the conditional attack fraction (large outbreaks only):
 n_large         mean           SE        CI_lo        CI_hi     CI width
------------------------------------------------------------------------------
      10      0.67852      0.00598      0.66680      0.69024      0.02344
      20      0.67870      0.00342      0.67199      0.68540      0.01342
      40      0.67993      0.00198      0.67605      0.68381      0.00776
      60      0.67957      0.00149      0.67666      0.68248      0.00583
     100      0.68100      0.00115      0.67875      0.68324      0.00449
     150      0.68144      0.00090      0.67969      0.68320      0.00351
     200      0.68170      0.00077      0.68019      0.68322      0.00304
     300      0.68191      0.00065      0.68064      0.68318      0.00255
     400      0.68162      0.00056      0.68053      0.68271      0.00218
     500      0.68154      0.00050      0.68057      0.68251      0.00194
     600      0.68147      0.00046      0.68057      0.68236      0.00179
     775      0.68113      0.00041      0.68033      0.68193      0.00160  <- all runs

Running estimate of P(large):
  n_runs     P(large)           SE     CI width
------------------------------------------------------------------------------
      10      0.60000      0.15492      0.60728
      20      0.55000      0.11124      0.43607
      40      0.55000      0.07866      0.30835
      60      0.56667      0.06397      0.25078
     100      0.54000      0.04984      0.19537
     150      0.54000      0.04069      0.15952
     200      0.54000      0.03524      0.13815
     300      0.56667      0.02861      0.11215
     400      0.53750      0.02493      0.09772
     500      0.52200      0.02234      0.08757
     600      0.52333      0.02039      0.07993
     800      0.51875      0.01767      0.06925
    1000      0.51400      0.01581      0.06196
    1200      0.52167      0.01442      0.05653
    1500      0.51667      0.01290      0.05058  <- all runs

Check that the error really falls as 1/sqrt(n): SE(n) * sqrt(n) should be flat
and equal to the sample standard deviation, 0.01134.
   n =   50   SE*sqrt(n) = 0.01194
   n =  100   SE*sqrt(n) = 0.01146
   n =  200   SE*sqrt(n) = 0.01095
   n =  400   SE*sqrt(n) = 0.01111

Distribution of final sizes at this configuration (all 1500 runs):
        0 -      1 :   432 runs (0.288)
        2 -      4 :   221 runs (0.147)
        5 -      9 :    46 runs (0.031)
       10 -     24 :    25 runs (0.017)
       25 -     49 :     1 runs (0.001)
     6000 -   6999 :   736 runs (0.491)
     7000 -   7999 :    39 runs (0.026)

The bimodality is the whole reason we quote outbreak probability and
conditional size as two separate numbers. Averaging across the gap between
the fizzles and the epidemics would produce a figure that describes neither.

-- Sensitivity of the answers to the large-outbreak cutoff -------------------
The cutoff is a judgement call. Here is what moves if we change it.

network       cutoff     P(large)  cond.attack     T_c meas
------------------------------------------------------------------------------
lattice           20       0.8450       0.9722      0.32275
lattice           50       0.8450       0.9722      0.29407
lattice          100       0.8450       0.9722      0.35107
lattice          200       0.8450       0.9722      0.35597
lattice          500       0.8450       0.9722      0.35569
ER                20       0.8367       0.9410      0.17083
ER                50       0.8367       0.9410      0.16921
ER               100       0.8367       0.9410      0.16791
ER               200       0.8367       0.9410      0.16116
ER               500       0.8367       0.9410      0.17817
WS                20       0.8333       0.9081      0.35577
WS                50       0.8217       0.9209      0.33035
WS               100       0.8200       0.9228      0.37408
WS               200       0.8200       0.9228      0.34141
WS               500       0.8200       0.9228      0.33998
BA                20       0.8283       0.9198      0.15000
BA                50       0.8283       0.9198      0.07245
BA               100       0.8283       0.9198      0.07281
BA               200       0.8283       0.9198      0.06427
BA               500       0.8283       0.9198      0.08593

==============================================================================
PART 7.  WHAT THE NUMBERS SAY, IN ONE PLACE
==============================================================================
At T = 0.50, identical disease, identical mean number of contacts:

network         attack frac       peak I/N  t_peak(steps)    fizzle frac
------------------------------------------------------------------------------
lattice              0.9722        0.06021          136.6         0.1550
ER                   0.9410        0.40697           20.2         0.1633
WS                   0.9228        0.14528           61.3         0.1800
BA                   0.9198        0.44060           12.8         0.1717
well-mixed           0.9405   (see part 2)              -              -

Spread across the four networks at T = 0.50:
   attack fraction : 0.9198 to 0.9722, ratio 1.06
   peak prevalence : 0.06021 to 0.44060, ratio 7.32
   time to peak    : 12.8 to 136.6 steps, ratio 10.69

Thresholds, lowest to highest:
   BA        T_c = 0.0728   [0.0507, 0.1198]   = 0.44 x the well-mixed value
   ER        T_c = 0.1679   [0.1610, 0.1705]   = 1.01 x the well-mixed value
   lattice   T_c = 0.3511   [0.3078, 0.3751]   = 2.11 x the well-mixed value
   WS        T_c = 0.3741   [0.3013, 0.3778]   = 2.24 x the well-mixed value
   spread from lowest to highest: factor 5.14

==============================================================================
Total runtime: 217.4 s
Seed: 20251215. Every number above is reproducible by rerunning this file.
==============================================================================