diff --git a/.github/workflows/docker-publish.yml b/.github/workflows/docker-publish.yml
index 8531c72..fa04e20 100644
--- a/.github/workflows/docker-publish.yml
+++ b/.github/workflows/docker-publish.yml
@@ -22,6 +22,10 @@ jobs:
       - name: Checkout
         uses: actions/checkout@v4
 
+      # QEMU emulates arm64 so the image can also be built for it.
+      - name: Set up QEMU
+        uses: docker/setup-qemu-action@v3
+
       # Buildx is needed for multi-arch builds and layer caching.
       - name: Set up Docker Buildx
         uses: docker/setup-buildx-action@v3
@@ -54,7 +58,6 @@ jobs:
           labels: ${{ steps.meta.outputs.labels }}
           cache-from: type=gha
           cache-to: type=gha,mode=max
-          # Publish for x86_64 Linux hosts. Add ,linux/arm64 here if you
-          # also want ARM (e.g. M-series Mac / Pi 4+) — it roughly doubles
-          # build time because there's no native arm64 runner on public GH.
-          platforms: linux/amd64
+          # x86_64 hosts and ARM hosts (e.g. Apple silicon Macs). The arm64
+          # build runs under QEMU, which roughly doubles the build time.
+          platforms: linux/amd64,linux/arm64
diff --git a/.github/workflows/docker.yml b/.github/workflows/docker.yml
deleted file mode 100644
index b551a85..0000000
--- a/.github/workflows/docker.yml
+++ /dev/null
@@ -1,52 +0,0 @@
-name: Docker
-
-on:
-  push:
-    branches: [master]
-    paths:
-      - 'app.py'
-      - 'src/**'
-      - 'templates/**'
-      - 'static/**'
-      - 'requirements.txt'
-      - 'requirements-solvers.txt'
-      - 'Dockerfile'
-
-env:
-  REGISTRY: ghcr.io
-  IMAGE_NAME: ${{ github.repository }}
-
-jobs:
-  build-and-push:
-    runs-on: ubuntu-22.04
-
-    permissions:
-      contents: read
-      packages: write
-
-    steps:
-      - uses: actions/checkout@v4
-
-      - name: Log in to GitHub Container Registry
-        uses: docker/login-action@v3
-        with:
-          registry: ${{ env.REGISTRY }}
-          username: ${{ github.actor }}
-          password: ${{ secrets.GITHUB_TOKEN }}
-
-      - name: Extract metadata
-        id: meta
-        uses: docker/metadata-action@v5
-        with:
-          images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }}
-          tags: |
-            type=raw,value=latest
-            type=sha,prefix=
-
-      - name: Build and push
-        uses: docker/build-push-action@v6
-        with:
-          context: .
-          push: true
-          tags: ${{ steps.meta.outputs.tags }}
-          labels: ${{ steps.meta.outputs.labels }}
diff --git a/DOCKER.md b/DOCKER.md
index ff49258..24f90a3 100644
--- a/DOCKER.md
+++ b/DOCKER.md
@@ -14,6 +14,8 @@ Run the pre-built image from the GitHub Container Registry (use a release tag su
 docker run -p 5000:5000 ghcr.io/kiguli/scen-opt:latest
 ```
 
+The image is built for both x86_64 and ARM (e.g. Apple silicon) machines.
+
 Or build the image yourself:
 
 ```bash
diff --git a/Dockerfile b/Dockerfile
index 234bb75..17b7418 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -30,4 +30,7 @@ USER appuser
 
 EXPOSE 5000
 
-CMD ["gunicorn", "--bind", "0.0.0.0:5000", "--workers", "2", "--timeout", "7200", "app:app"]
+# Default worker count. Override with -e GUNICORN_CMD_ARGS="--workers 4".
+ENV GUNICORN_CMD_ARGS="--workers 2"
+
+CMD ["gunicorn", "--bind", "0.0.0.0:5000", "--timeout", "7200", "app:app"]
diff --git a/README.md b/README.md
index 8bd1731..1e37b04 100644
--- a/README.md
+++ b/README.md
@@ -10,7 +10,7 @@ The tool is implemented as a Python Flask web application with a modern JavaScri
 
 ## Installation
 
-**Prerequisites:** Python 3.12+ and pip.
+**Prerequisites:** Python 3.10 or newer (3.12 recommended) and pip.
 
 ```bash
 # Clone the repository
@@ -26,7 +26,7 @@ python3 app.py
 
 The server starts at `http://127.0.0.1:5000`. Open this URL in a browser to access the web interface.
 
-**Optional:** For best performance on SDP problems, install [MOSEK](https://www.mosek.com/) (free academic license available). See [MOSEK Support](#mosek-support) for license configuration.
+**MOSEK:** [MOSEK](https://www.mosek.com/) is the recommended solver: it is the most stable for these problems and was used for all published results. Its Python package is installed by `requirements.txt`, but it needs a licence (free for academics). See [MOSEK Support](#mosek-support) for license configuration.
 
 ## Quick Start
 
@@ -94,7 +94,7 @@ where the scenario constraints are Linear Matrix Inequalities (LMIs) parameteriz
 | Scenarios | delta_i | Sampled uncertainty realizations (N samples) |
 | Slack penalty | rho | Penalty on constraint relaxation (rho = 0 means hard scenario constraints) |
 | Regularization | tau | Regularization strength toward reference point x_ref |
-| Norm type | p | Norm for regularization (1, 2, inf, fro, nuc) |
+| Norm type | p | Norm for regularization (any p ≥ 1, inf, or fro) |
 | Confidence | beta | Confidence parameter for risk bounds (e.g. beta = 0.01 for 99% confidence) |
 
 ### Risk Quantification
@@ -170,7 +170,6 @@ Each benchmark follows a standardized layout:
 │   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── scenarios.csv   Sampled uncertainty realizations
 │   ├── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
-│   ├── c.csv           Objective vector
 │   └── ...             Additional problem-specific data
 └── results/
     ├── metrics.json    Solver output (cost, risk bounds, complexity)
@@ -183,13 +182,13 @@ Each benchmark follows a standardized layout:
 ```bash
 cd benchmarks/<benchmark_name>
 
-# Regenerate scenarios (optional, data already provided)
+# Regenerate scenarios (optional, data already provided; some generators need extra packages, listed at the top of the script)
 python generate.py
 
-# Solve the program
+# Solve the program (MOSEK by default; e.g. --solver CLARABEL without a MOSEK licence)
 python run.py
 
-# Generate the paper figure
+# Generate the paper figure (needs matplotlib: pip install matplotlib)
 python plot.py
 ```
 
@@ -269,7 +268,7 @@ Identifies the support list of a solved scenario program: the constraints whose
 
 - **License upload:** When MOSEK is selected as the solver, a modal prompts you to upload your `mosek.lic` file. The license is **not** stored on the server.
 - **Browser caching:** After the first upload, the license is cached in your browser's session storage. Subsequent solves skip the modal automatically. A cache indicator and "Clear" button appear below the solver dropdown.
-- **Local license:** If MOSEK is installed locally with a system license (e.g., at `~/mosek/mosek.lic`), click "Skip" in the modal to use it directly.
+- **Local license:** If MOSEK is installed locally with a system license (e.g., at `~/mosek/mosek.lic`), click "I have a local license" in the modal to use it directly.
 - **Concurrent users:** Each MOSEK solve with an uploaded license runs in an isolated subprocess with the license written to an ephemeral `/tmp` directory (RAM-backed, auto-purged). Multiple users can solve simultaneously without conflicts.
 
 ## Docker
diff --git a/app.py b/app.py
index 9849f82..b4b8a8b 100644
--- a/app.py
+++ b/app.py
@@ -7,16 +7,9 @@
 import subprocess
 import sys
 import tempfile
-import time
 import uuid
-from src.LP import solve_lp
-from src.QP import solve_qp
-from src.SDP import solve_sdp
-from src.Risk import quantify_risk
 from src.Miscellaneous import get_solvers
-from src.parsing import (generate_matrix_function, generate_matrix,
-                         generate_tensor_function, generate_tensor,
-                         generate_numeric_A_b, generate_numeric_F)
+from src.solve_form import solve_form
 
 app = Flask(__name__)
 
@@ -129,8 +122,17 @@ def download_mat():
     if not data:
         return jsonify({"error": "No data provided"}), 400
 
+    def _none_to_nan(v):
+        # A failed run has no result (None), which savemat cannot store.
+        if v is None:
+            return np.nan
+        if isinstance(v, list):
+            return [_none_to_nan(item) for item in v]
+        return v
+
     mat_dict = {}
     for key, value in data.items():
+        value = _none_to_nan(value)
         try:
             mat_dict[key] = np.array(value)
         except (ValueError, TypeError):
@@ -155,7 +157,10 @@ def solve():
     )
 
     if not has_license:
-        return _solve_inner()
+        try:
+            return _solve_inner()
+        except (ValueError, TypeError) as error:
+            return jsonify({"error": str(error)}), 400
 
     # ---- MOSEK subprocess path ----
     # Use a platform-neutral temp directory (tempfile.mkdtemp → %TEMP% on
@@ -215,69 +220,11 @@ def solve():
         shutil.rmtree(tmp_dir, ignore_errors=True)
 
 
-def _parse_theta_bar(raw):
-    """Parse the x̄ (reference point) form field.
-
-    The field is populated by the matrix-grid modal as a JSON-encoded 2-D
-    array (e.g. ``[["10"]]`` or ``[[0], [1]]``). Older inputs or an empty
-    field may arrive as an empty/comma-separated string. This helper
-    accepts any of those and returns a 1-D numpy float array, or 0.0 when
-    the field is unset (backward compat with the previous default).
-    """
-    if raw is None:
-        return 0.0
-    s = raw.strip()
-    if s in ('', 'null', '""', '[]'):
-        return 0.0
-    # Try JSON first (the standard path for matrix-modal entries).
-    try:
-        parsed = json_module.loads(s)
-        return np.asarray(parsed, dtype=float).reshape(-1)
-    except (ValueError, TypeError):
-        pass
-    # Fallback: legacy comma-separated scalars (e.g. "0, 0").
-    try:
-        return np.array([float(x) for x in s.split(',') if x.strip() != ''])
-    except ValueError as e:
-        raise ValueError(
-            f"Could not parse x̄ (theta_bar) value {raw!r}: expected a JSON "
-            f"array (e.g. [[10]]) or comma-separated numbers. ({e})"
-        )
-
-
-def compute_base_cost(active_tab, x, c, Q):
-    """Compute c'x (LP) or c'x + ½x'Qx (QP/SDP) from the optimal x.
-
-    Regularization (τ·‖x − x̄‖) and relaxation (ρ·Σζ_i) penalty terms are
-    numerical tools for informing the solution; they are not part of the
-    underlying problem's objective and must not leak into reported cost.
-    """
-    if x is None:
-        return None
-    x_arr = np.asarray(x, dtype=float).ravel()
-    if x_arr.size == 0:
-        return None
-    c_arr = np.asarray(c, dtype=float).ravel() if c is not None and getattr(c, 'size', 0) else None
-    if c_arr is None or c_arr.size != x_arr.size:
-        return None
-    base = float(c_arr @ x_arr)
-    if active_tab in ('qp-tab', 'sdp-tab') and Q is not None and getattr(Q, 'size', 0):
-        Q_arr = np.asarray(Q, dtype=float)
-        if Q_arr.shape == (x_arr.size, x_arr.size):
-            base += 0.5 * float(x_arr @ Q_arr @ x_arr)
-    return base
-
-
 def _solve_inner():
-    option = request.form.get('option')
-    active_tab = request.form.get('active_tab')
-
     # Handle scenarios file if present
     scenarios = None
     if 'file' in request.files and request.files['file'].filename != '':
-        file = request.files['file']
-        filename = file.filename.lower()
-        parsed = parse_uploaded_file(file)
+        parsed = parse_uploaded_file(request.files['file'])
         if isinstance(parsed, dict):
             # MAT files return a dict; use the first array-like value
             for k, v in parsed.items():
@@ -289,220 +236,7 @@ def _solve_inner():
         else:
             scenarios = np.array(parsed, dtype=float)
 
-    mode = request.form.get('mode', 'symbolic')
-    n_x = int(request.form.get('n_x', 0)) if request.form.get('n_x') else 0
-    rows_A = int(request.form.get('rows_A', 0)) if request.form.get('rows_A') else 0
-    lmi_size = int(request.form.get('lmi_size', 0)) if request.form.get('lmi_size') else 0
-
-    A_d = None
-    b_d = None
-    F_d = None
-
-    if mode == 'numeric':
-        if scenarios is None:
-            raise ValueError("Numeric mode requires an uploaded scenario data file.")
-        if active_tab in ('lp-tab', 'qp-tab'):
-            if n_x <= 0 or rows_A <= 0:
-                raise ValueError("Numeric mode requires valid decision variable count (n_x) and soft constraint rows (rows_A).")
-            expected_cols = rows_A * n_x + rows_A
-            if scenarios.shape[1] != expected_cols:
-                raise ValueError(
-                    f"Each scenario row must have {expected_cols} elements "
-                    f"(rows_A*n_x + rows_A = {rows_A}*{n_x} + {rows_A}), "
-                    f"but got {scenarios.shape[1]}.")
-            A_d, b_d = generate_numeric_A_b(rows_A, n_x)
-        elif active_tab == 'sdp-tab':
-            if n_x <= 0 or lmi_size <= 0:
-                raise ValueError("Numeric mode requires valid decision variable count (n_x) and LMI matrix size.")
-            expected_cols = (n_x + 1) * lmi_size ** 2
-            if scenarios.shape[1] != expected_cols:
-                raise ValueError(
-                    f"Each scenario row must have {expected_cols} elements "
-                    f"((n_x+1)*lmi_size² = {n_x + 1}*{lmi_size}² = {expected_cols}), "
-                    f"but got {scenarios.shape[1]}.")
-            F_d = generate_numeric_F(lmi_size, n_x)
-    else:
-        if request.form.get('A_d'):
-            A_d = generate_matrix_function(request.form.get('A_d'))
-        else:
-            if (active_tab == 'lp-tab') or (active_tab == 'qp-tab'):
-                raise ValueError("\\(A(\\delta)\\) is ill-defined")
-        if request.form.get('b_d'):
-            b_d = generate_matrix_function(request.form.get('b_d'))
-        else:
-            if (active_tab == 'lp-tab') or (active_tab == 'qp-tab'):
-                raise ValueError("\\(b(\\delta)\\) is ill-defined")
-        if request.form.get('F_d'):
-            F_d = generate_tensor_function(request.form.get('F_d'))
-        else:
-            if (active_tab == 'sdp-tab'):
-                raise ValueError("\\(F_j(\\delta)\\) is ill-defined")
-    if request.form.get('A_da'):
-        A_da = generate_tensor_function(request.form.get('A_da'))
-    if request.form.get('b_da'):
-        b_da = generate_tensor_function(request.form.get('b_da'))
-
-    # Hard LP/QP constraints use G, h (form fields renamed from legacy A, b).
-    G = generate_matrix(request.form.get('G')) if request.form.get('G') else np.array([])
-    h_vec = generate_matrix(request.form.get('h')) if request.form.get('h') else np.array([])
-    c = generate_matrix(request.form.get('c')) if request.form.get('c') else np.array([])
-    Q = generate_matrix(request.form.get('Q')) if request.form.get('Q') else np.array([])
-    C = generate_matrix(request.form.get('C')) if request.form.get('C') else np.array([])
-    # Hard SDP LMI collection uses E (form field renamed from legacy F).
-    E = generate_tensor(request.form.get('E')) if request.form.get('E') else {}
-    A_a = generate_tensor(request.form.get('A_a')) if request.form.get('A_a') else {}
-    b_a = generate_matrix(request.form.get('b_a')) if request.form.get('b_a') else np.array([])
-
-    conf = float(request.form.get('confidence')) if request.form.get('confidence') else 0.0
-    # Get values from parameter boxes
-    form_data = request.form.to_dict()
-
-    rhos = np.array([float(x) for x in request.form.get('rho', 0).split(',')]) if request.form.get('rho') else np.array(
-        [0.0])
-    taus = np.array([float(x) for x in request.form.get('tau', 0).split(',')]) if request.form.get('tau') else np.array(
-        [0.0])
-    theta_bar = _parse_theta_bar(request.form.get('theta_bar'))
-    p_raw = request.form.get('p', '').strip().strip("'\"")
-    if not p_raw:
-        p = 2
-    elif p_raw.lower() in ('inf', 'infinity', 'np.inf'):
-        p = np.inf
-    elif p_raw.lower() == 'fro':
-        p = 'fro'
-    else:
-        try:
-            p = float(p_raw)
-        except ValueError:
-            raise ValueError(
-                f"p (norm order) must be a single number, 'inf', or 'fro'. "
-                f"You entered {p_raw!r} — comma-separated lists are not "
-                f"supported for p (sweep over τ or ρ instead)."
-            )
-
-    # Initialize lists to collect results for each run
-    optimal_x_list = []
-    optimal_s_list = []
-    optimal_cost_list = []
-    N_list = []
-    active_list = []
-    constraints_list = []
-    risk_list = []
-    e_list = []
-    degeneracy_list = []
-    rho_list = []
-    tau_list = []
-    solve_time_list = []
-    risk_time_list = []
-
-    # Overall confidence after q sweep attempts is 1 − β·q (union bound).
-    # Pass the user-entered β directly to the per-attempt risk computation;
-    # only the *reported* confidence folds in q.
-    n_attempts = len(taus) * len(rhos)
-    conf_reported = conf * n_attempts
-
-    # The relaxation formulation must keep its slack variables ζ_i even when
-    # the user has set ρ = 0 (treats that case as genuinely unbounded rather
-    # than silently coercing it into the robust LP). When `option` is not in
-    # the form payload (older clients / direct API callers), fall through to
-    # the legacy "infer slack from ρ" behaviour by leaving include_slack=None
-    # so each solver default kicks in.
-    selected_option = form_data.get('option')
-    if selected_option:
-        include_slack = selected_option in ('relaxation', 'regularization-relaxation')
-    else:
-        include_slack = None
-
-    for j in range(len(taus)):
-        tau = taus[j]
-        for i in range(len(rhos)):
-            rho = rhos[i]
-            solver = form_data.get('solver', 'CLARABEL')
-
-            optimal_x = np.array([])
-            optimal_s = np.array([])
-            optimal_cost = None
-            N = 0
-            complexity = []
-            constraints = []
-            degeneracy = False
-            risk = np.array([])
-            e = "None"
-            solve_time = 0.0
-            risk_time = 0.0
-
-            try:
-                t0 = time.perf_counter()
-                if active_tab == 'lp-tab':
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_lp(
-                        scenarios, A_d, b_d, G, h_vec, c, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack)
-                elif active_tab == 'qp-tab':
-                    Q = generate_matrix(request.form.get('Q')) if request.form.get('Q') else np.array([])
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_qp(
-                        scenarios, A_d, b_d, G, h_vec, c, Q, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack)
-                elif active_tab == 'sdp-tab':
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_sdp(
-                        scenarios, F_d, E, c, Q, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack)
-                solve_time = time.perf_counter() - t0
-
-                # Replace full-objective cost (which includes τ·‖x−x̄‖_p and
-                # ρ·Σζ_i penalty terms used only as solver guidance) with the
-                # base cost: c'x for LP, c'x + ½x'Qx for QP/SDP.
-                optimal_cost = compute_base_cost(active_tab, optimal_x, c, Q)
-
-                t1 = time.perf_counter()
-                risk = np.array(quantify_risk(complexity, N, conf))
-                risk_time = time.perf_counter() - t1
-            except Exception as error:
-                e = str(error)
-                # Shorten verbose MOSEK license errors
-                if 'license' in e.lower() and ('mosek' in e.lower() or 'rescode' in e.lower()):
-                    e = "MOSEK license cannot be located or is incorrect."
-                print(e)
-
-            # Append results for this run
-            optimal_x_list.append(optimal_x.tolist() if hasattr(optimal_x, 'tolist') else optimal_x)
-            optimal_s_list.append(optimal_s.tolist() if hasattr(optimal_s, 'tolist') else optimal_s)
-            optimal_cost_list.append(optimal_cost)
-            N_list.append(N)
-            active_list.append(complexity)
-            e_list.append(e)
-            constraints_list.append(len(constraints))
-            risk_list.append(risk.tolist() if hasattr(risk, 'tolist') else risk)
-            degeneracy_list.append(degeneracy)
-            rho_list.append(rho)
-            tau_list.append(tau)
-            solve_time_list.append(round(solve_time, 4))
-            risk_time_list.append(round(risk_time, 4))
-
-    print(tau_list)
-    print(rho_list)
-
-    # Prepare the result dictionary with lists for each parameter
-    result = {
-        "form_data": form_data,
-        "optimal_x": optimal_x_list,
-        "optimal_s": optimal_s_list,
-        "optimal_cost": optimal_cost_list,
-        "num_deltas": N_list,
-        "tot_con": constraints_list,
-        "active_con": active_list,
-        "risk": risk_list,
-        "conf": 1 - conf_reported,
-        "tau_": tau_list,
-        "rho_": rho_list,
-        "errorcode": e_list,
-        "degeneracy": degeneracy_list,
-        "solve_time": solve_time_list,
-        "risk_time": risk_time_list,
-        "n_attempts": n_attempts,
-    }
-
-    print(result)
-
-    return jsonify(result)
+    return jsonify(solve_form(request.form.to_dict(), scenarios))
 
 
 if __name__ == '__main__':
diff --git a/benchmarks/LP_growth_bound_12d/README.md b/benchmarks/LP_growth_bound_12d/README.md
index c577866..ed3ef60 100644
--- a/benchmarks/LP_growth_bound_12d/README.md
+++ b/benchmarks/LP_growth_bound_12d/README.md
@@ -20,10 +20,10 @@ This is a Linear Program with 12 decision variables: 9 entries of the growth mat
 
 The LP solves with N = 3127 sampled trajectories. The optimal growth matrix and bias are:
 
-    M* = [[-1.000, 0.000, 0.005], [0.000, -1.000, 0.009], [0.000, 0.000, -1.000]]
-    b* = [-1.000, -1.000, -1.000]
+    M* = [[1.000, 0.000, 0.005], [0.000, 1.000, 0.009], [0.000, 0.000, 1.000]]
+    b* = [1.000, 1.000, 1.000]
 
-With pre-selected bias γ = 0.067 and β = 10⁻⁶, the support set size is k = 6 and degeneracy was detected, so no lower bound on risk can be certified. The upper bound is ε̄ = 0.0101.
+With pre-selected bias γ = 0.067 and β = 10⁻⁶, the support set size is k = 6 and degeneracy was detected, so no lower bound on risk can be certified. The upper bound is ε̄ = 0.0103.
 
 ## Files
 
@@ -35,15 +35,12 @@ LP_growth_bound_12d/
 ├── run.py              Solve the LP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json         One-shot program definition (JSON)
-│   ├── program_symbolic.json          One-shot program definition (MATLAB)
-│   ├── A_d.csv                Scenario-dependent constraint matrix (3 × 12 expressions)
-│   ├── b_d.csv                Scenario-dependent RHS (3 × 1 expressions)
-│   ├── c.csv                  Objective vector (12 × 1)
-│   ├── G.csv                  Hard constraint matrix (6 × 12)
-│   ├── h.csv                  Hard constraint RHS (6 × 1)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── growth_bound.csv       Full dataset (3127 × 6)
-│   └── growth_bound_mini.csv  Mini dataset (99 × 6)
+│   ├── growth_bound_mini.csv  Mini dataset (99 × 6)
+│   ├── scenarios.csv          Full dataset for the web interface (same as growth_bound.csv)
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (growth matrix, bias, risk bounds)
     ├── solution.csv    Raw solution vector
@@ -57,9 +54,10 @@ LP_growth_bound_12d/
 # Regenerate trajectory data
 python generate.py
 
-# Solve the LP (use --mini for faster testing)
+# Solve the LP (use --mini for faster testing; MOSEK by default, use --solver to pick another, e.g. CLARABEL)
 python run.py
 python run.py --mini
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -71,25 +69,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **LP** tab, formulation: **Robust**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
 3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: -2.9862583119393014
+- Optimal cost: 3.0137485232272665
 - Complexity k: 6
-- Risk bounds: [0.0, 0.010092109023134854]
+- Risk bounds: [0.0, 0.010266487573906823] (degeneracy detected, so the lower bound is not certified)
diff --git a/benchmarks/LP_growth_bound_12d/parameters.txt b/benchmarks/LP_growth_bound_12d/parameters.txt
index 46f9db5..0cd6527 100644
--- a/benchmarks/LP_growth_bound_12d/parameters.txt
+++ b/benchmarks/LP_growth_bound_12d/parameters.txt
@@ -1,3 +1,4 @@
 rho = 0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/LP_growth_bound_12d/plot.py b/benchmarks/LP_growth_bound_12d/plot.py
index 66b4168..1646760 100644
--- a/benchmarks/LP_growth_bound_12d/plot.py
+++ b/benchmarks/LP_growth_bound_12d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the growth bound benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/LP_growth_bound_12d/run.py b/benchmarks/LP_growth_bound_12d/run.py
index b217b5e..e2fa1ec 100644
--- a/benchmarks/LP_growth_bound_12d/run.py
+++ b/benchmarks/LP_growth_bound_12d/run.py
@@ -3,7 +3,7 @@
 Solve the Growth Bound benchmark using the scenario approach LP solver.
 
 Usage:
-    python run.py [--mini]
+    python run.py [--mini] [--solver SOLVER]
 """
 
 import sys
@@ -67,6 +67,8 @@ def load_vector(filepath):
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument('--mini', action='store_true', help='Use mini dataset (100 samples)')
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
     args = parser.parse_args()
 
     print("=" * 65)
@@ -113,7 +115,7 @@ def main():
     print()
 
     # Solve LP
-    print("Solving LP with MOSEK...")
+    print(f"Solving LP with {args.solver}...")
     try:
         x, zeta, cost, N, k, constraints, degeneracy = solve_lp(
             deltas=scenarios,
@@ -126,7 +128,7 @@ def main():
             x_ref=np.zeros((n_vars, 1)),
             rho=rho,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
         print(f"Solved successfully")
diff --git a/benchmarks/LP_half_width_2d/README.md b/benchmarks/LP_half_width_2d/README.md
index e96bceb..dd925b9 100644
--- a/benchmarks/LP_half_width_2d/README.md
+++ b/benchmarks/LP_half_width_2d/README.md
@@ -16,7 +16,7 @@ No regularisation (τ = 0) or slack penalty (ρ = 0).
 
 ![Smallest Enclosing Interval](results/half_width.png)
 
-The LP solves with N = 100 data points. The optimal interval is [0.002, 0.999] with center 0.501 and half-width 0.498. The complexity is k = 2 (the two extreme data points define the interval). With β = 10⁻⁶, the risk bounds are ε̲ = 0 and ε̄ = 0.207.
+The LP solves with N = 100 data points. The optimal interval is [0.002, 0.999] with center 0.501 and half-width 0.498. The complexity is k = 2 (the two extreme data points define the interval). With β = 10⁻⁶, the risk bounds are ε̲ = 0 and ε̄ = 0.209.
 
 ## Files
 
@@ -27,12 +27,10 @@ LP_half_width_2d/
 ├── run.py              Solve the LP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Scenario-dependent constraint matrix (2 × 2)
-│   ├── b_d.csv         Scenario-dependent RHS (2 × 1)
-│   ├── c.csv           Objective vector (2 × 1)
-│   └── scenarios.csv   100 data points sampled from [0, 1]
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv          100 data points sampled from [0, 1]
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (center, half-width, risk bounds)
     ├── solution.csv    Raw solution vector
@@ -43,8 +41,9 @@ LP_half_width_2d/
 ## Usage
 
 ```bash
-# Solve the LP
+# Solve the LP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -56,17 +55,17 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **LP** tab, formulation: **Robust**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
 3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
@@ -75,4 +74,4 @@ python plot.py
 
 - Optimal cost: 0.498309968223239
 - Complexity k: 2
-- Risk bounds: [0.0, 0.20729497825494037]
+- Risk bounds: [0.0, 0.20852406037505716]
diff --git a/benchmarks/LP_half_width_2d/parameters.txt b/benchmarks/LP_half_width_2d/parameters.txt
index 46f9db5..0cd6527 100644
--- a/benchmarks/LP_half_width_2d/parameters.txt
+++ b/benchmarks/LP_half_width_2d/parameters.txt
@@ -1,3 +1,4 @@
 rho = 0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/LP_half_width_2d/plot.py b/benchmarks/LP_half_width_2d/plot.py
index 047b6d2..7793e12 100644
--- a/benchmarks/LP_half_width_2d/plot.py
+++ b/benchmarks/LP_half_width_2d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the smallest enclosing interval benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/LP_half_width_2d/run.py b/benchmarks/LP_half_width_2d/run.py
index 07a3ee3..c079501 100644
--- a/benchmarks/LP_half_width_2d/run.py
+++ b/benchmarks/LP_half_width_2d/run.py
@@ -3,12 +3,13 @@
 Solve the Smallest Enclosing Interval benchmark using the scenario approach LP solver.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 # Add project root to path for imports
@@ -52,6 +53,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: Smallest Enclosing Interval (LP)")
     print("=" * 65)
@@ -89,7 +95,7 @@ def main():
     beta = 1.0 - params.get('confidence', 0.999999)
 
     # Solve LP
-    print("Solving LP with MOSEK...")
+    print(f"Solving LP with {args.solver}...")
     try:
         x, zeta, cost, N, k, constraints, degeneracy = solve_lp(
             deltas=scenarios,
@@ -102,7 +108,7 @@ def main():
             x_ref=np.zeros((n_vars, 1)),
             rho=rho,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
         print(f"Solved successfully")
diff --git a/benchmarks/LP_inventory_11d/README.md b/benchmarks/LP_inventory_11d/README.md
index 51e0146..3eb2e22 100644
--- a/benchmarks/LP_inventory_11d/README.md
+++ b/benchmarks/LP_inventory_11d/README.md
@@ -26,7 +26,7 @@ Each scenario delta_i in R^15 encodes uncertain yields (delta[0:5]), demands (de
 
 ![Inventory Results](results/inventory.png)
 
-The LP solves in 9.07 seconds with N = 500 scenarios. The full budget of $12,000 is used. The penalty rho = 100 implies a cost per case of lost sales; the expected daily loss from stockouts is only $66.20 (about 0.5% of the ordering budget). The LP had complexity k = 7, no degeneracy, and risk bounds [0.000, 0.066] at 99.9999% confidence (beta = 10^{-6}).
+The LP solves in 9.07 seconds with N = 500 scenarios. The full budget of $12,000 is used. The penalty rho = 100 implies a cost per case of lost sales; the expected daily loss from stockouts is only $66.20 (about 0.5% of the ordering budget). The LP had complexity k = 7, no degeneracy, and risk bounds [0.000, 0.067] at 99.9999% confidence (beta = 10^{-6}).
 
 | Product | Order (cases) | Mean Demand | Service Level (%) |
 |---------|--------------|-------------|-------------------|
@@ -48,14 +48,10 @@ LP_inventory_11d/
 ├── run.py              Solve the LP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Augmented scenario constraint matrix (6 x 11)
-│   ├── b_d.csv         Scenario-dependent RHS vector (6 x 1)
-│   ├── c.csv           Augmented cost vector (11 x 1, includes rho)
-│   ├── G.csv           Augmented hard constraint matrix (17 x 11)
-│   ├── h.csv           Hard constraint RHS (17 x 1)
-│   └── scenarios.csv   500 x 15 uncertainty scenarios
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv          500 x 15 uncertainty scenarios
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (cost, risk bounds, etc.)
     ├── solution.csv    Raw solution vector
@@ -69,8 +65,9 @@ LP_inventory_11d/
 # Regenerate scenarios (optional)
 python generate.py --n_scenarios 500 --seed 42
 
-# Solve the LP
+# Solve the LP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -82,25 +79,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
-1. Select the **LP** tab, formulation: **Robust + Relaxation**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
-3. Set parameters: rho = 100.0, tau = 0, confidence (beta) = 1e-06
+1. Select the **LP** tab, formulation: **Robust** (the penalty rho = 100 is already in c)
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
+3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 45099.672322171755
+- Optimal cost: 45099.672322171806
 - Complexity k: 7
-- Risk bounds: [0.0, 0.06550756473687944]
+- Risk bounds: [0.0, 0.0666183782973094]
diff --git a/benchmarks/LP_inventory_11d/data/program_numeric.json b/benchmarks/LP_inventory_11d/data/program_numeric.json
index a1911fc..dac31da 100644
--- a/benchmarks/LP_inventory_11d/data/program_numeric.json
+++ b/benchmarks/LP_inventory_11d/data/program_numeric.json
@@ -311,7 +311,7 @@
       "0"
     ]
   ],
-  "rho": "100.0",
+  "rho": "0",
   "tau": "0",
   "confidence": "1e-06",
   "mode": "numeric",
diff --git a/benchmarks/LP_inventory_11d/data/program_symbolic.json b/benchmarks/LP_inventory_11d/data/program_symbolic.json
index fa4146b..d244bdf 100644
--- a/benchmarks/LP_inventory_11d/data/program_symbolic.json
+++ b/benchmarks/LP_inventory_11d/data/program_symbolic.json
@@ -411,7 +411,7 @@
       "0"
     ]
   ],
-  "rho": "100.0",
+  "rho": "0",
   "tau": "0",
   "confidence": "1e-06",
   "mode": "symbolic",
diff --git a/benchmarks/LP_inventory_11d/parameters.txt b/benchmarks/LP_inventory_11d/parameters.txt
index 8b38886..e64c9a9 100644
--- a/benchmarks/LP_inventory_11d/parameters.txt
+++ b/benchmarks/LP_inventory_11d/parameters.txt
@@ -16,4 +16,5 @@
 
 rho = 100.0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/LP_inventory_11d/plot.py b/benchmarks/LP_inventory_11d/plot.py
index 300b6f7..37b0768 100644
--- a/benchmarks/LP_inventory_11d/plot.py
+++ b/benchmarks/LP_inventory_11d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the inventory benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/LP_inventory_11d/run.py b/benchmarks/LP_inventory_11d/run.py
index c47342c..2d4463b 100644
--- a/benchmarks/LP_inventory_11d/run.py
+++ b/benchmarks/LP_inventory_11d/run.py
@@ -15,12 +15,13 @@
 subject to budget, warehouse capacity, and supplier constraints.
 
 Usage:
-    python test_and_visualize.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 # Add project root to path for imports
@@ -97,6 +98,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: Fresh Produce Distribution (LP)")
     print("Regional Distributor Ordering Under Uncertainty")
@@ -146,7 +152,7 @@ def main():
     # Data files contain the augmented formulation x_aug = [q; zeta] with
     # rho embedded in c, A_d augmented with -I slack columns, and G/h
     # extended with zeta >= 0 constraints. See generate.py for details.
-    print("Solving LP with MOSEK...")
+    print(f"Solving LP with {args.solver}...")
 
     try:
         x_full, _, cost, N_out, k, constraints, degeneracy = solve_lp(
@@ -160,7 +166,7 @@ def main():
             x_ref=np.zeros((n_aug, 1)),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         x = x_full[:n_vars]
         zeta = x_full[n_vars:]
@@ -198,9 +204,9 @@ def main():
     print(f"  - Service penalty: ${slack_cost:.2f}")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 65)
 
     print()
diff --git a/benchmarks/LP_portfolio_cvar_13d/README.md b/benchmarks/LP_portfolio_cvar_13d/README.md
index 1d55043..68ab781 100644
--- a/benchmarks/LP_portfolio_cvar_13d/README.md
+++ b/benchmarks/LP_portfolio_cvar_13d/README.md
@@ -27,7 +27,7 @@ Each scenario delta_i in R^12 encodes 8 asset returns (delta[0:8]), a market str
 
 ![CVaR Portfolio Results](results/portfolio_cvar.png)
 
-The LP solves with N = 1,255 market scenarios (one per historical trading day). The optimizer tilts heavily toward alternatives: GLD (23.2%) and VNQ (21.8%) dominate, while equities sit at the 30% regulatory minimum and fixed income at 25%. This defensive posture minimizes tail risk — the CVaR(95%) is only 1.21% daily loss. The LP had complexity k = 3, no degeneracy, and risk bounds [0.0000, 0.0200].
+The LP solves with N = 1,255 market scenarios (one per historical trading day). The optimizer tilts heavily toward alternatives: GLD (23.9%) and VNQ (21.1%) dominate, while equities sit at the 30% regulatory minimum and fixed income at 25%. This defensive posture minimizes tail risk — the CVaR(95%) is only 1.20% daily loss. The LP had complexity k = 3, no degeneracy, and risk bounds [0.0000, 0.0203].
 
 ## Files
 
@@ -39,14 +39,10 @@ LP_portfolio_cvar_13d/
 ├── run.py              Solve the LP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Augmented scenario constraint matrix (4 x 13)
-│   ├── b_d.csv         Scenario-dependent RHS vector (4 x 1, all zeros)
-│   ├── c.csv           Augmented objective vector (13 x 1, includes rho)
-│   ├── G.csv           Augmented hard constraint matrix (27 x 13)
-│   ├── h.csv           Hard constraint RHS (27 x 1)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── scenarios.csv   1255 x 12 uncertainty scenarios
+│   ├── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 │   ├── assets.csv      Asset metadata
 │   ├── historical_prices.csv   5-year ETF price history
 │   └── historical_returns.csv  Daily returns
@@ -63,8 +59,9 @@ LP_portfolio_cvar_13d/
 # Regenerate scenarios (requires yfinance and internet)
 python generate.py --seed 42
 
-# Solve the LP
+# Solve the LP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -76,25 +73,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
-1. Select the **LP** tab, formulation: **Robust + Relaxation**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
-3. Set parameters: rho = 0.015936, tau = 0, confidence (beta) = 1e-06
+1. Select the **LP** tab, formulation: **Robust** (the slack penalty is already in c)
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
+3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 0.05133721274593786
+- Optimal cost: 0.05133720473689373
 - Complexity k: 3
-- Risk bounds: [0.0, 0.02002552248652564]
+- Risk bounds: [0.0, 0.020289284798068057]
diff --git a/benchmarks/LP_portfolio_cvar_13d/data/program_numeric.json b/benchmarks/LP_portfolio_cvar_13d/data/program_numeric.json
index dce9264..05d6a88 100644
--- a/benchmarks/LP_portfolio_cvar_13d/data/program_numeric.json
+++ b/benchmarks/LP_portfolio_cvar_13d/data/program_numeric.json
@@ -29,16 +29,16 @@
       "1"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ]
   ],
   "G": [
@@ -531,7 +531,7 @@
       "0"
     ]
   ],
-  "rho": "0.040000",
+  "rho": "0",
   "tau": "0",
   "confidence": "1e-06",
   "mode": "numeric",
diff --git a/benchmarks/LP_portfolio_cvar_13d/data/program_symbolic.json b/benchmarks/LP_portfolio_cvar_13d/data/program_symbolic.json
index 728c342..8df7a75 100644
--- a/benchmarks/LP_portfolio_cvar_13d/data/program_symbolic.json
+++ b/benchmarks/LP_portfolio_cvar_13d/data/program_symbolic.json
@@ -105,16 +105,16 @@
       "1"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ],
     [
-      "0.04"
+      "0.015936"
     ]
   ],
   "G": [
@@ -607,7 +607,7 @@
       "0"
     ]
   ],
-  "rho": "0.040000",
+  "rho": "0",
   "tau": "0",
   "confidence": "1e-06",
   "mode": "symbolic",
diff --git a/benchmarks/LP_portfolio_cvar_13d/generate.py b/benchmarks/LP_portfolio_cvar_13d/generate.py
index 47a3816..0221c34 100644
--- a/benchmarks/LP_portfolio_cvar_13d/generate.py
+++ b/benchmarks/LP_portfolio_cvar_13d/generate.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires yfinance (not in requirements.txt): pip install yfinance
 """
 Golden State Teachers' Pension Fund CVaR Portfolio Benchmark
 
@@ -187,7 +188,7 @@ def save_constraint_files(data_dir, mean_returns, n_scenarios):
     With rho = 1/((1-beta)*N), minimizing alpha + rho*sum(zeta) approximates CVaR.
     """
     params = load_parameters()
-    rho = params.get('rho', 0.04)
+    rho = params.get('rho', 0.015936)
 
     n_vars = N_ASSETS + 1  # 8 weights + 1 alpha (VaR threshold)
 
diff --git a/benchmarks/LP_portfolio_cvar_13d/parameters.txt b/benchmarks/LP_portfolio_cvar_13d/parameters.txt
index 2ff2dcc..a678220 100644
--- a/benchmarks/LP_portfolio_cvar_13d/parameters.txt
+++ b/benchmarks/LP_portfolio_cvar_13d/parameters.txt
@@ -21,4 +21,5 @@
 
 rho = 0.015936
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/LP_portfolio_cvar_13d/plot.py b/benchmarks/LP_portfolio_cvar_13d/plot.py
index 7dcdbbf..237e5d5 100644
--- a/benchmarks/LP_portfolio_cvar_13d/plot.py
+++ b/benchmarks/LP_portfolio_cvar_13d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the CVaR portfolio benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/LP_portfolio_cvar_13d/results/metrics.json b/benchmarks/LP_portfolio_cvar_13d/results/metrics.json
index f3ef9b3..4d37efa 100644
--- a/benchmarks/LP_portfolio_cvar_13d/results/metrics.json
+++ b/benchmarks/LP_portfolio_cvar_13d/results/metrics.json
@@ -16,7 +16,7 @@
   "alternatives_allocation": 0.449999855590255,
   "mean_return": 0.0003031247620507457,
   "std_return": 0.005559309001023381,
-  "objective_value": 0.05335645682574842,
+  "objective_value": 0.05133720473689373,
   "N": 1255,
   "k": 3,
   "eps_lower": 0.0,
diff --git a/benchmarks/LP_portfolio_cvar_13d/run.py b/benchmarks/LP_portfolio_cvar_13d/run.py
index 2bd92f5..a7a11d4 100644
--- a/benchmarks/LP_portfolio_cvar_13d/run.py
+++ b/benchmarks/LP_portfolio_cvar_13d/run.py
@@ -3,12 +3,13 @@
 Solve the CVaR Portfolio benchmark using the scenario approach LP solver.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 import pandas as pd
 
@@ -92,6 +93,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: Pension Fund CVaR Portfolio Optimization (LP)")
     print("=" * 65)
@@ -117,7 +123,7 @@ def main():
                 key, val = line.split('=', 1)
                 params[key.strip()] = float(val.split('#')[0].strip())
 
-    rho = params.get('rho', 0.04)
+    rho = params.get('rho', 0.015936)
     tau = params.get('tau', 0.0)
     beta = 1.0 - params.get('confidence', 0.999999)
 
@@ -138,7 +144,7 @@ def main():
     # Data files contain the augmented formulation x_aug = [x; zeta] with
     # rho embedded in c, A_d augmented with -I slack columns, and G/h
     # extended with zeta >= 0 constraints. See generate.py for details.
-    print("Solving LP with MOSEK...")
+    print(f"Solving LP with {args.solver}...")
 
     try:
         x_full, _, cost, N_out, k, constraints, degeneracy = solve_lp(
@@ -152,7 +158,7 @@ def main():
             x_ref=np.zeros((n_aug, 1)),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         x = x_full[:n_vars]
         zeta = x_full[n_vars:]
@@ -196,9 +202,9 @@ def main():
     print(f"CVaR (95%): {cvar_actual*100:.4f}%")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 65)
 
     print()
diff --git a/benchmarks/LP_power_dispatch_120d/README.md b/benchmarks/LP_power_dispatch_120d/README.md
index c2be8ad..113a5a7 100644
--- a/benchmarks/LP_power_dispatch_120d/README.md
+++ b/benchmarks/LP_power_dispatch_120d/README.md
@@ -22,7 +22,7 @@ Each scenario delta_i in R^48 encodes uncertain hourly wind generation (delta[0:
 
 ![Power Dispatch Results](results/power_dispatch.png)
 
-The LP solves with N = 150 scenarios. Total cost is $831,357 comprising $537,670 in generation cost and $293,687 in power imbalance penalty. Coal (cheapest at $30/MWh) provides baseload at 284-500 MW. Gas 1 ($40/MWh) ramps to cover evening peak, reaching 300 MW at hour 18. Gas 2 ($50/MWh) sits at its 75 MW minimum throughout. The LP had complexity k = 34, no degeneracy, and risk bounds [0.079, 0.443] at 99.9999% confidence (beta = 10^{-6}).
+The LP solves with N = 150 scenarios. Total cost is $831,357 comprising $537,670 in generation cost and $293,687 in power imbalance penalty. Coal (cheapest at $30/MWh) provides baseload at 284-500 MW. Gas 1 ($40/MWh) ramps to cover evening peak, reaching 300 MW at hour 18. Gas 2 ($50/MWh) sits at its 75 MW minimum throughout. The LP had complexity k = 34, no degeneracy, and risk bounds [0.079, 0.450] at 99.9999% confidence (beta = 10^{-6}).
 
 ## Files
 
@@ -34,16 +34,12 @@ LP_power_dispatch_120d/
 ├── run.py              Solve the LP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Augmented scenario constraint matrix (48 x 120)
-│   ├── b_d.csv         Scenario-dependent RHS vector (48 x 1)
-│   ├── c.csv           Augmented cost vector (120 x 1, includes rho)
-│   ├── G.csv           Augmented hard constraint matrix (330 x 120)
-│   ├── h.csv           Hard constraint RHS (330 x 1)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── demand.csv      Hourly demand profile (24 rows)
 │   ├── generators.txt  Generator parameters
-│   └── scenarios.csv   150 x 48 renewable generation scenarios
+│   ├── scenarios.csv   150 x 48 renewable generation scenarios
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (cost, dispatch, risk bounds)
     ├── solution.csv    Raw solution vector
@@ -57,8 +53,9 @@ LP_power_dispatch_120d/
 # Regenerate scenarios (optional)
 python generate.py --n_scenarios 150 --seed 42
 
-# Solve the LP
+# Solve the LP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -70,25 +67,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **LP** tab, formulation: **Robust**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
 3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 831357.4274826923
+- Optimal cost: 831357.4274827515
 - Complexity k: 34
-- Risk bounds: [0.07906891773454847, 0.44250483819283526]
+- Risk bounds: [0.0790688063794126, 0.44961507555718216]
diff --git a/benchmarks/LP_power_dispatch_120d/parameters.txt b/benchmarks/LP_power_dispatch_120d/parameters.txt
index dc72be3..be6ca54 100644
--- a/benchmarks/LP_power_dispatch_120d/parameters.txt
+++ b/benchmarks/LP_power_dispatch_120d/parameters.txt
@@ -18,4 +18,5 @@
 # The penalty coefficient rho=100 is embedded in c_aug, not passed to the solver
 rho = 0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/LP_power_dispatch_120d/plot.py b/benchmarks/LP_power_dispatch_120d/plot.py
index b8eeaa9..dcd630a 100644
--- a/benchmarks/LP_power_dispatch_120d/plot.py
+++ b/benchmarks/LP_power_dispatch_120d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the power dispatch benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/LP_power_dispatch_120d/run.py b/benchmarks/LP_power_dispatch_120d/run.py
index 3222d73..af5db37 100644
--- a/benchmarks/LP_power_dispatch_120d/run.py
+++ b/benchmarks/LP_power_dispatch_120d/run.py
@@ -3,12 +3,13 @@
 Solve the Economic Dispatch benchmark using the scenario approach LP solver.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 # Add project root to path for imports
@@ -88,6 +89,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: Economic Dispatch with Renewable Uncertainty (LP)")
     print("=" * 65)
@@ -136,7 +142,7 @@ def main():
     # Data files contain the augmented formulation x_aug = [P; zeta] with
     # rho embedded in c, A_d augmented with -I slack columns, and G/h
     # extended with zeta >= 0 constraints. See generate.py for details.
-    print("Solving LP with MOSEK...")
+    print(f"Solving LP with {args.solver}...")
 
     try:
         x_full, _, cost, N_out, k, constraints, degeneracy = solve_lp(
@@ -150,7 +156,7 @@ def main():
             x_ref=np.zeros((n_aug, 1)),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         x = x_full[:n_vars]
         zeta = x_full[n_vars:]
@@ -190,9 +196,9 @@ def main():
     print(f"  - Imbalance penalty: ${slack_cost:.2f}")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 65)
 
     print()
diff --git a/benchmarks/QP_Iris_minimal_3d/README.md b/benchmarks/QP_Iris_minimal_3d/README.md
index 32d4079..738563c 100644
--- a/benchmarks/QP_Iris_minimal_3d/README.md
+++ b/benchmarks/QP_Iris_minimal_3d/README.md
@@ -16,7 +16,7 @@ A(δᵢ) = [-δᵢ(1)·δᵢ(3), -δᵢ(2)·δᵢ(3), -δᵢ(3)], b(δᵢ) = 1
 
 ![Iris SVM Classification](results/iris_svm.png)
 
-The QP solves with N = 150 data points achieving 100% training accuracy. The optimal classifier is 3.79 = 1.29x₁ + 0.82x₂ with margin 1.304. The complexity is k = 2 (two support vectors). With β = 10⁻⁶, the risk bounds are ε̲ = 0 and ε̄ = 0.143.
+The QP solves with N = 150 data points achieving 100% training accuracy. The optimal classifier is 3.79 = 1.29x₁ + 0.82x₂ with margin 1.304. The complexity is k = 2 (two support vectors). With β = 10⁻⁶, the risk bounds are ε̲ = 0 and ε̄ = 0.144.
 
 ## Files
 
@@ -26,26 +26,27 @@ QP_Iris_minimal_3d/
 ├── parameters.txt      Solver parameters (rho, tau, confidence)
 ├── run.py              Solve the QP and print results
 ├── plot.py             Generate the paper figure
+├── plot_data_only.py   Generate a figure of the data only (no classifier)
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Scenario-dependent constraint (1 × 3 expression)
-│   ├── b_d.csv         Scenario-dependent RHS (constant = 1)
-│   ├── c.csv           Linear objective vector (3 × 1)
-│   ├── Q.csv           Quadratic objective matrix (3 × 3)
-│   └── scenarios.csv   150 samples [petal_length, petal_width, label]
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv   150 samples [petal_length, petal_width, label]
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (hyperplane, margin, risk bounds)
     ├── solution.csv    Raw solution vector [w1, w2, b]
     ├── iris_svm.png    Paper figure (300 dpi)
-    └── iris_svm.pdf    Paper figure (vector)
+    ├── iris_svm.pdf    Paper figure (vector)
+    ├── iris_data_only.png  Data-only figure (300 dpi)
+    └── iris_data_only.pdf  Data-only figure (vector)
 ```
 
 ## Usage
 
 ```bash
-# Solve the QP
+# Solve the QP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -57,24 +58,24 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **QP** tab, formulation: **Robust**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **Q**: `data/Q.csv`
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **Q**: `Q`
 3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 1.176470626281622
+- Optimal cost: 1.1764706262816014
 - Complexity k: 2
-- Risk bounds: [0.0, 0.14300981911209715]
+- Risk bounds: [0.0, 0.14409952527843412]
diff --git a/benchmarks/QP_Iris_minimal_3d/parameters.txt b/benchmarks/QP_Iris_minimal_3d/parameters.txt
index 46f9db5..0cd6527 100644
--- a/benchmarks/QP_Iris_minimal_3d/parameters.txt
+++ b/benchmarks/QP_Iris_minimal_3d/parameters.txt
@@ -1,3 +1,4 @@
 rho = 0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/QP_Iris_minimal_3d/plot.py b/benchmarks/QP_Iris_minimal_3d/plot.py
index 8e95afe..15351b8 100644
--- a/benchmarks/QP_Iris_minimal_3d/plot.py
+++ b/benchmarks/QP_Iris_minimal_3d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the Iris SVM classification benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/QP_Iris_minimal_3d/plot_data_only.py b/benchmarks/QP_Iris_minimal_3d/plot_data_only.py
index ee1bbf7..9b7a810 100644
--- a/benchmarks/QP_Iris_minimal_3d/plot_data_only.py
+++ b/benchmarks/QP_Iris_minimal_3d/plot_data_only.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a paper-ready figure showing only the Iris data (no solution overlay)."""
 
 import os
diff --git a/benchmarks/QP_Iris_minimal_3d/run.py b/benchmarks/QP_Iris_minimal_3d/run.py
index 2de6ec1..afeeaa4 100644
--- a/benchmarks/QP_Iris_minimal_3d/run.py
+++ b/benchmarks/QP_Iris_minimal_3d/run.py
@@ -3,12 +3,13 @@
 Solve the Iris SVM Classification benchmark using the scenario approach QP solver.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 # Add project root to path for imports
@@ -57,6 +58,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: Iris SVM Classification (QP)")
     print("=" * 65)
@@ -107,7 +113,7 @@ def b_d(delta):
     print()
 
     # Solve QP
-    print("Solving QP with MOSEK...")
+    print(f"Solving QP with {args.solver}...")
     try:
         x, zeta, cost, N, k, constraints, degeneracy = solve_qp(
             deltas=scenarios,
@@ -121,7 +127,7 @@ def b_d(delta):
             x_ref=np.zeros((n_vars, 1)),
             rho=rho,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
         print(f"Solved successfully")
diff --git a/benchmarks/QP_cbf_robot_navigation_118d/README.md b/benchmarks/QP_cbf_robot_navigation_118d/README.md
index 8c0dab8..bac6902 100644
--- a/benchmarks/QP_cbf_robot_navigation_118d/README.md
+++ b/benchmarks/QP_cbf_robot_navigation_118d/README.md
@@ -24,7 +24,7 @@ This is a Quadratic Program with 118 decision variables: 80 states (position and
 
 ![QP Robot Navigation](results/cbf_navigation.png)
 
-The QP solves with N = 500 scenarios. The optimal trajectory has path length 9.79 m over 8.0 s, arcing smoothly above the wall with a minimum clearance of 0.19 m from the safety boundary. The complexity is k = 1 (one support constraint). With beta = 10^{-6}, the risk bounds are eps_lower = 0 and eps_upper = 0.0400.
+The QP solves with N = 500 scenarios. The optimal trajectory has path length 9.79 m over 8.0 s, arcing smoothly above the wall with a minimum clearance of 0.19 m from the safety boundary. The complexity is k = 1 (one support constraint). With beta = 10^{-6}, the risk bounds are eps_lower = 0 and eps_upper = 0.0403.
 
 ## Files
 
@@ -35,21 +35,18 @@ QP_cbf_robot_navigation_118d/
 ├── generate.py         Generate benchmark data (scenarios + QP matrices)
 ├── run.py              Solve the QP and print results
 ├── plot.py             Generate the paper figure
+├── animate.py          Generate the trajectory animation GIF
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
-│   ├── A_d.csv         Scenario-dependent constraint matrix (3 x 118)
-│   ├── b_d.csv         Scenario-dependent RHS (3 expressions with delta)
-│   ├── c.csv           Linear objective vector (118 x 1)
-│   ├── Q.csv           Quadratic objective matrix (118 x 118)
-│   ├── G.csv           Hard constraint matrix (320 x 118)
-│   ├── h.csv           Hard constraint RHS (320 x 1)
-│   └── scenarios.csv   500 samples (1D wall face perturbations)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv   500 samples (1D wall face perturbations)
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (trajectory, clearances, risk bounds)
     ├── solution.csv    Raw solution vector (118 x 1)
     ├── cbf_navigation.png   Paper figure (300 dpi)
-    └── cbf_navigation.pdf   Paper figure (vector)
+    ├── cbf_navigation.pdf   Paper figure (vector)
+    └── trajectory_animation.gif  Trajectory animation
 ```
 
 ## Usage
@@ -58,8 +55,9 @@ QP_cbf_robot_navigation_118d/
 # Generate benchmark data (optional, data already provided)
 python generate.py [--n_scenarios 500] [--sigma 0.05]
 
-# Solve the QP
+# Solve the QP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -71,26 +69,26 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **QP** tab, formulation: **Robust**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
-   - **Q**: `data/Q.csv`
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
+   - **Q**: `Q`
 3. Set parameters: rho = 0, tau = 0, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: -13777.807262403046
+- Optimal cost: -13777.807262404083
 - Complexity k: 1
-- Risk bounds: [0.0, 0.04000231437711046]
+- Risk bounds: [0.0, 0.04027646044222639]
diff --git a/benchmarks/QP_cbf_robot_navigation_118d/animate.py b/benchmarks/QP_cbf_robot_navigation_118d/animate.py
index 4e08dc4..cb555fe 100644
--- a/benchmarks/QP_cbf_robot_navigation_118d/animate.py
+++ b/benchmarks/QP_cbf_robot_navigation_118d/animate.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a trajectory animation GIF for the robot navigation benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/QP_cbf_robot_navigation_118d/parameters.txt b/benchmarks/QP_cbf_robot_navigation_118d/parameters.txt
index b66c637..424ddc0 100644
--- a/benchmarks/QP_cbf_robot_navigation_118d/parameters.txt
+++ b/benchmarks/QP_cbf_robot_navigation_118d/parameters.txt
@@ -11,4 +11,5 @@
 
 rho = 0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/QP_cbf_robot_navigation_118d/plot.py b/benchmarks/QP_cbf_robot_navigation_118d/plot.py
index afd58f1..9b85421 100644
--- a/benchmarks/QP_cbf_robot_navigation_118d/plot.py
+++ b/benchmarks/QP_cbf_robot_navigation_118d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the robot navigation benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/QP_cbf_robot_navigation_118d/run.py b/benchmarks/QP_cbf_robot_navigation_118d/run.py
index 6829b61..98ac99b 100644
--- a/benchmarks/QP_cbf_robot_navigation_118d/run.py
+++ b/benchmarks/QP_cbf_robot_navigation_118d/run.py
@@ -3,12 +3,13 @@
 Solve the CBF Robot Navigation benchmark using the scenario approach QP solver.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 # Add project root to path for imports
@@ -64,6 +65,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 65)
     print("BENCHMARK: QP Robot Navigation")
     print("=" * 65)
@@ -103,7 +109,7 @@ def main():
     print()
 
     # Solve QP
-    print("Solving QP with MOSEK...")
+    print(f"Solving QP with {args.solver}...")
     try:
         x, zeta, cost, N, k, constraints, degeneracy = solve_qp(
             deltas=scenarios,
@@ -117,12 +123,12 @@ def main():
             x_ref=np.zeros((n_vars, 1)),
             rho=rho,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
         print(f"Solved successfully")
     except Exception as e:
-        print(f"MOSEK failed: {e}")
+        print(f"{args.solver} failed: {e}")
         print("Trying CLARABEL...")
         try:
             x, zeta, cost, N, k, constraints, degeneracy = solve_qp(
diff --git a/benchmarks/QP_radiation_therapy_50d/README.md b/benchmarks/QP_radiation_therapy_50d/README.md
index 9157910..0b1f1b5 100644
--- a/benchmarks/QP_radiation_therapy_50d/README.md
+++ b/benchmarks/QP_radiation_therapy_50d/README.md
@@ -30,7 +30,7 @@ The optimal plan achieves a tumour D95 = 145.12 Gy (exceeding the 137.75 Gy mini
 
 Panel (a) shows the dose-volume histogram (DVH). The tumour curve drops steeply beyond the prescription dose, confirming adequate target coverage, while the long tail reflects the hot spots inherent to brachytherapy (high doses near dwell positions are physically unavoidable and clinically acceptable). Both OAR curves fall well within their respective dose constraints.
 
-Panel (b) shows the scenario approach risk bounds. The complexity is k = 6 with no degeneracy, yielding risk bounds [0.000, 0.149] at 99.9999% confidence (beta = 10^{-6}). This guarantees that the dose constraints will be satisfied for a new random catheter placement with probability at least 85.1%. The low complexity relative to d = 50 decision variables indicates that only a small number of catheter-shift scenarios are critical to the solution geometry.
+Panel (b) shows the scenario approach risk bounds. The complexity is k = 5 and degeneracy was detected, so no lower bound on risk can be certified; the upper risk bound is 0.141 at 99.9999% confidence (beta = 10^{-6}). This guarantees that the dose constraints will be satisfied for a new random catheter placement with probability at least 85.9%. The low complexity relative to d = 50 decision variables indicates that only a small number of catheter-shift scenarios are critical to the solution geometry.
 
 Panel (c) displays the dwell-position intensity profile. Most channels carry high intensity (~150 Gy.s, the upper bound), with a few positions at low or zero intensity. This reflects the geometry of the prostate target relative to the catheter array: dwell positions near the centre of the PTV deliver maximum intensity, while those far from the target or near OARs are suppressed.
 
@@ -44,17 +44,13 @@ QP_radiation_therapy_50d/
 ├── run.py              Solve the QP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json  One-shot program definition (JSON)
-│   ├── program_symbolic.json   One-shot program definition (MATLAB)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── anatomy.txt     Structure definitions (voxel counts, dose limits)
 │   ├── scenarios.csv   200 x 3 catheter shift scenarios
-│   ├── A_d.csv         Affine constraint coefficients (80 x 50, with delta expressions)
-│   ├── b_d.csv         Constraint RHS values (80 x 1, with delta expressions)
-│   ├── Q.csv           Quadratic objective matrix (50 x 50)
-│   ├── c.csv           Linear objective vector (50 x 1)
-│   ├── G.csv           Hard constraint matrix (100 x 50)
-│   ├── h.csv           Hard constraint RHS (100 x 1)
-│   └── D_nominal.csv   Nominal dose-influence matrix (100 x 50)
+│   ├── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
+│   ├── D_nominal.csv   Nominal dose-influence matrix (100 x 50)
+│   └── Prostate_BT_01.mat  TROTS source data used by generate.py
 └── results/
     ├── metrics.json                              Solver output (cost, risk bounds, doses)
     ├── solution.csv                              Raw solution vector
@@ -68,8 +64,9 @@ QP_radiation_therapy_50d/
 # Generate scenarios and constraint matrices from TROTS data (requires mat73)
 python generate.py
 
-# Solve the QP
+# Solve the QP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -81,26 +78,26 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
-1. Select the **QP** tab, formulation: **Robust + Regularization**
-2. Upload or enter each matrix:
-   - **A(delta)**: `data/A_d.csv`
-   - **b(delta)**: `data/b_d.csv`
-   - **c**: `data/c.csv`
-   - **G**: `data/G.csv`
-   - **h**: `data/h.csv`
-   - **Q**: `data/Q.csv`
+1. Select the **QP** tab, formulation: **Regularization**
+2. Enter each matrix with its **Edit** button, using the matching field of `data/program_symbolic.json`:
+   - **A(delta)**: `A_d`
+   - **b(delta)**: `b_d`
+   - **c**: `c`
+   - **G**: `G`
+   - **h**: `h`
+   - **Q**: `Q`
 3. Set parameters: rho = 0, tau = 0.1, confidence (beta) = 1e-06
 4. Upload `data/scenarios.csv` in the Scenarios box
 5. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 4698.836116109508
-- Complexity k: 6
-- Risk bounds: [0.0, 0.14861067387333604]
+- Optimal cost: 4698.836115903002 (objective including tau ||x||_2, as in `results/metrics.json`; the web interface reports the cost without this term, about 4626.85)
+- Complexity k: 5
+- Risk bounds: [0.0, 0.1413759133502026] (degeneracy detected, so the lower bound is not certified)
diff --git a/benchmarks/QP_radiation_therapy_50d/generate.py b/benchmarks/QP_radiation_therapy_50d/generate.py
index a0d9db1..9855a52 100644
--- a/benchmarks/QP_radiation_therapy_50d/generate.py
+++ b/benchmarks/QP_radiation_therapy_50d/generate.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires mat73 (not in requirements.txt): pip install mat73
 """
 Prostate Brachytherapy Benchmark Data Generator (TROTS Dataset)
 
diff --git a/benchmarks/QP_radiation_therapy_50d/parameters.txt b/benchmarks/QP_radiation_therapy_50d/parameters.txt
index cd4f041..4159def 100644
--- a/benchmarks/QP_radiation_therapy_50d/parameters.txt
+++ b/benchmarks/QP_radiation_therapy_50d/parameters.txt
@@ -12,4 +12,5 @@
 #
 rho = 0
 tau = 0.1
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/QP_radiation_therapy_50d/plot.py b/benchmarks/QP_radiation_therapy_50d/plot.py
index 8b1378f..9942aaf 100644
--- a/benchmarks/QP_radiation_therapy_50d/plot.py
+++ b/benchmarks/QP_radiation_therapy_50d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the prostate brachytherapy benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/QP_radiation_therapy_50d/run.py b/benchmarks/QP_radiation_therapy_50d/run.py
index 6dad747..6aca4f9 100644
--- a/benchmarks/QP_radiation_therapy_50d/run.py
+++ b/benchmarks/QP_radiation_therapy_50d/run.py
@@ -9,12 +9,13 @@
 Clinical Context: Prostate IMRT under patient setup uncertainty (±5mm in 3D).
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import re
 import numpy as np
 import pandas as pd
@@ -157,6 +158,11 @@ def load_anatomy(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 70)
     print("BENCHMARK: Prostate IMRT Radiation Therapy (QP)")
     print("Robust Treatment Under Patient Setup Uncertainty")
@@ -212,7 +218,7 @@ def main():
     print()
 
     # ===== SOLVE USING solve_qp (robust, ρ=0, τ>0) =====
-    print("Solving QP with MOSEK...")
+    print(f"Solving QP with {args.solver}...")
 
     try:
         x, zeta, cost, N_out, k, constraints, degeneracy = solve_qp(
@@ -227,7 +233,7 @@ def main():
             x_ref=np.zeros((n_vars, 1)),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
         print(f"Solved successfully")
@@ -261,9 +267,9 @@ def main():
     print(f"Objective Value: {cost:.4f}")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 70)
 
     print()
diff --git a/benchmarks/SDP_LPV_stability_3d/README.md b/benchmarks/SDP_LPV_stability_3d/README.md
index 8814d4f..66fa70d 100644
--- a/benchmarks/SDP_LPV_stability_3d/README.md
+++ b/benchmarks/SDP_LPV_stability_3d/README.md
@@ -25,33 +25,33 @@ where F_j encode the Lyapunov equation A(delta)'P + PA(delta). A hard constraint
 
 The SDP solves with N = 100 scenarios. The optimal Lyapunov matrix is:
 
-    P = [[10.68, 1.31],
-         [1.31,  9.00]]
+    P = [[10.67, 1.30],
+         [1.30,  9.01]]
 
-The solver found complexity k = 1, meaning one scenario constraint is active at the optimal solution. This occurs at the upper boundary of the parameter range (delta near 1.0), where the system is closest to instability. With 99.9999% confidence (beta = 10^{-6}), the risk bounds are [0.000, 0.185], meaning at most 18.5% of parameter values could violate the certified stability condition. There is no degeneracy.
+The solver found complexity k = 1, meaning one scenario constraint is active at the optimal solution. This occurs at the upper boundary of the parameter range (delta near 1.0), where the system is closest to instability. With 99.9999% confidence (beta = 10^{-6}), the risk bounds are [0.000, 0.186], meaning at most 18.6% of parameter values could violate the certified stability condition. There is no degeneracy.
 
 ## Files
 
 ```
-LPV_stability_3d/
+SDP_LPV_stability_3d/
 ├── README.md           This file
 ├── parameters.txt      Solver parameters (rho, tau, confidence)
 ├── generate.py         Generate scheduling parameter scenarios
 ├── run.py              Solve the SDP and print results
 ├── plot.py             Generate the paper figure
+├── plot_data_only.py   Generate a figure of the trajectories only (no Lyapunov ellipse)
 ├── data/
-│   ├── program_symbolic.json         One-shot program definition (JSON)
-│   ├── program_symbolic.json          One-shot program definition (MATLAB)
-│   ├── F_0.csv ... F_3.csv   Scenario-dependent LMI matrices
-│   ├── E_0.csv ... E_3.csv   Hard constraint (P > 0) matrices
-│   ├── c.csv                  Linear objective vector
-│   ├── Q.csv                  Quadratic objective matrix
-│   └── scenarios.csv          100 scheduling parameter samples
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv          100 scheduling parameter samples
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (cost, risk bounds, etc.)
     ├── solution.csv    Raw solution vector
     ├── lpv_stability.png   Paper figure (300 dpi)
-    └── lpv_stability.pdf   Paper figure (vector)
+    ├── lpv_stability.pdf   Paper figure (vector)
+    ├── lpv_trajectories_only.png   Trajectories-only figure (300 dpi)
+    └── lpv_trajectories_only.pdf   Trajectories-only figure (vector)
 ```
 
 ## Usage
@@ -60,8 +60,9 @@ LPV_stability_3d/
 # Generate scenarios (optional)
 python generate.py
 
-# Solve the SDP
+# Solve the SDP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -73,32 +74,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
-1. Select the **SDP** tab, formulation: **Robust + Relaxation**
-2. Upload matrices using one of these approaches:
-
-   **Option A — One-shot LMI upload:**
-   - Click **Edit F(delta)** -> upload `data/program_symbolic.json` (the F_d section) as a JSON file, or enter the F_d dict
-   - Click **Edit E** -> upload the E matrices similarly
-
-   **Option B — Individual matrix entry:**
-   - Click **Edit F(delta)** -> set n = 3 -> click each F_i button and upload `data/F_0.csv` through `data/F_3.csv`
-   - Click **Edit E** -> set n = 3 -> click each E_i button and upload `data/E_0.csv` through `data/E_3.csv`
-
-3. Upload or enter:
-   - **c**: `data/c.csv`
-   - **Q**: `data/Q.csv`
+1. Select the **SDP** tab, formulation: **Relaxation**
+2. Enter the matrices with their **Edit** buttons, using the matching fields of `data/program_symbolic.json`:
+   - **F(delta)**: `F_d` (matrices F_0 ... F_3)
+   - **E**: `E` (matrices E_0 ... E_3)
+3. Enter:
+   - **c**: `c`
+   - **Q**: `Q`
 4. Set parameters: rho = 1.0, tau = 0, confidence (beta) = 1e-06
 5. Upload `data/scenarios.csv` in the Scenarios box
 6. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: -9.843942363041855
+- Optimal cost: -9.843942363753102
 - Complexity k: 1
-- Risk bounds: [0.0, 0.185321964349132]
+- Risk bounds: [0.0, 0.18583116631314622]
diff --git a/benchmarks/SDP_LPV_stability_3d/parameters.txt b/benchmarks/SDP_LPV_stability_3d/parameters.txt
index f429cb6..1a960b8 100644
--- a/benchmarks/SDP_LPV_stability_3d/parameters.txt
+++ b/benchmarks/SDP_LPV_stability_3d/parameters.txt
@@ -10,4 +10,5 @@
 
 rho = 1.0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/SDP_LPV_stability_3d/plot.py b/benchmarks/SDP_LPV_stability_3d/plot.py
index c9a52a7..1346e9d 100644
--- a/benchmarks/SDP_LPV_stability_3d/plot.py
+++ b/benchmarks/SDP_LPV_stability_3d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the LPV Stability 3D benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/SDP_LPV_stability_3d/plot_data_only.py b/benchmarks/SDP_LPV_stability_3d/plot_data_only.py
index 1ebe2f3..524e4a0 100644
--- a/benchmarks/SDP_LPV_stability_3d/plot_data_only.py
+++ b/benchmarks/SDP_LPV_stability_3d/plot_data_only.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a paper-ready figure showing only the LPV trajectories (no Lyapunov ellipse)."""
 
 import os
diff --git a/benchmarks/SDP_LPV_stability_3d/run.py b/benchmarks/SDP_LPV_stability_3d/run.py
index 6f89325..8a0212e 100644
--- a/benchmarks/SDP_LPV_stability_3d/run.py
+++ b/benchmarks/SDP_LPV_stability_3d/run.py
@@ -6,12 +6,13 @@
 across all sampled parameter values delta in [-0.22, 1].
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..'))
@@ -66,6 +67,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 60)
     print("BENCHMARK: LPV_stability_3d (SDP)")
     print("Common Lyapunov Function for LPV System")
@@ -103,18 +109,18 @@ def main():
     print()
 
     # Solve
-    print("Solving SDP with MOSEK...")
+    print(f"Solving SDP with {args.solver}...")
     rho_value = 1.0
 
     try:
         x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
             deltas=scenarios, F_d=F_d, E=E, c=c, Q=Q,
             tau=0.0, x_ref=np.zeros(n_vars), rho=rho_value,
-            norm_type=2, solver='MOSEK'
+            norm_type=2, solver=args.solver
         )
-        print(f"Solved with MOSEK")
+        print(f"Solved with {args.solver}")
     except Exception as e:
-        print(f"MOSEK failed: {e}")
+        print(f"{args.solver} failed: {e}")
         print("Attempting with SCS solver...")
         try:
             x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -137,9 +143,9 @@ def main():
     print(f"Optimal Cost: {cost:.6f}")
     print(f"Complexity (k): {k}")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 60)
 
     print()
diff --git a/benchmarks/SDP_LPV_stability_3d_numeric/data/program_numeric.json b/benchmarks/SDP_LPV_stability_3d_numeric/data/program_numeric.json
deleted file mode 100644
index ac1f091..0000000
--- a/benchmarks/SDP_LPV_stability_3d_numeric/data/program_numeric.json
+++ /dev/null
@@ -1,80 +0,0 @@
-{
-  "type": "SDP",
-  "mode": "numeric",
-  "c": [
-    [
-      "-1"
-    ],
-    [
-      "0"
-    ],
-    [
-      "-1"
-    ]
-  ],
-  "Q": [
-    [
-      "0.1",
-      "0",
-      "0"
-    ],
-    [
-      "0",
-      "0.1",
-      "0"
-    ],
-    [
-      "0",
-      "0",
-      "0.1"
-    ]
-  ],
-  "E": {
-    "0": [
-      [
-        "0",
-        "0"
-      ],
-      [
-        "0",
-        "0"
-      ]
-    ],
-    "1": [
-      [
-        "-1",
-        "0"
-      ],
-      [
-        "0",
-        "0"
-      ]
-    ],
-    "2": [
-      [
-        "0",
-        "-1"
-      ],
-      [
-        "-1",
-        "0"
-      ]
-    ],
-    "3": [
-      [
-        "0",
-        "0"
-      ],
-      [
-        "0",
-        "-1"
-      ]
-    ]
-  },
-  "n_x": 3,
-  "lmi_size": 2,
-  "rho": "1.0",
-  "tau": "0",
-  "confidence": "1e-06",
-  "lmi_e_size": 2
-}
\ No newline at end of file
diff --git a/benchmarks/SDP_LPV_stability_3d_numeric/data/program_symbolic.json b/benchmarks/SDP_LPV_stability_3d_numeric/data/program_symbolic.json
deleted file mode 100644
index 178be87..0000000
--- a/benchmarks/SDP_LPV_stability_3d_numeric/data/program_symbolic.json
+++ /dev/null
@@ -1,80 +0,0 @@
-{
-  "type": "SDP",
-  "mode": "symbolic",
-  "c": [
-    [
-      "-1"
-    ],
-    [
-      "0"
-    ],
-    [
-      "-1"
-    ]
-  ],
-  "Q": [
-    [
-      "0.1",
-      "0",
-      "0"
-    ],
-    [
-      "0",
-      "0.1",
-      "0"
-    ],
-    [
-      "0",
-      "0",
-      "0.1"
-    ]
-  ],
-  "E": {
-    "0": [
-      [
-        "0",
-        "0"
-      ],
-      [
-        "0",
-        "0"
-      ]
-    ],
-    "1": [
-      [
-        "-1",
-        "0"
-      ],
-      [
-        "0",
-        "0"
-      ]
-    ],
-    "2": [
-      [
-        "0",
-        "-1"
-      ],
-      [
-        "-1",
-        "0"
-      ]
-    ],
-    "3": [
-      [
-        "0",
-        "0"
-      ],
-      [
-        "0",
-        "-1"
-      ]
-    ]
-  },
-  "n_x": 3,
-  "lmi_size": 2,
-  "rho": "1.0",
-  "tau": "0",
-  "confidence": "1e-06",
-  "lmi_e_size": 2
-}
\ No newline at end of file
diff --git a/benchmarks/SDP_LPV_stability_3d_numeric/data/scenarios.csv b/benchmarks/SDP_LPV_stability_3d_numeric/data/scenarios.csv
deleted file mode 100644
index 6520317..0000000
--- a/benchmarks/SDP_LPV_stability_3d_numeric/data/scenarios.csv
+++ /dev/null
@@ -1,100 +0,0 @@
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.112119463901388272e+00,-1.018686577316898045e+00,-1.018686577316898045e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.056059731950694136e+00,-2.056059731950694136e+00,-2.037373154633796091e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.101785268712619725e+00,-1.016964211452103362e+00,-1.016964211452103362e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.050892634356309863e+00,-2.050892634356309863e+00,-2.033928422904206723e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.099266417871516133e+00,-1.016544402978586170e+00,-1.016544402978586170e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.049633208935758066e+00,-2.049633208935758066e+00,-2.033088805957172340e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.093053472563251916e+00,-1.015508912093875393e+00,-1.015508912093875393e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.046526736281625958e+00,-2.046526736281625958e+00,-2.031017824187750787e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.091382360719866718e+00,-1.015230393453311120e+00,-1.015230393453311120e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.045691180359933359e+00,-2.045691180359933359e+00,-2.030460786906622239e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.082877396977558249e+00,-1.013812899496259634e+00,-1.013812899496259634e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.041438698488779124e+00,-2.041438698488779124e+00,-2.027625798992519268e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.077805892378880692e+00,-1.012967648729813375e+00,-1.012967648729813375e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.038902946189440346e+00,-2.038902946189440346e+00,-2.025935297459626749e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.070825213281538346e+00,-1.011804202213589798e+00,-1.011804202213589798e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.035412606640769173e+00,-2.035412606640769173e+00,-2.023608404427179597e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.066600794579191103e+00,-1.011100132429865184e+00,-1.011100132429865184e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.033300397289595551e+00,-2.033300397289595551e+00,-2.022200264859730368e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.065406924413758283e+00,-1.010901154068959640e+00,-1.010901154068959640e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.032703462206879141e+00,-2.032703462206879141e+00,-2.021802308137919280e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.060168364526724361e+00,-1.010028060754454060e+00,-1.010028060754454060e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.030084182263362180e+00,-2.030084182263362180e+00,-2.020056121508908120e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.059936635518561765e+00,-1.009989439253093702e+00,-1.009989439253093702e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.029968317759280882e+00,-2.029968317759280882e+00,-2.019978878506187403e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.044691298553543923e+00,-1.007448549758924061e+00,-1.007448549758924061e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.022345649276771962e+00,-2.022345649276771962e+00,-2.014897099517848122e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.038227421969654785e+00,-1.006371236994942464e+00,-1.006371236994942464e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.019113710984827392e+00,-2.019113710984827392e+00,-2.012742473989884928e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.028796208339347373e+00,-1.004799368056557896e+00,-1.004799368056557896e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.014398104169673687e+00,-2.014398104169673687e+00,-2.009598736113115791e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.021557424694345073e+00,-1.003592904115724105e+00,-1.003592904115724105e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.010778712347172537e+00,-2.010778712347172537e+00,-2.007185808231448210e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.016981874661724738e+00,-1.002830312443620864e+00,-1.002830312443620864e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.008490937330862369e+00,-2.008490937330862369e+00,-2.005660624887241728e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.009784483346217954e+00,-1.001630747224369733e+00,-1.001630747224369733e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.004892241673108977e+00,-2.004892241673108977e+00,-2.003261494448739466e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.001643885532265799e+00,-1.000273980922044448e+00,-1.000273980922044448e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.000821942766132899e+00,-2.000821942766132899e+00,-2.000547961844088896e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-4.000301897762891024e+00,-1.000050316293815245e+00,-1.000050316293815245e+00,2.000000000000000000e+00,0.000000000000000000e+00,-2.000150948881445512e+00,-2.000150948881445512e+00,-2.000100632587630489e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.999696176606913234e+00,-9.999493627678188723e-01,-9.999493627678188723e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.999848088303456617e+00,-1.999848088303456617e+00,-1.999898725535637745e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.990284873555570933e+00,-9.983808122592618961e-01,-9.983808122592618961e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.995142436777785466e+00,-1.995142436777785466e+00,-1.996761624518523792e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.988665699950435073e+00,-9.981109499917392158e-01,-9.981109499917392158e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.994332849975217536e+00,-1.994332849975217536e+00,-1.996221899983478432e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.981284718068835549e+00,-9.968807863448059248e-01,-9.968807863448059248e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.990642359034417774e+00,-1.990642359034417774e+00,-1.993761572689611850e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.959597372777234003e+00,-9.932662287962057412e-01,-9.932662287962057412e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.979798686388617002e+00,-1.979798686388617002e+00,-1.986532457592411482e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.957530291616739859e+00,-9.929217152694566062e-01,-9.929217152694566062e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.978765145808369930e+00,-1.978765145808369930e+00,-1.985843430538913212e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.947227841127700643e+00,-9.912046401879500701e-01,-9.912046401879500701e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.973613920563850321e+00,-1.973613920563850321e+00,-1.982409280375900140e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.946295188037403534e+00,-9.910491980062339223e-01,-9.910491980062339223e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.973147594018701767e+00,-1.973147594018701767e+00,-1.982098396012467845e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.942238504692537227e+00,-9.903730841154227971e-01,-9.903730841154227971e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.971119252346268613e+00,-1.971119252346268613e+00,-1.980746168230845594e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.941438450081389266e+00,-9.902397416802315444e-01,-9.902397416802315444e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.970719225040694633e+00,-1.970719225040694633e+00,-1.980479483360463089e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.937748798774406200e+00,-9.896247997957343667e-01,-9.896247997957343667e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.968874399387203100e+00,-1.968874399387203100e+00,-1.979249599591468733e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.934987861979514445e+00,-9.891646436632524075e-01,-9.891646436632524075e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.967493930989757223e+00,-1.967493930989757223e+00,-1.978329287326504815e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.932403097319645546e+00,-9.887338495532742577e-01,-9.887338495532742577e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.966201548659822773e+00,-1.966201548659822773e+00,-1.977467699106548515e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.920677188953607395e+00,-9.867795314922679362e-01,-9.867795314922679362e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.960338594476803697e+00,-1.960338594476803697e+00,-1.973559062984535872e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.918386043444387568e+00,-9.863976739073979649e-01,-9.863976739073979649e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.959193021722193784e+00,-1.959193021722193784e+00,-1.972795347814795930e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.911122289839947364e+00,-9.851870483066578199e-01,-9.851870483066578199e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.955561144919973682e+00,-1.955561144919973682e+00,-1.970374096613315640e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.882451814953700175e+00,-9.804086358256167699e-01,-9.804086358256167699e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.941225907476850088e+00,-1.941225907476850088e+00,-1.960817271651233540e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.882235534418983391e+00,-9.803725890698306022e-01,-9.803725890698306022e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.941117767209491696e+00,-1.941117767209491696e+00,-1.960745178139661204e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.880728579646800291e+00,-9.801214299411333819e-01,-9.801214299411333819e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.940364289823400146e+00,-1.940364289823400146e+00,-1.960242859882266764e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.880334135125258488e+00,-9.800556891875430443e-01,-9.800556891875430443e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.940167067562629244e+00,-1.940167067562629244e+00,-1.960111378375086089e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.869351375544245020e+00,-9.782252292573740959e-01,-9.782252292573740959e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.934675687772122510e+00,-1.934675687772122510e+00,-1.956450458514748192e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.864078892605931514e+00,-9.773464821009885117e-01,-9.773464821009885117e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.932039446302965757e+00,-1.932039446302965757e+00,-1.954692964201977023e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.860577962523235929e+00,-9.767629937538726548e-01,-9.767629937538726548e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.930288981261617964e+00,-1.930288981261617964e+00,-1.953525987507745310e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.857147345116424209e+00,-9.761912241860707384e-01,-9.761912241860707384e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.928573672558212104e+00,-1.928573672558212104e+00,-1.952382448372141477e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.845989690128260197e+00,-9.743316150213767735e-01,-9.743316150213767735e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.922994845064130098e+00,-1.922994845064130098e+00,-1.948663230042753547e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.839572139990564459e+00,-9.732620233317607061e-01,-9.732620233317607061e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.919786069995282229e+00,-1.919786069995282229e+00,-1.946524046663521412e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.835355881347024987e+00,-9.725593135578375348e-01,-9.725593135578375348e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.917677940673512493e+00,-1.917677940673512493e+00,-1.945118627115675070e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.830284674414752732e+00,-9.717141124024587517e-01,-9.717141124024587517e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.915142337207376366e+00,-1.915142337207376366e+00,-1.943428224804917503e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.819036109862361528e+00,-9.698393516437269213e-01,-9.698393516437269213e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.909518054931180764e+00,-1.909518054931180764e+00,-1.939678703287453843e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.815919392219989437e+00,-9.693198987033315728e-01,-9.693198987033315728e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.907959696109994718e+00,-1.907959696109994718e+00,-1.938639797406663146e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.814404611277317780e+00,-9.690674352128862967e-01,-9.690674352128862967e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.907202305638658890e+00,-1.907202305638658890e+00,-1.938134870425772593e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.812278326683153473e+00,-9.687130544471922455e-01,-9.687130544471922455e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.906139163341576737e+00,-1.906139163341576737e+00,-1.937426108894384491e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.802728301223997409e+00,-9.671213835373329015e-01,-9.671213835373329015e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.901364150611998705e+00,-1.901364150611998705e+00,-1.934242767074665803e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.801470313725866745e+00,-9.669117189543110502e-01,-9.669117189543110502e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.900735156862933373e+00,-1.900735156862933373e+00,-1.933823437908622100e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.795039876687112734e+00,-9.658399794478520484e-01,-9.658399794478520484e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.897519938343556367e+00,-1.897519938343556367e+00,-1.931679958895704097e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.792746698677678108e+00,-9.654577831129463883e-01,-9.654577831129463883e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.896373349338839054e+00,-1.896373349338839054e+00,-1.930915566225892777e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.792409306119812751e+00,-9.654015510199687178e-01,-9.654015510199687178e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.896204653059906375e+00,-1.896204653059906375e+00,-1.930803102039937436e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.791711274005299170e+00,-9.652852123342164914e-01,-9.652852123342164914e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.895855637002649585e+00,-1.895855637002649585e+00,-1.930570424668432983e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.782796656829930715e+00,-9.637994428049884155e-01,-9.637994428049884155e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.891398328414965357e+00,-1.891398328414965357e+00,-1.927598885609976831e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.752262165531660809e+00,-9.587103609219435052e-01,-9.587103609219435052e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.876131082765830405e+00,-1.876131082765830405e+00,-1.917420721843887010e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.743087247934401063e+00,-9.571812079890668068e-01,-9.571812079890668068e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.871543623967200531e+00,-1.871543623967200531e+00,-1.914362415978133614e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.741616496238377998e+00,-9.569360827063964070e-01,-9.569360827063964070e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.870808248119188999e+00,-1.870808248119188999e+00,-1.913872165412792814e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.725315841746747569e+00,-9.542193069577912246e-01,-9.542193069577912246e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.862657920873373785e+00,-1.862657920873373785e+00,-1.908438613915582449e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.724722574509128137e+00,-9.541204290848547265e-01,-9.541204290848547265e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.862361287254564068e+00,-1.862361287254564068e+00,-1.908240858169709453e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.722102470834789845e+00,-9.536837451391316778e-01,-9.536837451391316778e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.861051235417394922e+00,-1.861051235417394922e+00,-1.907367490278263356e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.700023368117618450e+00,-9.500038946862697786e-01,-9.500038946862697786e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.850011684058809225e+00,-1.850011684058809225e+00,-1.900007789372539557e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.695088112288900106e+00,-9.491813520481500177e-01,-9.491813520481500177e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.847544056144450053e+00,-1.847544056144450053e+00,-1.898362704096300035e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.693244430561250624e+00,-9.488740717602084374e-01,-9.488740717602084374e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.846622215280625312e+00,-1.846622215280625312e+00,-1.897748143520416875e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.657962662677138432e+00,-9.429937771128563684e-01,-9.429937771128563684e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.828981331338569216e+00,-1.828981331338569216e+00,-1.885987554225712737e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.654007641064829226e+00,-9.423346068441381673e-01,-9.423346068441381673e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.827003820532414613e+00,-1.827003820532414613e+00,-1.884669213688276335e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.635756223692909295e+00,-9.392927039488181418e-01,-9.392927039488181418e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.817878111846454647e+00,-1.817878111846454647e+00,-1.878585407897636284e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.631425403061260937e+00,-9.385709005102101932e-01,-9.385709005102101932e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.815712701530630468e+00,-1.815712701530630468e+00,-1.877141801020420386e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.614141268838645171e+00,-9.356902114731074915e-01,-9.356902114731074915e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.807070634419322586e+00,-1.807070634419322586e+00,-1.871380422946214983e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.604331661160563804e+00,-9.340552768600939304e-01,-9.340552768600939304e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.802165830580281902e+00,-1.802165830580281902e+00,-1.868110553720187861e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.599802385648055658e+00,-9.333003976080092023e-01,-9.333003976080092023e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.799901192824027829e+00,-1.799901192824027829e+00,-1.866600795216018405e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.592526933402571565e+00,-9.320878222337619645e-01,-9.320878222337619645e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.796263466701285783e+00,-1.796263466701285783e+00,-1.864175644467523929e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.585246223407735577e+00,-9.308743705679559666e-01,-9.308743705679559666e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.792623111703867789e+00,-1.792623111703867789e+00,-1.861748741135911933e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.574639666690437423e+00,-9.291066111150729778e-01,-9.291066111150729778e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.787319833345218711e+00,-1.787319833345218711e+00,-1.858213222230145956e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.571562283697898188e+00,-9.285937139496497350e-01,-9.285937139496497350e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.785781141848949094e+00,-1.785781141848949094e+00,-1.857187427899299470e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.563377270533726993e+00,-9.272295450889544988e-01,-9.272295450889544988e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.781688635266863496e+00,-1.781688635266863496e+00,-1.854459090177908998e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.556867922180706909e+00,-9.261446536967844478e-01,-9.261446536967844478e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.778433961090353455e+00,-1.778433961090353455e+00,-1.852289307393568896e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.555056047007615572e+00,-9.258426745012693360e-01,-9.258426745012693360e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.777528023503807786e+00,-1.777528023503807786e+00,-1.851685349002538672e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.551527144481729792e+00,-9.252545240802882986e-01,-9.252545240802882986e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.775763572240864896e+00,-1.775763572240864896e+00,-1.850509048160576597e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.540907673409620759e+00,-9.234846122349368303e-01,-9.234846122349368303e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.770453836704810380e+00,-1.770453836704810380e+00,-1.846969224469873661e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.523762859051509100e+00,-9.206271431752515166e-01,-9.206271431752515166e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.761881429525754550e+00,-1.761881429525754550e+00,-1.841254286350503033e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.515674551524414770e+00,-9.192790919207357581e-01,-9.192790919207357581e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.757837275762207385e+00,-1.757837275762207385e+00,-1.838558183841471516e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.509548359518279703e+00,-9.182580599197133209e-01,-9.182580599197133209e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.754774179759139852e+00,-1.754774179759139852e+00,-1.836516119839426642e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.504072691806095818e+00,-9.173454486343493031e-01,-9.173454486343493031e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.752036345903047909e+00,-1.752036345903047909e+00,-1.834690897268698606e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.491707515567259890e+00,-9.152845859278766483e-01,-9.152845859278766483e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.745853757783629945e+00,-1.745853757783629945e+00,-1.830569171855753297e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.477807679631721882e+00,-9.129679466052869063e-01,-9.129679466052869063e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.738903839815860941e+00,-1.738903839815860941e+00,-1.825935893210573813e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.471002427707633231e+00,-9.118337379512722052e-01,-9.118337379512722052e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.735501213853816616e+00,-1.735501213853816616e+00,-1.823667475902544410e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.470065832835994790e+00,-9.116776388059990577e-01,-9.116776388059990577e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.735032916417997395e+00,-1.735032916417997395e+00,-1.823355277611998115e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.467925616591222227e+00,-9.113209360985370378e-01,-9.113209360985370378e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.733962808295611113e+00,-1.733962808295611113e+00,-1.822641872197074076e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.448421792588201917e+00,-9.080702987647003566e-01,-9.080702987647003566e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.724210896294100959e+00,-1.724210896294100959e+00,-1.816140597529400713e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.440997473465970735e+00,-9.068329122443283818e-01,-9.068329122443283818e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.720498736732985368e+00,-1.720498736732985368e+00,-1.813665824488656764e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.439173303098520762e+00,-9.065288838497533863e-01,-9.065288838497533863e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.719586651549260381e+00,-1.719586651549260381e+00,-1.813057767699506773e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.436833241632538893e+00,-9.061388736054231119e-01,-9.061388736054231119e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.718416620816269447e+00,-1.718416620816269447e+00,-1.812277747210846224e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.434185455985974578e+00,-9.056975759976624296e-01,-9.056975759976624296e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.717092727992987289e+00,-1.717092727992987289e+00,-1.811395151995324859e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.421367941493065779e+00,-9.035613235821775557e-01,-9.035613235821775557e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.710683970746532889e+00,-1.710683970746532889e+00,-1.807122647164355111e+00
-0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,1.000000000000000000e+00,1.000000000000000000e+00,0.000000000000000000e+00,-3.401255409798296103e+00,-9.002092349663826099e-01,-9.002092349663826099e-01,2.000000000000000000e+00,0.000000000000000000e+00,-1.700627704899148052e+00,-1.700627704899148052e+00,-1.800418469932765220e+00
diff --git a/benchmarks/SDP_Quadratic_Stability_6d/README.md b/benchmarks/SDP_Quadratic_Stability_6d/README.md
index ff20630..af911d7 100644
--- a/benchmarks/SDP_Quadratic_Stability_6d/README.md
+++ b/benchmarks/SDP_Quadratic_Stability_6d/README.md
@@ -43,20 +43,17 @@ with trace(P) = 0.033. The solver found complexity k = 1 with no degeneracy, mea
 ## Files
 
 ```
-Quadratic_Stability_6d/
+SDP_Quadratic_Stability_6d/
 ├── README.md           This file
 ├── parameters.txt      Solver parameters (rho, tau, confidence)
 ├── generate.py         Generate scenarios and constraint matrices
 ├── run.py              Solve the SDP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json         One-shot program definition (JSON)
-│   ├── program_symbolic.json          One-shot program definition (MATLAB)
-│   ├── F_0.csv ... F_6.csv   Scenario-dependent LMI matrices
-│   ├── E_0.csv ... E_6.csv   Hard constraint (P >= 0.01*I) matrices
-│   ├── c.csv                  Linear objective vector (trace)
-│   ├── Q.csv                  Quadratic regularisation matrix
-│   └── scenarios.csv          500 x 2 parameter samples
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
+│   ├── scenarios.csv          500 x 2 parameter samples
+│   └── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 └── results/
     ├── metrics.json    Solver output (cost, risk bounds, etc.)
     ├── solution.csv    Raw solution vector
@@ -70,8 +67,9 @@ Quadratic_Stability_6d/
 # Generate scenarios and constraint matrices (optional)
 python generate.py
 
-# Solve the SDP
+# Solve the SDP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -83,32 +81,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **SDP** tab, formulation: **Robust**
-2. Upload matrices using one of these approaches:
-
-   **Option A — One-shot LMI upload:**
-   - Click **Edit F(delta)** -> upload `data/program_symbolic.json` (the F_d section) as a JSON file, or enter the F_d dict
-   - Click **Edit E** -> upload the E matrices similarly
-
-   **Option B — Individual matrix entry:**
-   - Click **Edit F(delta)** -> set n = 6 -> click each F_i button and upload `data/F_0.csv` through `data/F_6.csv`
-   - Click **Edit E** -> set n = 6 -> click each E_i button and upload `data/E_0.csv` through `data/E_6.csv`
-
-3. Upload or enter:
-   - **c**: `data/c.csv`
-   - **Q**: `data/Q.csv`
+2. Enter the matrices with their **Edit** buttons, using the matching fields of `data/program_symbolic.json`:
+   - **F(delta)**: `F_d` (matrices F_0 ... F_6)
+   - **E**: `E` (matrices E_0 ... E_6)
+3. Enter:
+   - **c**: `c`
+   - **Q**: `Q`
 4. Set parameters: rho = 0.0, tau = 0, confidence (beta) = 1e-06
 5. Upload `data/scenarios.csv` in the Scenarios box
 6. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 0.032858958578525506
+- Optimal cost: 0.03285895857852664
 - Complexity k: 1
-- Risk bounds: [0.0, 0.04000231437711046]
+- Risk bounds: [0.0, 0.04027646044222639]
diff --git a/benchmarks/SDP_Quadratic_Stability_6d/parameters.txt b/benchmarks/SDP_Quadratic_Stability_6d/parameters.txt
index 7d5f3c9..a96848a 100644
--- a/benchmarks/SDP_Quadratic_Stability_6d/parameters.txt
+++ b/benchmarks/SDP_Quadratic_Stability_6d/parameters.txt
@@ -10,4 +10,5 @@
 
 rho = 0.0
 tau = 0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/SDP_Quadratic_Stability_6d/plot.py b/benchmarks/SDP_Quadratic_Stability_6d/plot.py
index cc7b6e2..78a0fe0 100644
--- a/benchmarks/SDP_Quadratic_Stability_6d/plot.py
+++ b/benchmarks/SDP_Quadratic_Stability_6d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the Quadratic Stability 6D benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/SDP_Quadratic_Stability_6d/run.py b/benchmarks/SDP_Quadratic_Stability_6d/run.py
index 3fea3c3..38eec9c 100644
--- a/benchmarks/SDP_Quadratic_Stability_6d/run.py
+++ b/benchmarks/SDP_Quadratic_Stability_6d/run.py
@@ -6,12 +6,13 @@
 oscillator LPV system with 2D parameter uncertainty.
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..'))
@@ -66,6 +67,11 @@ def load_vector(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 60)
     print("BENCHMARK: Quadratic_Stability_6d (SDP)")
     print("6D Quadratic Stability with 2D Parameter Uncertainty")
@@ -101,17 +107,17 @@ def main():
     print()
 
     # Solve
-    print("Solving SDP with MOSEK...")
+    print(f"Solving SDP with {args.solver}...")
 
     try:
         x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
             deltas=scenarios, F_d=F_d, E=E_mats, c=c, Q=Q,
             tau=0.0, x_ref=np.zeros(n_vars), rho=0.0,
-            norm_type=2, solver='MOSEK'
+            norm_type=2, solver=args.solver
         )
-        print(f"Solved with MOSEK")
+        print(f"Solved with {args.solver}")
     except Exception as e:
-        print(f"MOSEK failed: {e}")
+        print(f"{args.solver} failed: {e}")
         print("Attempting with SCS solver...")
         try:
             x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -134,9 +140,9 @@ def main():
     print(f"Optimal Cost: {cost:.6f}")
     print(f"Complexity (k): {k}")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 60)
 
     print()
diff --git a/benchmarks/SDP_covariance_wine_15d/README.md b/benchmarks/SDP_covariance_wine_15d/README.md
index 1c0a782..9fd9e0a 100644
--- a/benchmarks/SDP_covariance_wine_15d/README.md
+++ b/benchmarks/SDP_covariance_wine_15d/README.md
@@ -22,7 +22,7 @@ so that (delta_i - mean)(delta_i - mean)' - Sigma <= 0 (PSD sense), i.e., Sigma
 
 ![Robust Covariance Wine](results/robust_covariance_wine.png)
 
-The optimal Sigma_hat has trace 75.71 and eigenvalues {6.20, 9.29, 12.59, 21.05, 26.58}, confirming positive definiteness. The Frobenius error ||Sigma_hat - Sigma_full||_F = 35.61 reflects the conservatism of the domination requirement: every individual point's outer product must be dominated, inflating the diagonal entries substantially. The complexity is k = 7, no degeneracy, and risk bounds [0.000, 0.175] at 99.9999% confidence.
+The optimal Sigma_hat has trace 75.71 and eigenvalues {6.20, 9.29, 12.59, 21.05, 26.58}, confirming positive definiteness. The Frobenius error ||Sigma_hat - Sigma_full||_F = 35.61 reflects the conservatism of the domination requirement: every individual point's outer product must be dominated, inflating the diagonal entries substantially. The complexity is k = 7, no degeneracy, and risk bounds [0.000, 0.178] at 99.9999% confidence.
 
 ## Files
 
@@ -34,21 +34,22 @@ SDP_covariance_wine_15d/
 ├── run.py              Solve the SDP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json         One-shot program definition (JSON)
-│   ├── program_symbolic.json          One-shot program definition (MATLAB)
+│   ├── program_symbolic.json  Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json   Upload-Program definition (numeric mode)
 │   ├── scenarios.csv          178 x 5 data point scenarios
+│   ├── scenarios_numeric.csv  Per-row-flattened matrices for numeric mode
 │   ├── basis_matrices.npy     15 x 5 x 5 symmetric basis matrices Pi_k
 │   ├── free_entries.csv       Mapping of free entries in symmetric matrix
-│   ├── c.csv                  Linear objective vector (trace)
-│   ├── Q.csv                  Quadratic objective matrix (zero)
-│   ├── E_0.csv ... E_15.csv   Hard constraint matrices (Sigma > 0)
 │   ├── data_standardized.csv  Standardized wine data
 │   ├── feature_names.txt      Feature labels
-│   └── cov_full.csv           Full-sample covariance for comparison
+│   ├── cov_full.csv           Full-sample covariance for comparison
+│   ├── entry_info.csv         Legacy entry table (not used by the current scripts)
+│   └── entry_matrices.npy     Legacy entry matrices (not used by the current scripts)
 └── results/
     ├── metrics.json                Solver output (cost, risk bounds, etc.)
     ├── solution.csv                Raw solution vector
     ├── Sigma_solution.csv          Reconstructed covariance matrix
+    ├── X_solution.csv              Legacy solution matrix (not written by the current run.py)
     ├── robust_covariance_wine.png  Paper figure (300 dpi)
     └── robust_covariance_wine.pdf  Paper figure (vector)
 ```
@@ -59,8 +60,9 @@ SDP_covariance_wine_15d/
 # Generate scenarios from Wine dataset (optional)
 python generate.py
 
-# Solve the SDP
+# Solve the SDP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -72,32 +74,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **SDP** tab, formulation: **Robust**
-2. Upload matrices using one of these approaches:
-
-   **Option A — One-shot LMI upload:**
-   - Click **Edit F(delta)** -> upload `data/program_symbolic.json` (the F_d section) as a JSON file, or enter the F_d dict
-   - Click **Edit E** -> upload the E matrices similarly
-
-   **Option B — Individual matrix entry:**
-   - Click **Edit F(delta)** -> set n = 15 -> click each F_i button and upload `data/F_0.csv` through `data/F_15.csv`
-   - Click **Edit E** -> set n = 15 -> click each E_i button and upload `data/E_0.csv` through `data/E_15.csv`
-
-3. Upload or enter:
-   - **c**: `data/c.csv`
-   - **Q**: `data/Q.csv`
+2. Enter the matrices with their **Edit** buttons, using the matching fields of `data/program_symbolic.json`:
+   - **F(delta)**: `F_d` (matrices F_0 ... F_15)
+   - **E**: `E` (matrices E_0 ... E_15)
+3. Enter:
+   - **c**: `c`
+   - **Q**: `Q`
 4. Set parameters: rho = 0.0, tau = 0.0, confidence (beta) = 1e-06
 5. Upload `data/scenarios.csv` in the Scenarios box
 6. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: 75.71438155818954
+- Optimal cost: 75.71438155818643
 - Complexity k: 7
-- Risk bounds: [0.0, 0.17538841818544354]
+- Risk bounds: [0.0, 0.17795097645963384]
diff --git a/benchmarks/SDP_covariance_wine_15d/generate.py b/benchmarks/SDP_covariance_wine_15d/generate.py
index c225580..4d441f4 100644
--- a/benchmarks/SDP_covariance_wine_15d/generate.py
+++ b/benchmarks/SDP_covariance_wine_15d/generate.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires scikit-learn (not in requirements.txt): pip install scikit-learn
 """
 Generate robust covariance estimation data from UCI Wine dataset for SDP formulation.
 
diff --git a/benchmarks/SDP_covariance_wine_15d/parameters.txt b/benchmarks/SDP_covariance_wine_15d/parameters.txt
index c703ff0..a4612f5 100644
--- a/benchmarks/SDP_covariance_wine_15d/parameters.txt
+++ b/benchmarks/SDP_covariance_wine_15d/parameters.txt
@@ -8,4 +8,5 @@
 #
 rho = 0.0
 tau = 0.0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/SDP_covariance_wine_15d/plot.py b/benchmarks/SDP_covariance_wine_15d/plot.py
index f40be3b..a3e5da2 100644
--- a/benchmarks/SDP_covariance_wine_15d/plot.py
+++ b/benchmarks/SDP_covariance_wine_15d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the wine covariance estimation benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/SDP_covariance_wine_15d/run.py b/benchmarks/SDP_covariance_wine_15d/run.py
index 971b65b..dc947e0 100644
--- a/benchmarks/SDP_covariance_wine_15d/run.py
+++ b/benchmarks/SDP_covariance_wine_15d/run.py
@@ -14,12 +14,13 @@
   Objective: minimize trace(Sigma)
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..'))
@@ -78,6 +79,11 @@ def load_parameters(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 60)
     print("BENCHMARK: Robust Covariance Estimation (SDP)")
     print("UCI Wine — Minimum-Trace Covariance Domination")
@@ -120,7 +126,7 @@ def main():
     print()
 
     # Solve
-    print("Solving SDP with MOSEK...")
+    print(f"Solving SDP with {args.solver}...")
 
     try:
         x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -133,12 +139,12 @@ def main():
             x_ref=np.zeros(n_vars),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
-        print(f"Solved with MOSEK")
+        print(f"Solved with {args.solver}")
     except Exception as e:
-        print(f"MOSEK failed: {e}")
+        print(f"{args.solver} failed: {e}")
         print("Attempting with SCS solver...")
         try:
             x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -191,9 +197,9 @@ def main():
     print(f"Objective Value: {cost:.4f}")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 60)
 
     print()
diff --git a/benchmarks/SDP_ellipsoid_enclosure_15d/README.md b/benchmarks/SDP_ellipsoid_enclosure_15d/README.md
index 59ac6cd..032a7e4 100644
--- a/benchmarks/SDP_ellipsoid_enclosure_15d/README.md
+++ b/benchmarks/SDP_ellipsoid_enclosure_15d/README.md
@@ -22,7 +22,7 @@ so that (delta_i - x_bar)' P (delta_i - x_bar) - 1 <= 0. The hard constraint enf
 
 ![Minimum Enclosing Ellipsoid](results/minimum_enclosing_ellipsoid.png)
 
-The optimal P has trace 44.59 and eigenvalues {~0, ~0, ~0, ~0, 44.59}, revealing that the standardized data effectively lies in a 1-dimensional subspace (the first 5 breast cancer features are highly correlated). Since all 569 data points are used as scenarios, the ellipsoid contains all of them (100% containment). The complexity is k = 5 (five boundary data points define the ellipsoid), no degeneracy, and risk bounds [0.000, 0.051] at 99.9999% confidence, guaranteeing that a new random data point will fall inside the ellipsoid with probability at least 94.9%.
+The optimal P has trace 44.59 and eigenvalues {~0, ~0, ~0, ~0, 44.59}, revealing that the standardized data effectively lies in a 1-dimensional subspace (the first 5 breast cancer features are highly correlated). Since all 569 data points are used as scenarios, the ellipsoid contains all of them (100% containment). The complexity is k = 5 (five boundary data points define the ellipsoid), no degeneracy, and risk bounds [0.000, 0.052] at 99.9999% confidence, guaranteeing that a new random data point will fall inside the ellipsoid with probability at least 94.8%.
 
 ## Files
 
@@ -34,16 +34,15 @@ SDP_ellipsoid_enclosure_15d/
 ├── run.py              Solve the SDP and print results
 ├── plot.py             Generate the paper figure
 ├── data/
-│   ├── program_symbolic.json             One-shot program definition (JSON)
-│   ├── program_symbolic.json              One-shot program definition (MATLAB)
+│   ├── program_symbolic.json      Upload-Program definition (symbolic mode)
+│   ├── program_numeric.json       Upload-Program definition (numeric mode)
 │   ├── scenarios.csv              569 x 5 data point scenarios
+│   ├── scenarios_numeric.csv      Per-row-flattened matrices for numeric mode
 │   ├── basis_matrices.npy         15 x 5 x 5 symmetric basis matrices Pi_k
 │   ├── free_entries.csv           Mapping of free entries in symmetric matrix
 │   ├── coeffs.csv                 Precomputed containment coefficients
 │   ├── center.csv                 Data center x_bar
-│   ├── c.csv                      Linear objective vector (-trace)
-│   ├── Q.csv                      Quadratic objective matrix (zero)
-│   ├── E_0.csv ... E_15.csv       Hard constraint matrices (P > 0)
+│   ├── cov_full.csv               Full-sample covariance
 │   ├── data_standardized.csv      Full standardized dataset (569 x 5)
 │   └── feature_names.txt          Feature labels
 └── results/
@@ -60,8 +59,9 @@ SDP_ellipsoid_enclosure_15d/
 # Generate scenarios from Breast Cancer dataset (optional)
 python generate.py
 
-# Solve the SDP
+# Solve the SDP (MOSEK by default; use --solver to pick another, e.g. CLARABEL)
 python run.py
+python run.py --solver CLARABEL
 
 # Generate paper figure
 python plot.py
@@ -73,32 +73,25 @@ python plot.py
 
 1. Start the web server: `python3 app.py`
 2. Click **"Upload Program"** button (next to LP/QP/SDP tabs)
-3. Upload `data/program_symbolic.json` or `data/program_symbolic.json`
-4. Upload `data/scenarios.csv` in the Scenarios box
+3. Upload `data/program_symbolic.json` (or `data/program_numeric.json` for numeric mode)
+4. Upload `data/scenarios.csv` (or `data/scenarios_numeric.csv` with `program_numeric.json`) in the Scenarios box
 5. Set solver to **MOSEK** and press **Solve**
 
 ### Manual Method
 
 1. Select the **SDP** tab, formulation: **Robust**
-2. Upload matrices using one of these approaches:
-
-   **Option A — One-shot LMI upload:**
-   - Click **Edit F(delta)** -> upload `data/program_symbolic.json` (the F_d section) as a JSON file, or enter the F_d dict
-   - Click **Edit E** -> upload the E matrices similarly
-
-   **Option B — Individual matrix entry:**
-   - Click **Edit F(delta)** -> set n = 15 -> click each F_i button and upload `data/F_0.csv` through `data/F_15.csv`
-   - Click **Edit E** -> set n = 15 -> click each E_i button and upload `data/E_0.csv` through `data/E_15.csv`
-
-3. Upload or enter:
-   - **c**: `data/c.csv`
-   - **Q**: `data/Q.csv`
+2. Enter the matrices with their **Edit** buttons, using the matching fields of `data/program_symbolic.json`:
+   - **F(delta)**: `F_d` (matrices F_0 ... F_15)
+   - **E**: `E` (matrices E_0 ... E_15)
+3. Enter:
+   - **c**: `c`
+   - **Q**: `Q`
 4. Set parameters: rho = 0.0, tau = 0.0, confidence (beta) = 1e-06
 5. Upload `data/scenarios.csv` in the Scenarios box
 6. Press **Solve**
 
 ### Expected Results
 
-- Optimal cost: -44.59438886598514
+- Optimal cost: -44.594388865949526
 - Complexity k: 5
-- Risk bounds: [0.0, 0.05102523978920701]
+- Risk bounds: [0.0, 0.05181163547131739]
diff --git a/benchmarks/SDP_ellipsoid_enclosure_15d/generate.py b/benchmarks/SDP_ellipsoid_enclosure_15d/generate.py
index e951f82..0a22aaa 100644
--- a/benchmarks/SDP_ellipsoid_enclosure_15d/generate.py
+++ b/benchmarks/SDP_ellipsoid_enclosure_15d/generate.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires scikit-learn (not in requirements.txt): pip install scikit-learn
 """
 Generate minimum enclosing ellipsoid data from Breast Cancer Wisconsin dataset for SDP.
 
diff --git a/benchmarks/SDP_ellipsoid_enclosure_15d/parameters.txt b/benchmarks/SDP_ellipsoid_enclosure_15d/parameters.txt
index 29c9dee..be1c092 100644
--- a/benchmarks/SDP_ellipsoid_enclosure_15d/parameters.txt
+++ b/benchmarks/SDP_ellipsoid_enclosure_15d/parameters.txt
@@ -8,4 +8,5 @@
 #
 rho = 0.0
 tau = 0.0
+# confidence is the confidence level 1 - beta (here beta = 1e-06); the program JSON's "confidence" holds beta
 confidence = 0.999999
diff --git a/benchmarks/SDP_ellipsoid_enclosure_15d/plot.py b/benchmarks/SDP_ellipsoid_enclosure_15d/plot.py
index d708324..13729ae 100644
--- a/benchmarks/SDP_ellipsoid_enclosure_15d/plot.py
+++ b/benchmarks/SDP_ellipsoid_enclosure_15d/plot.py
@@ -1,4 +1,5 @@
 #!/usr/bin/env python3
+# Requires matplotlib (not in requirements.txt): pip install matplotlib
 """Generate a clean, paper-ready figure for the minimum enclosing ellipsoid benchmark."""
 
 import sys, os, json
diff --git a/benchmarks/SDP_ellipsoid_enclosure_15d/run.py b/benchmarks/SDP_ellipsoid_enclosure_15d/run.py
index 17f1836..7d81602 100644
--- a/benchmarks/SDP_ellipsoid_enclosure_15d/run.py
+++ b/benchmarks/SDP_ellipsoid_enclosure_15d/run.py
@@ -13,12 +13,13 @@
   Objective: maximize trace(P) = minimize -trace(P)  (tightest ellipsoid)
 
 Usage:
-    python run.py
+    python run.py [--solver SOLVER]
 """
 
 import sys
 import os
 import json
+import argparse
 import numpy as np
 
 sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..', '..'))
@@ -77,6 +78,11 @@ def load_parameters(filepath):
 
 
 def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--solver', default='MOSEK',
+                        help='CVXPY solver, e.g. MOSEK or CLARABEL (default: MOSEK)')
+    args = parser.parse_args()
+
     print("=" * 60)
     print("BENCHMARK: Minimum Enclosing Ellipsoid (SDP)")
     print("Breast Cancer Wisconsin — Tightest Containment Ellipsoid")
@@ -120,7 +126,7 @@ def main():
     print()
 
     # Solve
-    print("Solving SDP with MOSEK...")
+    print(f"Solving SDP with {args.solver}...")
 
     try:
         x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -133,12 +139,12 @@ def main():
             x_ref=np.zeros(n_vars),
             rho=0.0,
             norm_type=2,
-            solver='MOSEK'
+            solver=args.solver
         )
         status = "SUCCESS"
-        print(f"Solved with MOSEK")
+        print(f"Solved with {args.solver}")
     except Exception as e:
-        print(f"MOSEK failed: {e}")
+        print(f"{args.solver} failed: {e}")
         print("Attempting with SCS solver...")
         try:
             x, zeta, cost, N_out, k, constraints, degeneracy = solve_sdp(
@@ -206,9 +212,9 @@ def main():
     print(f"Objective Value: {cost:.4f} (= -trace(P))")
     print(f"Complexity (k): {k} support constraints")
     if degeneracy:
-        print(f"Risk Bounds (99%): [unreliable, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [unreliable, {eps_upper:.4f}]")
     else:
-        print(f"Risk Bounds (99%): [{eps_lower:.4f}, {eps_upper:.4f}]")
+        print(f"Risk Bounds ({(1 - beta) * 100:g}%): [{eps_lower:.4f}, {eps_upper:.4f}]")
     print("-" * 60)
 
     print()
diff --git a/screenshots/scrn_1_tool.png b/screenshots/scrn_1_tool.png
index c2d4f1b..d7fbd5e 100644
Binary files a/screenshots/scrn_1_tool.png and b/screenshots/scrn_1_tool.png differ
diff --git a/screenshots/scrn_2_results.png b/screenshots/scrn_2_results.png
index 85fd6f9..372b7ba 100644
Binary files a/screenshots/scrn_2_results.png and b/screenshots/scrn_2_results.png differ
diff --git a/screenshots/scrn_3_tau_graph.png b/screenshots/scrn_3_tau_graph.png
index 3340a2d..d605197 100644
Binary files a/screenshots/scrn_3_tau_graph.png and b/screenshots/scrn_3_tau_graph.png differ
diff --git a/screenshots/scrn_4_rho_graph.png b/screenshots/scrn_4_rho_graph.png
index fcb2479..9bf2867 100644
Binary files a/screenshots/scrn_4_rho_graph.png and b/screenshots/scrn_4_rho_graph.png differ
diff --git a/src/LP.py b/src/LP.py
index 6b1a5a1..390835c 100644
--- a/src/LP.py
+++ b/src/LP.py
@@ -1,6 +1,6 @@
 import cvxpy as cp
 import numpy as np
-from src.Miscellaneous import get_support
+from src.Miscellaneous import get_support, check_scenarios, regularization_term
 
 def 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):
     r"""Solve a linear program with scenario constraints via the scenario approach.
@@ -43,10 +43,14 @@ def solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
     rho : float, optional
         Penalty on slack variables. When ``0.0`` the scenario constraints are
         hard. Default is ``0.0``.
-    norm_type : int, optional
-        Norm type for the regularization term (1, 2, etc.). Default is ``2``.
+    norm_type : int, float or str, optional
+        Order p of the vector norm in the regularization term: any number
+        ``p >= 1``, ``'inf'`` or ``'fro'`` (Euclidean). Default is ``2``.
     solver : str or None, optional
         CVXPY solver name. ``None`` for automatic selection.
+    include_slack : bool or None, optional
+        Whether to add the slack variables ζ_i (the relaxation formulation).
+        ``None`` adds them only when ``rho != 0``.
 
     Returns
     -------
@@ -68,14 +72,12 @@ def solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
     Raises
     ------
     ValueError
-        If the solver does not reach an optimal or near-optimal status.
+        If ``deltas`` is empty, or the solver does not reach an optimal or
+        near-optimal status.
     """
+    N = check_scenarios(deltas)  # Number of scenarios
     n = A_d(deltas[0]).shape[1]  # Number of variables
     m_scenario = A_d(deltas[0]).shape[0]  # Number of scenario constraints
-    try:
-        N = deltas.shape[0]  # Number of scenarios
-    except IndexError:
-        raise ValueError("The input `deltas` must have at least one row.")
 
     # Variables
     x = cp.Variable((n, 1))
@@ -94,12 +96,14 @@ def solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
 
     constraints = []
     for i in range(N):
-        constraints.append(A_d(deltas[i]) @ x + b_d(deltas[i]) <= zeta[i])  # Per-scenario slack
+        # b(δ) as a column: a 1-D (m,) array would otherwise be broadcast to (m, m).
+        b_i = np.asarray(b_d(deltas[i]), dtype=float).reshape(-1, 1)
+        constraints.append(A_d(deltas[i]) @ x + b_i <= zeta[i])  # Per-scenario slack
 
     # Hard constraints are kept separate from the scenario constraints: they
     # are always enforced but never candidates for the support list.
     if not (G.size == 0 or h.size == 0):
-        non_risk_constraints = [G @ x + h <= 0]
+        non_risk_constraints = [np.atleast_2d(G) @ x + np.asarray(h, dtype=float).reshape(-1, 1) <= 0]
     else:
         non_risk_constraints = []
 
@@ -119,7 +123,7 @@ def solve_lp(deltas, A_d, b_d, G, h, c, tau=0.0, x_ref=np.array([0.0]), rho=0.0,
         )
 
     # Objective Function
-    objective = cp.Minimize(c.T @ x + tau * cp.norm(x - x_ref, norm_type) + rho * cp.sum(zeta))
+    objective = cp.Minimize(c.T @ x + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))
 
     # Solve the problem
     prob = cp.Problem(objective, constraints + non_risk_constraints)
diff --git a/src/Miscellaneous.py b/src/Miscellaneous.py
index 9edc53e..ffac9b9 100644
--- a/src/Miscellaneous.py
+++ b/src/Miscellaneous.py
@@ -4,20 +4,38 @@
 
 SOLVER_CAPABILITIES = {
     'CLARABEL': ['LP', 'QP', 'SDP'], 'SCS': ['LP', 'QP', 'SDP'],
-    # OSQP and SCIPY are removed from the LP list — both technically claim
-    # LP support but fail on the scenario LPs we build (slack reformulation
-    # produces degenerate QPs OSQP can't handle; SCIPY's wrapper hits
-    # NotImplemented paths for our equality reductions).
-    'OSQP': ['QP'], 'ECOS': ['LP', 'QP', 'SDP'],
+    # OSQP is kept off the LP list: it fails on the scenario LPs we build
+    # (the slack reformulation produces degenerate QPs OSQP can't handle).
+    'OSQP': ['QP'], 'ECOS': ['LP', 'QP'],
     'CVXOPT': ['LP', 'QP', 'SDP'], 'GLOP': ['LP'], 'GLPK': ['LP'],
-    'GLPK_MI': ['LP'], 'SCIPY': ['QP'], 'HIGHS': ['LP', 'QP'],
-    'SCIP': ['LP', 'QP', 'SDP'], 'CBC': ['LP', 'QP'],
+    'GLPK_MI': ['LP'], 'SCIPY': ['LP'], 'HIGHS': ['LP', 'QP'],
+    'SCIP': ['LP', 'QP'], 'CBC': ['LP'],
     'DAQP': ['LP', 'QP'], 'PIQP': ['LP', 'QP'],
-    'PROXQP': ['LP', 'QP'], 'QPALM': ['LP', 'QP'],
-    'MOSEK': ['LP', 'QP', 'SDP'], 'GUROBI': ['LP', 'QP', 'SDP'],
-    'CPLEX': ['LP', 'QP', 'SDP'], 'SDPA': ['SDP'],
+    'PROXQP': ['LP', 'QP'],
+    'MOSEK': ['LP', 'QP', 'SDP'], 'GUROBI': ['LP', 'QP'],
+    'CPLEX': ['LP', 'QP'], 'SDPA': ['SDP'],
 }
 
+
+def _cvxpy_capabilities(name):
+    """Program types a solver missing from SOLVER_CAPABILITIES can handle, from CVXPY's own metadata.
+
+    LP needs any conic or QP interface, QP needs a QP interface or second-order
+    cones, and SDP needs positive semidefinite cones.
+    """
+    from cvxpy.constraints import PSD, SOC
+    from cvxpy.reductions.solvers.defines import SOLVER_MAP_CONIC, SOLVER_MAP_QP
+    conic, qp = SOLVER_MAP_CONIC.get(name), SOLVER_MAP_QP.get(name)
+    supported = conic.SUPPORTED_CONSTRAINTS if conic else []
+    types = []
+    if conic or qp:
+        types.append('LP')
+    if qp or SOC in supported:
+        types.append('QP')
+    if PSD in supported:
+        types.append('SDP')
+    return types
+
 def get_solvers():
     """Retrieve installed CVXPY solvers and their supported problem types.
 
@@ -32,8 +50,47 @@ def get_solvers():
     `CVXPY solver list <https://www.cvxpy.org/tutorial/solvers/index.html#choosing-a-solver>`_
     """
     installed = cp.installed_solvers()
-    solver_dict = {s: SOLVER_CAPABILITIES.get(s, ['LP', 'QP', 'SDP']) for s in installed}
-    return solver_dict
+    solver_dict = {s: SOLVER_CAPABILITIES[s] if s in SOLVER_CAPABILITIES else _cvxpy_capabilities(s)
+                   for s in installed}
+    return {s: types for s, types in solver_dict.items() if types}
+
+def check_scenarios(deltas):
+    """Return the number of scenarios N, raising ValueError when there are none."""
+    if deltas is None or np.size(deltas) == 0:
+        raise ValueError("The input `deltas` must have at least one row.")
+    return len(deltas)
+
+
+def check_psd(Q, tol=1e-9):
+    """Return Q symmetrised, after checking it is symmetric positive semidefinite.
+
+    Both checks allow a small relative tolerance, so rank-deficient matrices
+    such as D.T @ D, whose zero eigenvalues come out as about -1e-16, pass.
+    """
+    Q = np.asarray(Q, dtype=float)
+    scale = max(1.0, float(np.abs(Q).max()))
+    assert np.allclose(Q, Q.T, rtol=0, atol=tol * scale), "Q must be symmetric"
+    Q = (Q + Q.T) / 2
+    eigvals = np.linalg.eigvalsh(Q)
+    assert eigvals.min() >= -tol * scale, f"Q must be positive semi-definite (eigenvalues: {eigvals})"
+    return Q
+
+
+def regularization_term(tau, x, x_ref, norm_type):
+    """Return tau * ||x - x_ref||_p, or 0 when tau = 0.
+
+    The term is left out when tau = 0, so a robust or relaxed LP stays a pure
+    LP: with the default p = 2 the norm would add a second-order cone, which
+    LP-only solvers cannot handle. The norm is the vector p-norm of
+    x - x_ref, so any p >= 1 works, and ``'fro'`` is the Euclidean norm.
+    """
+    if tau == 0:
+        return 0
+    diff = x - x_ref
+    if len(diff.shape) == 2:
+        diff = diff[:, 0]  # (n, 1) column -> vector
+    return tau * cp.norm(diff, 2 if norm_type == 'fro' else norm_type)
+
 
 def get_norm_types():
     """Retrieve available norm types for the regularization term.
@@ -41,16 +98,19 @@ def get_norm_types():
     Returns
     -------
     list
-        Supported CVXPY norm types: ``1`` (L1), ``2`` (L2), ``'inf'``
-        (L-infinity), ``'fro'`` (Frobenius), ``'nuc'`` (nuclear).
+        Common norm types: ``1`` (L1), ``2`` (L2), ``'inf'`` (L-infinity)
+        and ``'fro'`` (the same as L2 for the vector ``x - x_ref``). Any
+        other number ``p >= 1`` also works.
     """
-    norm_types = [1, 2, "inf", "fro", "nuc"]
+    norm_types = [1, 2, "inf", "fro"]
     return norm_types
 
 def load_file(file_path):
     """Read a data file and convert it to a NumPy array.
 
-    Supports CSV, TXT, XLSX, and JSON formats.
+    Supports CSV, TXT, XLSX, and JSON formats. A ``.txt`` file may be
+    separated by whitespace or by commas; a ``.json`` file holds a list of
+    rows (one per scenario) or a table in pandas' default JSON format.
 
     Parameters
     ----------
@@ -61,22 +121,21 @@ def load_file(file_path):
     Returns
     -------
     numpy.ndarray or None
-        Array containing the file data, or ``None`` if an error occurred.
-
-    Raises
-    ------
-    ValueError
-        If the file extension is not supported.
+        Array containing the file data, or ``None`` if the file cannot be
+        read or has an unsupported extension (the error is printed).
     """
     try:
         if file_path.endswith('.csv'):
             data = pd.read_csv(file_path,header=None).values  # Read CSV using pandas and convert to NumPy
         elif file_path.endswith('.txt'):
-            data = np.loadtxt(file_path)  # Read TXT using NumPy (default to expect floats)
+            try:
+                data = np.loadtxt(file_path)  # whitespace-separated numbers
+            except ValueError:
+                data = np.loadtxt(file_path, delimiter=',')  # comma-separated numbers
         elif file_path.endswith('.xlsx'):
             data = pd.read_excel(file_path, engine='openpyxl',header=None).values  # Read Excel file and convert to NumPy
         elif file_path.endswith('.json'):
-            data = pd.read_json(file_path, orient='record').values  # Read JSON and convert to NumPy
+            data = pd.read_json(file_path).values  # a list of rows, or pandas' column format
         else:
             raise ValueError("Unsupported file format")
 
diff --git a/src/QP.py b/src/QP.py
index 36020d6..8c9c4bd 100644
--- a/src/QP.py
+++ b/src/QP.py
@@ -1,6 +1,6 @@
 import cvxpy as cp
 import numpy as np
-from src.Miscellaneous import get_support
+from src.Miscellaneous import get_support, check_scenarios, check_psd, regularization_term
 
 def 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):
     r"""Solve a quadratic program with scenario constraints via the scenario approach.
@@ -46,10 +46,14 @@ def solve_qp(deltas, A_d, b_d, G, h, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0
     rho : float, optional
         Penalty on slack variables. When ``0.0`` the scenario constraints are
         hard. Default is ``0.0``.
-    norm_type : int, optional
-        Norm type for the regularization term (1, 2, etc.). Default is ``2``.
+    norm_type : int, float or str, optional
+        Order p of the vector norm in the regularization term: any number
+        ``p >= 1``, ``'inf'`` or ``'fro'`` (Euclidean). Default is ``2``.
     solver : str or None, optional
         CVXPY solver name. ``None`` for automatic selection.
+    include_slack : bool or None, optional
+        Whether to add the slack variables ζ_i (the relaxation formulation).
+        ``None`` adds them only when ``rho != 0``.
 
     Returns
     -------
@@ -75,16 +79,11 @@ def solve_qp(deltas, A_d, b_d, G, h, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0
     AssertionError
         If *Q* is not positive semidefinite or not symmetric.
     """
-    # Check Q is positive semi-definite and symmetric.
-    eigvals = np.linalg.eigvals(Q)
-    assert np.all(eigvals >= 0), f"Q must be positive semi-definite (eigenvalues: {eigvals})"
-    assert (Q == Q.T).all(), "Q must be symmetric"
+    # Check Q is symmetric positive semi-definite (up to a small tolerance).
+    Q = check_psd(Q)
+    N = check_scenarios(deltas)  # Number of scenarios
     n = A_d(deltas[0]).shape[1]  # Number of variables
     m_scenario = A_d(deltas[0]).shape[0]  # Number of scenario constraints
-    try:
-        N = deltas.shape[0]  # Number of scenarios
-    except IndexError:
-        raise ValueError("The input `deltas` must have at least one row.")
 
     # Variables
     x = cp.Variable((n, 1))
@@ -99,12 +98,14 @@ def solve_qp(deltas, A_d, b_d, G, h, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0
         zeta = np.zeros(N)
     constraints = []
     for i in range(N):
-        constraints.append(A_d(deltas[i]) @ x + b_d(deltas[i]) <= zeta[i])  # Per-scenario slack
+        # b(δ) as a column: a 1-D (m,) array would otherwise be broadcast to (m, m).
+        b_i = np.asarray(b_d(deltas[i]), dtype=float).reshape(-1, 1)
+        constraints.append(A_d(deltas[i]) @ x + b_i <= zeta[i])  # Per-scenario slack
 
     # Hard constraints are kept separate from the scenario constraints: they
     # are always enforced but never candidates for the support list.
     if not (G.size == 0 or h.size == 0):
-        non_risk_constraints = [G @ x + h <= 0]
+        non_risk_constraints = [np.atleast_2d(G) @ x + np.asarray(h, dtype=float).reshape(-1, 1) <= 0]
     else:
         non_risk_constraints = []
 
@@ -122,7 +123,7 @@ def solve_qp(deltas, A_d, b_d, G, h, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0
         )
 
     # Objective Function
-    objective = cp.Minimize((1/2)*cp.quad_form(x, cp.psd_wrap(Q)) + c.T @ x + tau * cp.norm(x - x_ref, norm_type) + rho * cp.sum(zeta))
+    objective = cp.Minimize((1/2)*cp.quad_form(x, cp.psd_wrap(Q)) + c.T @ x + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))
 
     # Solve the problem
     prob = cp.Problem(objective, constraints + non_risk_constraints)
diff --git a/src/Risk.py b/src/Risk.py
index acc340b..ee630cb 100644
--- a/src/Risk.py
+++ b/src/Risk.py
@@ -1,10 +1,10 @@
-import jax
 # The bisection compares two betainc terms whose difference is ~1e-8 near the
 # upper root; in JAX's default float32 that difference is pure noise, which
 # flips the comparison and can leave epsU understated by up to ~0.01 at
 # beta=1e-6 (the high-confidence regime this tool targets). float64 makes the
-# root agree with a float64 SciPy reference to ~1e-10.
-jax.config.update("jax_enable_x64", True)
+# root agree with a float64 SciPy reference to ~1e-10. It is switched on only
+# inside quantify_risk (enable_x64), not for the whole Python process.
+from jax.experimental import enable_x64
 from jax.scipy.special import betainc
 
 def quantify_risk(k,N,beta):
@@ -36,12 +36,30 @@ def quantify_risk(k,N,beta):
     epsU : float
         Upper bound on the constraint violation probability.
 
+    Raises
+    ------
+    ValueError
+        If ``N < 1``, ``k`` is not between 0 and ``N``, or ``beta`` is not
+        strictly between 0 and 1.
+
     Notes
     -----
     The incomplete beta function ``betainc`` used here follows the JAX/SciPy
     argument convention, which differs from MATLAB (the first argument appears
     last in JAX).
     """
+    if int(N) != N or N < 1:
+        raise ValueError(f"N must be a positive integer (got {N}).")
+    if int(k) != k or not 0 <= k <= N:
+        raise ValueError(f"k must be an integer between 0 and N = {N} (got {k}).")
+    if not 0 < beta < 1:
+        raise ValueError(f"beta must lie strictly between 0 and 1 (got {beta}).")
+    with enable_x64():
+        return _risk_bounds(int(k), int(N), beta)
+
+
+def _risk_bounds(k, N, beta):
+    """Bisection for the lower and upper risk bounds (see quantify_risk)."""
     t1 = 0.0
     t2 = k/N
     threshold = 1e-10
@@ -61,7 +79,7 @@ def quantify_risk(k,N,beta):
     #print("epsL = ", epsL)
 
     if (k==N):
-        epsU = 1 #set upper bound
+        epsU = 1.0 #set upper bound
     else:
         t1 = k/N
         t2 = 1
@@ -78,13 +96,4 @@ def quantify_risk(k,N,beta):
         epsU = t2 #set upper bound
     #print("epsU = ", epsU)
     # Return results
-    return epsL, epsU
-
-def quantify_conf():
-    """
-        To create...
-    """
-
-
-    # Return results
-    return 0.0
\ No newline at end of file
+    return float(epsL), float(epsU)
\ No newline at end of file
diff --git a/src/SDP.py b/src/SDP.py
index ebb2543..f6df8cb 100644
--- a/src/SDP.py
+++ b/src/SDP.py
@@ -1,6 +1,6 @@
 import cvxpy as cp
 import numpy as np
-from src.Miscellaneous import get_support
+from src.Miscellaneous import get_support, check_scenarios, check_psd, regularization_term
 
 
 def 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):
@@ -27,14 +27,16 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
         Function mapping a scenario to a dict of symmetric matrices keyed by
         variable index (``'0'`` for the constant term, ``'1'`` for
         :math:`x_1`, etc.):  ``F_d(delta) -> {'0': F0, '1': F1, ...}``.
+        Integer keys (``0, 1, ...``) work too.
     E : dict or None
         Dict of symmetric matrices for the hard LMI constraint, same key
         format as *F_d* output. Pass ``None`` if there are no hard constraints.
     c : numpy.ndarray
         Linear objective coefficient vector with shape ``(n,)``.
-    Q : numpy.ndarray
+    Q : numpy.ndarray or None
         Quadratic cost matrix with shape ``(n, n)``. Must be positive
-        semidefinite and symmetric.
+        semidefinite and symmetric. ``None`` or an empty array means no
+        quadratic term.
     tau : float, optional
         Regularization strength toward *x_ref*. Default is ``0.0``.
     x_ref : numpy.ndarray, optional
@@ -42,10 +44,14 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
     rho : float, optional
         Penalty on slack variables. When ``0.0`` the scenario constraints are
         hard. Default is ``0.0``.
-    norm_type : int, optional
-        Norm type for the regularization term (1, 2, etc.). Default is ``2``.
+    norm_type : int, float or str, optional
+        Order p of the vector norm in the regularization term: any number
+        ``p >= 1``, ``'inf'`` or ``'fro'`` (Euclidean). Default is ``2``.
     solver : str or None, optional
         CVXPY solver name. ``None`` for automatic selection.
+    include_slack : bool or None, optional
+        Whether to add the slack variables ζ_i (the relaxation formulation).
+        ``None`` adds them only when ``rho != 0``.
 
     Returns
     -------
@@ -72,18 +78,13 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
         If *Q* or any :math:`F_j(\delta)` matrix is not symmetric.
     """
 
-    # Check Q is positive semi-definite and symmetric.
-    eigvals = np.linalg.eigvals(Q)
-    assert np.all(eigvals >= 0), f"Q must be positive semi-definite (eigenvalues: {eigvals})"
-    assert (Q == Q.T).all(), "Q must be symmetric"
-    n = Q.shape[1]
+    n = np.asarray(c).reshape(-1).size  # Number of variables
+    # Check Q is symmetric positive semi-definite (up to a small tolerance);
+    # an empty Q means there is no quadratic term.
+    Q = np.zeros((n, n)) if Q is None or np.size(Q) == 0 else check_psd(Q)
+    N = check_scenarios(deltas)  # Number of scenarios
     m = list(F_d(deltas[0]).values())[0].shape[0]
 
-    try:
-        N = deltas.shape[0]  # Number of scenarios
-    except IndexError:
-        raise ValueError("The input `deltas` must have at least one row.")
-
     # Variables
     x = cp.Variable(n)
     # See solve_lp for the rationale: respect include_slack independently of
@@ -103,7 +104,7 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
         for k, Fk in F_dict.items():
             #assert np.all(np.linalg.eigvals(Fk) >= 0), "\\(F_j(\\delta)\\) need to be positive semi-definite and symmetric"
             assert (Fk == Fk.T).all(), "\\(F_j(\\delta)\\) need to be symmetric"
-            if k == '0':
+            if int(k) == 0:  # keys may be '0', '1', ... or 0, 1, ...
                 term = Fk
             else:
                 xk = x[int(k)-1]
@@ -116,7 +117,7 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
         for k, Fk in E.items():
             #assert np.all(np.linalg.eigvals(Fk) >= 0), "\\(F_j\\) need to be positive semi-definite and symmetric"
             assert (Fk == Fk.T).all(), "\\(F_j\\) needs to be symmetric"
-            if k == '0':
+            if int(k) == 0:  # keys may be '0', '1', ... or 0, 1, ...
                 term = Fk
             else:
                 xk = x[int(k) - 1]
@@ -129,8 +130,21 @@ def solve_sdp(deltas, F_d, E, c, Q, tau=0.0, x_ref=np.array([0.0]), rho=0.0, nor
     else:
         non_risk_constraints = []
 
+    # x is a vector here, so x_ref must be a scalar or a length-n vector: a
+    # (n, 1) column would be broadcast against x into an (n, n) matrix.
+    _xr = np.asarray(x_ref, dtype=float).reshape(-1)
+    if _xr.size == 1:
+        x_ref = float(_xr[0])
+    elif _xr.size == n:
+        x_ref = _xr
+    else:
+        raise ValueError(
+            f"x̄ has {_xr.size} entries but d = {n}. Provide either a scalar "
+            f"(broadcast across all coordinates) or a length-{n} vector."
+        )
+
     # Objective Function
-    objective = cp.Minimize((1 / 2) * cp.quad_form(x, Q) + c.T @ x + tau * cp.norm(x - x_ref, norm_type) + rho * cp.sum(zeta))
+    objective = cp.Minimize((1 / 2) * cp.quad_form(x, cp.psd_wrap(Q)) + c.T @ x + regularization_term(tau, x, x_ref, norm_type) + rho * cp.sum(zeta))
 
     # Solve the problem
     prob = cp.Problem(objective, constraints + non_risk_constraints)
diff --git a/src/mosek_solve.py b/src/mosek_solve.py
index c318e60..3682c6b 100644
--- a/src/mosek_solve.py
+++ b/src/mosek_solve.py
@@ -5,8 +5,8 @@
 reads the license path on first import. This enables concurrent
 users to each supply their own license without conflicts.
 
-Usage:
-    echo '<json>' | python3 src/mosek_solve.py /tmp/mosek_<uuid>/mosek.lic
+Usage (from the repository root):
+    echo '<json>' | python -m src.mosek_solve /tmp/mosek_<uuid>/mosek.lic
 """
 import os
 import sys
@@ -16,22 +16,11 @@
     os.environ["MOSEKLM_LICENSE_FILE"] = sys.argv[1]
 
 import json
-import time
 import numpy as np
 
-# Now safe to import solve functions (which import CVXPY → MOSEK)
-from src.LP import solve_lp
-from src.QP import solve_qp
-from src.SDP import solve_sdp
-from src.Risk import quantify_risk
-from src.parsing import (
-    generate_matrix_function,
-    generate_matrix,
-    generate_tensor_function,
-    generate_tensor,
-    generate_numeric_A_b,
-    generate_numeric_F,
-)
+# Now safe to import the solve code (which imports CVXPY → MOSEK). It is the
+# same code the app uses, so both paths report results identically.
+from src.solve_form import solve_form
 
 
 def _safe_tolist(v):
@@ -52,204 +41,9 @@ def main():
     input_data = json.loads(sys.stdin.read())
     form = input_data["form"]
     scenarios_raw = input_data.get("scenarios")
-
     scenarios = np.array(scenarios_raw, dtype=float) if scenarios_raw is not None else None
 
-    active_tab = form.get("active_tab")
-
-    # ---- build matrix functions / constants from raw form strings ----
-    A_d = b_d = F_d = A_da = b_da = None
-
-    mode = form.get("mode", "symbolic")
-    n_x = int(form.get("n_x", 0)) if form.get("n_x") else 0
-    rows_A = int(form.get("rows_A", 0)) if form.get("rows_A") else 0
-    lmi_size = int(form.get("lmi_size", 0)) if form.get("lmi_size") else 0
-
-    if mode == "numeric":
-        if scenarios is None:
-            raise ValueError("Numeric mode requires an uploaded scenario data file.")
-        if active_tab in ("lp-tab", "qp-tab"):
-            if n_x <= 0 or rows_A <= 0:
-                raise ValueError("Numeric mode requires valid n_x and rows_A.")
-            expected_cols = rows_A * n_x + rows_A
-            if scenarios.shape[1] != expected_cols:
-                raise ValueError(
-                    f"Each scenario row must have {expected_cols} elements "
-                    f"(rows_A*n_x + rows_A = {rows_A}*{n_x} + {rows_A}), "
-                    f"but got {scenarios.shape[1]}.")
-            A_d, b_d = generate_numeric_A_b(rows_A, n_x)
-        elif active_tab == "sdp-tab":
-            if n_x <= 0 or lmi_size <= 0:
-                raise ValueError("Numeric mode requires valid n_x and lmi_size.")
-            expected_cols = (n_x + 1) * lmi_size ** 2
-            if scenarios.shape[1] != expected_cols:
-                raise ValueError(
-                    f"Each scenario row must have {expected_cols} elements "
-                    f"((n_x+1)*lmi_size² = {n_x + 1}*{lmi_size}² = {expected_cols}), "
-                    f"but got {scenarios.shape[1]}.")
-            F_d = generate_numeric_F(lmi_size, n_x)
-    else:
-        if form.get("A_d"):
-            A_d = generate_matrix_function(form["A_d"])
-        elif active_tab in ("lp-tab", "qp-tab"):
-            raise ValueError("\\(A(\\delta)\\) is ill-defined")
-
-        if form.get("b_d"):
-            b_d = generate_matrix_function(form["b_d"])
-        elif active_tab in ("lp-tab", "qp-tab"):
-            raise ValueError("\\(b(\\delta)\\) is ill-defined")
-
-        if form.get("F_d"):
-            F_d = generate_tensor_function(form["F_d"])
-        elif active_tab == "sdp-tab":
-            raise ValueError("\\(F_j(\\delta)\\) is ill-defined")
-
-    if form.get("A_da"):
-        A_da = generate_tensor_function(form["A_da"])
-    if form.get("b_da"):
-        b_da = generate_tensor_function(form["b_da"])
-
-    # Hard LP/QP constraints are G, h (form fields renamed from legacy A, b).
-    G = generate_matrix(form["G"]) if form.get("G") else np.array([])
-    h_vec = generate_matrix(form["h"]) if form.get("h") else np.array([])
-    c = generate_matrix(form["c"]) if form.get("c") else np.array([])
-    Q = generate_matrix(form["Q"]) if form.get("Q") else np.array([])
-    # Hard SDP LMI collection is E (form field renamed from legacy F).
-    E = generate_tensor(form["E"]) if form.get("E") else {}
-
-    conf = float(form["confidence"]) if form.get("confidence") else 0.0
-
-    rhos = (
-        np.array([float(x) for x in form["rho"].split(",")])
-        if form.get("rho")
-        else np.array([0.0])
-    )
-    taus = (
-        np.array([float(x) for x in form["tau"].split(",")])
-        if form.get("tau")
-        else np.array([0.0])
-    )
-    # x̄ is usually sent as a JSON-encoded 2-D array by the matrix-grid modal
-    # (e.g. "[[10]]"); accept legacy comma-separated scalars as a fallback.
-    _raw_theta = form.get("theta_bar", "")
-    if not _raw_theta:
-        theta_bar = 0.0
-    else:
-        _s = _raw_theta.strip()
-        try:
-            theta_bar = np.asarray(json.loads(_s), dtype=float).reshape(-1)
-        except (ValueError, TypeError):
-            theta_bar = np.array([float(x) for x in _s.split(",") if x.strip() != ""])
-    p_raw = form.get("p", "").strip().strip("'\"")
-    if not p_raw:
-        p = 2
-    elif p_raw.lower() in ("inf", "infinity", "np.inf"):
-        p = np.inf
-    elif p_raw.lower() == "fro":
-        p = "fro"
-    else:
-        try:
-            p = float(p_raw)
-        except ValueError:
-            raise ValueError(
-                f"p (norm order) must be a single number, 'inf', or 'fro'. "
-                f"You entered {p_raw!r} — comma-separated lists are not "
-                f"supported for p (sweep over τ or ρ instead)."
-            )
-
-    # ---- solve loop (mirrors _solve_inner in app.py) ----
-    conf = conf / (len(taus) * len(rhos))
-
-    optimal_x_list, optimal_s_list, optimal_cost_list = [], [], []
-    N_list, active_list, constraints_list = [], [], []
-    risk_list, e_list, degeneracy_list = [], [], []
-    rho_list, tau_list = [], []
-    solve_time_list, risk_time_list = [], []
-
-    # Honour relaxation slack independently of ρ=0 (see solve_lp). Leave as
-    # None when no explicit option is provided so the solver defaults apply.
-    selected_option = form.get("option")
-    if selected_option:
-        include_slack = selected_option in ("relaxation", "regularization-relaxation")
-    else:
-        include_slack = None
-
-    for tau in taus:
-        for rho in rhos:
-            solver = "MOSEK"
-
-            optimal_x = np.array([])
-            optimal_s = np.array([])
-            optimal_cost = None
-            N = 0
-            complexity = []
-            constraints = []
-            degeneracy = False
-            risk = np.array([])
-            e = "None"
-            solve_time = 0.0
-            risk_time = 0.0
-
-            try:
-                t0 = time.perf_counter()
-                if active_tab == "lp-tab":
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_lp(
-                        scenarios, A_d, b_d, G, h_vec, c, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack,
-                    )
-                elif active_tab == "qp-tab":
-                    Q_inner = generate_matrix(form["Q"]) if form.get("Q") else np.array([])
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_qp(
-                        scenarios, A_d, b_d, G, h_vec, c, Q_inner, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack,
-                    )
-                elif active_tab == "sdp-tab":
-                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_sdp(
-                        scenarios, F_d, E, c, Q, tau, theta_bar, rho, p, solver,
-                        include_slack=include_slack,
-                    )
-                solve_time = time.perf_counter() - t0
-
-                t1 = time.perf_counter()
-                risk = np.array(quantify_risk(complexity, N, conf))
-                risk_time = time.perf_counter() - t1
-            except Exception as error:
-                e = str(error)
-                if "license" in e.lower() and ("mosek" in e.lower() or "rescode" in e.lower()):
-                    e = "MOSEK license cannot be located or is incorrect."
-
-            optimal_x_list.append(_safe_tolist(optimal_x))
-            optimal_s_list.append(_safe_tolist(optimal_s))
-            optimal_cost_list.append(_safe_tolist(optimal_cost))
-            N_list.append(N)
-            active_list.append(_safe_tolist(complexity))
-            e_list.append(e)
-            constraints_list.append(len(constraints))
-            risk_list.append(_safe_tolist(risk))
-            degeneracy_list.append(degeneracy)
-            rho_list.append(_safe_tolist(rho))
-            tau_list.append(_safe_tolist(tau))
-            solve_time_list.append(round(solve_time, 4))
-            risk_time_list.append(round(risk_time, 4))
-
-    result = {
-        "form_data": form,
-        "optimal_x": optimal_x_list,
-        "optimal_s": optimal_s_list,
-        "optimal_cost": optimal_cost_list,
-        "num_deltas": N_list,
-        "tot_con": constraints_list,
-        "active_con": active_list,
-        "risk": risk_list,
-        "conf": 1 - conf,
-        "tau_": tau_list,
-        "rho_": rho_list,
-        "errorcode": e_list,
-        "degeneracy": degeneracy_list,
-        "solve_time": solve_time_list,
-        "risk_time": risk_time_list,
-    }
-
+    result = solve_form(form, scenarios, solver="MOSEK")
     json.dump(result, real_stdout, default=_safe_tolist)
 
 
diff --git a/src/parsing.py b/src/parsing.py
index 50e81f1..f373779 100644
--- a/src/parsing.py
+++ b/src/parsing.py
@@ -5,7 +5,8 @@
 
 # Untrusted matrix/tensor cells arrive from the public web form, so they are
 # evaluated by a restricted AST walker instead of eval(): only numeric
-# literals, + - * / // % **, unary +/-, the name `delta` (indexed), and
+# literals, + - * / // % ** (and ^ as a power, as in SymPy), unary +/-, the
+# name `delta` (indexed), the functions and constants below, and
 # `math.<fn>(...)` are permitted. Anything else raises ValueError, which
 # closes the remote-code-execution path that eval() with an open namespace
 # would otherwise expose.
@@ -16,18 +17,32 @@
 }
 _UNARYOPS = {ast.UAdd: operator.pos, ast.USub: operator.neg}
 
+# SymPy-style names, so cells can say sin(delta[0]) instead of math.sin(delta[0]).
+_FUNCTIONS = {
+    'sin': math.sin, 'cos': math.cos, 'tan': math.tan,
+    'asin': math.asin, 'acos': math.acos, 'atan': math.atan, 'atan2': math.atan2,
+    'sinh': math.sinh, 'cosh': math.cosh, 'tanh': math.tanh,
+    'exp': math.exp, 'log': math.log, 'ln': math.log, 'log10': math.log10, 'log2': math.log2,
+    'sqrt': math.sqrt, 'abs': abs, 'Abs': abs,
+    'floor': math.floor, 'ceiling': math.ceil, 'ceil': math.ceil,
+    'min': min, 'Min': min, 'max': max, 'Max': max,
+}
+_CONSTANTS = {'pi': math.pi, 'e': math.e, 'E': math.e}
+
 
 def safe_eval(expr, delta=None):
     """Safely evaluate a scalar arithmetic expression over `delta` and `math`.
 
     Parameters
     ----------
-    expr : str
+    expr : str or number
         Expression such as ``"delta[0]"``, ``"-delta[3]"``, ``"0.5"``,
-        or ``"math.sin(delta[1])"``.
-    delta : sequence or None
-        The uncertainty vector referenced by ``delta[i]``; ``None`` for
-        constant expressions.
+        ``"sin(delta[1])"``, ``"delta[0]^2"`` or ``"math.sin(delta[1])"``.
+        A plain number (e.g. from a JSON or MAT file) is returned as a float.
+    delta : sequence, scalar or None
+        The uncertainty vector referenced by ``delta[i]`` (a scalar is
+        treated as a vector with one component); ``None`` for constant
+        expressions.
 
     Returns
     -------
@@ -54,9 +69,19 @@ def ev(node):
                 if delta is None:
                     raise ValueError("delta is not available in this expression")
                 return delta
+            if node.id in _CONSTANTS:
+                return _CONSTANTS[node.id]
+            if node.id in _FUNCTIONS:
+                return _FUNCTIONS[node.id]
             raise ValueError(f"disallowed name: {node.id}")
         if isinstance(node, ast.Subscript):
-            return ev(node.value)[int(ev(node.slice))]
+            base, index = ev(node.value), int(ev(node.slice))
+            if np.ndim(base) == 0:  # a scalar scenario has one component
+                base = [base]
+            try:
+                return base[index]
+            except IndexError:
+                raise ValueError(f"delta[{index}] does not exist: the scenario has {len(base)} component(s)")
         if isinstance(node, ast.Attribute):
             if isinstance(node.value, ast.Name) and node.value.id == "math" \
                     and not node.attr.startswith("_"):
@@ -69,7 +94,17 @@ def ev(node):
             return func(*[ev(a) for a in node.args])
         raise ValueError(f"disallowed expression element: {type(node).__name__}")
 
-    return ev(ast.parse(expr, mode="eval"))
+    if isinstance(expr, (int, float, np.integer, np.floating)) and not isinstance(expr, bool):
+        return float(expr)
+    # ^ means power, as in SymPy. It is replaced before parsing so that it gets
+    # the precedence of ** (Python's ^ is XOR, which binds more loosely than +).
+    try:
+        value = ev(ast.parse(expr.replace("^", "**"), mode="eval"))
+    except TypeError as error:  # e.g. a function applied to the wrong arguments
+        raise ValueError(f"invalid expression {expr!r}: {error}")
+    if isinstance(value, bool) or not isinstance(value, (int, float, np.integer, np.floating)):
+        raise ValueError(f"expression {expr!r} does not evaluate to a number")
+    return value
 
 
 def generate_matrix_function(expr_matrix_str):
@@ -114,7 +149,7 @@ def generate_tensor(expr_matrix_str):
     """Evaluate a dict of matrix expression strings (no delta dependency) into numpy arrays."""
     expr_dict = ast.literal_eval(expr_matrix_str)
     return {
-        key: np.array([[float(cell) for cell in row] for row in expr_dict[key]])
+        key: np.array([[safe_eval(cell) for cell in row] for row in expr_dict[key]])
         for key in expr_dict
     }
 
diff --git a/src/solve_form.py b/src/solve_form.py
new file mode 100644
index 0000000..4dc08be
--- /dev/null
+++ b/src/solve_form.py
@@ -0,0 +1,318 @@
+"""Solve a scenario program described by the web form's fields.
+
+Used by the Flask app (``app.py``) and by the MOSEK licence subprocess
+(``src/mosek_solve.py``), so that both paths compute and report results in
+exactly the same way.
+"""
+import json
+import time
+
+import numpy as np
+
+from src.LP import solve_lp
+from src.QP import solve_qp
+from src.SDP import solve_sdp
+from src.Risk import quantify_risk
+from src.parsing import (generate_matrix_function, generate_matrix,
+                         generate_tensor_function, generate_tensor,
+                         generate_numeric_A_b, generate_numeric_F)
+
+
+def parse_theta_bar(raw):
+    """Parse the x̄ (reference point) form field.
+
+    The field is populated by the matrix-grid modal as a JSON-encoded 2-D
+    array (e.g. ``[["10"]]`` or ``[[0], [1]]``). Older inputs or an empty
+    field may arrive as an empty/comma-separated string. This helper
+    accepts any of those and returns a 1-D numpy float array, or 0.0 when
+    the field is unset (backward compat with the previous default).
+    """
+    if raw is None:
+        return 0.0
+    s = raw.strip()
+    if s in ('', 'null', '""', '[]'):
+        return 0.0
+    # Try JSON first (the standard path for matrix-modal entries).
+    try:
+        parsed = json.loads(s)
+        return np.asarray(parsed, dtype=float).reshape(-1)
+    except (ValueError, TypeError):
+        pass
+    # Fallback: legacy comma-separated scalars (e.g. "0, 0").
+    try:
+        return np.array([float(x) for x in s.split(',') if x.strip() != ''])
+    except ValueError as e:
+        raise ValueError(
+            f"Could not parse x̄ (theta_bar) value {raw!r}: expected a JSON "
+            f"array (e.g. [[10]]) or comma-separated numbers. ({e})"
+        )
+
+
+def compute_base_cost(active_tab, x, c, Q):
+    """Compute c'x (LP) or c'x + ½x'Qx (QP/SDP) from the optimal x.
+
+    Regularization (τ·‖x − x̄‖) and relaxation (ρ·Σζ_i) penalty terms are
+    numerical tools for informing the solution; they are not part of the
+    underlying problem's objective and must not leak into reported cost.
+    """
+    if x is None:
+        return None
+    x_arr = np.asarray(x, dtype=float).ravel()
+    if x_arr.size == 0:
+        return None
+    c_arr = np.asarray(c, dtype=float).ravel() if c is not None and getattr(c, 'size', 0) else None
+    if c_arr is None or c_arr.size != x_arr.size:
+        return None
+    base = float(c_arr @ x_arr)
+    if active_tab in ('qp-tab', 'sdp-tab') and Q is not None and getattr(Q, 'size', 0):
+        Q_arr = np.asarray(Q, dtype=float)
+        if Q_arr.shape == (x_arr.size, x_arr.size):
+            base += 0.5 * float(x_arr @ Q_arr @ x_arr)
+    return base
+
+
+def parse_confidence(raw):
+    """Parse the confidence parameter β, which must lie strictly between 0 and 1."""
+    try:
+        beta = float(raw)
+    except (TypeError, ValueError):
+        raise ValueError("Enter the confidence parameter β, a number between 0 and 1 (e.g. 1e-06).")
+    if not 0 < beta < 1:
+        raise ValueError(f"The confidence parameter β must lie strictly between 0 and 1 (got {raw}).")
+    return beta
+
+
+def parse_norm(raw):
+    """Parse the norm order p: a number p ≥ 1, ``inf`` or ``fro``. Empty means p = 2."""
+    p_raw = (raw or '').strip().strip("'\"")
+    if not p_raw:
+        return 2
+    if p_raw.lower() in ('inf', 'infinity', 'np.inf'):
+        return np.inf
+    if p_raw.lower() == 'fro':
+        return 'fro'
+    try:
+        p = float(p_raw)
+    except ValueError:
+        raise ValueError(
+            f"p (norm order) must be a single number, 'inf', or 'fro'. "
+            f"You entered {p_raw!r} — comma-separated lists are not "
+            f"supported for p (sweep over τ or ρ instead)."
+        )
+    if p < 1:
+        raise ValueError(f"p (norm order) must be at least 1 (got {p_raw}).")
+    return p
+
+
+def solve_form(form, scenarios, solver=None):
+    """Solve every (τ, ρ) combination described by the web form.
+
+    Parameters
+    ----------
+    form : dict
+        The web form's fields as strings (``active_tab``, ``mode``, the
+        matrices as JSON strings, ``rho``, ``tau``, ``confidence``, ...).
+    scenarios : numpy.ndarray or None
+        The uploaded scenarios, one row per scenario.
+    solver : str or None
+        Solver to use. ``None`` uses the form's ``solver`` field
+        (default ``CLARABEL``); the MOSEK subprocess passes ``"MOSEK"``.
+
+    Returns
+    -------
+    dict
+        The result lists shown in the results table, one entry per run.
+    """
+    active_tab = form.get('active_tab')
+
+    mode = form.get('mode', 'symbolic')
+    n_x = int(form.get('n_x', 0)) if form.get('n_x') else 0
+    rows_A = int(form.get('rows_A', 0)) if form.get('rows_A') else 0
+    lmi_size = int(form.get('lmi_size', 0)) if form.get('lmi_size') else 0
+
+    A_d = None
+    b_d = None
+    F_d = None
+
+    if mode == 'numeric':
+        if scenarios is None:
+            raise ValueError("Numeric mode requires an uploaded scenario data file.")
+        if active_tab in ('lp-tab', 'qp-tab'):
+            if n_x <= 0 or rows_A <= 0:
+                raise ValueError("Numeric mode requires valid decision variable count (n_x) and soft constraint rows (rows_A).")
+            expected_cols = rows_A * n_x + rows_A
+            if scenarios.shape[1] != expected_cols:
+                raise ValueError(
+                    f"Each scenario row must have {expected_cols} elements "
+                    f"(rows_A*n_x + rows_A = {rows_A}*{n_x} + {rows_A}), "
+                    f"but got {scenarios.shape[1]}.")
+            A_d, b_d = generate_numeric_A_b(rows_A, n_x)
+        elif active_tab == 'sdp-tab':
+            if n_x <= 0 or lmi_size <= 0:
+                raise ValueError("Numeric mode requires valid decision variable count (n_x) and LMI matrix size.")
+            expected_cols = (n_x + 1) * lmi_size ** 2
+            if scenarios.shape[1] != expected_cols:
+                raise ValueError(
+                    f"Each scenario row must have {expected_cols} elements "
+                    f"((n_x+1)*lmi_size² = {n_x + 1}*{lmi_size}² = {expected_cols}), "
+                    f"but got {scenarios.shape[1]}.")
+            F_d = generate_numeric_F(lmi_size, n_x)
+    else:
+        if form.get('A_d'):
+            A_d = generate_matrix_function(form.get('A_d'))
+        elif active_tab in ('lp-tab', 'qp-tab'):
+            raise ValueError("\\(A(\\delta)\\) is ill-defined")
+        if form.get('b_d'):
+            b_d = generate_matrix_function(form.get('b_d'))
+        elif active_tab in ('lp-tab', 'qp-tab'):
+            raise ValueError("\\(b(\\delta)\\) is ill-defined")
+        if form.get('F_d'):
+            F_d = generate_tensor_function(form.get('F_d'))
+        elif active_tab == 'sdp-tab':
+            raise ValueError("\\(F_j(\\delta)\\) is ill-defined")
+
+    # Hard LP/QP constraints use G, h (form fields renamed from legacy A, b).
+    G = generate_matrix(form.get('G')) if form.get('G') else np.array([])
+    h_vec = generate_matrix(form.get('h')) if form.get('h') else np.array([])
+    c = generate_matrix(form.get('c')) if form.get('c') else np.array([])
+    Q = generate_matrix(form.get('Q')) if form.get('Q') else np.array([])
+    # Hard SDP LMI collection uses E (form field renamed from legacy F).
+    E = generate_tensor(form.get('E')) if form.get('E') else {}
+
+    conf = parse_confidence(form.get('confidence'))
+
+    rhos = np.array([float(x) for x in form.get('rho').split(',')]) if form.get('rho') else np.array([0.0])
+    taus = np.array([float(x) for x in form.get('tau').split(',')]) if form.get('tau') else np.array([0.0])
+
+    # The selected formulation decides which parameters apply. The fields of
+    # the other formulations are hidden in the form but still submitted, so
+    # they are ignored here rather than silently applied or swept over.
+    selected_option = form.get('option')
+    if selected_option in ('robust', 'regularization'):
+        rhos = np.array([0.0])
+    if selected_option in ('robust', 'relaxation'):
+        taus = np.array([0.0])
+
+    theta_bar = parse_theta_bar(form.get('theta_bar'))
+    p = parse_norm(form.get('p'))
+
+    # Hard constraints counted in "Total Constraints": the rows of G for
+    # LP/QP, or the single hard LMI for SDP.
+    if active_tab == 'sdp-tab':
+        n_hard = 1 if E else 0
+    else:
+        n_hard = G.shape[0] if (G.size and h_vec.size) else 0
+
+    # Initialize lists to collect results for each run
+    optimal_x_list = []
+    optimal_s_list = []
+    optimal_cost_list = []
+    N_list = []
+    active_list = []
+    constraints_list = []
+    risk_list = []
+    e_list = []
+    degeneracy_list = []
+    rho_list = []
+    tau_list = []
+    solve_time_list = []
+    risk_time_list = []
+
+    # Overall confidence after q sweep attempts is 1 − β·q (union bound).
+    # Pass the user-entered β directly to the per-attempt risk computation;
+    # only the *reported* confidence folds in q.
+    n_attempts = len(taus) * len(rhos)
+    conf_reported = conf * n_attempts
+
+    # The relaxation formulation must keep its slack variables ζ_i even when
+    # the user has set ρ = 0 (treats that case as genuinely unbounded rather
+    # than silently coercing it into the robust LP). When `option` is not in
+    # the form payload (older clients / direct API callers), fall through to
+    # the legacy "infer slack from ρ" behaviour by leaving include_slack=None
+    # so each solver default kicks in.
+    if selected_option:
+        include_slack = selected_option in ('relaxation', 'regularization-relaxation')
+    else:
+        include_slack = None
+
+    run_solver = solver or form.get('solver', 'CLARABEL')
+
+    for tau in taus:
+        for rho in rhos:
+            optimal_x = np.array([])
+            optimal_s = np.array([])
+            optimal_cost = None
+            N = 0
+            complexity = []
+            constraints = []
+            degeneracy = False
+            risk = np.array([])
+            e = "None"
+            solve_time = 0.0
+            risk_time = 0.0
+
+            try:
+                t0 = time.perf_counter()
+                if active_tab == 'lp-tab':
+                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_lp(
+                        scenarios, A_d, b_d, G, h_vec, c, tau, theta_bar, rho, p, run_solver,
+                        include_slack=include_slack)
+                elif active_tab == 'qp-tab':
+                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_qp(
+                        scenarios, A_d, b_d, G, h_vec, c, Q, tau, theta_bar, rho, p, run_solver,
+                        include_slack=include_slack)
+                elif active_tab == 'sdp-tab':
+                    optimal_x, optimal_s, optimal_cost, N, complexity, constraints, degeneracy = solve_sdp(
+                        scenarios, F_d, E, c, Q, tau, theta_bar, rho, p, run_solver,
+                        include_slack=include_slack)
+                solve_time = time.perf_counter() - t0
+
+                # Replace full-objective cost (which includes τ·‖x−x̄‖_p and
+                # ρ·Σζ_i penalty terms used only as solver guidance) with the
+                # base cost: c'x for LP, c'x + ½x'Qx for QP/SDP.
+                optimal_cost = compute_base_cost(active_tab, optimal_x, c, Q)
+
+                t1 = time.perf_counter()
+                risk = np.array(quantify_risk(complexity, N, conf))
+                risk_time = time.perf_counter() - t1
+            except Exception as error:
+                e = str(error)
+                # Shorten verbose MOSEK license errors
+                if 'license' in e.lower() and ('mosek' in e.lower() or 'rescode' in e.lower()):
+                    e = "MOSEK license cannot be located or is incorrect."
+                print(e)
+
+            # Append results for this run
+            optimal_x_list.append(optimal_x.tolist() if hasattr(optimal_x, 'tolist') else optimal_x)
+            optimal_s_list.append(optimal_s.tolist() if hasattr(optimal_s, 'tolist') else optimal_s)
+            optimal_cost_list.append(optimal_cost)
+            N_list.append(N)
+            active_list.append(complexity)
+            e_list.append(e)
+            constraints_list.append(len(constraints) + n_hard if len(constraints) else 0)
+            risk_list.append(risk.tolist() if hasattr(risk, 'tolist') else risk)
+            degeneracy_list.append(degeneracy)
+            rho_list.append(float(rho))
+            tau_list.append(float(tau))
+            solve_time_list.append(round(solve_time, 4))
+            risk_time_list.append(round(risk_time, 4))
+
+    # Prepare the result dictionary with lists for each parameter
+    return {
+        "form_data": dict(form),
+        "optimal_x": optimal_x_list,
+        "optimal_s": optimal_s_list,
+        "optimal_cost": optimal_cost_list,
+        "num_deltas": N_list,
+        "tot_con": constraints_list,
+        "active_con": active_list,
+        "risk": risk_list,
+        "conf": 1 - conf_reported,
+        "tau_": tau_list,
+        "rho_": rho_list,
+        "errorcode": e_list,
+        "degeneracy": degeneracy_list,
+        "solve_time": solve_time_list,
+        "risk_time": risk_time_list,
+        "n_attempts": n_attempts,
+    }
diff --git a/static/js/main.js b/static/js/main.js
index b35150b..0141106 100644
--- a/static/js/main.js
+++ b/static/js/main.js
@@ -232,9 +232,9 @@ function openMatrixModal(tab, matrix) {
 
 function getDynamicTip(matrix) {
     if (matrix === 'A_d') {
-        return "Manually enter the matrix. Cells accept SymPy expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, **.";
+        return "Manually enter the matrix. Cells accept expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, ** (or ^) and functions such as <code>sin()</code>, <code>exp()</code>, <code>sqrt()</code>.";
     } else if (matrix === 'b_d') {
-        return "Manually enter the vector. Cells accept SymPy expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, **.";
+        return "Manually enter the vector. Cells accept expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, ** (or ^) and functions such as <code>sin()</code>, <code>exp()</code>, <code>sqrt()</code>.";
     } else if (matrix === 'Q') {
         return "Manually enter the matrix. \\(Q\\) must be <b>symmetric</b> and <b>positive semidefinite</b> (\\(Q \\succeq 0\\)).";
     } else if (matrix === 'c') {
@@ -244,11 +244,11 @@ function getDynamicTip(matrix) {
     } else if (matrix === 'h') {
         return "Manually enter the hard constraint vector.";
     } else if (matrix === 'F_d') {
-        return "Manually enter the matrices. Each \\(F_j(\\delta)\\) must be <b>symmetric</b>. Cells accept SymPy expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, **.";
+        return "Manually enter the matrices. Each \\(F_j(\\delta)\\) must be <b>symmetric</b>. Cells accept expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, ** (or ^) and functions such as <code>sin()</code>, <code>exp()</code>, <code>sqrt()</code>.";
     } else if (matrix === 'E') {
         return "Manually enter the matrices. Each \\(E_j\\) must be <b>symmetric</b>.";
     } else if (matrix === 'A_da') {
-        return "Manually enter the matrices. Each \\(A_j(\\delta)\\) must be <b>symmetric</b>. Cells accept SymPy expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, **.";
+        return "Manually enter the matrices. Each \\(A_j(\\delta)\\) must be <b>symmetric</b>. Cells accept expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, ** (or ^) and functions such as <code>sin()</code>, <code>exp()</code>, <code>sqrt()</code>.";
     } else if (matrix === 'A_a') {
         return "Manually enter the matrices. Each \\(G_j\\) must be <b>symmetric</b>.";
     } else if (matrix === 'C') {
@@ -258,7 +258,7 @@ function getDynamicTip(matrix) {
     } else if (matrix === 'theta-bar') {
         return "Manually enter the reference vector.";
     } else if (matrix === 'b_da') {
-        return "Manually enter the vectors. Cells accept SymPy expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, **.";
+        return "Manually enter the vectors. Cells accept expressions with <code>delta[k]</code> (delta[0] = 1st component of vector \\(\\delta\\), delta[1] = 2nd component of vector \\(\\delta\\), ...) including +, -, *, /, ** (or ^) and functions such as <code>sin()</code>, <code>exp()</code>, <code>sqrt()</code>.";
     } else if (matrix === 'b_a') {
         return "Manually enter the vectors.";
     }
@@ -294,9 +294,9 @@ function updateFormatReference(matrix) {
     // ── Cell content rules ──
     let html = '<h6><b>Cell content:</b></h6>';
     if (isDelta) {
-        html += '<p>This matrix is <b>delta-capable</b> — cells accept <b>numeric values</b> <i>or</i> SymPy expressions that reference ' +
+        html += '<p>This matrix is <b>delta-capable</b> — cells accept <b>numeric values</b> <i>or</i> expressions that reference ' +
                 '<code>delta[k]</code> (<code>delta[0]</code> = 1st component of \\(\\delta\\), <code>delta[1]</code> = 2nd component of \\(\\delta\\), \\(\\ldots\\)). ' +
-                'Operators <code>+ - * / **</code> and functions like <code>sin()</code>, <code>exp()</code> are allowed. ' +
+                'Operators <code>+ - * / **</code> and functions such as <code>sin()</code>, <code>cos()</code>, <code>exp()</code>, <code>log()</code>, <code>sqrt()</code>, <code>abs()</code> and the constants <code>pi</code>, <code>e</code> are allowed (<code>^</code> also means power). ' +
                 'Example cells: <code>delta[0] + 1</code>, <code>sin(delta[1])</code>, <code>2.5</code>.</p>';
         html += '<p><b>Mode note:</b> this matrix is only editable when the <b>Symbolic</b> input mode is selected. ' +
                 'In <b>Numeric</b> mode, \\(A(\\delta_i)\\) and \\(b(\\delta_i)\\) come directly from the uploaded scenario data ' +
@@ -372,9 +372,9 @@ function updateSDPFormatReference(matrix) {
 
     html += '<h6><b>Cell content:</b></h6>';
     if (isSoft) {
-        html += '<p>Each matrix is <b>delta-capable</b> — cells accept <b>numeric values</b> <i>or</i> SymPy expressions that reference ' +
+        html += '<p>Each matrix is <b>delta-capable</b> — cells accept <b>numeric values</b> <i>or</i> expressions that reference ' +
                 '<code>delta[k]</code> (<code>delta[0]</code> = 1st component of \\(\\delta\\), <code>delta[1]</code> = 2nd component of \\(\\delta\\), \\(\\ldots\\)). ' +
-                'Operators <code>+ - * / **</code> and functions like <code>sin()</code>, <code>exp()</code> are allowed. ' +
+                'Operators <code>+ - * / **</code> and functions such as <code>sin()</code>, <code>cos()</code>, <code>exp()</code>, <code>log()</code>, <code>sqrt()</code>, <code>abs()</code> and the constants <code>pi</code>, <code>e</code> are allowed (<code>^</code> also means power). ' +
                 'Every matrix must be <b>symmetric</b>.</p>';
         html += '<p><b>Mode note:</b> scenario LMI matrices \\(F_j(\\delta)\\) are only editable when the <b>Symbolic</b> input mode is selected. ' +
                 'In <b>Numeric</b> mode, \\(F_{0,i}, F_{1,i}, \\ldots, F_{d,i}\\) come directly from the uploaded scenario data ' +
@@ -431,7 +431,7 @@ function isCurrentMatrixDeltaCapable() {
 /**
  * Validate a single cell's text against the delta-capable rule:
  *   • numeric-only → must parse as a finite number
- *   • delta-capable → either a finite number, or a SymPy-ish expression
+ *   • delta-capable → either a finite number, or an expression in delta[k]
  * Returns '' if valid, otherwise a short reason.
  */
 function validateCellText(text, deltaCapable) {
@@ -1218,6 +1218,16 @@ function solve() {
     const check = validateSolveInputs();
     if (!check.valid) return;
 
+    const beta = Number((document.getElementById('confidence') || {}).value);
+    if (!(beta > 0 && beta < 1)) {
+        const errEl = document.getElementById('solve-error');
+        if (errEl) {
+            errEl.textContent = 'Enter the confidence parameter β, a number between 0 and 1 (e.g. 1e-06).';
+            errEl.style.display = 'block';
+        }
+        return;
+    }
+
     const selectedSolver = document.getElementById('solver').value;
 
     if (selectedSolver === 'MOSEK') {
@@ -1295,6 +1305,15 @@ function executeSolve(mosekLicenseFile) {
 
     document.getElementById('result-box').innerHTML = '<p>Solve button pressed, processing results...</p>';
 
+    // On the hosted app, suggest a local install once a solve runs for a while.
+    const isLocalHost = ['localhost', '127.0.0.1', '[::1]', ''].includes(window.location.hostname);
+    const slowSolveTimer = isLocalHost ? null : setTimeout(function() {
+        const box = document.getElementById('result-box');
+        if (box) box.insertAdjacentHTML('beforeend',
+            '<div class="alert alert-info mt-2">This problem is taking a while: problems of this size are better solved locally. ' +
+            'See the <a href="https://github.com/Kiguli/Scen-Opt/wiki/Local_Install" target="_blank" rel="noopener">local install guide</a>.</div>');
+    }, 30000);
+
     fetch(activeForm.action, {
         method: activeForm.method,
         body: formData
@@ -1302,7 +1321,12 @@ function executeSolve(mosekLicenseFile) {
         .then(response => {
             if (!response.ok) {
                 return response.text().then(text => {
-                    throw new Error('Server error (' + response.status + '): ' + text);
+                    let message = 'Server error (' + response.status + '): ' + text;
+                    try {
+                        const body = JSON.parse(text);
+                        if (body && body.error) message = body.error;
+                    } catch (e) { /* not JSON: keep the raw text */ }
+                    throw new Error(message);
                 });
             }
             return response.json();
@@ -1379,7 +1403,11 @@ function executeSolve(mosekLicenseFile) {
             $('#result-box .result-help-btn').popover({ trigger: 'hover', html: true });
 
             // Prepare risk lower and upper arrays
-            const riskLower = Array.isArray(data.risk) ? data.risk.map(r => Array.isArray(r) ? r[0] : r) : [data.risk];
+            // A degenerate run's lower bound is not certified, so it is left off the graphs (null = gap).
+            const degenerateRuns = Array.isArray(data.degeneracy) ? data.degeneracy : [data.degeneracy];
+            const riskLower = Array.isArray(data.risk)
+                ? data.risk.map((r, i) => degenerateRuns[i] ? null : (Array.isArray(r) ? r[0] : r))
+                : [degenerateRuns[0] ? null : data.risk];
             const riskUpper = Array.isArray(data.risk) ? data.risk.map(r => Array.isArray(r) ? r[1] : r) : [data.risk];
 
             // Assume data.rho_ and data.tau_ are arrays of equal length, and each (rho, tau) pair is unique
@@ -1581,7 +1609,8 @@ function executeSolve(mosekLicenseFile) {
             console.error('Error:', error);
             document.getElementById('result-box').innerHTML =
                 '<div class="alert alert-danger"><strong>Error:</strong> ' + error.message + '</div>';
-        });
+        })
+        .finally(() => clearTimeout(slowSolveTimer));
 }
 
 function generateResultTable(data) {
@@ -1646,12 +1675,18 @@ function generateResultTable(data) {
             if (opt === 'regularization-relaxation') return 'Confidence \\(1 - \\beta \\cdot n_{\\tau} \\cdot n_{\\rho}\\)';
             return 'Confidence \\(1 - \\beta\\)';
         })(), values: data.conf},
-        {label: 'Risk Bounds &epsilon;', values: data.risk},
+        {label: `Risk Bounds &epsilon; <button type="button" class="btn btn-link p-0 ml-1 dim-help-btn result-help-btn"
+                data-toggle="popover" data-trigger="hover" data-placement="right" data-html="true"
+                data-content="Lower and upper bounds on the risk, each holding with the confidence above. The upper bound always holds. The lower bound needs non-degeneracy: when degeneracy is detected it is not certified and is shown as N/A."><span class="dim-help-icon">?</span></button>`, values: data.risk},
         {label: 'Degeneracy Detected?', values: data.degeneracy},
         {label: 'Complexity (support list size)', values: data.active_con},
         {label: 'Number of data samples', values: data.num_deltas},
-        {label: 'Total Constraints', values: data.tot_con},
-        {label: 'Optimization Time (s)', values: data.solve_time},
+        {label: `Total Constraints <button type="button" class="btn btn-link p-0 ml-1 dim-help-btn result-help-btn"
+                data-toggle="popover" data-trigger="hover" data-placement="right" data-html="true"
+                data-content="The \\(N\\) scenario constraints plus the hard constraints (the rows of \\(G\\), or the hard LMI)."><span class="dim-help-icon">?</span></button>`, values: data.tot_con},
+        {label: `Optimization Time (s) <button type="button" class="btn btn-link p-0 ml-1 dim-help-btn result-help-btn"
+                data-toggle="popover" data-trigger="hover" data-placement="right" data-html="true"
+                data-content="Time to solve the program and compute its support list (which re-solves the program)."><span class="dim-help-icon">?</span></button>`, values: data.solve_time},
         {label: 'Risk Computation Time (s)', values: data.risk_time},
         {label: '<i>Error</i>', values: data.errorcode}
     ];
@@ -1690,14 +1725,18 @@ function generateResultTable(data) {
             }
 
             if (row.label.includes('Risk Bounds') && Array.isArray(val) && val.length === 2) {
-                val = `<strong>[${fmt(val[0])}, ${fmt(val[1])}]</strong>`;
+                const degenerate = getValue(data.degeneracy, i) === true;
+                val = `<strong>[${degenerate ? 'N/A' : fmt(val[0])}, ${fmt(val[1])}]</strong>`;
             } else if (row.label.includes('Confidence')) {
                 // Clamp at 0 — over-aggressive sweeps can push 1 − β·q below
                 // zero, which is nonsensical to display. Reuse fmt() so a
                 // tiny but positive confidence still prints in scientific
                 // form rather than rounding to 0.0000.
-                const v = Number(val);
-                val = Number.isFinite(v) ? fmt(Math.max(0, v)) : val;
+                // Show enough digits that 1 − β is not rounded to 1.0000.
+                const v = Math.max(0, Number(val));
+                val = Number.isFinite(v)
+                    ? (v > 0.9999 && v < 1 ? String(Number(v.toPrecision(12))) : fmt(v))
+                    : val;
             } else if (row.label.includes('Optimal Cost')) {
                 val = `<strong>${fmt(val)}</strong>`;
             } else if (row.label.includes('Error')) {
diff --git a/take_screenshots.py b/take_screenshots.py
index 7121a36..1bf85f5 100644
--- a/take_screenshots.py
+++ b/take_screenshots.py
@@ -1,12 +1,24 @@
 #!/usr/bin/env python3
-"""Take screenshots of the web interface for documentation."""
+"""Take screenshots of the web interface for documentation.
+
+Requires playwright (not in requirements.txt):
+    pip install playwright && playwright install chromium
+
+Start the app first (python app.py), then run from the repository root:
+    python take_screenshots.py [URL]     # default URL: http://127.0.0.1:5000/
+The screenshots are written to screenshots/.
+"""
 import asyncio
+import os
+import sys
 from playwright.async_api import async_playwright
 
-SCREENSHOTS_DIR = "/Users/ben/Documents/Tools/ScenarioApproachTool/screenshots"
+ROOT = os.path.dirname(os.path.abspath(__file__))
+SCREENSHOTS_DIR = os.path.join(ROOT, "screenshots")
+URL = sys.argv[1] if len(sys.argv) > 1 else "http://127.0.0.1:5000/"
+BENCHMARK_DATA = os.path.join(ROOT, "benchmarks", "LP_half_width_2d", "data")
 
 async def main():
-    import os
     os.makedirs(SCREENSHOTS_DIR, exist_ok=True)
 
     async with async_playwright() as p:
@@ -16,17 +28,17 @@ async def main():
         # Handle alert dialogs automatically
         page.on("dialog", lambda d: asyncio.ensure_future(d.accept()))
 
-        await page.goto("http://127.0.0.1:5000/")
+        await page.goto(URL)
         await page.wait_for_load_state("networkidle")
         await asyncio.sleep(2)
 
-        # --- Load benchmark via Detect Program modal ---
+        # --- Load the benchmark via the Upload Program dialog ---
         await page.click('button[data-target="#detectProgramModal"]')
         await asyncio.sleep(1)
 
         # Upload the benchmark JSON
         await page.set_input_files('#detect-program-file',
-            "benchmarks/LP_half_width_2d/data/benchmark.json")
+            os.path.join(BENCHMARK_DATA, "program_symbolic.json"))
         await asyncio.sleep(0.5)
 
         # Click "Load" button
@@ -35,16 +47,17 @@ async def main():
 
         # Set multiple rho and tau values for interesting graphs
         # First select the general formulation to show rho+tau fields
-        await page.select_option('#lp-options', 'robust-regularization-relaxation')
+        await page.select_option('#lp-options', 'regularization-relaxation')
         await asyncio.sleep(0.5)
 
-        await page.fill('#lp-rho', '0,0.25,0.5,0.75,1')
+        # (rho = 0 is left out: relaxation with rho = 0 is unbounded.)
+        await page.fill('#lp-rho', '0.25,0.5,0.75,1')
         await page.fill('#lp-tau', '0,0.01,0.1,1')
         await asyncio.sleep(0.3)
 
         # Upload scenarios file
         await page.set_input_files('#file',
-            "benchmarks/LP_half_width_2d/data/scenarios.csv")
+            os.path.join(BENCHMARK_DATA, "scenarios.csv"))
         await asyncio.sleep(0.5)
 
         # Screenshot 1: Tool (everything above results, including Solve button)
@@ -57,7 +70,7 @@ async def main():
         result_box = await result_card.bounding_box()
 
         await page.screenshot(
-            path=f"{SCREENSHOTS_DIR}/1_tool.png",
+            path=f"{SCREENSHOTS_DIR}/scrn_1_tool.png",
             clip={"x": 0, "y": 0, "width": 1400, "height": result_box["y"] - 5},
             full_page=True
         )
@@ -76,21 +89,21 @@ async def main():
         result_card_el = page.locator('.card.mt-4:has(.card-title:has-text("Result"))').first
         await result_card_el.scroll_into_view_if_needed()
         await asyncio.sleep(0.5)
-        await result_card_el.screenshot(path=f"{SCREENSHOTS_DIR}/2_results.png")
+        await result_card_el.screenshot(path=f"{SCREENSHOTS_DIR}/scrn_2_results.png")
         print("Screenshot 2: Results table - done")
 
         # Screenshot 3: tau graph tab
         await page.click('a#result-graph2-tab')
         await asyncio.sleep(1.5)
         tau_pane = page.locator('#result-graph2')
-        await tau_pane.screenshot(path=f"{SCREENSHOTS_DIR}/3_tau_graph.png")
+        await tau_pane.screenshot(path=f"{SCREENSHOTS_DIR}/scrn_3_tau_graph.png")
         print("Screenshot 3: Tau graph - done")
 
         # Screenshot 4: rho graph tab
         await page.click('a#result-graph-tab')
         await asyncio.sleep(1.5)
         rho_pane = page.locator('#result-graph')
-        await rho_pane.screenshot(path=f"{SCREENSHOTS_DIR}/4_rho_graph.png")
+        await rho_pane.screenshot(path=f"{SCREENSHOTS_DIR}/scrn_4_rho_graph.png")
         print("Screenshot 4: Rho graph - done")
 
         await browser.close()
diff --git a/tests/test_app.py b/tests/test_app.py
index acfd1c0..8d89c87 100644
--- a/tests/test_app.py
+++ b/tests/test_app.py
@@ -141,3 +141,64 @@ def test_download_mat_route(client):
 
     assert resp.status_code == 200
     assert "matlab" in resp.content_type or "octet-stream" in resp.content_type
+
+
+def _interval_form(**overrides):
+    """Form fields for the smallest-enclosing-interval LP, as the browser sends them."""
+    csv_bytes = "\n".join(str(v) for v in [0.1, 0.4, 0.9]).encode()
+    data = {
+        "active_tab": "lp-tab",
+        "A_d": '[["-1", "-1"], ["1", "-1"]]',
+        "b_d": '[["delta[0]"], ["-delta[0]"]]',
+        "c": '[["0"], ["1"]]',
+        "rho": "0",
+        "tau": "0",
+        "option": "robust",
+        "solver": "CLARABEL",
+        "confidence": "1e-06",
+        "file": (io.BytesIO(csv_bytes), "scenarios.csv"),
+    }
+    data.update(overrides)
+    return data
+
+
+def test_solve_accepts_plain_numbers(client):
+    """Matrices with plain numbers (not strings) solve like their string form."""
+    resp = client.post("/solve", data=_interval_form(c="[[0], [1]]", A_d='[[-1, -1], [1, -1]]'),
+                       content_type="multipart/form-data")
+    assert resp.status_code == 200
+    np.testing.assert_allclose(np.ravel(resp.get_json()["optimal_x"][0]), [0.5, 0.4], atol=1e-6)
+
+
+@pytest.mark.parametrize("beta", ["", "0", "1", "abc"])
+def test_solve_rejects_invalid_beta(client, beta):
+    """A missing or out-of-range β is reported as an error, not treated as 0."""
+    resp = client.post("/solve", data=_interval_form(confidence=beta), content_type="multipart/form-data")
+    assert resp.status_code == 400
+    assert "confidence parameter" in resp.get_json()["error"]
+
+
+def test_solve_robust_ignores_hidden_parameters(client):
+    """With the robust option, hidden τ and ρ fields are neither applied nor swept."""
+    resp = client.post("/solve", data=_interval_form(tau="5", rho="0.1, 1"), content_type="multipart/form-data")
+    data = resp.get_json()
+    assert data["n_attempts"] == 1 and data["tau_"] == [0.0] and data["rho_"] == [0.0]
+    np.testing.assert_allclose(np.ravel(data["optimal_x"][0]), [0.5, 0.4], atol=1e-6)
+
+
+def test_solve_sweep_confidence_and_total_constraints(client):
+    """A sweep reports 1 − β·n_runs; Total Constraints counts the hard constraints too."""
+    resp = client.post("/solve", data=_interval_form(option="relaxation", rho="0.5, 1", confidence="0.01",
+                                                     G='[["0", "-1"]]', h='[["0"]]'),
+                       content_type="multipart/form-data")
+    data = resp.get_json()
+    assert data["n_attempts"] == 2
+    assert abs(data["conf"] - 0.98) < 1e-12
+    assert data["tot_con"] == [4, 4]  # 3 scenario constraints + 1 row of G
+
+
+def test_download_mat_with_a_failed_run(client):
+    """A run without results (None) no longer breaks the .MAT download."""
+    resp = client.post("/download-mat", json={"optimal_cost": [0.5, None], "risk": [[0.0, 0.2], []]})
+    assert resp.status_code == 200
+    assert resp.data[:6] == b"MATLAB"
diff --git a/tests/test_lp.py b/tests/test_lp.py
index 1693116..acbf560 100644
--- a/tests/test_lp.py
+++ b/tests/test_lp.py
@@ -131,3 +131,52 @@ def test_lp_infeasible():
             c=np.array([[-5], [-3]]),
             tau=0.0, x_ref=np.array([0, 0]), rho=0.0, norm_type=2, solver="SCS",
         )
+
+
+# ── Interval LP (smallest enclosing interval) used by the regression tests below ──
+_INTERVAL = np.array([0.1, 0.4, 0.9])
+
+
+def _interval_A(delta):
+    return np.array([[-1.0, -1.0], [1.0, -1.0]])
+
+
+def test_lp_empty_deltas_raises_value_error():
+    """No scenarios gives a clear ValueError (not an IndexError)."""
+    with pytest.raises(ValueError, match="at least one row"):
+        solve_lp(np.empty((0, 1)), _interval_A, lambda d: np.array([[d[0]], [-d[0]]]),
+                 np.array([]), np.array([]), np.array([[0.0], [1.0]]), solver="CLARABEL")
+
+
+def test_lp_one_dimensional_b_is_not_broadcast():
+    """A 1-D b(δ) of shape (m,) gives the same solution as a (m, 1) column."""
+    deltas = _INTERVAL.reshape(-1, 1)
+    c = np.array([[0.0], [1.0]])
+    x_col, *_ = solve_lp(deltas, _interval_A, lambda d: np.array([[d[0]], [-d[0]]]),
+                         np.array([]), np.array([]), c, solver="CLARABEL")
+    x_1d, *_ = solve_lp(deltas, _interval_A, lambda d: np.array([d[0], -d[0]]),
+                        np.array([]), np.array([]), c, solver="CLARABEL")
+    np.testing.assert_allclose(x_1d, x_col, atol=1e-6)
+    np.testing.assert_allclose(x_col.ravel(), [0.5, 0.4], atol=1e-6)
+
+
+def test_lp_without_regularization_is_a_pure_lp():
+    """With τ = 0 no norm term is added, so an LP-only solver can solve it."""
+    import cvxpy as cp
+    if "SCIPY" not in cp.installed_solvers():
+        pytest.skip("SCIPY not installed")
+    x, zeta, cost, N, k, cons, degenerate = solve_lp(
+        _INTERVAL.reshape(-1, 1), _interval_A, lambda d: np.array([[d[0]], [-d[0]]]),
+        np.array([]), np.array([]), np.array([[0.0], [1.0]]), solver="SCIPY")
+    np.testing.assert_allclose(x.ravel(), [0.5, 0.4], atol=1e-6)
+    assert k == 2 and not degenerate
+
+
+@pytest.mark.parametrize("p", [1, 1.5, 2, 3, "inf", "fro"])
+def test_lp_regularization_accepts_any_vector_norm(p):
+    """The regularization term is a vector norm, so any p >= 1, inf and fro work."""
+    x, *_ = solve_lp(_INTERVAL.reshape(-1, 1), _interval_A, lambda d: np.array([[d[0]], [-d[0]]]),
+                     np.array([]), np.array([]), np.array([[0.0], [1.0]]),
+                     tau=0.01, x_ref=np.array([0.5, 0.4]), norm_type=float("inf") if p == "inf" else p,
+                     solver="CLARABEL")
+    np.testing.assert_allclose(x.ravel(), [0.5, 0.4], atol=1e-5)
diff --git a/tests/test_miscellaneous.py b/tests/test_miscellaneous.py
index 572427b..ba68901 100644
--- a/tests/test_miscellaneous.py
+++ b/tests/test_miscellaneous.py
@@ -1,4 +1,4 @@
-from src.Miscellaneous import get_solvers, get_norm_types
+from src.Miscellaneous import get_solvers, get_norm_types, SOLVER_CAPABILITIES, _cvxpy_capabilities
 
 
 def test_get_solvers():
@@ -20,4 +20,17 @@ def test_get_solvers_capabilities():
 def test_get_norm_types():
     """get_norm_types returns the expected list of norms."""
     norms = get_norm_types()
-    assert norms == [1, 2, "inf", "fro", "nuc"]
+    assert norms == [1, 2, "inf", "fro"]
+
+
+def test_solver_capabilities_match_cvxpy():
+    """Table entries never claim a program type CVXPY says the solver cannot handle."""
+    for name, types in SOLVER_CAPABILITIES.items():
+        assert set(types) <= set(_cvxpy_capabilities(name)), name
+
+
+def test_unlisted_solvers_use_cvxpy_metadata():
+    """A solver missing from the table gets its types from CVXPY, not all three."""
+    assert _cvxpy_capabilities("CUOPT") == ["LP"]
+    assert _cvxpy_capabilities("COPT") == ["LP", "QP", "SDP"]
+    assert _cvxpy_capabilities("NOT_A_SOLVER") == []
diff --git a/tests/test_numeric_mode.py b/tests/test_numeric_mode.py
index 9620422..73794d7 100644
--- a/tests/test_numeric_mode.py
+++ b/tests/test_numeric_mode.py
@@ -387,9 +387,9 @@ def test_detect_program_sdp(client):
 def test_detect_program_sdp_numeric(client):
     """Detect Program with numeric SDP benchmark loads and matches symbolic results."""
     # Load numeric benchmark
-    with open("benchmarks/SDP_LPV_stability_3d_numeric/data/program_symbolic.json") as f:
+    with open("benchmarks/SDP_LPV_stability_3d/data/program_numeric.json") as f:
         num_bm = json.load(f)
-    num_scenarios = np.loadtxt("benchmarks/SDP_LPV_stability_3d_numeric/data/scenarios.csv",
+    num_scenarios = np.loadtxt("benchmarks/SDP_LPV_stability_3d/data/scenarios_numeric.csv",
                                 delimiter=",", ndmin=2)
     num_csv = _to_csv_bytes(num_scenarios)
 
diff --git a/tests/test_parsing.py b/tests/test_parsing.py
index a3a6568..6f0fa09 100644
--- a/tests/test_parsing.py
+++ b/tests/test_parsing.py
@@ -53,3 +53,43 @@ def test_generate_tensor():
     )
     np.testing.assert_array_equal(result["0"], np.array([[1, 0], [0, 1]]))
     np.testing.assert_array_equal(result["1"], np.array([[2, 0], [0, 3]]))
+
+
+def test_numeric_entries_are_accepted():
+    """Plain numbers (from JSON or MAT files) work like their string form."""
+    np.testing.assert_allclose(generate_matrix('[[0], [1]]'), [[0.0], [1.0]])
+    f = generate_matrix_function('[[-1, 2.5], ["delta[0]", 0]]')
+    np.testing.assert_allclose(f(np.array([3.0])), [[-1.0, 2.5], [3.0, 0.0]])
+    t = generate_tensor_function('{"0": [[1, "delta[0]"]]}')
+    np.testing.assert_allclose(t(np.array([2.0]))["0"], [[1.0, 2.0]])
+
+
+def test_constant_tensor_accepts_expressions():
+    """E matrices accept constant expressions like '1/2', as other matrices do."""
+    E = generate_tensor('{"0": [["1/2", "0"], ["0", "-1"]]}')
+    np.testing.assert_allclose(E["0"], [[0.5, 0.0], [0.0, -1.0]])
+
+
+def test_scalar_scenario_is_indexable():
+    """A scalar scenario (from a 1-D scenarios array) has one component, delta[0]."""
+    f = generate_matrix_function('[["delta[0]", "-delta[0]"]]')
+    np.testing.assert_allclose(f(0.25), [[0.25, -0.25]])
+
+
+def test_sympy_style_names():
+    """SymPy-style functions, constants and ^ work without the math. prefix."""
+    from src.parsing import safe_eval
+    d = np.array([0.5, 2.0])
+    assert math.isclose(safe_eval("sin(delta[0])", d), math.sin(0.5))
+    assert math.isclose(safe_eval("delta[1]^2 + sqrt(delta[1])", d), 4 + math.sqrt(2))
+    assert math.isclose(safe_eval("Abs(-pi) + ln(E)", d), math.pi + 1)
+    assert math.isclose(safe_eval("math.exp(delta[0])", d), math.exp(0.5))  # still accepted
+
+
+def test_disallowed_names_are_rejected():
+    """Names outside the whitelist are still rejected."""
+    import pytest
+    from src.parsing import safe_eval
+    for bad in ["__import__('os')", "open('x')", "eval('1')", "sin"]:
+        with pytest.raises(ValueError):
+            safe_eval(bad, np.array([1.0]))
diff --git a/tests/test_qp.py b/tests/test_qp.py
index 1149ea3..a38b1be 100644
--- a/tests/test_qp.py
+++ b/tests/test_qp.py
@@ -102,3 +102,13 @@ def test_qp_non_symmetric_q():
             c=C_VEC, Q=Q_bad,
             tau=0.0, x_ref=np.array([0, 0]), rho=0.0, norm_type=2, solver="SCS",
         )
+
+
+def test_qp_rank_deficient_q_is_accepted():
+    """A PSD Q whose zero eigenvalues come out as about -1e-16 is accepted."""
+    D = np.array([[1.0, 2.0, 3.0], [2.0, 4.0, 6.0]])
+    Q = D.T @ D  # rank 1: two eigenvalues are zero up to rounding
+    scenarios = np.array([[0.0], [0.5]])
+    x, *_ = solve_qp(scenarios, lambda d: np.array([[-1.0, 0.0, 0.0]]), lambda d: np.array([[d[0]]]),
+                     np.array([]), np.array([]), np.zeros((3, 1)), Q, solver="CLARABEL")
+    assert x is not None
diff --git a/tests/test_risk.py b/tests/test_risk.py
index 0900ec4..01cb4fe 100644
--- a/tests/test_risk.py
+++ b/tests/test_risk.py
@@ -1,3 +1,6 @@
+import numpy as np
+import pytest
+
 from src.Risk import quantify_risk
 
 
@@ -39,3 +42,19 @@ def test_risk_bounds_high_beta():
     _, epsU_low_beta = quantify_risk(k=10, N=4000, beta=1e-9)
     _, epsU_high_beta = quantify_risk(k=10, N=4000, beta=1e-3)
     assert epsU_high_beta < epsU_low_beta
+
+
+def test_risk_invalid_inputs_raise():
+    """k > N, N = 0 and β outside (0, 1) raise ValueError."""
+    for k, N, beta in [(3, 2, 1e-6), (0, 0, 1e-6), (1, 10, 0.0), (1, 10, 1.0)]:
+        with pytest.raises(ValueError):
+            quantify_risk(k, N, beta)
+
+
+def test_risk_returns_floats_and_keeps_x64_local():
+    """Bounds are plain floats, and importing Risk does not switch JAX to 64-bit globally."""
+    import jax
+    epsL, epsU = quantify_risk(5, 5, 1e-6)
+    assert isinstance(epsL, float) and isinstance(epsU, float) and epsU == 1.0
+    assert not jax.config.jax_enable_x64
+    np.testing.assert_allclose(quantify_risk(2, 100, 1e-6)[1], 0.20852406, atol=1e-8)
diff --git a/tests/test_sdp.py b/tests/test_sdp.py
index 6a0c221..20086bb 100644
--- a/tests/test_sdp.py
+++ b/tests/test_sdp.py
@@ -102,3 +102,46 @@ def test_sdp_non_psd_q():
             Q=Q_bad,
             tau=0.0, x_ref=np.zeros(2), rho=0.0, norm_type=2, solver="SCS",
         )
+
+
+def _lpv_problem():
+    """LPV stability example (Example SDP1): x = (p11, p12, p22)."""
+    def F_d(delta):
+        d = delta[0]
+        return {'0': np.zeros((2, 2)),
+                '1': np.array([[0.0, 1.0], [1.0, 0.0]]),
+                '2': np.array([[-4 + 0.6 * d, -1 + 0.1 * d], [-1 + 0.1 * d, 2.0]]),
+                '3': np.array([[0.0, 0.3 * d - 2], [0.3 * d - 2, -2 + 0.2 * d]])}
+    E = {'0': np.zeros((2, 2)), '1': -np.array([[1.0, 0], [0, 0]]),
+         '2': -np.array([[0, 1.0], [1.0, 0]]), '3': -np.array([[0, 0], [0, 1.0]])}
+    deltas = np.linspace(-0.2, 1.0, 7).reshape(-1, 1)
+    c = np.array([[-1.0], [0.0], [-1.0]])
+    return deltas, F_d, E, c, 0.1 * np.eye(3)
+
+
+def test_sdp_integer_keys_match_string_keys():
+    """Matrix dicts keyed 0, 1, ... give the same solution as keys '0', '1', ..."""
+    deltas, F_d, E, c, Q = _lpv_problem()
+    x_str, *_ = solve_sdp(deltas, F_d, E, c, Q, rho=1.0, solver="CLARABEL")
+    F_int = lambda delta: {int(k): v for k, v in F_d(delta).items()}
+    E_int = {int(k): v for k, v in E.items()}
+    x_int, *_ = solve_sdp(deltas, F_int, E_int, c, Q, rho=1.0, solver="CLARABEL")
+    np.testing.assert_allclose(x_int, x_str, atol=1e-5)
+
+
+def test_sdp_empty_q_means_no_quadratic_term():
+    """An empty Q is the same as a zero matrix."""
+    deltas, F_d, E, c, _ = _lpv_problem()
+    E_bounded = dict(E, **{'0': 0.0 * np.eye(2)})
+    # Without Q, the regularization (τ = 2 > √2) is what keeps max trace(P) bounded.
+    x_empty, *_ = solve_sdp(deltas, F_d, E_bounded, c, np.array([]), tau=2.0, rho=1.0, solver="CLARABEL")
+    x_zero, *_ = solve_sdp(deltas, F_d, E_bounded, c, np.zeros((3, 3)), tau=2.0, rho=1.0, solver="CLARABEL")
+    np.testing.assert_allclose(x_empty, x_zero, atol=1e-5)
+
+
+def test_sdp_column_x_ref_matches_vector_x_ref():
+    """x_ref given as a (d, 1) column gives the same solution as a 1-D vector."""
+    deltas, F_d, E, c, Q = _lpv_problem()
+    x_vec, *_ = solve_sdp(deltas, F_d, E, c, Q, tau=0.5, x_ref=np.array([1.0, 0.0, 1.0]), rho=1.0, solver="CLARABEL")
+    x_col, *_ = solve_sdp(deltas, F_d, E, c, Q, tau=0.5, x_ref=np.array([[1.0], [0.0], [1.0]]), rho=1.0, solver="CLARABEL")
+    np.testing.assert_allclose(x_col, x_vec, atol=1e-5)