API_Example

First Example: Python Calls

The web interface is a thin layer over a few Python functions in src/. You can call these functions directly from your own scripts or notebooks, with no browser and no Flask server. This page solves the smallest-enclosing-interval LP this way. It then lists the available solvers and shows how to rerun a benchmark.

The examples were run with Python 3.10.11, CVXPY 1.7.2, CLARABEL 0.11.1 and JAX 0.6.2, installed from requirements.txt. Scen-Opt needs Python 3.10 or newer; 3.12 is recommended. The outputs below are the real outputs of those runs. The examples use the open-source solver CLARABEL so that they run without a licence. If you have a MOSEK licence, pass solver="MOSEK": MOSEK is the recommended solver and was used for all published results.

How to import Scen-Opt

Scen-Opt is not a pip-installable package: setup.py is empty, and pyproject.toml only configures pytest. To use it as a library:

1. Clone the repository and install the dependencies, as on Local Install:

   git clone https://github.com/Kiguli/Scen-Opt.git
   cd Scen-Opt
   pip install -r requirements.txt

2. Import the modules under the package name src, with the repository root on the Python path. Any one of these works:

   • put your script in the repository root and run it from there;
   • run it from elsewhere with PYTHONPATH=/path/to/Scen-Opt;
   • add the root to sys.path at the top of the script, as each benchmark's run.py does.

The functions you will use most:

Function	Module	Purpose
solve_lp	src.LP	Scenario LP: solve, then compute the support list and degeneracy flag
solve_qp	src.QP	Scenario QP (adds $\tfrac12 x^\top Q x$)
solve_sdp	src.SDP	Scenario SDP with LMI constraints
quantify_risk	src.Risk	Risk bounds $(\underline{\varepsilon}, \overline{\varepsilon})$ from $k$, $N$, $\beta$
get_solvers	src.Miscellaneous	Installed CVXPY solvers and the program types each one is offered for
load_file	src.Miscellaneous	Read a .csv, .txt, .xlsx or .json file into a NumPy array
generate_matrix_function, generate_numeric_A_b, ...	src.parsing	Build $A(\delta)$, $b(\delta)$, $F_j(\delta)$ callables from strings or numeric rows
load_symbolic_program, load_scenarios_csv	benchmarks._loader	Load a program_symbolic.json and a scenarios CSV

Two things to know about these modules:

• quantify_risk computes in JAX 64-bit mode, but it switches that mode on only while it runs. Other JAX code in your process keeps JAX's default (32-bit).
• The support-list routine can print short status messages to stdout (see below).

The example problem

Given $N = 100$ points $\delta_i \in [0,1]$, find the narrowest interval $[o - w, o + w]$ that contains them all. With $x = [o, w]^\top$, this is the robust LP

$$\min_{x} \ c^\top x \quad \text{s.t.} \quad A(\delta_i) x + b(\delta_i) \le 0, \quad i = 1, \ldots, N,$$

where

$$c = \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \qquad A(\delta) = \begin{bmatrix} -1 & -1 \\ 1 & -1 \end{bmatrix}, \qquad b(\delta) = \begin{bmatrix} \delta \\ -\delta \end{bmatrix}.$$

The two rows say $o - w \le \delta_i$ and $\delta_i \le o + w$. The data are in benchmarks/LP_half_width_2d/data/scenarios.csv: 100 rows, one number per row. There are no hard constraints, $\tau = \rho = 0$ and $\beta = 10^{-6}$.

Complete script

Save this as half_width_example.py in the repository root:

import numpy as np

from src.LP import solve_lp
from src.Risk import quantify_risk
from src.Miscellaneous import load_file

# 1. Scenarios: an (N, q) array, one sampled delta per row. Here N = 100, q = 1.
deltas = load_file("benchmarks/LP_half_width_2d/data/scenarios.csv")

# 2. Scenario constraint A(delta) x + b(delta) <= 0, written as functions of one row.
def A_d(delta):
    return np.array([[-1.0, -1.0],
                     [ 1.0, -1.0]])              # (m_s, d) = (2, 2)

def b_d(delta):
    return np.array([[ delta[0]],
                     [-delta[0]]])               # (m_s, 1) = (2, 1)

# 3. Objective c^T x with x = [centre, half-width]; no hard constraints.
c = np.array([[0.0],
              [1.0]])                            # (d, 1)
G = np.array([])
h = np.array([])

# 4. Solve the robust LP (tau = rho = 0 by default).
x, zeta, cost, N, k, constraints, degeneracy = solve_lp(
    deltas, A_d, b_d, G, h, c, solver="CLARABEL")

# 5. Risk bounds at confidence 1 - beta.
beta = 1e-6
eps_lower, eps_upper = quantify_risk(k, N, beta)

print("x*         =", x.ravel())
print("cost       =", cost)
print("N          =", N)
print("k          =", k)
print("degenerate =", degeneracy)
print(f"risk       = [{eps_lower:.4f}, {eps_upper:.4f}] with confidence 1 - {beta:g}")

Run it from the repository root:

python half_width_example.py

Output:

x*         = [0.5006347  0.49830997]
cost       = 0.49830996824508517
N          = 100
k          = 2
degenerate = False
risk       = [0.0000, 0.2085] with confidence 1 - 1e-06

Reading the result:

• The optimal interval has centre $0.5006$ and half-width $0.4983$, which is $[0.0023, 0.9989]$.
• Two scenarios support it: the smallest and the largest point, so $k = 2$.
• With confidence $1 - 10^{-6}$, the probability that a new point falls outside the interval is at most $0.2085$.
• The problem is not degenerate, so the lower bound $0$ is also certified.

See Scenario Theory for what these numbers mean.

The support-list computation (get_support) prints a status line when something needs explaining; this run needed none. Three messages can appear:

• "Discarded ... spuriously screened constraint(s) ..." means an interior-point solver left a tiny dual value on a constraint that is strictly slack, and Scen-Opt recognized it as solver noise. The result is still non-degenerate. (The sweep below shows this message.)
• "PRUNE removed ... flagging degeneracy" means degeneracy was detected.
• "Dual-screened candidates did not reproduce the optimum ..." means the screen failed, and the degeneracy flag is raised as a precaution.

Inputs and outputs of solve_lp

solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
         norm_type=2, solver=None, include_slack=None)

Argument	What to pass
deltas	(N, q) NumPy array, one scenario per row. load_file returns this shape for a CSV. An empty array raises ValueError.
A_d, b_d	Callables that take one row of deltas and return arrays of shape (m_s, d) and (m_s, 1). A 1-D b_d result of length m_s is reshaped to a column.
G, h	Hard constraints $Gx + h \le 0$, shapes (n, d) and (n, 1) (a 1-D h is reshaped to a column). Pass np.array([]) for none.
c	Column vector of shape (d, 1).
tau, x_ref, norm_type	Regularization $\tau \lVert x - \bar{x} \rVert_p$, added only when tau > 0. x_ref may be a scalar (applied to every coordinate) or a length-d vector. norm_type is the order of the vector norm: any number $p \ge 1$, np.inf (or "inf"), or "fro" (the same as 2).
rho, include_slack	Relaxation $\rho \sum_i \zeta_i$. include_slack=None creates slacks exactly when rho != 0. Pass include_slack=True to force the relaxation formulation.
solver	CVXPY solver name as a string, or None to let CVXPY choose (see below).

Return value	Meaning
x	Optimal $x$, shape (d, 1).
zeta	Optimal slacks, length N. All zeros when there are no slacks.
cost	Optimal value of the full objective, including $\tau \lVert x - \bar{x} \rVert_p$ and $\rho \sum_i \zeta_i$. The web interface reports only $c^\top x$ (LP) or $c^\top x + \tfrac12 x^\top Q x$ (QP/SDP). Compute that yourself if you need it, e.g. (c.T @ x).item().
N	Number of scenarios.
complexity	$k$, the size of the support list that was found.
constraints	List of the N CVXPY scenario constraints, one per scenario. This is all of them, not the support list.
degeneracy	True if degeneracy was detected. Then eps_lower from quantify_risk is not certified.

If the solver does not return optimal or optimal_inaccurate, solve_lp raises ValueError.

quantify_risk(k, N, beta) returns the pair (eps_lower, eps_upper) as Python floats. It raises ValueError unless $N \ge 1$, $0 \le k \le N$ and $0 < \beta < 1$. It does not look at the degeneracy flag, so check degeneracy yourself before using eps_lower.

Other ways to build $A(\delta)$ and $b(\delta)$

From strings, the same way the web interface's symbolic mode does. A restricted evaluator accepts numbers, + - * / // % ** (and ^ for power), delta[i], common functions and constants such as sin, exp, sqrt, abs, pi and e, and math.<function>(...); the full list is on Miscellaneous.py. Entries may also be plain numbers:

from src.parsing import generate_matrix_function

A_d = generate_matrix_function("[[-1, -1], [1, -1]]")
b_d = generate_matrix_function("[['delta[0]'], ['-delta[0]']]")

From a benchmark's program file. benchmarks/_loader.py parses program_symbolic.json into ready-to-use arrays and callables:

from benchmarks._loader import load_symbolic_program, load_scenarios_csv
from src.LP import solve_lp
from src.Risk import quantify_risk

data_dir = "benchmarks/LP_half_width_2d/data"
prog = load_symbolic_program(data_dir)        # parses program_symbolic.json
deltas = load_scenarios_csv(data_dir)          # reads scenarios.csv -> (N, q) array

x, zeta, cost, N, k, constraints, degenerate = solve_lp(
    deltas, prog["A_d"], prog["b_d"], prog["G"], prog["h"], prog["c"],
    tau=prog["tau"], rho=prog["rho"], solver="CLARABEL")
print(prog["type"], prog["option"], prog["confidence"])
print(x.ravel(), k, quantify_risk(k, N, prog["confidence"]))

LP robust 1e-06
[0.5006347  0.49830997] 2 (0.0, 0.20852406037505716)

(The support-list status line is omitted here and below.) In program_symbolic.json, the confidence field holds $\beta$ itself.

From pre-evaluated numeric rows, as in the web interface's numeric mode. Each row of scenarios_numeric.csv is $[A_i \text{ flattened row-wise} \mid b_i]$:

import numpy as np
from src.LP import solve_lp
from src.Risk import quantify_risk
from src.Miscellaneous import load_file
from src.parsing import generate_numeric_A_b

rows = load_file("benchmarks/LP_half_width_2d/data/scenarios_numeric.csv")   # N x (m*d + m)
A_d, b_d = generate_numeric_A_b(rows_A=2, n_x=2)
c = np.array([[0.0], [1.0]])
x, zeta, cost, N, k, constraints, degenerate = solve_lp(
    rows, A_d, b_d, np.array([]), np.array([]), c, solver="CLARABEL")
print(rows.shape, x.ravel(), k, quantify_risk(k, N, 1e-6))

(100, 6) [0.5006347  0.49830997] 2 (0.0, 0.20852406037505716)

Choosing a solver

get_solvers() lists the CVXPY solvers installed on your machine. For each one it gives the program types Scen-Opt offers it for; the web interface uses the same list to fill its solver menu.

from src.Miscellaneous import get_solvers

solvers = get_solvers()
print(solvers)
print("LP-capable:", [name for name, kinds in solvers.items() if "LP" in kinds])

With the packages from requirements.txt, the output was:

{'CLARABEL': ['LP', 'QP', 'SDP'], 'MOSEK': ['LP', 'QP', 'SDP'], 'OSQP': ['QP'], 'SCIPY': ['LP'], 'SCS': ['LP', 'QP', 'SDP']}
LP-capable: ['CLARABEL', 'MOSEK', 'SCIPY', 'SCS']

requirements-solvers.txt lists further open-source solvers (CVXOPT, ECOS, PIQP, PROXQP, HiGHS, DAQP). For an installed solver that is missing from Scen-Opt's capability table, the program types are taken from CVXPY's own solver metadata (see Miscellaneous.py).

Pass the chosen name as solver="...". Some practical points:

• solver=None lets CVXPY choose. requirements.txt installs the Mosek package, and CVXPY may then pick MOSEK. Without a valid licence, the solve fails with a licence error; the example above failed this way with an expired licence. Name the solver explicitly to avoid this.
• Accuracy matters for $k$. The support-list routine compares optimal values to a relative tolerance of $10^{-6}$, so low-accuracy solvers can inflate $k$. On this example, SCS, SCIPY and CLARABEL all return $k = 2$. On the 200-point example of LP.py, however, the first-order solver SCS returned $k = 31$ with the degeneracy flag raised, while CLARABEL returned $k = 2$. Prefer an interior-point solver such as MOSEK or CLARABEL.
• LP-only and QP-only solvers. Without regularization ($\tau = 0$) a scenario LP is a pure LP, so LP solvers such as SCIPY can solve it. The regularization term with $p \ne 1, \infty$ adds a cone that SCIPY and OSQP cannot handle.
• MOSEK needs a licence file (by default ~/mosek/mosek.lic). The benchmark run.py scripts use MOSEK unless you pass --solver; see below.

Relaxation and parameter sweeps

The web interface's comma-separated $\rho$ and $\tau$ lists are just a loop over solve_lp. Here is the Relaxation formulation of the same example for four values of $\rho$:

import numpy as np

from src.LP import solve_lp
from src.Risk import quantify_risk
from src.Miscellaneous import load_file

deltas = load_file("benchmarks/LP_half_width_2d/data/scenarios.csv")
A_d = lambda d: np.array([[-1.0, -1.0], [1.0, -1.0]])
b_d = lambda d: np.array([[d[0]], [-d[0]]])
c = np.array([[0.0], [1.0]])
G, h = np.array([]), np.array([])

rhos = [0.05, 0.1, 0.25, 1.0]
beta = 1e-6
for rho in rhos:
    x, zeta, cost, N, k, _, degeneracy = solve_lp(
        deltas, A_d, b_d, G, h, c, rho=rho, include_slack=True, solver="CLARABEL")
    lo, up = quantify_risk(k, N, beta)
    violated = int(np.sum(zeta > 1e-7))
    print(f"rho={rho:<5} c'x={(c.T @ x).item():.4f}  violated={violated:<3} k={k:<3} "
          f"risk=[{lo:.4f}, {up:.4f}]")

print("overall confidence:", 1 - beta * len(rhos))

rho=0.05  c'x=0.4205  violated=20  k=20  risk=[0.0433, 0.4719]
rho=0.1   c'x=0.4644  violated=10  k=10  risk=[0.0000, 0.3432]
Discarded 1 spuriously screened constraint(s) (strictly slack, near-threshold duals); remaining support scenarios are irreducible: non-degenerate.
rho=0.25  c'x=0.4911  violated=4   k=4   risk=[0.0000, 0.2476]
Discarded 1 spuriously screened constraint(s) (strictly slack, near-threshold duals); remaining support scenarios are irreducible: non-degenerate.
rho=1.0   c'x=0.4983  violated=0   k=2   risk=[0.0000, 0.2085]
overall confidence: 0.999996

Some notes on this output:

• Violated scenarios always belong to the support list, so a cheaper interval (small $\rho$) comes with a larger $k$ and a higher certified risk.
• The last line applies the union bound: all four certificates hold together with confidence $1 - 4\beta$. The web interface does this bookkeeping for you; in your own loops you must do it yourself.
• With include_slack=True and rho=0, this problem is unbounded and solve_lp raises ValueError (status unbounded).
• These values of $\rho$ are of the form $0.5/n$, for which several intervals are optimal (see First Example: Web Interface). The optimal objective is unique, but $c^\top x$ can differ in the last digit between solvers or versions.

See Scenario Theory for the background.

QP and SDP

The other two solvers follow the same pattern and return the same seven values:

solve_qp(deltas, A_d, b_d, G, h, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
         norm_type=2, solver=None, include_slack=None)

solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
          norm_type=2, solver=None, include_slack=None)

• Q must be symmetric positive semidefinite, up to a small tolerance (so rank-deficient matrices such as $D^\top D$ are accepted); otherwise an AssertionError is raised.
• For solve_sdp:
  • F_d(delta) returns a dict of symmetric matrices keyed '0' (constant term), '1', ..., 'd'; integer keys 0, 1, ..., d work too. The scenario constraint is $F_0(\delta_i) + \sum_j x_j F_j(\delta_i) \preceq \zeta_i I$.
  • E is a dict in the same format for the hard LMI. Pass None or {} for none.
  • c is a vector of shape (d,) or (d, 1); it sets $d$.
  • Q may be None or an empty array when there is no quadratic term.

See Quadratic Programs and Semidefinite Programs for worked examples.

Reproducing a benchmark

Each folder in benchmarks/ has a run.py that solves the program and a plot.py that draws its figure (see In-built Graphics). Most also have a generate.py that regenerates the data, which is already shipped. For the interval example:

cd benchmarks/LP_half_width_2d
python run.py

What run.py does:

• It adds the repository root to sys.path itself, so you can run it from its own folder.
• It loads data/program_symbolic.json through benchmarks/_loader.py, and loads data/scenarios.csv and parameters.txt. In parameters.txt, confidence is $1 - \beta$ (here 0.999999).
• It prints the results and overwrites results/metrics.json and results/solution.csv.

Every run.py takes --solver, with MOSEK as the default because it produced the published results. Without a valid MOSEK licence the default run stops with a message such as Solver failed: rescode.err_license_expired(1001): The license has expired. To use an open-source solver, pass it on the command line:

python run.py --solver CLARABEL

The same works from the repository root, e.g. python benchmarks/LP_half_width_2d/run.py --solver CLARABEL. With CLARABEL the interval benchmark prints (after a header and "Solving LP with CLARABEL..."):

-----------------------------------------------------------------
                      OPTIMIZATION RESULTS                       
-----------------------------------------------------------------
Status: SUCCESS
Data Points (N): 100
Optimal Half-width: 0.498310
Complexity (k): 2 support constraints
Risk Bounds (99.9999%): [0.0000, 0.2085]
-----------------------------------------------------------------

OPTIMAL INTERVAL:
  Center:     0.500635
  Half-width: 0.498310
  Interval:   [0.002325, 0.998945]
  All points contained: True
  Data range: [0.002325, 0.998945]
  Support scenarios: 2 points at boundaries

These numbers match the shipped results/metrics.json and the benchmark table in the README. When degeneracy is detected, run.py prints the bounds as [unreliable, <upper>].

Then run python plot.py to redraw the figure from results/metrics.json. This needs matplotlib, which is not in requirements.txt, so install it first with pip install matplotlib.