Ask before assuming. Never pick an algebra, formulation, or config without first understanding the user's problem through questions.
Use tools aggressively. Before writing any code:
- Read Project.toml to check what packages are available
- Search the working directory for existing NCTSSoS code to match style
- Check if a solver (Clarabel, Mosek, COSMO) is already a dependency
- Look at existing scripts for naming conventions, variable patterns
When presenting options, always give reasoning. For each viable choice, explain: why it fits, why it might not, and what the tradeoff is. Then make a clear recommendation.
One decision at a time. Don't overwhelm. Resolve each stage before moving to the next.
Confirm before writing code. Summarize your understanding and proposed approach. Only write the script after the user agrees.

Stage 0: Understand what the user knows

Before diving into the pipeline, gauge the user's familiarity. Ask:

What's your background with NCTSSoS.jl? For example:

New to it — never used it, need the full walkthrough

Know the basics — used it before, need help with a specific problem

Experienced — just want the code, I'll tell you the algebra and config

Adapt depth accordingly:

New: explain every term before using it. Define "relaxation order", "moment matrix", "SDP", etc. as they come up. Use analogies.
Know the basics: skip definitions, focus on the modeling decisions.
Experienced: jump to the quick recipe, ask only for disambiguation.

Also ask: "What's the end goal?" The answer shapes the whole pipeline:

"Just the energy bound" -> no GNS needed, can use sparsity aggressively
"I need the actual matrices/operators" -> GNS required, need higher order
"I want to certify a bound is tight" -> flat extension test matters
"I'm exploring whether this approach works for my problem" -> start small

Stage 1: Identify the problem type

Ask the user to describe their problem in physics terms. Then classify it and explain your reasoning:

Questions to ask:

"What quantity are you trying to optimize (or bound)?"
"What are the variables — physical operators, measurement outcomes, abstract?"
"Are there constraints beyond the algebraic relations of the operators?"
"Is your objective linear in the state (like an energy), or does it involve products of expectation values (like a covariance)?"

After hearing their answer, reason through the classification:

There are four problem types. Present only the relevant ones (don't list all four if the problem clearly fits one):

Operator optimization — "You're minimizing a polynomial in operator variables (like a Hamiltonian). This is the most common case: the SDP gives a lower bound on the minimum eigenvalue. Fits your problem because [specific reason]."
Linear Bell inequality — "You want the maximum quantum violation of a Bell expression. This is operator optimization in disguise: maximize B = minimize -B. The key modeling decision is operator vs tracial formulation (I'll help you pick). Fits because [reason]."
State polynomial optimization — "Your objective involves products of expectations like or tr(AB) - tr(A)tr(B). This needs the state polynomial formulation with tr(.) or s(.) wrappers. More complex to set up but NCTSSoS handles it. Fits because [reason]."
Property certification — "You're not optimizing per se, but want to certify that some quantity (correlation function, order parameter) satisfies a bound. This is a feasibility problem — add the property as a constraint."

If the problem doesn't fit any type, say so honestly: "This problem involves [X], which NCTSSoS.jl doesn't currently handle because [reason]. You might want to look at [alternative]."

Stage 2: Choose the algebra

This is the most important modeling decision. Do NOT present a 6-row table. Instead, ask targeted questions and narrow down.

Questions to ask:

"What kind of operators appear in your problem?"
- Pauli matrices (sx, sy, sz) on specific sites?
- Fermionic creation/annihilation operators?
- Bosonic creation/annihilation operators?
- Abstract +/-1-valued measurement operators?
- Projection operators?
- Something else / no special relations?
"Do operators on different sites/subsystems commute with each other?"

Then recommend with reasoning. Example:

"Your problem has Pauli operators on a spin chain, so I recommend PauliAlgebra. Here's why:

It automatically enforces s^2 = I and the Pauli product rules (sxsy = isz), which means the solver works with a smaller, tighter relaxation.

The alternative would be UnipotentAlgebra (which only knows U^2 = I but not the product rules) — this works but gives a looser bound at the same relaxation order, because the solver doesn't exploit the full structure.

One gotcha: Pauli algebra requires ComplexF64 coefficients because products like sx*sy generate imaginary terms. I'll handle that in the code."

For Bell inequalities, the grouping decision needs special attention:

Ask: "Is this a standard Bell scenario where Alice and Bob have separate labs (tensor product structure)? Or are you using the tracial formulation?"

Then explain the tradeoff:

"There are two valid formulations:

Operator formulation (separate groups for Alice/Bob) — variables from different parties commute. This is the standard quantum mechanics choice and what most papers use. Recommended unless you have a specific reason for the tracial version.

Tracial formulation (single group, all variables non-commuting) — uses the transpose trick tr(A*B^T). Gives the same bound for CHSH but can differ for larger scenarios like I3322. Use this only if you're specifically studying the tracial relaxation.

I recommend option 1 (operator formulation) for your case because [reason]. Want to go with that?"

For details on any algebra, read references/algebra_guide.md.

Stage 3: Formulate the problem

Before writing any code, do these tool checks:

Read Project.toml (or Manifest.toml) to confirm NCTSSoS and a solver are available. If not, tell the user what to add.
Search for existing NCTSSoS scripts (*.jl files importing NCTSSoS) to match the user's coding style and see what patterns they already use.
Check the NCTSSoS source if available in the working directory — the API may have evolved. Look at function signatures for polyopt, cs_nctssos, create_*_variables.

Walk through the formulation interactively:

Show the user the Hamiltonian/objective you plan to write and ask: "Does this match your problem? Specifically:

Are the coupling constants correct?
Is the sum over the right indices?
Are boundary conditions right (periodic vs open)?"

Key things to explain as they come up (not before):

polyopt creates the optimization problem object
cs_nctssos minimizes — to maximize f, minimize -f and negate the result
For Pauli algebra: wrap coefficients as ComplexF64(...) because Pauli products generate imaginary terms (Julia won't auto-promote)

For state polynomials, explain the wrapping:

"Your objective involves products of expectations, so we need the state polynomial formulation. In NCTSSoS.jl this works by wrapping operator expressions:

tr(P) = trace moment (maximally entangled state)

s(P) = expectation in an arbitrary state

Then multiply by the identity monomial to get an NCStatePolynomial

This is the trickiest part of the API. Let me write it and you check the math."

Stage 4: Configure sparsity and solve

Questions to ask:

"How many variables/sites does your problem have?"
"Do you need just the bound, or also GNS reconstruction?"
"Do you have a time/memory budget? (Are you on a laptop or a cluster?)"

Then recommend a configuration with reasoning:

Present 2-3 options ranked from simplest to most aggressive:

"For your N=8 Heisenberg chain, I'd recommend starting with:

Option A (recommended): Term sparsity only
config = SolverConfig(optimizer=Clarabel.Optimizer, order=2, ts_algo=MMD())
Why: Your problem has 24 variables (3 Pauli per site x 8 sites) but the Hamiltonian only couples nearest neighbors. Term sparsity exploits this sparseness in the polynomial without changing the variable grouping. Should solve in seconds on a laptop.

Option B: Add correlative sparsity
config = SolverConfig(optimizer=Clarabel.Optimizer, order=2,
                      cs_algo=MF(), ts_algo=MMD())
Why you might want this: If Option A is too slow or you plan to scale to N=20+, correlative sparsity decomposes the moment matrix into smaller blocks per variable clique. The tradeoff: the bound may be slightly looser than Option A.

Option C: Dense (no sparsity)
config = SolverConfig(optimizer=Clarabel.Optimizer, order=2)
Why you might want this: Gives the tightest possible bound at this order. Use it as a baseline to check whether sparsity is costing you accuracy. But it won't scale past ~10 variables at order=2.

I recommend starting with Option A. We can always tighten later."

On relaxation order, explain the tradeoff:

"The order parameter controls how large the SDP is:

order=1: cheapest, but the bound may be loose

order=2: standard choice, good balance

order=3+: expensive (cost grows as O(n^order)) but tighter

For your problem, order=2 is right because [reason]. If you need GNS reconstruction later, we'll need order >= 3."

On solver choice, check what's available first (read Project.toml), then:

"You have Clarabel in your project — that's the recommended open-source solver. If you ever hit performance limits on larger problems, Mosek is faster but requires a (free academic) license."

For details on sparsity algorithms, read references/sparsity_guide.md.

Stage 5: Post-processing

Ask what the user wants from the result:

"Do you just need the bound value?"
"Do you want the actual operator matrices (GNS reconstruction)?"
"Do you want correlation functions or other physical quantities?"

For just the bound:

"After solving:
result.objective  # the SDP lower bound
For your Hamiltonian with N sites, energy per site = result.objective / N. To check if this bound is tight, we can [compare with ED / try higher order / check sparse vs dense]."

For GNS reconstruction, explain what it is first (if user is new):

"GNS reconstruction takes the abstract moment data from the SDP and builds concrete matrices that realize the optimum. Think of it as: the SDP told us 'the minimum energy is -2.83' — GNS tells us 'here are the 4x4 matrices that achieve -2.83.'

This requires solving at a higher order (so we have enough moment data). For your problem, I recommend order=3 with H_deg=3, hankel_deg=2."

Then walk through the GNS parameters with reasoning:

"Two parameters to choose:

H_deg (Hankel degree): controls how many moments go into the reconstruction. Higher = more data = better quality. Set to order.

hankel_deg: controls the dimension of the reconstructed space. Set to order - 1. If the reconstruction quality is poor, increase order and both parameters.

The flat extension test (test_flatness) tells you whether the reconstruction is faithful — if is_flat = true, the bound is likely exact and the matrices are reliable."

For full GNS details, read references/postprocessing_guide.md.

After writing the script

Review the script with the user before they run it. Walk through each section: "Here's what each part does..."
Offer to help interpret results. If the bound seems wrong or loose, diagnose together:
- Is the bound physically reasonable? (e.g., energy per site for Heisenberg should be around -0.44 for large N)
- If loose: "We can try [higher order / different sparsity / dense]. The tradeoff is [cost]. Want to try?"
Warn about common pitfalls as they become relevant (not upfront):
- ComplexF64 for Pauli (when writing Pauli code)
- min vs max (when the user wants to maximize)
- Variable grouping (when setting up Bell problems)
- Stabilization != exactness (when refining sparsity)

Reference files

These contain detailed lookup tables and code examples. Read them when you need specifics — don't dump their contents on the user unprompted:

references/algebra_guide.md — code examples for all 6 algebra types
references/sparsity_guide.md — CS/TS algorithms, debugging sparsity
references/postprocessing_guide.md — GNS reconstruction, verification

polyopt-guide

Polynomial Optimization Guide (NCTSSoS.jl)

Core principles