Documentation

PYCSE

Module containing useful scientific and engineering functions.

• Linear regression

• Nonlinear regression.

• Differential equation solvers.

See http://kitchingroup.cheme.cmu.edu/pycse

Copyright 2025, John Kitchin (see accompanying license files for details).

pycse.PYCSE.Rsquared(y, Y)

Return R^2, or coefficient of determination.

y is a 1d array of observations. Y is a 1d array of predictions from a model.

Returns:

The R^2 value for the fit.

pycse.PYCSE.bic(x, y, model, popt)

Compute the Bayesian information criterion (BIC).

Parameters:

modelfunction(x, …) returns prediction for y

poptoptimal parameters

yarray, known y-values

Returns:

BICfloat

https://en.wikipedia.org/wiki/Bayesian_information_criterion#Gaussian_special_case

pycse.PYCSE.ivp(f, tspan, y0, *args, **kwargs)

Solve an ODE initial value problem.

This provides some convenience defaults that I think are better than solve_ivp.

Parameters:

ffunction

callable y’(x, y) = f(x, y)

tspanarray

The x points you want the solution at. The first and last points are used in

tspan in solve_ivp.

y0array

Initial conditions

*argstype

extra arguments passed to f, i.e. f(t, y, *args). They are passed to

solve_ivp as ``args=args``. Do not also pass ``args=`` as a keyword.

**kwargstype arbitrary kwargs to pass to solve_ivp.

max_step is set to be the min diff of tspan. dense_output is set to True.

t_eval is set to the array specified in tspan.

Returns:

solution from solve_ivp

pycse.PYCSE.lbic(X, y, popt)

Compute the Bayesian information criterion for a linear model.

Returns:

BICfloat

pycse.PYCSE.nlinfit(model, x, y, p0, alpha=0.05, **kwargs)

Nonlinear regression with confidence intervals.

Parameters:

modelfunction f(x, p0, p1, …) = y

xarray of the independent data

yarray of the dependent data

p0array of the initial guess of the parameters

alpha100*(1 - alpha) is the confidence interval

i.e. alpha = 0.05 is 95% confidence

kwargs are passed to curve_fit.

Returns:

[p, pint, SE]

p is an array of the fitted parameters pint is an array of confidence intervals SE is an array of standard errors for the parameters.

pycse.PYCSE.nlpredict(X, y, model, popt, xnew, loss=None, alpha=0.05, ub=1e-05, ef=1.05)

Prediction error for a nonlinear fit.

Parameters:

Xarray-like

Independent variable data used for fitting

yarray-like

Dependent variable data used for fitting

modelcallable

Model function with signature model(x, …)

poptarray-like

Optimized parameters from fitting (e.g., from curve_fit or nlinfit)

xnewarray-like

x-values to predict at

losscallable, optional

Loss function that was minimized during fitting, with signature loss(*params). If None (default), assumes scipy.optimize.curve_fit was used and automatically constructs the correct loss function: loss = 0.5 * sum((y - model(X, *p))**2). This is the ½SSE convention used by scipy’s least_squares optimizer. If you used a different optimizer or loss function, provide it explicitly. The loss must be proportional to the sum of squared residuals (e.g. ½SSE, SSE or MSE); the parameter covariance 2 * loss / (n - p) * H^-1 (H is the Hessian of the loss) is then independent of the proportionality constant.

alphafloat, optional

100*(1 - alpha) is the prediction interval level (default: 0.05 for 95%)

ubfloat, optional

Smallest allowed eigenvalue of the Hessian of the loss (absolute, default 1e-5). If the smallest eigenvalue is at least ub, no regularization is applied. Otherwise H + eps * I is inverted, with eps = max(ub, ub - ef * lambda_min).

effloat, optional

Eigenvalue factor used when regularizing (default: 1.05, see ub)

This function uses numdifftools for the Hessian and Jacobian.

Returns:

yarray

Predicted values at xnew

yintarray

Prediction intervals (for a new observation) at alpha confidence level, shape (n, 2); yint[:, 0] is the lower bound and yint[:, 1] the upper bound. Same layout as predict.

searray

Standard error of the fitted mean at xnew, sqrt(diag(J cov J^T)). This is NOT the prediction standard error used for yint, which is sqrt(se**2 + sigma**2). Note predict returns the prediction standard error.

Notes

The default loss function (½SSE) matches the convention used by scipy.optimize.curve_fit, which internally minimizes 0.5 * sum(residuals**2). If you provide a custom loss function, ensure it uses the same convention as your fitting procedure. The noise variance sigma**2 is estimated from the residuals as SSE / (n - p).

pycse.PYCSE.polyfit(x, y, deg, alpha=0.05, *args, **kwargs)

Least squares polynomial fit with parameter confidence intervals.

Parameters:

xarray_like, shape (M,)

x-coordinates of the M sample points (x[i], y[i]).

yarray_like, shape (M,) or (M, K)

y-coordinates of the sample points. Several data sets of sample points sharing the same x-coordinates can be fitted at once by passing in a 2D-array that contains one dataset per column.

degint

Degree of the fitting polynomial

*args and **kwargs are passed to regress.

Returns:

[b, bint, se]

b is a vector of the fitted parameters

bint is a 2D array of confidence intervals

se is an array of standard error for each parameter.

pycse.PYCSE.polyval(p, x, X, y, alpha=0.05, ub=1e-05, ef=1.05)

Evaluate the polynomial p at x with prediction intervals.

Parameters:

p: parameters from pycse.polyfit

x: array_like, shape (M,)

x-coordinates to evaluate the polynomial at.

X: array_like, shape (N,)

the original x-data that p was fitted from.

y: array-like, shape (N,)

the original y-data that p was fitted from.

alpha100*(1 - alpha) is the prediction interval level, 95% = 0.05

ubsmallest allowed eigenvalue of X^T X. The matrix is regularized only

if its smallest eigenvalue is below ub (see predict).

efeigenvalue factor used when regularizing (see predict).

Returns:

y, yint, pred_se

ythe predicted values

yint: prediction interval (for a new observation), shape (M, 2);

yint[:, 0] is the lower bound and yint[:, 1] the upper bound.

pred_se: standard error of prediction for a new observation,

sqrt(sigma^2 + se_mean^2), shape (M,).

pycse.PYCSE.predict(X, y, pars, XX, alpha=0.05, ub=1e-05, ef=1.05)

Prediction interval for linear regression.

Based on the delta method.

Parameters:

Xknown x-value array, one row for each y-point

yknown y-value array

parsfitted parameters

XXx-value array to make predictions for

alpha100*(1 - alpha) is the prediction interval level, 95% = 0.05

ubsmallest allowed eigenvalue of X^T X (absolute, default 1e-5). If the

smallest eigenvalue of X^T X is at least ub, no regularization is applied and the exact textbook prediction interval is returned. Otherwise (near-singular X^T X), X^T X + eps * I is inverted instead, with eps = max(ub, ub - ef * lambda_min).

efeigenvalue factor used when regularizing (see ub).

See https://en.wikipedia.org/wiki/Prediction_interval#Unknown_mean,_unknown_variance

Returns:

y, yint, pred_se

ythe predicted values

yint: prediction interval (for a new observation), shape (n, 2) for a

single output (yint[:, 0] is the lower bound, yint[:, 1] the upper bound), or (n, k, 2) for k outputs. Same layout as nlpredict.

pred_se: standard error of prediction for a new observation,

sqrt(sigma^2 + se_mean^2), where se_mean is the standard error of the fitted mean. Note nlpredict returns se_mean instead.

pycse.PYCSE.regress(A, y, alpha=0.05, *args, **kwargs)

Linear least squares regression with confidence intervals.

Solve the matrix equation (A p = y) for p.

The confidence intervals account for sample size using a student T multiplier.

This code is derived from the descriptions at http://www.weibull.com/DOEWeb/confidence_intervals_in_multiple_linear_regression.htm and http://www.weibull.com/DOEWeb/estimating_regression_models_using_least_squares.htm

Parameters:

Aa matrix of function values in columns, e.g.

A = np.column_stack([T**0, T**1, T**2, T**3, T**4])

ya vector of values you want to fit

alpha100*(1 - alpha) confidence level

args and kwargs are passed to np.linalg.lstsq

Returns:

[b, bint, se]

b is a vector of the fitted parameters

bint is an array of confidence intervals. The ith row is for the ith parameter.

se is an array of standard error for each parameter.

If alpha is None, bint and se are None.

pycse.utils

Provides utility functions in pycse.

1. Fuzzy comparisons for float numbers.

2. An ignore exception context manager

3. A handy function to read a google sheet.

pycse.utils.feq(x, y, epsilon=np.float64(2.220446049250313e-16), rtol=0)

Fuzzy equals.

x == y with tolerance epsilon + rtol * max(|x|, |y|).

epsilon is an absolute tolerance, rtol is an optional relative tolerance. Works elementwise on arrays; returns a bool for scalar inputs.

pycse.utils.fge(x, y, epsilon=np.float64(2.220446049250313e-16), rtol=0)

Fuzzy greater than or equal to.

x >= y with tolerance epsilon + rtol * max(|x|, |y|). Works elementwise on arrays; returns a bool for scalar inputs.

pycse.utils.fgt(x, y, epsilon=np.float64(2.220446049250313e-16), rtol=0)

Fuzzy greater than.

x > y with tolerance epsilon + rtol * max(|x|, |y|). Works elementwise on arrays; returns a bool for scalar inputs.

pycse.utils.fle(x, y, epsilon=np.float64(2.220446049250313e-16), rtol=0)

Fuzzy less than or equal to.

x <= y with tolerance epsilon + rtol * max(|x|, |y|). Works elementwise on arrays; returns a bool for scalar inputs.

pycse.utils.flt(x, y, epsilon=np.float64(2.220446049250313e-16), rtol=0)

Fuzzy less than.

x < y with tolerance epsilon + rtol * max(|x|, |y|). Works elementwise on arrays; returns a bool for scalar inputs.

pycse.utils.ignore_exception(*exceptions)

Context manager to ignore EXCEPTIONS raised in its body.

A message is printed when an exception is caught.

>>> with ignore_exception(ZeroDivisionError):
...     print(1/0)
caught division by zero

pycse.utils.read_gsheet(url, *args, **kwargs)

Return a dataframe for the Google Sheet at url.

args and kwargs are passed to pd.read_csv The url should be viewable by anyone with the link.

pycse.plotly

Module for using plotly with orgmode.

This monkey-patches go.Figure.show to provide a png image for org-mode, and an html file that is saved that you can click on in org-mode to see the interactive version.

Example:

import pycse.plotly import plotly.graph_objects as go

# Sample data x = [0, 1, 2, 3, 4] y = [0, 1, 4, 9, 16]

# Create a figure fig = go.Figure(

data=go.Scatter(x=x, y=y, mode=’lines+markers’, name=’y = x^2’), layout=go.Layout(

title=’Simple Line Plot’, xaxis_title=’x’, yaxis_title=’y’,

)

)

# Show the plot fig.show()

pycse.plotly.myshow(self, *args, **kwargs)

Make a PNG image to display for plotly.

pycse.hashcache

hashcache - a class decorator for persistent, file/hash-based cache.

I found some features of joblib were unsuitable for how I want to use a cache.

1. The “file” Python thinks the function is in is used to save the results in joblib, which leads to repeated runs if you run the same code in Python, notebook or stdin, and means the cache is not portable to other machines, and maybe not even in time since temp directories and kernel parameters are involved. I could not figure out how to change those in joblib.

2. joblib uses the function source code in the hash, so inconsequential changes like whitespace, docstrings and comments change the hash.

This library aims to provide a simpler version of what I wish joblib did for me.

Results are cached based on a hash of the function name, bytecode, constants (e.g. numeric literals, recursively including nested functions), referenced names, arg values and kwarg values. I use joblib.hash for this. The values of global variables used by the function are NOT part of the hash, so changing a global will not invalidate the cache.

The cache location is set as a class attribute:

HashCache.cache = ‘./cache’

HashCache - stores joblib.dump pickle strings in files named by hash

SqlCache - stores orjson serialized data in a sqlite3 database by hash key

JsonCache - stores orjson serialized data in json files, compatible with maggma

SqlCache and JsonCache tag numpy arrays, numpy scalars and tuples when they are serialized, so a cache hit returns the same types as the miss that stored it. Other objects must be JSON serializable by orjson (or by the default function you set on the class).

This is still alpha, proof of concept code. Test it a lot for your use case. The API is not stable, and subject to change.

Pros:

1. File-based cache which means many functions can run in parallel reading and writing, and you are limited only by file io speeds, and disk space.

2. semi-portability. The cache could be synced across machines, and caches can be merged with little risk of conflict.

3. No server is required. Everything is done at the OS level.

4. Extendability. You can define your own functions for loading and dumping data.

Cons:

1. hashes are fragile and not robust. They are fragile with respect to any changes in how byte-code is made, or via mutable arguments, etc. The hashes are not robust to system level changes like library versions, or global variables. The only advantage of hashes is you can compute them.

2. File-based cache which means if you generate thousands of files, it can be slow to delete them. Although it should be fast to access the results since you access them directly by path, it will not be fast to iterate over all the results, e.g. if you want to implement some kind of search or reporting.

3. No server. You have to roll your own update strategy if you run things on multiple machines that should all cache to a common location.

Changelog

[2023-09-23 Sat] Changed hash signature (breaking change). It is too difficult to figure out how to capture global state, and the use of internal variable names is not consistent with using the bytecode to be insensitive to unimportant variable name changes.

Pulled out some functions for loading and dumping data. This is a precursor to enabling other backends like lmdb or sqlite instead of files. You can then simply provide new functions for this.

[2024-06-18 Tue] Changed from function to class decorator (breaking change).

[2026-09-26 Sat] The hash now includes the code constants and names (breaking change, old caches will not be found). Previously a change like x * 2 -> x * 3 did not change the hash, and stale results were returned.

class pycse.hashcache.HashCache(function)

Class decorator to cache using hashes and pickle (via joblib).

Data is stored in directories named by the hash.

Methods

__call__(*args, **kwargs)	Code to run around self.function.
dump(**kwargs)	Dump KWARGS to the cache.
dump_data(hsh, data)	Dump DATA into HSH.
get_hash(args, kwargs)	Get a hash for running FUNC(ARGS, KWARGS).
get_hashpath(hsh)	Return path to file for HSH.
get_standardized_args(args, kwargs)	Return a standardized dictionary of kwargs for func(args, kwargs).
load(hsh[, cache])	Load saved variables from HSH.
load_data(hsh)	Load data for HSH.

static dump(**kwargs)

Dump KWARGS to the cache.

Returns a hash string for future lookup.

cache is a special kwarg that is not saved. It defaults to HashCache.cache, the same directory the decorator uses.

dump_data(hsh, data)

Dump DATA into HSH.

get_hash(args, kwargs)

Get a hash for running FUNC(ARGS, KWARGS).

This is the most critical feature of hashcache as it provides a key to store and look up results later. You should think carefully before changing this function, it breaks past caches.

FUNC should be as pure as reasonable. This hash is insensitive to the values of global variables the function uses; if you change a global, the cached result will still be returned.

The hash is on the function name, bytecode, constants (recursively including nested code objects), referenced names, and a standardized kwargs including defaults. The docstring is not included. We use bytecode because it is insensitive to things like whitespace, comments, docstrings, and variable name changes that don’t affect results. It is assumed that two functions with the same name and bytecode will evaluate to the same result. However, this makes the hash fragile to changes in Python version that affect bytecode.

get_hashpath(hsh)

Return path to file for HSH.

get_standardized_args(args, kwargs)

Return a standardized dictionary of kwargs for func(args, kwargs).

This dictionary includes default values, even if they were not called.

static load(hsh, cache=None)

Load saved variables from HSH.

CACHE defaults to HashCache.cache, the same directory the decorator uses.

load_data(hsh)

Load data for HSH.

HSH is a string for the hash associated with the data you want.

Returns success, data. If it succeeds, success with be True. If the data does not exist yet, sucess will be False, and data will be None.

class pycse.hashcache.JsonCache(function)

Json-based cache.

This is compatible with maggma.

Attributes:

default

Methods

__call__(*args, **kwargs)	Code to run around self.function.
dump(**kwargs)	Dump KWARGS to the cache.
dump_data(hsh, data)	Dump DATA into HSH.
get_hash(args, kwargs)	Get a hash for running FUNC(ARGS, KWARGS).
get_hashpath(hsh)	Return path to file for HSH.
get_standardized_args(args, kwargs)	Return a standardized dictionary of kwargs for func(args, kwargs).
load(hsh)	Load data from HSH.
load_data(hsh)	Load data from hsh.

static dump(**kwargs)

Dump KWARGS to the cache.

Returns a hash string for future lookup.

dump_data(hsh, data)

Dump DATA into HSH.

static load(hsh)

Load data from HSH.

load_data(hsh)

Load data from hsh.

class pycse.hashcache.SqlCache(function)

Class decorator to cache using orjson and sqlite.

Data is stored in a sqlite database as json.

Attributes:

default

Methods

__call__(*args, **kwargs)	Code to run around self.function.
dump(**kwargs)	Dump KWARGS to the cache.
dump_data(hsh, data)	Dump DATA into HSH.
get_hash(args, kwargs)	Get a hash for running FUNC(ARGS, KWARGS).
get_hashpath(hsh)	Return path to file for HSH.
get_standardized_args(args, kwargs)	Return a standardized dictionary of kwargs for func(args, kwargs).
load(hsh)	Load data from HSH.
load_data(hsh)	Load data for HSH.
search(query, *args)	Run a sql QUERY with args.

static dump(**kwargs)

Dump KWARGS to the cache.

Returns a hash string for future lookup.

dump_data(hsh, data)

Dump DATA into HSH.

DATA must be serializable to json.

static load(hsh)

Load data from HSH.

load_data(hsh)

Load data for HSH.

HSH is a string for the hash associated with the data you want.

Returns success, data. If it succeeds, success with be True. If the data does not exist yet, sucess will be False, and data will be None.

static search(query, *args)

Run a sql QUERY with args.

args are substituted in ? placeholders in the query.

This is just a light wrapper on con.execute.

pycse.hashcache.hashcache(*args, **kwargs)

Raise an exception if the old hashcache decorator is used.

pycse.sklearn.mip

Mixed-integer programming (MIP) export of trained piecewise-linear models.

A trained model built from affine maps, ReLUs, piecewise-linear splines and tree splits is an exactly piecewise-linear function, so it can be embedded in an optimization problem as mixed-integer linear constraints. This module holds the solver-side machinery once, and each estimator only describes its trained weights:

• Affine, ReLU and PiecewiseLinear are layer specs holding numpy arrays.

• MIPNetwork is a list of layers plus optional uncertainty (StdSpec).

• TreeModel is a list of Leaf boxes, each with an affine prediction.

• propagate_bounds computes interval bounds for every layer, which give the big-M constants.

• add_to_pyomo writes a MIPNetwork or TreeModel into a Pyomo block.

• MIPExportMixin gives an estimator to_pyomo() from its mip_network().

The ReLU y = max(0, a) with bounds L <= a <= U (L < 0 < U) is encoded with one binary d:

y >= 0,  y >= a,  y <= a - L (1 - d),  y <= U d

Units that are always active (L >= 0) or always inactive (U <= 0) need no binary. Tight input bounds give tight big-M constants and faster solves.

formulation="convex" drops the two upper inequalities, so the model is a pure LP (no binaries). It is exact for input-convex networks (every weight after the first ReLU nonnegative) when the output is minimized or bounded above: then y >= f(x) holds everywhere and y = f(x) at the optimum.

Uncertainty, when a model has it, has the form std = sqrt(||A z + c||^2 + d) where z is the output of some layer (the last-layer features or the ensemble members). This is a second-order cone, so optimizing e.g. mean + k * std is a mixed-integer SOCP that solvers such as Gurobi handle directly.

Example:

import pyomo.environ as pyo

m = pyo.ConcreteModel()
m.T = pyo.Var(bounds=(300, 500))
m.P = pyo.Var(bounds=(1, 10))
m.surrogate = pyo.Block()
model.to_pyomo(block=m.surrogate, inputs=[m.T, m.P])
m.obj = pyo.Objective(expr=m.surrogate.y[0], sense=pyo.maximize)
pyo.SolverFactory("highs").solve(m)

Pyomo is an optional dependency (pip install pycse[mip]); this module only imports it when a model is built.

class pycse.sklearn.mip.Affine(W: ndarray, b: ndarray, Wx: ndarray = None)

Affine layer z_out = z_in @ W + x @ Wx + b.

Args:

W: Weight matrix, shape (n_in, n_out). b: Bias vector, shape (n_out,). Wx: Optional skip connection from the network inputs x, shape

(n_inputs, n_out), as in input-convex and monotonic networks.

Attributes:

Wx

n_in

n_out

class pycse.sklearn.mip.Elementwise1D(n_in: int, index: ndarray, funcs: list, slope: ndarray = None, intercept: ndarray = None, spacing: ndarray = None, curvature: ndarray = None)

Features z_out[j] = f_j(z_in[index[j]]) of one input each.

Linear features (funcs[j] is None) are slope[j] * z + intercept[j] and are exported exactly. Nonlinear features (e.g. the sin/cos of a periodic model) are interpolated in the MIP on a grid with spacing at most spacing[j] over the bounds of their input, which add_to_pyomo knows only when it is called; see MIPNetwork.encoded. With |f_j''| <= curvature[j] the interpolation error of feature j is at most curvature[j] * h**2 / 8 for grid spacing h.

Args:

n_in: Number of inputs. index: Input of each feature, shape (n_out,). funcs: Per feature, a vectorized numpy function, or None for a linear feature. slope, intercept: Linear features’ coefficients (ignored for nonlinear ones). spacing: Maximum grid spacing per feature (nonlinear features only). curvature: Bound on |f''| per feature, for max_error.

Attributes:

curvature

intercept

n_out

slope

spacing

Methods

__call__(z)	Call self as a function.
grid(i, lo, hi)	The interpolation grid of input i on [lo, hi].
max_error(lo, hi)	Bound on each feature's interpolation error for inputs in [lo, hi].
to_piecewise(lo, hi)	The PiecewiseLinear layer the MIP uses for inputs in [lo, hi].

grid(i, lo, hi)

The interpolation grid of input i on [lo, hi].

max_error(lo, hi)

Bound on each feature’s interpolation error for inputs in [lo, hi].

to_piecewise(lo, hi)

The PiecewiseLinear layer the MIP uses for inputs in [lo, hi].

class pycse.sklearn.mip.Leaf(lower: ndarray, upper: ndarray, w: ndarray, c: float)

A box lower < x <= upper on which the prediction is x @ w + c.

Bounds may be infinite. This matches sklearn trees, which send x <= threshold to the left child.

class pycse.sklearn.mip.MIPExportMixin

Adds to_pyomo() to an estimator that implements mip_network().

mip_network() must return the trained model as a MIPNetwork (or TreeModel) in the original (unscaled) input and output units, and raise ValueError when the configuration is not piecewise-linear (e.g. a smooth activation). An estimator can also define _mip_default_input_bounds() to supply bounds when the caller gives none.

Methods

to_pyomo([input_bounds, block, inputs, ...])

Export the trained model as Pyomo MILP constraints.

mip_network

to_pyomo(input_bounds=None, block=None, inputs=None, include_std=False, formulation='bigm', **options)

Export the trained model as Pyomo MILP constraints.

Args:

input_bounds: (lower, upper) per input feature. Required unless inputs

are bounded Pyomo vars. Tighter bounds give faster solves.

block: Pyomo block to populate (e.g. m.surrogate = pyo.Block()).

A new ConcreteModel is created and returned if None.

inputs: Existing Pyomo vars to use as the model inputs. include_std: Also add block.std (models with uncertainty only):

True/"cone" for a second-order cone, exact when the objective penalizes the std; "exact" for a nonconvex equality, needed when it rewards the std. See pycse.sklearn.mip.add_to_pyomo.

formulation: "bigm" (exact MILP) or "convex" (an LP for

input-convex networks whose output is minimized or bounded above).

**options: Model-specific options passed to mip_network(), e.g.

split_gap for LeafModelRegressor or points_per_period for JAXPeriodicRegressor.

Returns:

The populated block. block.x[i] are the inputs, block.y[j] the predictions, and block.variance[j] the predictive variance when the model has uncertainty.

class pycse.sklearn.mip.MIPNetwork(layers: list, std: list = None)

A trained model as a sequence of layers, with optional uncertainty.

Args:

layers: List of Affine, ReLU and PiecewiseLinear layers; the first

must be Affine or PiecewiseLinear.

std: Optional list of StdSpec, one per output.

Attributes:

n_inputs

n_outputs

std

Methods

encoded(input_bounds)	The network exactly as the MIP represents it for these input bounds.
forward(X[, return_all])	Evaluate the network with numpy.
max_feature_error(input_bounds)	Bounds on the interpolation error of each Elementwise1D feature.
predict_std(X)	Evaluate the std of each output with numpy, shape (n_samples, n_outputs).
scale_inputs(mean, scale)	Network on raw inputs x for one trained on (x - mean) / scale.
scale_outputs(mean, scale)	Network predicting y * scale + mean (and std times |scale|).

encoded(input_bounds)

The network exactly as the MIP represents it for these input bounds.

Elementwise1D layers are replaced by their interpolating PiecewiseLinear layers; every other layer is unchanged. Compare encoded(bounds).forward(X) with forward(X) to see the approximation.

forward(X, return_all=False)

Evaluate the network with numpy.

Args:

X: Inputs, shape (n_samples, n_inputs). return_all: If True, also return the output of every layer.

Returns:

Outputs, shape (n_samples, n_outputs), and the list of layer outputs if return_all.

max_feature_error(input_bounds)

Bounds on the interpolation error of each Elementwise1D feature.

Returns a list with one array per Elementwise1D layer (empty when the network is exact).

predict_std(X)

Evaluate the std of each output with numpy, shape (n_samples, n_outputs).

scale_inputs(mean, scale)

Network on raw inputs x for one trained on (x - mean) / scale.

scale_outputs(mean, scale)

Network predicting y * scale + mean (and std times |scale|).

class pycse.sklearn.mip.PiecewiseLinear(breakpoints: list, values: list, W: ndarray = None, b: ndarray = None)

Sum of 1D piecewise-linear functions: z_out[j] = sum_i f_ij(z_in[i]) + z_in @ W + b.

f_ij interpolates values[i][:, j] linearly between breakpoints[i] and is constant outside them, as in a KAN layer with linear splines and a clipped grid. All functions of one input share its breakpoints, so the MIP needs one set of segment binaries per input, not per edge. Inputs may have different numbers of breakpoints.

Args:

breakpoints: Increasing breakpoints for each input: an (n_in, K) array, or a

list of n_in 1D arrays.

values: Function values at the breakpoints: an (n_in, K, n_out) array, or a

list of (K_i, n_out) arrays.

W: Optional linear term, shape (n_in, n_out). b: Optional bias, shape (n_out,).

Attributes:

W

b

n_in

n_out

Methods

__call__(z)	Call self as a function.
segments(i, lo, hi)	Breakpoints of input i restricted to [lo, hi], and the values there.

segments(i, lo, hi)

Breakpoints of input i restricted to [lo, hi], and the values there.

class pycse.sklearn.mip.ReLU(concave: ndarray = None)

Elementwise max(0, z).

Args:

concave: Optional boolean mask. Masked units compute ``-max(0, -z) =

min(0, z)`` instead, the concave counterpart used by monotonic networks.

Attributes:

concave

Methods

sign(n)

+1 for ReLU units and -1 for concave units: z_out = s * max(0, s * z).

sign(n)

+1 for ReLU units and -1 for concave units: z_out = s * max(0, s * z).

class pycse.sklearn.mip.StdSpec(tap: int, A: ndarray, c: ndarray = None, d: float = 0.0)

Predictive standard deviation sqrt(||A z + c||^2 + d).

Args:

tap: Index of the layer whose output is z (-1 means the inputs). A: Matrix, shape (k, len(z)). c: Offset, shape (k,). Defaults to zeros. d: Constant variance term (>= 0).

Attributes:

c

Methods

__call__(z)

Evaluate the std for a batch z of shape (n, len(z)).

class pycse.sklearn.mip.SymbolicModel(n_inputs: int, outputs: object, variance: object = None)

A model given by explicit formulas, exported as exact nonlinear expressions.

This is not a MILP: the Pyomo block holds y[j] == f_j(x) with the original nonlinear functions (division, exp, log, …), so it needs a nonlinear solver, e.g. Gurobi (pyo.SolverFactory("gurobi_direct_minlp")), BARON, SCIP or Ipopt (local only).

Args:

n_inputs: Number of inputs. outputs: outputs(x, ops) returns a list with one expression per output,

built from the list of inputs x and the functions in ops.

variance: Optional variance(x, ops) returning one predictive variance

expression per output.

Attributes:

std

variance

Methods

forward(X)	Evaluate with numpy, shape (n_samples, n_outputs).
predict_std(X)	Evaluate the std of each output with numpy, shape (n_samples, n_outputs).

forward(X)

Evaluate with numpy, shape (n_samples, n_outputs).

predict_std(X)

Evaluate the std of each output with numpy, shape (n_samples, n_outputs).

class pycse.sklearn.mip.TreeModel(leaves: list, split_gap: float = 0.0)

A piecewise-affine model on a partition of the input space into boxes.

Args:

leaves: The Leaf boxes. split_gap: In the MIP, a leaf’s strict lower bounds x > lower become

x >= lower + split_gap. With the default 0 a point exactly on a split may take either neighboring leaf’s value (the tree takes the left one); a positive gap rules that out, at the cost of excluding inputs within the gap above each split. It must be well above the solver’s feasibility and integrality tolerances (about 1e-6), e.g. 1e-4.

Attributes:

n_inputs

std

Methods

forward(X)

Evaluate with numpy, shape (n_samples, 1).

forward(X)

Evaluate with numpy, shape (n_samples, 1).

pycse.sklearn.mip.add_to_pyomo(model, block=None, inputs=None, input_bounds=None, include_std=False, formulation='bigm')

Write a MIPNetwork, TreeModel or SymbolicModel into a Pyomo block.

Networks and trees become MILP constraints; a SymbolicModel becomes exact nonlinear constraints (see its docstring).

Args:

model: The MIPNetwork, TreeModel or SymbolicModel to encode. block: Pyomo block (or ConcreteModel) to populate. A new ConcreteModel is

created if None.

inputs: Optional list of existing Pyomo vars to use as the inputs, e.g. the

decision variables of your own model. If None, block.x is created.

input_bounds: (lower, upper) per input. If None, the bounds of inputs are

used. Required either way: they set the big-M constants.

include_std: Add block.std[j] >= 0 for models with uncertainty

(block.variance is always added). True or "cone" adds the second-order cone block.std[j] >= sqrt(block.variance[j]): convex, and exact when the objective penalizes the std (minimize y + k * std or maximize y - k * std). "exact" adds block.std[j]**2 == block.variance[j], a nonconvex quadratic needed when the objective rewards the std (e.g. maximize y + k * std for exploration). Both need a solver for quadratic constraints, e.g. Gurobi.

formulation: "bigm" (exact MILP) or "convex" (LP without binaries;

exact for input-convex networks whose output is minimized or bounded above, see check_convex). The variance is not exported with "convex", since the hidden units are then only bounded below.

Returns:

The populated block, with block.x (the inputs) and block.y (the outputs) indexed by position.

pycse.sklearn.mip.affine_from_predict(predict, x0, step, n_check=20, rtol=1e-06, seed=0)

Recover f(x) = x @ w + c from a model that is affine, by probing predict.

The slope along each input comes from a finite difference of size step; the result is then checked at n_check random points so that a nonlinear model (e.g. a pipeline with PolynomialFeatures) raises instead of being linearized.

Args:

predict: Function mapping an (n, n_inputs) array to n predictions. x0: Center point for the probes, shape (n_inputs,). step: Probe step per input, shape (n_inputs,) (positive).

Returns:

(w, c).

Raises:

ValueError: If predict is not affine.

pycse.sklearn.mip.check_convex(network, tol=1e-12)

Raise ValueError unless formulation="convex" is exact for network.

That needs plain ReLUs (no concave units, no splines) and nonnegative weights on every hidden unit after the first ReLU, which makes the output convex and nondecreasing in every ReLU output.

pycse.sklearn.mip.eval_expression(expression, names, values, ops)

Evaluate a formula string such as "2.1*(x0*x1) - exp(x1)".

Args:

expression: The formula, using names and the functions in ops. names: Variable names, e.g. ["x0", "x1"]. values: One value per name: numpy arrays or Pyomo variables/expressions. ops: NUMPY_OPS or pyomo_ops().

pycse.sklearn.mip.llpr_std(cov, alpha_squared, zeta_squared, tap, y_scale=1.0)

StdSpec``s for LLPR: ``std_j = y_scale_j sqrt(alpha_j f^T (C + zeta_j I)^-1 f).

Args:

cov: Feature covariance F^T F, shape (n_feat, n_feat). alpha_squared, zeta_squared: Scalars or one value per output. tap: Layer index of the features f. y_scale: Output scale(s) when the targets were standardized.

Returns:

A list of StdSpec, one per output.

pycse.sklearn.mip.propagate_bounds(layers, input_bounds, network_input_bounds=None)

Interval bounds on the output of every layer.

Args:

layers: List of Affine / ReLU / PiecewiseLinear layers

(Elementwise1D layers are bounded through their interpolation).

input_bounds: Array-like of (lower, upper) per input, shape (n_inputs, 2). network_input_bounds: (lower, upper) arrays of the network inputs, for skip

connections when layers is not the whole network.

Returns:

List of (lower, upper) array pairs, one per layer.

pycse.sklearn.mip.pyomo_ops()

The functions of NUMPY_OPS as Pyomo expression builders.

pycse.sklearn.mip.std_from_covariance(cov, tap, scale=1.0, offset=None, d=0.0)

StdSpec for std = scale * sqrt((z - offset) @ cov @ (z - offset) + d).

Factors the (symmetric positive semidefinite) cov with an eigendecomposition so that (z - offset) @ cov @ (z - offset) = ||A (z - offset)||^2.