Value types

Complex numbers

To fit complex-valued expressions, pass feature and target arrays with a complex NumPy dtype.

The code constructs a complex cosine target and restricts the search to +, -, *, and cos:

python

import numpy as np

from pysr import PySRRegressor

X = np.random.randn(100, 1) + 1j * np.random.randn(100, 1)
y = (1 + 2j) * np.cos(X[:, 0] * (0.5 - 0.2j))

model = PySRRegressor(
    binary_operators=["+", "-", "*"], unary_operators=["cos"], niterations=100,
)

model.fit(X, y)

The fitted constants belong to the complex-valued expression domain. Use model.sympy() to inspect the selected equation in SymPy form:

python

model.sympy()

Pass an equation-row index to predict; -1 selects the final row of equations_:

python

model.predict(X, -1)

Julia packages and types

PySR delegates its symbolic-regression search to the pure Julia package SymbolicRegression.jl. A custom operator can call functionality from another Julia package through that backend.

Load the package that supplies the prime lookup before defining the operator. The target in this example is

y = p_{3 x + 1} - 5,

where $p_{i}$ denotes the $i$ th prime and $x$ is the input feature. The lookup comes from Primes.jl.

First, import the handle to PySR's Julia runtime:

python

from pysr import jl

jl is the juliacall runtime object, and jl.seval evaluates Julia source from Python. Use it to add Primes.jl to the active PySR environment:

python

from pysr import jl

jl.seval("""
import Pkg
Pkg.add("Primes")
""")

Adding a package and importing its namespace are separate steps. Import Primes before defining p:

python

from pysr import jl

jl.seval("import Primes")

With the package available, define the unary operator used by the search:

python

from pysr import jl

jl.seval("""
function p(i::T) where T
    if (0.5 < i < 1000)
        return T(Primes.prime(round(Int, i)))
    else
        return T(NaN)
    end
end
""")

The operator rounds inputs within 0.5 < i < 1000 to prime indices and preserves the input type T. Outside that domain it returns T(NaN) to reject invalid evaluations.

Call the Julia function through jl.p to build the prime values used by the synthetic data:

python

from pysr import jl

primes = {i: jl.p(i*1.0) for i in range(1, 999)}

Use the lookup table to generate noisy samples of the target:

python

import numpy as np

X = np.random.randint(0, 100, 100)[:, None]
y = [primes[3*X[i, 0] + 1] - 5 + np.random.randn()*0.001 for i in range(100)]

Register p with the regular PySR operator interface. The sympy_p subclass and extra_sympy_mappings give SymPy a symbolic form for the custom function:

python

from pysr import PySRRegressor
import sympy

class sympy_p(sympy.Function):
    pass

model = PySRRegressor(
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["p"],
    niterations=100,
    extra_sympy_mappings={"p": sympy_p}
)

Fit the model on the generated pairs:

python

model.fit(X, y)

Since p rounds its argument, different real offsets can produce the same integer prime index, so the offset in a fitted expression need not equal 1 exactly. Inspect the selected equation in SymPy with:

python

model.sympy()

Custom value types

TypeSpec declares a custom Julia value that can flow through an expression tree. Vectors, tensors, strings, structs, and other Julia value types can serve as payloads when their behavior is described by the type's hooks, operators, and loss.

A model configured with type_spec=TypeSpec(...) needs three pieces:

Define hooks that sample and change values. For continuously optimized constants, scalar_constants exposes scalar coordinates and with_scalar_constants rebuilds the value. A custom mutate hook handles changes that need a discrete or structural operation.
Choose exactly one of elementwise_loss, loss_function, or loss_function_expression. The selected loss must accept the custom values and return a real-valued score. Elementwise-loss return types are inferred; a full objective requires an explicit loss_type.
Supply operators through operators={...}, grouped by arity. Each operator must accept and return the declared value type so the type-stability check succeeds. Add an explicit return annotation such as ::Vec2 when Julia cannot infer it.

The examples below combine these pieces for vector values, strings as discrete constants, and payloads with several shapes.

Vector-valued expression trees

This example searches for an expression whose values are two-dimensional vectors:

y = rotate 90 (x_{1}) + 2 x_{2} + [\begin{matrix} 0.5 \\ - 1.0 \end{matrix}] .

The 128-by-2 DataFrame has two vector-valued features per row, and each entry of y is one target vector. The table dimensions describe samples and features, while each cell of the object array holds a single logical vector value.

python

import numpy as np
import pandas as pd

from pysr import PySRRegressor, TypeSpec

rng = np.random.default_rng(0)
x1 = [rng.normal(size=2) for _ in range(128)]
x2 = [rng.normal(size=2) for _ in range(128)]
X = pd.DataFrame({"x1": x1, "x2": x2})
y = np.empty(128, dtype=object)
offset = np.array([0.5, -1.0])
y[:] = [np.array([-a[1], a[0]]) + 2 * b + offset for a, b in zip(x1, x2)]

Declare Vec2 with one Vector{Float64} field. sample creates a random payload of length two, while scalar_constants and with_scalar_constants expose and rebuild the entries used by continuous constant optimization:

python

from pysr import TypeSpec

type_spec = TypeSpec(
    "Vec2",
    fields={"data": "Vector{Float64}"},
    sample="rng -> Vec2(randn(rng, 2))",
    scalar_constants="value -> value.data",
    with_scalar_constants=(
        "(value, scalar_constants) -> Vec2(scalar_constants)"
    ),
)

PySR derives how values are initialized, mutated, validated, counted, packed, and unpacked from these hooks. Add init or mutate for specialized search behavior, or string for custom display, when the derived defaults are insufficient.

TypeSpec generates the Julia value type from that declaration. Its equivalent Julia struct has one vector field:

julia

struct Vec2
    data::Vector{Float64}
end

Define operators that consume and return Vec2 values:

python

operators = {
    1: [
        "rotate90(a::Vec2) = Vec2([-a.data[2], a.data[1]])",
        "double(a::Vec2) = Vec2(2a.data)",
    ],
    2: ["add_vectors(a::Vec2, b::Vec2) = Vec2(a.data + b.data)"],
}

The elementwise loss compares vector payloads by summing squared differences, which gives the real-valued score the search requires:

python

from pysr import PySRRegressor

model = PySRRegressor(
    type_spec=type_spec,
    operators=operators,
    elementwise_loss="vector_loss(a::Vec2, b::Vec2) = sum(abs2, a.data - b.data)",
    niterations=40,
    populations=4,
    maxsize=10,
)

model.fit(X, y)
print(model.equations_)

The search space can represent the target as add_vectors(add_vectors(rotate90(x1), double(x2)), Vec2([0.5, -1.0])). The vector constant carries both offset components, which the scalar hook pair exposes for joint optimization.

String-valued expressions and discrete constants

This example treats each string as a value and chooses a separator from a finite set during mutation. It lowercases the first input, uppercases the second, and joins them with the selected separator:

python

import numpy as np
import pandas as pd

from pysr import PySRRegressor, TypeSpec

X = pd.DataFrame(
    {
        "first": ["Py", "symbolic", "hello", "left"],
        "second": ["SR", "regression", "world", "right"],
    }
)
y = np.array(
    [f"{a.lower()}-{b.upper()}" for a, b in X.itertuples(index=False)],
    dtype=object,
)

type_spec = TypeSpec(
    "StringValue",
    fields={"data": "String"},
    sample='rng -> StringValue(rand(rng, ("", "-", "_")))',
    mutate="""
    mutate_string(rng, value, temperature) = StringValue(rand(rng, ("", "-", "_")))
    """,
)

model = PySRRegressor(
    type_spec=type_spec,
    operators={
        1: [
            "string_lowercase(x::StringValue) = StringValue(lowercase(x.data))",
            "string_uppercase(x::StringValue) = StringValue(uppercase(x.data))",
        ],
        2: [
            "string_concat(a::StringValue, b::StringValue) = StringValue(a.data * b.data)"
        ],
    },
    elementwise_loss="""
    string_loss(a::StringValue, b::StringValue) = Float64(Base.editdistance(a.data, b.data))
    """,
    niterations=40,
)

model.fit(X, y)
print(model.equations_)

The sample and mutate hooks choose discrete separator constants. Without a scalar_constants and with_scalar_constants pair, BFGS leaves them unchanged. Edit distance provides the real-valued loss for string outputs.

Advanced: recovering a neural network with tensor constants

A TypeSpec can place scalar, vector, and matrix payloads in the same Julia value type. Its scalar-constant hooks flatten each payload for BFGS and rebuild its original shape.

The advanced example searches for this two-layer neural network:

y = W_{2} relu (W_{1} x + b_{1}) + b_{2}

The targets are vectors. safe_matmul and safe_add return NaN for incompatible shapes, so expressions using those helpers receive an invalid value instead of a dimension error:

python

import numpy as np
import pandas as pd

from pysr import PySRRegressor, TypeSpec

preamble = """
const NNPayload = Union{Float64, Vector{Float64}, Matrix{Float64}}

safe_matmul(a::Matrix{Float64}, b::Vector{Float64}) =
    size(a, 2) == length(b) ? a * b : NaN
safe_matmul(::NNPayload, ::NNPayload) = NaN

safe_add(a::Float64, b::Float64) = a + b
safe_add(a::T, b::T) where {T<:Union{Vector{Float64}, Matrix{Float64}}} =
    size(a) == size(b) ? a + b : NaN
safe_add(::NNPayload, ::NNPayload) = NaN

# Constants sample a random rank, generating only the payload that was chosen:
function random_nn_payload(rng)
    rank = rand(rng, 0:2)
    rank == 0 ? randn(rng) : rank == 1 ? randn(rng, 2) : randn(rng, 2, 2)
end
"""

type_spec = TypeSpec(
    "NNValue",
    fields={"data": "NNPayload"},
    sample="rng -> NNValue(random_nn_payload(rng))",
    # Mutations usually perturb every scalar in the payload, but occasionally
    # resample a fresh rank:
    mutate="""
    (rng, value, temperature) -> if rand(rng) < 0.1
        NNValue(random_nn_payload(rng))
    else
        NNValue(value.data .+ temperature .* randn(rng, size(value.data)...))
    end
    """,
    scalar_constants="""
    function scalar_constants(value)
        return value.data isa Float64 ? [value.data] : vec(value.data)
    end
    """,
    with_scalar_constants="""
    function with_scalar_constants(value, scalar_constants)
        data = value.data isa Float64 ? scalar_constants[1] :
            reshape(collect(scalar_constants), size(value.data))
        return NNValue(data)
    end
    """,
    preamble=preamble,
)

preamble makes the payload union and helpers available before the generated type is defined. Mutation can change a constant's rank; continuous optimization preserves its shape through the scalar hook pair.

Generate vector-valued targets from fixed weights and biases:

python

import numpy as np
import pandas as pd

rng = np.random.default_rng(0)
x_values = rng.normal(size=(64, 2))
W1 = np.array([[1.2, -0.7], [0.5, 1.1]])
b1 = np.array([0.3, -0.2])
W2 = np.array([[0.8, -1.0], [1.3, 0.4]])
b2 = np.array([-0.4, 0.2])
y_values = (W2 @ np.maximum(x_values @ W1.T + b1, 0).T).T + b2

X = pd.DataFrame({"x": list(x_values)})
y = pd.Series(list(y_values), dtype=object)

Register the neural-network operations and score matching vector outputs with mean squared error. The loss penalizes other output shapes:

python

from pysr import PySRRegressor

model = PySRRegressor(
    type_spec=type_spec,
    operators={
        1: ["nn_relu(a) = NNValue(max.(a.data, 0.0))"],
        2: [
            "nn_matmul(a, b) = NNValue(safe_matmul(a.data, b.data))",
            "nn_add(a, b) = NNValue(safe_add(a.data, b.data))",
        ],
    },
    elementwise_loss="""
    function nn_mse(a, b)::Float64
        valid = a.data isa Vector && b.data isa Vector && size(a.data) == size(b.data)
        return valid ? sum(abs2, a.data .- b.data) / length(a.data) : 1.0e6
    end
    """,
    niterations=100,
    populations=4,
    maxsize=11,
)

model.fit(X, y)
print(model.equations_)

A recovered expression can have the form nn_matmul(W2, nn_add(b, nn_relu(nn_matmul(W1, nn_add(x, c))))). The inner constant c represents the first bias through $b_{1} = W_{1} c$ , and the outer constant b represents the second through $b_{2} = W_{2} b$ . Each displayed constant is an NNValue containing its fitted scalar, vector, or matrix payload.

The same custom-value interface can seed a search with guesses. In the vector example above, Vec2 is the generated constructor from its TypeSpec, add_vectors is the declared binary operator, and x1 is the first named feature column. A matching guess is guesses=["add_vectors(x1, Vec2([1.0, 2.0]))"].

Value types ​

Complex numbers ​

Julia packages and types ​

Custom value types ​

Vector-valued expression trees ​

Value types

Complex numbers

Julia packages and types

Custom value types

Vector-valued expression trees