Anchor Sampling AS

Pronunciation note: say AS like the French as 🂡, meaning the ace of a deck, or spell it A, S. Please do not pronounce it like “ass”.

Anchor Sampling (AS) is a location-anchored variant of QuickSampling. It preserves the same path-based simulation workflow, neighborhood template logic, continuous/categorical handling, multivariate support, kernels, and random or explicit simulation paths. The core conceptual change is the sampling constraint:

• In QS, once a neighborhood match is found, the sampled center can come from any admissible location in the training image(s).
• In AS, when simulating a cell at coordinate (i,j,...), the sampled center value must come from that exact same coordinate across an aligned stack of training images.

This makes AS suitable when you have many realizations, scenarios, priors, or historical fields defined on the same grid geometry and you want neighborhood-based selection without losing spatial anchoring.

Parameters for AS

Usage: [sim,index,time,finalprogress,jobid] = g2s(flag1,value1, flag2,value2, ...) Outputs: sim = simulation, index = selected TI id for each simulated value, time = simulation time, finalprogress = final progression of the simulation (normally 100), jobid = job ID.

Flag	Description	Mandatory
-a	Use as, AS, or AnchorSampling.	✔
-ti	Training images provided as an aligned stack. All TIs must have exactly the same spatial size as the simulation grid. In multivariate cases, the last dimension still represents the variables. NaN TI values are ignored.	✔
-di	Destination image (simulation grid). NaN values are simulated. Non-NaN values are conditioning data. Its spatial size must match every TI.	✔
-dt	Data type. 0 → continuous and 1 → categorical. Provide one value per variable.	✔
-k	Number of best anchored TI candidates retained before sampling.	✔
-n	Number of closest neighbors to consider. As in QS, this can be scalar or variable-specific.	✔
-ki	Optional kernel image defining the neighborhood footprint and weights.
-cnorm / -cn	Continuous mismatch norm power(s). Use one strictly positive value to apply the same p to all continuous variables, or one value per continuous variable. Continuous mismatch uses the rooted Lp form pow(sum(pow(abs(diff),p)),1/p) (kernel-weighted internally). Default is 2 (L2).
-sp	Optional simulation path. Default is a random path. As in QS, explicit path images are supported.
-j	Parallel execution settings. Same meaning as in QS.
-s	Random seed value.
-mi	Optional prior or mask image with one value per TI at each spatial location. In Python and MATLAB this is typically a stack with trailing TI axis (spatial..., n_ti). Mismatch is computed first, then -mi is applied as an admissibility mask before top-k: non-finite (NaN/inf) or non-positive mask values are excluded. The retained top-k candidates are finally sampled using -mi as relative weights.

Advanced parameters

Flag	Description	Mandatory
-kii	Per-pixel kernel-selection image, same semantics as QS.
-ni	Per-pixel number of neighbors.
-kvi	Per-pixel number of retained best candidates.
-cti	Treat each TI as periodic over each dimension when evaluating TI-side neighborhoods.
-csim	Create a periodic simulation over each dimension.
-fs	Full simulation with distinct paths per variable, same semantics as QS.

Important differences from QS

• The index image exported by AS stores the selected TI id, not a flattened TI pixel address.
• -ii is intentionally not supported in AS. Use -mi instead when you want to constrain or weight which TI can be selected at a given location.
• -adsim is not supported in AS. AS requires all TIs and the simulation grid to share the same geometry.
• AS is CPU-first. GPU flags are accepted for interface compatibility, but the anchored matcher currently runs on CPU.

Examples

The following examples assume the G2S server is running.

Minimal Python example

#This code requires the G2S server to be running

from argparse import ArgumentParser
from pathlib import Path
import sys

import matplotlib.pyplot as plt
import numpy as np


def load_g2s():
    try:
        from g2s import g2s

        return g2s
    except ImportError as first_error:
        repo_root = Path(__file__).resolve().parents[2]
        python_build = repo_root / "build" / "python-build"
        if python_build.exists():
            sys.path.insert(0, str(python_build))
            try:
                from g2s import g2s

                return g2s
            except ImportError as second_error:
                raise ImportError(
                    "Unable to import the Python g2s client. "
                    "Build/install the Python interface and ensure pyzmq is available."
                ) from second_error
        raise ImportError(
            "Unable to import the Python g2s client. "
            "Build/install the Python interface first."
        ) from first_error


def build_synthetic_case(height=100, width=120, num_ti=7, keep_ratio=0.05, seed=4):
    y, x = np.mgrid[0:height, 0:width].astype(np.float32)
    x_norm = (x - 0.5 * (width - 1)) / (0.5 * width)
    y_norm = (y - 0.5 * (height - 1)) / (0.5 * height)

    tis = []
    for ti_id in range(num_ti):
        stretch_x = 0.9 + 0.10 * ti_id
        stretch_y = 1.15 - 0.05 * ti_id
        radius = np.sqrt((x_norm / stretch_x) ** 2 + (y_norm / stretch_y) ** 2)
        disk = 1.0 / (1.0 + np.exp((radius - (0.24 + 0.055 * ti_id)) / 0.035))
        stripes = 0.5 + 0.5 * np.sin(6.0 * x_norm + 3.5 * y_norm + 0.55 * ti_id)
        plume = 0.5 + 0.5 * np.cos(4.0 * radius - 0.35 * ti_id + 0.8 * x_norm)
        tis.append((0.60 * disk + 0.25 * stripes + 0.15 * plume).astype(np.float32))
    tis = np.stack(tis, axis=0)

    x01 = x / max(width - 1, 1)
    y01 = y / max(height - 1, 1)
    truth_index = np.clip(
        np.rint(
            (num_ti - 1)
            * (0.10 + 0.78 * x01 + 0.08 * np.sin(2.0 * np.pi * y01) - 0.04 * np.cos(2.5 * np.pi * x01 * y01))
        ),
        0,
        num_ti - 1,
    ).astype(np.int32)

    row = np.arange(height)[:, None]
    col = np.arange(width)[None, :]
    truth = tis[truth_index, row, col].astype(np.float32)

    rng = np.random.default_rng(seed)
    conditioning = np.full_like(truth, np.nan, dtype=np.float32)
    known = rng.random(truth.shape) < keep_ratio
    known[:4, :] = True
    known[-4:, :] = True
    known[:, :4] = True
    known[:, -4:] = True
    conditioning[known] = truth[known]

    mask = np.full((height, width, num_ti), np.nan, dtype=np.float32)
    rr, cc = np.indices(truth_index.shape)
    mask[rr, cc, truth_index] = 1.0

    return tis, truth, truth_index, conditioning, mask


def run_as(g2s, tis, conditioning, mask, k, n, seed, jobs):
    args = ["-a", "as"]
    for ti in tis:
        args.extend(["-ti", ti.astype(np.float32, copy=False)])
    args.extend(
        [
            "-di",
            conditioning.astype(np.float32, copy=False),
            "-mi",
            mask.astype(np.float32, copy=False),
            "-dt",
            np.array([0.0], dtype=np.float32),
            "-k",
            int(k),
            "-n",
            int(n),
            "-s",
            int(seed),
            "-j",
            int(jobs),
            "-silent",
        ]
    )
    simulation, selected_ti, *_ = g2s(*args)
    return simulation, selected_ti.astype(np.int64, copy=False)


def rmse(reference, estimate, valid_mask):
    if not np.any(valid_mask):
        return float("nan")
    delta = estimate[valid_mask] - reference[valid_mask]
    return float(np.sqrt(np.mean(delta * delta)))


def show_minimal_figure(truth, conditioning, simulation, selected_ti):
    conditioning_display = np.nan_to_num(conditioning, nan=-0.2)

    fig, axes = plt.subplots(1, 4, figsize=(12, 3), sharex=True, sharey=True)
    fig.suptitle("Anchor Sampling: minimal synthetic demo")

    for ax, image, title in zip(
        axes,
        [truth, conditioning_display, simulation, selected_ti],
        ["Ground truth", "Conditioning", "AS result", "Selected TI id"],
    ):
        ax.imshow(image, cmap="viridis")
        ax.set_title(title)
        ax.axis("off")
    plt.tight_layout()
    plt.show()


def main():
    parser = ArgumentParser(description="Minimal Anchor Sampling synthetic example.")
    parser.add_argument("--height", type=int, default=100)
    parser.add_argument("--width", type=int, default=120)
    parser.add_argument("--num-ti", type=int, default=7)
    parser.add_argument("--keep-ratio", type=float, default=0.05)
    parser.add_argument("--k", type=int, default=None, help="defaults to num-ti")
    parser.add_argument("--n", type=int, default=24)
    parser.add_argument("--seed", type=int, default=17)
    parser.add_argument("--jobs", type=int, default=1)
    args = parser.parse_args()

    g2s = load_g2s()
    tis, truth, truth_index, conditioning, mask = build_synthetic_case(
        height=args.height,
        width=args.width,
        num_ti=args.num_ti,
        keep_ratio=args.keep_ratio,
        seed=args.seed,
    )
    simulation, selected_ti = run_as(
        g2s,
        tis,
        conditioning,
        mask,
        k=args.num_ti if args.k is None else args.k,
        n=args.n,
        seed=args.seed,
        jobs=args.jobs,
    )

    unknown = np.isnan(conditioning)
    agreement = float(np.mean(selected_ti[unknown] == truth_index[unknown]))
    print(f"Unknown-cell RMSE: {rmse(truth, simulation, unknown):.4f}")
    print(f"Unknown-cell TI-index agreement: {agreement:.4f}")

    show_minimal_figure(truth, conditioning, simulation, selected_ti)


if __name__ == "__main__":
    main()

%This code requires the G2S server to be running

% Pronunciation note:
% say AS like the French "as" (ace of a deck), or spell it "A, S".

ti1=ones(21,21);
ti2=2*ones(21,21);
ti1(11,11)=5;
ti2(11,11)=9;

conditioning=nan(21,21,1);
conditioning(10,11)=1;
conditioning(11,10)=1;
conditioning(11,12)=1;
conditioning(12,11)=1;

path=100+reshape(0:numel(conditioning(:,:,1))-1,21,21);
path(11,11)=0;

mask=zeros(21,21,2);
mask(:,:,1)=1;
mask(:,:,2)=0.2;

[simulation,selected_ti]=g2s('-a','as', ...
	'-ti',{ti1 ti2}, ...
	'-di',conditioning, ...
	'-sp',path, ...
	'-mi',mask, ...
	'-dt',[0], ...
	'-k',2, ...
	'-n',8, ...
	'-j',0.5);

figure(1); clf
subplot(1,5,1); imagesc(ti1); axis image off; title('TI 1');
subplot(1,5,2); imagesc(ti2); axis image off; title('TI 2');
subplot(1,5,3); imagesc(conditioning(:,:,1)); axis image off; title('Conditioning');
subplot(1,5,4); imagesc(simulation(:,:,1)); axis image off; title('AS result');
subplot(1,5,5); imagesc(selected_ti(:,:,1)); axis image off; title('Selected TI id');

Larger diagnostic example in Python

For a richer diagnostic script with a rectangular synthetic case, -mi edge cases, exported TI-id checks, and a live figure window, see:

• example/python/AnchorSamplingSyntheticExperiment.py

When to use AS instead of QS

Use AS when:

• your training data are available as a large stack of aligned realizations or scenarios
• the same coordinate across TIs has a stable physical meaning
• you want neighborhood matching to choose which TI to use at a location, rather than which location to sample from

Use QS when:

• stationarity in TI coordinates is acceptable
• the center value should be allowed to come from any TI location
• you need workflows such as TI-index maps through -ii or augmented-dimensional simulation

Notes

• AS keeps the same general workflow as QS, including categorical, continuous, multivariate, kernel-based, and path-image driven simulations.
• In masked runs, AS computes mismatch first, applies -mi admissibility masking before top-k ranking (dropping non-finite and non-positive entries), then uses -mi as relative weights inside the retained set.
• For -cnorm / -cn, if you provide more than one value, the number of values must match the number of continuous variables in -dt.