User guide

Everything below is runnable as-is; outputs shown are from real runs.

Install

pip install beamfeat            # core: numpy + scikit-learn only
pip install "beamfeat[units]"   # + pint, for dimensional analysis

See Installation for running the test suite or reproducing the benchmarks.

Sixty seconds to a vetted equation

import numpy as np
from beamfeat import BeamFeatRegressor

rng = np.random.default_rng(0)
X = rng.uniform(1, 6, (400, 4))
y = X[:, 0] * X[:, 1] + rng.normal(0, 0.05, 400)

model = BeamFeatRegressor(max_depth=2, beam_width=25, random_state=0).fit(X, y)
print(model.equation())
print(model.fdr_controlled_)

y = 0.9974*(x0 * x1) + 0.0264
True

Two lines matter. equation() is the fitted model itself — evaluate it on raw feature values and you reproduce predict(). fdr_controlled_ states whether the features carry the false-discovery-rate guarantee; check it before treating them as statistically vetted.

Column names flow into formulas

Fit on a DataFrame and formulas use your names:

import pandas as pd

df = pd.DataFrame({"mass": rng.uniform(1, 5, 300),
                   "vol":  rng.uniform(1, 5, 300)})
target = df["mass"] / df["vol"] + rng.normal(0, 0.02, 300)
print(BeamFeatRegressor(random_state=0).fit(df, target).formulas())

['(mass / vol)']

Dimensional analysis

Give units as pint quantities or plain strings; dimensionally invalid expressions (metres plus kilograms) are rejected before any numerical work — so x0 + x1 (kg plus m) is never even scored:

y_phys = X[:, 0] * X[:, 1] / X[:, 2]          # kg·m/s
model = BeamFeatRegressor(
    units={"x0": "kg", "x1": "m", "x2": "s"}, random_state=0
).fit(X[:, :3], y_phys)
print(model.formulas())

['((x1 / x2) * x0)']

The recovered form, kg·m/s, is exactly the target's dimension — and the dimensionally invalid spellings never consumed a beam slot.

Audit what selection did

Every candidate's exact p-value and q-value is kept:

for row in model.selection_report_[:3]:
    print(row["formula"], round(row["q_value"], 4), row["kept"])

Features are kept only if they pass FDR screening at target_fdr (default 0.1, Benjamini–Yekutieli) on a held-out split, then survive a parsimony pass that keeps the compact predictive subset.

Honest failure, by default

On a target with no structure, nothing passes — and the model says so instead of returning junk:

noise = rng.standard_normal(400)
model = BeamFeatRegressor(random_state=3).fit(X, noise)
# NoDiscoveriesWarning: ... Returning no constructed features ...
print(model.equation())

y = 0.0353  (no feature passed selection)

Prefer the old behaviour? on_no_discoveries="fallback" keeps the unfiltered search output (flagged fdr_controlled_=False); "raise" raises.

Classification

from beamfeat import BeamFeatClassifier

labels = (X[:, 0] * X[:, 1] > X[:, 2] * X[:, 3]).astype(int)
clf = BeamFeatClassifier(max_depth=2, random_state=0).fit(X, labels)
print(clf.equation())          # log-odds form; boundary = zero level set

Real data with missing values or categoricals

beamfeat never imputes or encodes silently — compose explicitly:

from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer

pipe = make_pipeline(SimpleImputer(strategy="median"),
                     BeamFeatRegressor(random_state=0))

When to reach for something else

If a tree model beats beamfeat by a wide margin, your signal is likely piecewise, not algebraic — use the tree. If you need constants fitted inside expressions (exp(-3.2*x)), use a symbolic regressor such as PySR. The guarantees page states every boundary with the measurement behind it.