Notebooks · lib/hf/sweep

Sweeping array-valued strategies, and blaming the one that dissents

Every stage of the Hartree-Fock code has more than one implementation: McMurchie-Davidson, Obara-Saika and Rys for the electron repulsion integrals, in-core and direct J/K builds, several orthogonalizers. In debug mode a stage runs all of them on the same input and requires them to agree. Agreement between algorithms that share no code is strong evidence each one is right, and a disagreement names the one that is wrong.

The drivers here do that for strategies that return arrays. Each member of a registry runs, its result is flattened into one row of a table, and the consensus is the elementwise median of the rows. The member farthest from the median at any entry is blamed, together with the flat (row-major) index of that entry. The median matters: if the consensus were member 0's result, a bad member 0 would put the blame on every good member instead.

Deviation is relative, abs(v - median) / (1 + abs(median)), so a tolerance reads as relative for large values and absolute for small ones.

Only the first ncmp entries are compared. An SCF result vector carries the energy and densities first and then the iteration count; two correct accelerators converge in different numbers of iterations, and comparing that tail would fail on correct code. Past ncmp the result carries member 0's values.

release_mode() runs member 0 only. The sw_check* variants never assert: they return 1 when the members agree and 0 when they do not, so a test can check that a planted bug is caught.

A Viv def is typed by the one registry signature it is called with (finding F4 of the contracts), so there is one driver per signature:

sweep3_arr   members (Arr, Arr, Arr) -> Arr            integrals, transforms, eigen
sweep5_arr   members (Arr, Arr, Arr, Arr, Arr) -> Arr  J/K
sweep3_num   members (Arr, Arr, Arr) -> number
sweepnx_arr  members (nmax, x) -> Arr                   Boys

N is the number of entries a member returns; Viv has no element-count builtin. rep receives [strategy, flat index, worst, nstrat, lo, hi]; lo and hi are reserved for the interval wave and stay 0.

The coverage counters at the end record which integral classes, Boys regions and SCF regimes each strategy has actually been exercised on, so a build can fail when a class was only ever checked by one strategy.

Used with import hf.sweep.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 51–746

SW_REP_STRAT = 0
SW_REP_INDEX = 1
SW_REP_WORST = 2
SW_REP_NSTRAT = 3
SW_REP_LO = 4
SW_REP_HI = 5
SW_REP_LEN = 6

SW_COV_NSTRAT = 16
SW_LMAX1 = 5
SW_NCLS_ERI = 625
SW_NXBIN = 8
SW_NCLS_BOYS = 136
SW_NCLS_SCF = 16

# the test tensor: SW_T_DIM^4 entries, standing in for an ERI block
SW_T_DIM = 3
SW_T_LEN = SW_T_DIM * SW_T_DIM * SW_T_DIM * SW_T_DIM
# the J/K test members take an options vector of lib/hf/opts's length (OPT_LEN),
# which this module does not use
SW_T_OPT_LEN = 32

#| A fresh report vector.
def sw_rep_new() {
    return zeros(SW_REP_LEN)
}

#| The flat indices `0 .. N-1`, for `take` and `put`.
def sw_flat_index(N) {
    ix = zeros(N)
    for i to N {
        ix[i] = i
    }
    return ix
}

#| Elementwise median of the rows of `T` (`nstrat` by `N`) into `C`, over the
#| first `ncmp` columns; later columns copy row `0` and are never compared.
#| Returns the worst relative deviation and records the member and flat index
#| that produced it in `rep`. With an even count the upper middle value is
#| the median, the same choice `sweep` in lang/viv/std/contracts.viv makes.
def sw_table_median(T: Arr, nstrat, N, ncmp, C: Arr, rep: Arr) {
    col = zeros(nstrat)
    worst = 0.0
    rep[SW_REP_STRAT] = 0
    rep[SW_REP_INDEX] = 0
    for i to N {
        if i < ncmp {
            for m to nstrat {
                col[m] = T[m, i]
            }
            sort(col)
            C[i] = col[floor(nstrat / 2.0)]
            for m to nstrat {
                dev = abs(T[m, i] - C[i]) / (1.0 + abs(C[i]))
                if dev > worst {
                    worst = dev
                    rep[SW_REP_STRAT] = m
                    rep[SW_REP_INDEX] = i
                }

SW_REP_STRAT = 0
SW_REP_INDEX = 1
SW_REP_WORST = 2
SW_REP_NSTRAT = 3
SW_REP_LO = 4
SW_REP_HI = 5
SW_REP_LEN = 6

SW_COV_NSTRAT = 16
SW_LMAX1 = 5
SW_NCLS_ERI = 625
SW_NXBIN = 8
SW_NCLS_BOYS = 136
SW_NCLS_SCF = 16

# the test tensor: SW_T_DIM^4 entries, standing in for an ERI block
SW_T_DIM = 3
SW_T_LEN = SW_T_DIM * SW_T_DIM * SW_T_DIM * SW_T_DIM
# the J/K test members take an options vector of lib/hf/opts's length (OPT_LEN),
# which this module does not use
SW_T_OPT_LEN = 32

#| A fresh report vector.
def sw_rep_new() {
    return zeros(SW_REP_LEN)
}

#| The flat indices `0 .. N-1`, for `take` and `put`.
def sw_flat_index(N) {
    ix = zeros(N)
    for i to N {
        ix[i] = i
    }
    return ix
}

#| Elementwise median of the rows of `T` (`nstrat` by `N`) into `C`, over the
#| first `ncmp` columns; later columns copy row `0` and are never compared.
#| Returns the worst relative deviation and records the member and flat index
#| that produced it in `rep`. With an even count the upper middle value is
#| the median, the same choice `sweep` in lang/viv/std/contracts.viv makes.
def sw_table_median(T: Arr, nstrat, N, ncmp, C: Arr, rep: Arr) {
    col = zeros(nstrat)
    worst = 0.0
    rep[SW_REP_STRAT] = 0
    rep[SW_REP_INDEX] = 0
    for i to N {
        if i < ncmp {
            for m to nstrat {
                col[m] = T[m, i]
            }
            sort(col)
            C[i] = col[floor(nstrat / 2.0)]
            for m to nstrat {
                dev = abs(T[m, i] - C[i]) / (1.0 + abs(C[i]))
                if dev > worst {
                    worst = dev
                    rep[SW_REP_STRAT] = m
                    rep[SW_REP_INDEX] = i
                }