example_run.py

#!/usr/bin/env python3
# pyVIDE -- pure-python re-implementation of the VIDE/ZOBOV void finder
# Copyright (C) 2026 Nico Schuster
#
# An independent implementation of the algorithms of VIDE
# (Copyright (C) 2010-2025 Guilhem Lavaux, 2011-2014 P. M. Sutter)
# and ZOBOV (Mark Neyrinck), reproducing VIDE's output exactly.
#
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; version 2 of the License.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, see <https://www.gnu.org/licenses/>.

"""A complete, runnable pyVIDE session -- edit ONE block and it is your run.

Run it as it stands and it uses the public reference catalog if you have
downloaded it (see REFERENCE_CATALOG below), else it invents a synthetic box, so
you can watch the whole pipeline work end to end in a few minutes::

    python3 example_run.py

Then replace the marked block near the top with your own catalog.  Everything
below that block is generic and does not need editing.

What it demonstrates, in order:

1. ``identify_voids`` with the options that actually matter, spelled out;
2. reading the catalog back at all three tiers with ``load_catalog``;
3. void properties as plain numpy arrays;
4. ``members`` -- the tracers of one void, and their properties;
5. ``filter_catalog`` -- reproducing a VIDE catalog variant;
6. ``rerun_watershed`` -- new weights without re-tessellating;
7. ``rebuild_catalog`` -- a different merging threshold in milliseconds;
8. ``identify_voids_from_field`` -- the same finder on a density grid.

Nothing here writes outside ``OUTPUT_DIR``.
"""

import os

import numpy as np

import pyvide

# ===========================================================================
#  EDIT THIS BLOCK -- everything below it is generic
# ===========================================================================

OUTPUT_DIR = "example_output"
RUN_NAME = "example"

#: Box length in Mpc/h.  A scalar means a cube; use (Lx, Ly, Lz) otherwise.
BOX_LEN = 200.0

#: Per-axis periodicity.  True / False / (True, True, False) / 'xy'.
PERIODIC = True

#: Which VIDE behaviour to reproduce.  'VIDE_halos' and 'VIDE_direct' are the
#: VIDE-compatible float32 paths; 'pyvide' is pyVIDE's own float64 mode.
#: There is no default on purpose -- see pyvide.explain('VIDE_mode').
VIDE_MODE = "pyvide"

#: The public reference catalog of the pyVIDE data record on Zenodo
#: (10.5281/zenodo.22736465): tracers/magneticum_box2b_hr_z0.25_m3e10.txt,
#: 1 214 928 subhalos in a 640 Mpc/h box, the sample every number in the
#: README's verification section refers to.  Download it and put the path
#: here; the script then uses it (with BOX_LEN 640 and VIDE_mode 'VIDE_halos',
#: as VIDE was run) whenever load_my_tracers() has not been edited.  Budget
#: 15 minutes on one core for the whole script; set NUM_THREADS below to use
#: more.
REFERENCE_CATALOG = "magneticum_box2b_hr_z0.25_m3e10.txt"
NUM_THREADS = 1


def load_my_tracers():
    """Return positions, and optionally velocities / weights / extra columns.

    REPLACE THE BODY OF THIS FUNCTION.  For example::

        table = np.loadtxt("my_halos.txt")
        return dict(
            positions=table[:, 1:4],                  # (N, 3) in Mpc/h
            velocities=table[:, 4:7],                  # (N, 3) in km/s, optional
            weights=table[:, 0],                       # (N,) positive, optional
            otherProperties={"mass": table[:, 0]},     # anything per tracer
        )

    If you do not, the script looks for the public reference catalog
    (REFERENCE_CATALOG above) and uses it when the file exists.  Failing both,
    the fallback below is a synthetic box: a uniform background with a few
    spherical underdensities carved into it, so the run finds something.
    """
    global BOX_LEN, VIDE_MODE
    if os.path.exists(REFERENCE_CATALOG):
        # columns: mass, x, y, z, vx, vy, vz -- the file VIDE's prepareInputs
        # was fed; BOX_LEN and VIDE_MODE follow the VIDE run it is compared with
        table = np.loadtxt(REFERENCE_CATALOG, dtype=np.float64)
        BOX_LEN = 640.0
        VIDE_MODE = "VIDE_halos"
        print("data: the public reference catalog", REFERENCE_CATALOG)
        print("      -> BOX_LEN {0} and VIDE_mode {1!r} adopted from the VIDE "
              "run of this sample; your settings above are overridden".format(
                  BOX_LEN, VIDE_MODE))
        return dict(
            positions=table[:, 1:4],
            velocities=table[:, 4:7],
            weights=None,
            otherProperties={"mass": table[:, 0]},
        )
    #
    print("data: a synthetic box (edit load_my_tracers, or download the "
          "reference catalog named by REFERENCE_CATALOG)")
    rng = np.random.RandomState(20260818)
    numTracers = 40000
    positions = rng.uniform(0.0, BOX_LEN, size=(numTracers, 3))
    #
    # carve four spherical holes so the catalog has recognisable voids
    holes = [(50.0, 50.0, 50.0, 28.0), (150.0, 60.0, 40.0, 24.0),
             (60.0, 150.0, 150.0, 30.0), (160.0, 160.0, 120.0, 22.0)]
    keep = np.ones(numTracers, dtype=bool)
    for x, y, z, radius in holes:
        offset = np.abs(positions - np.array([x, y, z]))
        offset = np.minimum(offset, BOX_LEN - offset)      # periodic distance
        distance = np.sqrt((offset ** 2).sum(axis=1))
        # keep 15 % of the tracers inside each hole, so it is under- not empty
        inside = distance < radius
        keep &= ~(inside & (rng.rand(numTracers) > 0.15))
    #
    positions = positions[keep]
    masses = 10.0 ** rng.uniform(11.0, 14.0, size=positions.shape[0])
    return dict(
        positions=positions,
        velocities=None,
        weights=None,
        otherProperties={"mass": masses},
    )


# ===========================================================================
#  GENERIC FROM HERE ON
# ===========================================================================


def banner(text):
    print("\n" + "=" * 74)
    print(text)
    print("=" * 74)


def main():
    data = load_my_tracers()
    positions = data["positions"]
    numTracers = positions.shape[0]
    #
    # BOX_LEN is a scalar or (Lx, Ly, Lz); the three lines of this script that
    # need the sides one by one go through boxSides, so both spellings work.
    # The finder itself takes BOX_LEN as it stands.
    boxSides = np.broadcast_to(np.asarray(BOX_LEN, dtype=float), (3,))
    boxVolume = float(boxSides[0]) * float(boxSides[1]) * float(boxSides[2])
    separation = (boxVolume / numTracers) ** (1.0 / 3.0)

    banner("1. find the voids")
    print("{0} tracers, mean separation {1:.2f} Mpc/h".format(
        numTracers, separation))
    #
    # buffer: VIDE hard-codes 0.1, which on a VIDE-sized run is ~10 mean
    # tracer separations -- but a buffer thinner than ~3 separations silently
    # gives the tracers near a sub-box face the wrong Voronoi volumes.  In
    # 'pyvide' the default (buffer=None) sizes it from the catalog:
    # max(10 x mean separation, 0.02), capped below the whole-catalog regime.
    # The resolved value lands in catalog.parameters['buffer'] and the log.
    # In 'VIDE_halos'/'VIDE_direct' the default stays VIDE's 0.1: size it yourself there.
    #
    catalog = pyvide.identify_voids(
        saveDir=OUTPUT_DIR,
        saveName=RUN_NAME,
        boxLen=BOX_LEN,
        positions=positions,
        velocities=data.get("velocities"),
        weights=data.get("weights"),
        otherProperties=data.get("otherProperties"),
        VIDE_mode=VIDE_MODE,
        periodicBox=PERIODIC,
        numDivisions=2,          # >= 2 on every periodic axis in the VIDE modes
        numThreads=NUM_THREADS,  # sub-boxes in parallel; bit-identical to 1
        buffer=None,             # adaptive in 'pyvide' (see above)
        mergingThreshold=1.0e-9,  # 1e-9 = no merging (VIDE's catalogs);
                                  # 0.2 = classical ZOBOV; 0 = one hierarchy
        minRadius=None,           # None = keep every basin (pyVIDE default)
        maxCentralDen=None,       # None = compute centralDen, flag nothing
        zones=True,               # also write the pre-merging basins
        saveIntermediate=True,    # cache the tessellation for step 6
        overwrite=True,
        verbose=True,
    )
    print("\n->", catalog)

    banner("2. read it back")
    loaded = pyvide.load_catalog(
        OUTPUT_DIR, RUN_NAME,
        loadTracers=True,        # per-tracer positions, volumes, zone ids
        loadMembership=True,     # void -> zone -> tracer lists
        loadZones=True,          # the pre-merging basins
    )
    print("columns:", ", ".join(loaded.columns()))
    print("statistics:")
    for key in sorted(loaded.statistics):
        print("   {0:<22}: {1}".format(key, loaded.statistics[key]))

    banner("3. void properties are plain numpy arrays")
    biggest = np.argsort(loaded.radius)[::-1][:5]
    print("{0:>6} {1:>8} {2:>8} {3:>26} {4:>8}".format(
        "voidID", "R [Mpc/h]", "numPart", "macrocenter [Mpc/h]", "ellip"))
    for row in biggest:
        print("{0:>6} {1:>9.2f} {2:>8} {3:>26} {4:>8.3f}".format(
            loaded.voidID[row], loaded.radius[row], loaded.numPart[row],
            np.array2string(loaded.macrocenter[row], precision=1),
            loaded.ellipticity[row]))

    banner("4. the tracers of one void")
    # NOTE the voidID, not the row number.  They differ after filtering and on
    # a walled box, so always go through the voidID column.
    target = int(loaded.voidID[biggest[0]])
    index = pyvide.members(loaded, target)
    print("void {0} has {1} member tracers".format(target, index.size))
    print("their mean Voronoi volume :", loaded.tracers["volume"][index].mean())
    if "mass" in loaded.tracers:
        print("their median halo mass    : {0:.3e}".format(
            np.median(loaded.tracers["mass"][index])))
    print("first five positions:\n", loaded.tracers["positions"][index[:5]])

    banner("5. reproduce a VIDE catalog variant")
    # pyVIDE keeps every density basin and flags; VIDE deleted as it went.
    # filter_catalog applies VIDE's filters, in VIDE's order, to a pyVIDE
    # catalog.  maxCentralDen is what the 'dencut' variants cut on -- VIDE fed
    # its single mergingThreshold in here, which is the coupling pyVIDE undoes.
    # The two variants without a central-density cut ignore the number, so
    # only the other two are given one.
    for variant, denCut in (("untrimmed_all", {}),
                            ("trimmed_nodencut_all", {}),
                            ("untrimmed_dencut_all", {"maxCentralDen": 1.0e-9}),
                            ("all", {"maxCentralDen": 1.0e-9})):
        subset = pyvide.filter_catalog(
            loaded, variant, minRadius=0.0, **denCut
        )
        print("{0:<24} {1:>6} of {2} voids".format(
            variant, subset.numVoids, loaded.numVoids))

    banner("6. new weights without re-tessellating (rerun_watershed)")
    weights = loaded.tracers.get("mass")
    if weights is None:
        weights = np.ones(loaded.tracers["volume"].size)
    weighted = pyvide.rerun_watershed(
        OUTPUT_DIR, RUN_NAME, newSaveName=RUN_NAME + "_weighted",
        weights=weights, overwrite=True, verbose=False,
    )
    print("mass-weighted run: {0} voids (unweighted: {1})".format(
        weighted.numVoids, loaded.numVoids))
    print("the tessellation was reused -- only the density field changed.")
    print("(the synthetic masses here are random and uncorrelated with")
    print(" position, so they act as noise and fragment the field; real mass")
    print(" or luminosity weights are spatially coherent and do not.)")

    banner("7. a different merging threshold (rebuild_catalog)")
    #   for a sweep, pass outputs="voids" here: the tracer file a rebuild
    #   writes is a copy of the one it read
    for threshold in (0.2, 0.0):
        rebuilt = pyvide.rebuild_catalog(
            OUTPUT_DIR, RUN_NAME, mergingThreshold=threshold,
            overwrite=True, verbose=False,
        )

        print("mergingThreshold={0:<5} {1:>5} voids, largest holds {2} zones, "
              "{3} top-level".format(
                  threshold, rebuilt.numVoids, rebuilt.numZones.max(),
                  int(rebuilt.isTopLevel.sum())))

    banner("8. the same finder on a density grid")
    grid = 32
    counts, _ = np.histogramdd(
        positions, bins=(grid, grid, grid),
        range=[(0.0, side) for side in boxSides],
    )
    #
    # SMOOTH THE FIELD.  A raw count grid is shot-noise dominated: almost
    # every voxel is a local minimum of its own, so the watershed returns
    # thousands of one-voxel "voids" and none of them mean anything.  Smoothing
    # on a scale of a few voxels is not cosmetic -- it is what makes the
    # density field a density field.  (Without the filter below this same run
    # returns ~11600 voids from 32768 voxels; with it, around two hundred.)
    try:
        from scipy.ndimage import gaussian_filter
        field = gaussian_filter(counts, sigma=1.5, mode="wrap") + 0.05
        print("field smoothed with a 1.5-voxel Gaussian "
              "({0:.1f} Mpc/h)".format(1.5 * float(boxSides.min()) / grid))
    except ImportError:                                   # pragma: no cover
        field = counts + 0.5
        print("scipy.ndimage unavailable -- using the raw counts, expect noise")
    #
    fieldCatalog = pyvide.identify_voids_from_field(
        OUTPUT_DIR, RUN_NAME + "_field", BOX_LEN, field,
        connectivity=6, periodicBox=PERIODIC, overwrite=True, verbose=False,
    )
    print("{0}^3 grid -> {1} voids".format(grid, fieldCatalog.numVoids))
    biggestCell = np.argsort(fieldCatalog.radius)[::-1][:3]
    for row in biggestCell:
        print("   R = {0:6.2f} Mpc/h, {1:6d} voxels, centre {2}".format(
            fieldCatalog.radius[row], fieldCatalog.numCells[row],
            np.array2string(fieldCatalog.macrocenter[row], precision=1)))

    banner("done")
    print("everything written to {0}/".format(os.path.abspath(OUTPUT_DIR)))
    print("try pyvide.explain('buffer') or pyvide.explain('VIDE_mode') next.")


if __name__ == "__main__":
    main()