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.