Columnar Programs

Rediscovers a causal program over a table of numeric columns. The search uses the vectorized primitive kit and the causal window primitives, and the fitness function masks the nan warmup that those windows produce.

See the Columnar Programs tutorial for a walkthrough of the pieces used here.

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization

COLUMNS = ["level", "flow", "noise"]


def make_columns(size=512):
    steps = numpy.linspace(0.0, 12.0, size)
    level = numpy.sin(steps) + 0.25 * steps
    flow = numpy.cos(steps * 0.7) * 2.0
    noise = numpy.abs(numpy.sin(steps * 3.1)) + 0.5
    return tuple(numpy.ascontiguousarray(c, dtype=numpy.float64) for c in (level, flow, noise))


def make_target(columns):
    # The program the search is expected to rediscover. It is causal,
    # so its first samples are nan and the fitness must ignore them.
    level, flow, _noise = columns
    return gp.rolling_mean(level, 4) - gp.delay(flow, 2)


def evaluate(individual, toolbox, columns, target):
    func = toolbox.compile(expr=individual)
    predicted = numpy.broadcast_to(numpy.asarray(func(*columns), dtype=numpy.float64), target.shape)
    valid = numpy.isfinite(predicted) & numpy.isfinite(target)
    if valid.sum() < target.size // 2:
        return (1.0e6,)  # too much warmup or too many holes to judge
    return (float(numpy.mean((predicted[valid] - target[valid]) ** 2)),)


def setup():
    columns = make_columns()
    target = make_target(columns)

    pset = gp.columnar_pset(COLUMNS, window=(2, 8))

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    gp.register_gp(
        toolbox,
        pset,
        individual=creator.Individual,
        min_depth=1,
        max_depth=3,
        height_limit=8,
    )
    toolbox.register("evaluate", evaluate, toolbox=toolbox, columns=columns, target=target)

    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", numpy.mean)
    stats.register("min", numpy.min)

    return toolbox, stats


def print_results(best_ind):
    print(f"\nBest program: {best_ind}")
    print(f"Mean squared error: {best_ind.fitness.values[0]:.6g}")


def main():
    toolbox, stats = setup()
    pop = toolbox.population(size=300)
    hof = tools.HallOfFame(1)
    args = {
        "toolbox": toolbox,
        "population": pop,
        "generations": 40,
        "cx_prob": 0.5,
        "mut_prob": 0.2,
        "hof": hof,
        "stats": stats,
        "verbose": True,  # prints stats
    }
    tools.ea_simple(**args)
    print_results(hof[0])


if __name__ == "__main__":
    main()

Batch evaluation

The same columnar contract, scored with evaluate_batch, interpret_tapes, and per-fold case_errors. Pair-window and time-series primitives are on the primitive set.

import numpy
from deap_er import Fitness, Toolbox, creator, gp, tools

tools.rng.seed(1234)  # disables randomization

COLUMNS = ["level", "flow"]
CASES = [(0, 256), (256, 512)]
N_ROWS = 512


def make_columns(size=N_ROWS):
    steps = numpy.linspace(0.0, 12.0, size)
    level = numpy.sin(steps) + 0.25 * steps
    flow = numpy.cos(steps * 0.7) * 2.0
    return tuple(numpy.ascontiguousarray(c, dtype=numpy.float64) for c in (level, flow))


def make_target(columns):
    level, flow = columns
    # Pair-window beta and a time-series rank — both causal.
    return gp.rolling_beta(level, flow, 4) + gp.ts_rank(level, 6)


def setup():
    columns = make_columns()
    target = make_target(columns)
    matrix = numpy.column_stack(columns)

    pset = gp.columnar_pset(COLUMNS, window=(2, 8), pair_windows=True, ts=True)

    creator.create_type("FitnessMin", Fitness, weights=(-1.0,))
    creator.create_type("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin)

    toolbox = Toolbox()
    gp.register_gp(
        toolbox,
        pset,
        individual=creator.Individual,
        min_depth=1,
        max_depth=3,
        height_limit=8,
        backend="opcode",
    )
    toolbox.register(
        "evaluate_batch",
        gp.evaluate_columnar,
        pset=pset,
        matrix=matrix,
        target=target,
        cases=CASES,
    )

    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", numpy.mean)
    stats.register("min", numpy.min)
    return toolbox, stats, pset, matrix, target


def print_results(best_ind, toolbox, matrix, target):
    func = toolbox.compile(expr=best_ind)
    predicted = func(*[matrix[:, i] for i in range(matrix.shape[1])])
    errors = tools.case_errors(predicted, target, CASES)
    print(f"\nBest program: {best_ind}")
    print(f"Case MSE: {tuple(round(e, 6) for e in errors)}")
    print(f"Mean case MSE: {float(numpy.mean(errors)):.6g}")


def main():
    toolbox, stats, _pset, matrix, target = setup()
    pop = toolbox.population(size=80)
    hof = tools.HallOfFame(1)
    tools.ea_simple(
        toolbox,
        pop,
        generations=15,
        cx_prob=0.5,
        mut_prob=0.2,
        hof=hof,
        stats=stats,
        verbose=True,
    )
    print_results(hof[0], toolbox, matrix, target)


if __name__ == "__main__":
    main()