np_app/service/server.py

import json
import mimetypes
import os
import re
import threading
import time
import traceback
from decimal import Decimal
from concurrent.futures import FIRST_COMPLETED, ThreadPoolExecutor, wait
from http import HTTPStatus
from http.server import BaseHTTPRequestHandler, ThreadingHTTPServer
from pathlib import Path
import subprocess
import tempfile
from urllib.parse import parse_qs, urlparse
from uuid import uuid4

import numpy
from nupack import (
    Complex,
    Complementarity,
    Domain,
    Diversity,
    EnergyMatch,
    DesignOptions,
    Library,
    Match,
    Model,
    Pattern,
    Similarity,
    SSM,
    SetSpec,
    Strand,
    TargetComplex,
    TargetStrand,
    TargetTube,
    Tube,
    Weights,
    Window,
    complex_analysis,
    complex_design,
    tube_analysis,
    tube_design,
)
from nupack import config as nupack_config
from split_strand_svg import render_split_strands_svg

try:
    import redis
except ImportError:
    redis = None


ROOT = Path(__file__).resolve().parent
INDEX_PATH = ROOT / "index.html"
FAVICON_PATH = ROOT / "favicon.svg"
GUIDE_PATH = ROOT / "design-guide.html"

HOST = os.environ.get("NP_HOST", "0.0.0.0")
PORT = int(os.environ.get("NP_PORT", "18765"))
RNA_PLOT_CMD = os.environ.get("RNA_PLOT_CMD", "RNAplot")
ENABLE_RNAPLOT = os.environ.get("ENABLE_RNAPLOT", "1") != "0"
STRUCTURE_PLOT_MODE = os.environ.get("NP_STRUCTURE_PLOT_MODE", "auto").strip().lower()
RUN_MODE = os.environ.get("NP_RUN_MODE", "server")
REDIS_URL = os.environ.get("NP_REDIS_URL", "").strip()
JOB_QUEUE_KEY = os.environ.get("NP_JOB_QUEUE_KEY", "np_replica:jobs")
JOB_RUNNING_KEY = os.environ.get("NP_JOB_RUNNING_KEY", "np_replica:jobs:running")
WORKER_CONCURRENCY = int(os.environ.get("NP_WORKER_CONCURRENCY", "2"))
PER_JOB_THREAD_LIMIT = int(os.environ.get("NP_PER_JOB_THREAD_LIMIT", "1"))
NUPACK_CACHE_GB = float(os.environ.get("NP_NUPACK_CACHE_GB", "2.0"))

UNIT_SCALE = {
    "M": 1.0,
    "mM": 1e-3,
    "uM": 1e-6,
    "nM": 1e-9,
    "pM": 1e-12,
}
IUPAC_CODES = "ACGTUWSMKRYBDHVN"
IUPAC_CONSTRAINT_TOKEN = re.compile(rf"[{IUPAC_CODES}](?:\d+)?")

VALID_COMPUTE = {"pfunc", "pairs", "mfe", "sample", "subopt", "ensemble_size"}
JOB_STORE = {}
JOB_LOCK = threading.Lock()
JOB_TTL_SECONDS = int(os.environ.get("NP_JOB_TTL_SECONDS", "3600"))
JOB_MAX_COUNT = int(os.environ.get("NP_JOB_MAX_COUNT", "64"))
SHARE_STORE = {}
SHARE_LOCK = threading.Lock()
SHARE_MAX_COUNT = int(os.environ.get("NP_SHARE_MAX_COUNT", "100"))
SHARE_MAX_BYTES = int(os.environ.get("NP_SHARE_MAX_BYTES", str(12 * 1024 * 1024)))
SHARE_KEY_PREFIX = os.environ.get("NP_SHARE_KEY_PREFIX", "np_replica:shares")
SHARE_INDEX_KEY = f"{SHARE_KEY_PREFIX}:index"
REDIS_CLIENT = None


def json_bytes(payload, status=HTTPStatus.OK):
    return (
        status,
        "application/json; charset=utf-8",
        json.dumps(payload, ensure_ascii=False, indent=2).encode("utf-8"),
    )


def log_event(message):
    print(f"[np-replica] {time.strftime('%Y-%m-%d %H:%M:%S')} {message}", flush=True)


def apply_thread_limits():
    # Keep native math libraries and NUPACK's own executor from using all CPU cores inside one job.
    thread_limit_int = max(1, PER_JOB_THREAD_LIMIT)
    thread_limit = str(thread_limit_int)
    for key in (
        "OMP_NUM_THREADS",
        "OPENBLAS_NUM_THREADS",
        "MKL_NUM_THREADS",
        "NUMEXPR_NUM_THREADS",
        "VECLIB_MAXIMUM_THREADS",
        "GOTO_NUM_THREADS",
    ):
        os.environ[key] = thread_limit
    nupack_config.threads = thread_limit_int
    nupack_config.cache = max(0.1, NUPACK_CACHE_GB)


def html_bytes(path):
    body = path.read_bytes()
    mime, _ = mimetypes.guess_type(str(path))
    return HTTPStatus.OK, mime or "text/html; charset=utf-8", body


def file_etag(path):
    stat = path.stat()
    return f'W/"{stat.st_mtime_ns:x}-{stat.st_size:x}"'


def redis_enabled():
    return bool(REDIS_URL and redis is not None)


def redis_client():
    global REDIS_CLIENT
    if not redis_enabled():
        return None
    if REDIS_CLIENT is None:
        REDIS_CLIENT = redis.Redis.from_url(REDIS_URL, decode_responses=True)
    return REDIS_CLIENT


def queue_size():
    if not redis_enabled():
        with JOB_LOCK:
            return sum(1 for job in JOB_STORE.values() if job.get("status") in {"queued", "running"})
    try:
        return int(redis_client().llen(JOB_QUEUE_KEY))
    except Exception:
        return None


def prune_redis_running_jobs():
    if not redis_enabled():
        return 0
    client = redis_client()
    removed = 0
    try:
        for job_id in client.smembers(JOB_RUNNING_KEY):
            raw = client.get(job_key(job_id))
            if raw is None:
                client.srem(JOB_RUNNING_KEY, job_id)
                removed += 1
                continue
            try:
                if json.loads(raw).get("status") != "running":
                    client.srem(JOB_RUNNING_KEY, job_id)
                    removed += 1
            except Exception:
                client.srem(JOB_RUNNING_KEY, job_id)
                removed += 1
    except Exception:
        return removed
    return removed


def job_stats():
    stats = {
        "queued": 0,
        "running": 0,
        "success": 0,
        "error": 0,
    }
    if redis_enabled():
        prune_redis_running_jobs()
        client = redis_client()
        stats["queued"] = int(client.llen(JOB_QUEUE_KEY))
        stats["running"] = int(client.scard(JOB_RUNNING_KEY))
        return stats

    with JOB_LOCK:
        prune_jobs()
        for job in JOB_STORE.values():
            status = job.get("status")
            if status in stats:
                stats[status] += 1
    return stats


def job_key(job_id):
    return f"np_replica:job:{job_id}"


def share_key(share_id):
    return f"{SHARE_KEY_PREFIX}:item:{share_id}"


def normalize_sequence(sequence):
    return "".join(sequence.upper().split())


def normalize_design_sequence(sequence):
    return re.sub(r"\s+", "", str(sequence or "").upper())


def is_valid_iupac_constraint(sequence):
    seq = normalize_design_sequence(sequence)
    if not seq:
        return False
    index = 0
    while index < len(seq):
        match = IUPAC_CONSTRAINT_TOKEN.match(seq, index)
        if not match:
            return False
        index = match.end()
    return True


def ensure_unit_interval_limits(lower, upper, label):
    if not (0 <= lower <= upper <= 1):
        raise ValueError(f"{label} limits must satisfy 0 <= lower <= upper <= 1.")


def unit_to_molar(value, unit):
    if unit not in UNIT_SCALE:
        raise ValueError(f"Unsupported concentration unit: {unit}")
    return float(value) * UNIT_SCALE[unit]


def build_model(model_input):
    return Model(
        material=model_input.get("material", "rna"),
        ensemble=model_input.get("ensemble", "stacking"),
        celsius=float(model_input.get("celsius", 37.0)),
        sodium=float(model_input.get("sodium", 1.0)),
        magnesium=float(model_input.get("magnesium", 0.0)),
    )


def build_model_summary(model_input):
    return {
        "material": model_input.get("material", "rna"),
        "ensemble": model_input.get("ensemble", "stacking"),
        "celsius": float(model_input.get("celsius", 37.0)),
        "sodium": float(model_input.get("sodium", 1.0)),
        "magnesium": float(model_input.get("magnesium", 0.0)),
    }


def parse_compute(payload):
    compute = payload.get("compute") or ["pfunc", "mfe"]
    compute = [item for item in compute if item in VALID_COMPUTE]
    if not compute:
        raise ValueError("At least one compute option is required.")
    return compute


def parse_options(payload):
    raw = payload.get("options") or {}
    options = {
        "num_sample": int(raw.get("num_sample", 20)),
        "energy_gap": float(raw.get("energy_gap", 1.0)),
        "sparsity_fraction": float(raw.get("sparsity_fraction", 1.0)),
        "sparsity_threshold": float(raw.get("sparsity_threshold", 0.0)),
        "single_mfe": bool(raw.get("single_mfe", False)),
        "result_limit": int(raw.get("result_limit", 25)),
        "pairs_preview_size": int(raw.get("pairs_preview_size", 24)),
    }
    if options["num_sample"] < 0 or options["num_sample"] > 1000:
        raise ValueError("num_sample must be between 0 and 1000.")
    if options["energy_gap"] < 0:
        raise ValueError("energy_gap must be non-negative.")
    if not 0 <= options["sparsity_fraction"] <= 1:
        raise ValueError("sparsity_fraction must be between 0 and 1.")
    if options["sparsity_threshold"] < 0:
        raise ValueError("sparsity_threshold must be non-negative.")
    if options["result_limit"] < 1:
        raise ValueError("result_limit must be at least 1.")
    if options["pairs_preview_size"] < 8 or options["pairs_preview_size"] > 128:
        raise ValueError("pairs_preview_size must be between 8 and 128.")
    return options


def parse_design_options(payload):
    raw = payload.get("design") or {}
    options = {
        "trials": int(raw.get("trials", 1)),
        "result_limit": int(raw.get("result_limit", 25)),
        "off_target_max_size": int(raw.get("off_target_max_size", payload.get("tube", {}).get("max_size", 2))),
        "stop_condition": float(raw.get("stop_condition", 0.02)),
        "seed": int(raw.get("seed", 0)),
        "wobble_mutations": bool(raw.get("wobble_mutations", False)),
        "max_time_seconds": int(raw.get("max_time_seconds", 0)),
    }
    if options["trials"] < 1 or options["trials"] > 8:
        raise ValueError("design trials must be between 1 and 8.")
    if options["result_limit"] < 1:
        raise ValueError("design result_limit must be at least 1.")
    if options["off_target_max_size"] < 1 or options["off_target_max_size"] > 8:
        raise ValueError("design off_target_max_size must be between 1 and 8.")
    if not 0 < options["stop_condition"] < 1:
        raise ValueError("design stop_condition must be between 0 and 1.")
    if options["max_time_seconds"] < 0:
        raise ValueError("design max_time_seconds must be non-negative.")
    return options


def build_strands(strand_payload):
    if not strand_payload:
        raise ValueError("At least one strand is required.")

    strand_map = {}
    ordered = []

    for row in strand_payload:
        name = (row.get("name") or "").strip()
        sequence = normalize_sequence(row.get("sequence") or "")
        if not name:
            raise ValueError("Every strand needs a name.")
        if not sequence:
            raise ValueError(f"Strand {name} is missing a sequence.")
        if name in strand_map:
            raise ValueError(f"Duplicate strand name: {name}")
        if not re.fullmatch(r"[ACGTUWSMKRYBDHVN]+", sequence):
            raise ValueError(f"Strand {name} contains unsupported characters.")
        concentration = float(row.get("concentration", 0))
        if concentration < 0:
            raise ValueError(f"Strand {name} concentration must be non-negative.")

        strand = Strand(sequence, name=name)
        strand_map[name] = strand
        ordered.append(
            {
                "name": name,
                "sequence": sequence,
                "concentration": concentration,
                "unit": row.get("unit", "uM"),
                "object": strand,
            }
        )

    return strand_map, ordered


def build_design_domains(domain_payload):
    domain_map = {}
    ordered = []
    for row in domain_payload:
        name = (row.get("name") or "").strip()
        sequence_constraint = normalize_design_sequence(row.get("sequence") or "")
        if not name:
            raise ValueError("Every design domain needs a name.")
        if not sequence_constraint:
            raise ValueError(f"Design domain {name} is missing a sequence constraint.")
        if name in domain_map:
            raise ValueError(f"Duplicate domain name: {name}")
        if not is_valid_iupac_constraint(sequence_constraint):
            raise ValueError(f"Design domain {name} contains unsupported constraint characters.")

        domain = Domain(sequence_constraint, name=name)
        domain_map[name] = domain
        ordered.append(
            {
                "name": name,
                "constraint": sequence_constraint,
                "object": domain,
            }
        )
    return domain_map, ordered


def parse_domain_composition(text, domain_map):
    tokens = [token.strip() for token in re.split(r"[\s,]+", str(text or "").strip()) if token.strip()]
    if not tokens:
        raise ValueError("Domain composition cannot be empty.")

    domains = []
    for token in tokens:
        complement = False
        domain_name = token
        if token.startswith("~"):
            complement = True
            domain_name = token[1:]
        elif token.endswith("*"):
            complement = True
            domain_name = token[:-1]

        if not domain_name:
            raise ValueError(f"Invalid domain token: {token}")
        if domain_name not in domain_map:
            raise ValueError(f"Unknown domain name in strand composition: {domain_name}")

        domain = domain_map[domain_name]
        domains.append(~domain if complement else domain)
    return domains


def build_design_strands(strand_payload, domain_map=None):
    if not strand_payload:
        raise ValueError("At least one design strand is required.")

    domain_map = domain_map or {}
    use_domain_composition = bool(domain_map)
    target_strand_map = {}
    ordered = []

    for row in strand_payload:
        name = (row.get("name") or "").strip()
        raw_definition = str(row.get("sequence") or "").strip()
        sequence_constraint = normalize_design_sequence(raw_definition)
        if not name:
            raise ValueError("Every design strand needs a name.")
        if not raw_definition:
            raise ValueError(f"Design strand {name} is missing a strand definition.")
        if name in target_strand_map:
            raise ValueError(f"Duplicate strand name: {name}")

        if use_domain_composition:
            try:
                strand_domains = parse_domain_composition(raw_definition, domain_map)
            except ValueError:
                if (" " in raw_definition) or ("," in raw_definition) or ("~" in raw_definition) or ("*" in raw_definition):
                    raise
                if not is_valid_iupac_constraint(sequence_constraint):
                    raise
                inline_domain = Domain(sequence_constraint, name=name)
                strand_domains = [inline_domain]
                target_strand = TargetStrand(strand_domains, name=name)
                constraint_kind = "sequence_constraint"
                constraint_value = sequence_constraint
            else:
                inline_domain = None
                target_strand = TargetStrand(strand_domains, name=name)
                constraint_kind = "domain_composition"
                constraint_value = raw_definition
        else:
            if not is_valid_iupac_constraint(sequence_constraint):
                raise ValueError(f"Design strand {name} contains unsupported constraint characters.")
            inline_domain = Domain(sequence_constraint, name=name)
            strand_domains = [inline_domain]
            target_strand = TargetStrand(strand_domains, name=name)
            constraint_kind = "sequence_constraint"
            constraint_value = sequence_constraint

        target_strand_map[name] = target_strand
        ordered.append(
            {
                "name": name,
                "definition": raw_definition,
                "constraint": constraint_value,
                "constraint_kind": constraint_kind,
                "object": target_strand,
                "domains": strand_domains,
                "domain": inline_domain,
            }
        )

    return target_strand_map, ordered


def resolve_scope_reference(expr, domain_map, strand_map, *, allow_global=False):
    text = str(expr or "").strip()
    if not text:
        if allow_global:
            return None
        raise ValueError("Constraint scope cannot be empty.")

    if text in strand_map:
        return strand_map[text]
    return parse_domain_composition(text, domain_map)


def resolve_domain_scope(expr, domain_map):
    text = str(expr or "").strip()
    if not text:
        raise ValueError("Domain scope cannot be empty.")
    return parse_domain_composition(text, domain_map)


def resolve_target_complex_scope(expr, target_complex_map):
    text = str(expr or "").strip()
    if not text:
        raise ValueError("Target complex scope cannot be empty.")
    tokens = [token.strip() for token in re.split(r"[\s,]+", text) if token.strip()]
    if not tokens:
        raise ValueError("Target complex scope cannot be empty.")
    missing = [token for token in tokens if token not in target_complex_map]
    if missing:
        raise ValueError(f"Unknown target complex name(s): {', '.join(missing)}")
    return [target_complex_map[token] for token in tokens]


def parse_sequence_list(text):
    if isinstance(text, list):
        values = text
    else:
        normalized = str(text or "").replace("|", "\n").replace(";", "\n")
        values = normalized.splitlines()
    output = [item.strip() for item in values if item and item.strip()]
    if not output:
        raise ValueError("Constraint sequence list cannot be empty.")
    for seq in output:
        if not re.fullmatch(rf"[{IUPAC_CODES}]+", seq.upper()):
            raise ValueError(f"Invalid sequence source: {seq!r}")
    return output


def parse_pattern_list(text):
    if isinstance(text, list):
        values = text
    else:
        values = re.split(r"[\n,]+", str(text or ""))
    output = [item.strip().upper() for item in values if item and item.strip()]
    if not output:
        raise ValueError("Pattern constraint requires at least one pattern.")
    for pattern in output:
        if not is_valid_iupac_constraint(pattern):
            raise ValueError(f"Invalid pattern: {pattern!r}")
    return output


def parse_catalog_list(text):
    raw_lines = [line.strip() for line in str(text or "").splitlines() if line.strip()]
    if not raw_lines:
        raise ValueError("Library constraint requires at least one catalog row.")
    catalog = []
    for line in raw_lines:
        library = [item.strip().upper() for item in re.split(r"[,|]+", line) if item.strip()]
        if not library:
            raise ValueError("Library constraint contains an empty catalog row.")
        for seq in library:
            if not re.fullmatch(rf"[{IUPAC_CODES}]+", seq):
                raise ValueError(f"Invalid library sequence: {seq!r}")
        catalog.append(library)
    return catalog


def build_hard_constraints(payload, domain_map, strand_map):
    constraints = []
    for index, row in enumerate(payload.get("hard_constraints") or [], start=1):
        constraint_type = (row.get("type") or "").strip().lower()
        if not constraint_type:
            continue

        try:
            def require_text(field, label):
                value = str(row.get(field) or "").strip()
                if not value:
                    raise ValueError(f"{label} is required.")
                return value

            if constraint_type == "match":
                constraints.append(
                    Match(
                        resolve_scope_reference(require_text("left", "Match left scope"), domain_map, strand_map),
                        resolve_scope_reference(require_text("right", "Match right scope"), domain_map, strand_map),
                    )
                )
            elif constraint_type == "complementarity":
                constraints.append(
                    Complementarity(
                        resolve_scope_reference(require_text("left", "Complementarity left scope"), domain_map, strand_map),
                        resolve_scope_reference(require_text("right", "Complementarity right scope"), domain_map, strand_map),
                        wobble_mutations=bool(row.get("wobble_mutations", False)),
                    )
                )
            elif constraint_type == "diversity":
                word = int(row.get("word", 4))
                types = int(row.get("types", 2))
                if word < 1:
                    raise ValueError("Diversity word must be >= 1.")
                if types < 1 or types > 4:
                    raise ValueError("Diversity types must be between 1 and 4.")
                kwargs = {
                    "word": word,
                    "types": types,
                }
                scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True)
                if scope is not None:
                    kwargs["scope"] = scope
                constraints.append(Diversity(**kwargs))
            elif constraint_type == "similarity":
                lower = float(row.get("min_fraction", 0.0))
                upper = float(row.get("max_fraction", 1.0))
                ensure_unit_interval_limits(lower, upper, "Similarity")
                reference = normalize_design_sequence(row.get("reference") or "")
                if not is_valid_iupac_constraint(reference):
                    raise ValueError("Similarity reference must be a valid IUPAC constraint.")
                constraints.append(
                    Similarity(
                        resolve_scope_reference(require_text("scope", "Similarity scope"), domain_map, strand_map),
                        reference,
                        limits=[lower, upper],
                    )
                )
            elif constraint_type == "window":
                constraints.append(
                    Window(
                        resolve_scope_reference(require_text("scope", "Window scope"), domain_map, strand_map),
                        sources=parse_sequence_list(row.get("sources") or ""),
                    )
                )
            elif constraint_type == "library":
                constraints.append(
                    Library(
                        resolve_scope_reference(require_text("scope", "Library scope"), domain_map, strand_map),
                        catalog=parse_catalog_list(row.get("catalog") or ""),
                    )
                )
            elif constraint_type == "pattern":
                kwargs = {"patterns": parse_pattern_list(row.get("patterns") or "")}
                scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True)
                if scope is not None:
                    kwargs["scope"] = scope
                constraints.append(Pattern(**kwargs))
            else:
                raise ValueError(f"Unsupported hard constraint type: {constraint_type}")
        except Exception as exc:
            raise ValueError(f"Invalid hard constraint #{index}: {exc}") from exc

    return constraints


def build_soft_constraints(payload, domain_map, strand_map, target_complex_map):
    constraints = []
    for index, row in enumerate(payload.get("soft_constraints") or [], start=1):
        constraint_type = (row.get("type") or "").strip().lower()
        if not constraint_type:
            continue

        try:
            def require_text(field, label):
                value = str(row.get(field) or "").strip()
                if not value:
                    raise ValueError(f"{label} is required.")
                return value

            if constraint_type == "pattern":
                weight = float(row.get("weight", 1.0))
                if weight < 0:
                    raise ValueError("Pattern weight must be non-negative.")
                kwargs = {
                    "patterns": parse_pattern_list(row.get("patterns") or ""),
                    "weight": weight,
                }
                scope = resolve_scope_reference(row.get("scope"), domain_map, strand_map, allow_global=True)
                if scope is not None:
                    kwargs["scope"] = scope
                constraints.append(Pattern(**kwargs))
            elif constraint_type == "similarity":
                lower = float(row.get("min_fraction", 0.0))
                upper = float(row.get("max_fraction", 1.0))
                ensure_unit_interval_limits(lower, upper, "Similarity")
                reference = normalize_design_sequence(row.get("reference") or "")
                if not is_valid_iupac_constraint(reference):
                    raise ValueError("Similarity reference must be a valid IUPAC constraint.")
                weight = float(row.get("weight", 1.0))
                if weight < 0:
                    raise ValueError("Similarity weight must be non-negative.")
                constraints.append(
                    Similarity(
                        resolve_scope_reference(require_text("scope", "Similarity scope"), domain_map, strand_map),
                        reference,
                        limits=[lower, upper],
                        weight=weight,
                    )
                )
            elif constraint_type == "ssm":
                word = int(row.get("word", 4))
                if word < 1:
                    raise ValueError("SSM word must be >= 1.")
                weight = float(row.get("weight", 1.0))
                if weight < 0:
                    raise ValueError("SSM weight must be non-negative.")
                kwargs = {
                    "word": word,
                    "weight": weight,
                }
                scope = resolve_target_complex_scope(row.get("scope"), target_complex_map) if row.get("scope") else None
                if scope is not None:
                    kwargs["scope"] = scope
                constraints.append(SSM(**kwargs))
            elif constraint_type == "energy_match":
                weight = float(row.get("weight", 1.0))
                if weight < 0:
                    raise ValueError("EnergyMatch weight must be non-negative.")
                kwargs = {
                    "domains": resolve_domain_scope(require_text("scope", "EnergyMatch domain scope"), domain_map),
                    "weight": weight,
                }
                if row.get("energy_ref") not in {None, ""}:
                    kwargs["energy_ref"] = float(row.get("energy_ref"))
                constraints.append(EnergyMatch(**kwargs))
            else:
                raise ValueError(f"Unsupported soft constraint type: {constraint_type}")
        except Exception as exc:
            raise ValueError(f"Invalid soft constraint #{index}: {exc}") from exc

    return constraints


def parse_design_complexes(payload, target_strand_map):
    target_rows = payload.get("design_complexes") or payload.get("design_targets") or []
    if not target_rows:
        raise ValueError("At least one design target complex is required.")

    targets = []
    target_complex_map = {}

    for idx, row in enumerate(target_rows, start=1):
        name = (row.get("name") or f"target_{idx}").strip() or f"target_{idx}"
        if name in target_complex_map:
            raise ValueError(f"Duplicate design target complex name: {name}")
        strand_text = (row.get("strands") or "").strip()
        structure = (row.get("structure") or "").strip()
        if not strand_text:
            raise ValueError(f"Design target complex {name} is missing strand membership.")
        if not structure:
            raise ValueError(f"Design target complex {name} is missing a target structure.")

        tokens = [token.strip() for token in strand_text.split("+") if token.strip()]
        if not tokens:
            raise ValueError(f"Design target complex {name} has an invalid strand list.")
        missing = [token for token in tokens if token not in target_strand_map]
        if missing:
            raise ValueError(
                f"Unknown strand name(s) in design target complex {name}: {', '.join(missing)}"
            )

        target_complex = TargetComplex(
            [target_strand_map[token] for token in tokens],
            structure,
            name=name,
        )
        target_payload = {
            "name": name,
            "strands": tokens,
            "structure": structure,
            "object": target_complex,
        }
        targets.append(target_payload)
        target_complex_map[name] = target_complex

    return targets, target_complex_map


def parse_design_tubes(payload, target_rows, target_complex_map, default_max_size):
    tube_rows = payload.get("design_tubes") or []
    if not tube_rows:
        legacy_targets = payload.get("design_targets") or []
        if legacy_targets:
            on_targets = []
            for row in legacy_targets:
                if row.get("concentration") in {None, ""}:
                    continue
                on_targets.append(
                    {
                        "complex": row.get("name"),
                        "concentration": row.get("concentration"),
                        "unit": row.get("unit", "uM"),
                    }
                )
            if on_targets:
                tube_rows = [{
                    "name": ((payload.get("tube") or {}).get("name") or "design_tube"),
                    "max_size": default_max_size,
                    "on_targets": on_targets,
                }]

    if not tube_rows:
        raise ValueError("Tube design requires at least one target tube.")

    tubes = []
    ordered_rows = []
    tube_name_set = set()
    for index, row in enumerate(tube_rows, start=1):
        tube_name = (row.get("name") or f"tube_{index}").strip() or f"tube_{index}"
        if tube_name in tube_name_set:
            raise ValueError(f"Duplicate target tube name: {tube_name}")
        tube_name_set.add(tube_name)
        max_size = int(row.get("max_size", default_max_size))
        if max_size < 1 or max_size > 8:
            raise ValueError(f"Target tube {tube_name} max_size must be between 1 and 8.")

        on_targets_input = row.get("on_targets") or []
        if not on_targets_input:
            raise ValueError(f"Target tube {tube_name} must include at least one on-target complex.")

        on_targets = {}
        serialized_on_targets = []
        seen_complex_names = set()
        for entry in on_targets_input:
            complex_name = (entry.get("complex") or "").strip()
            if complex_name not in target_complex_map:
                raise ValueError(f"Unknown target complex {complex_name} in tube {tube_name}.")
            if complex_name in seen_complex_names:
                raise ValueError(f"Duplicate on-target complex {complex_name} in tube {tube_name}.")
            seen_complex_names.add(complex_name)
            concentration = float(entry.get("concentration", 0))
            unit = entry.get("unit", "uM")
            if concentration <= 0:
                raise ValueError(f"Target tube {tube_name} requires positive target concentrations.")
            concentration_M = unit_to_molar(concentration, unit)
            complex_obj = target_complex_map[complex_name]
            on_targets[complex_obj] = concentration_M
            serialized_on_targets.append(
                {
                    "complex": complex_name,
                    "concentration": concentration,
                    "unit": unit,
                    "target_concentration_M": concentration_M,
                }
            )

        tube = TargetTube(
            on_targets=on_targets,
            off_targets=SetSpec(max_size=max_size),
            name=tube_name,
        )
        ordered_rows.append(
            {
                "name": tube_name,
                "max_size": max_size,
                "on_targets": serialized_on_targets,
                "object": tube,
            }
        )
        tubes.append(tube)

    return ordered_rows, tubes


def validate_design_object_names(design_domains, design_strands, target_rows, tube_rows):
    name_map = {}
    for kind, rows in (
        ("domain", design_domains),
        ("strand", design_strands),
        ("target complex", target_rows),
        ("target tube", tube_rows),
    ):
        for row in rows:
            name = (row.get("name") or "").strip()
            if not name:
                continue
            prev = name_map.get(name)
            if prev and prev != kind:
                raise ValueError(
                    f"Object name {name!r} is reused across {prev} and {kind}. "
                    "Design object names must be unique in a job."
                )
            name_map[name] = kind


def build_defect_weights(payload, design_mode, design_domains, design_strands, target_rows, tube_rows):
    rows = payload.get("defect_weights") or []
    if not rows:
        return None

    domain_index = {item["name"]: item["object"] for item in design_domains}
    strand_index = {item["name"]: item["object"] for item in design_strands}
    complex_index = {item["name"]: item["object"] for item in target_rows}
    tube_index = {item["name"]: item["object"] for item in tube_rows}
    weight_base = [item["object"] for item in tube_rows] if design_mode == "tube" else [item["object"] for item in target_rows]
    weights = Weights(weight_base)

    for index, row in enumerate(rows, start=1):
        kind = (row.get("kind") or "").strip().lower()
        value = float(row.get("weight", 1.0))
        if value < 0:
            raise ValueError(f"Invalid defect weight #{index}: weight must be non-negative")
        name = (row.get("name") or "").strip()
        try:
            if kind == "global":
                weights[:] = value
            elif kind == "domain":
                weights[domain_index[name]] = value
            elif kind == "strand":
                weights[:, strand_index[name]] = value
            elif kind == "complex":
                weights[:, :, complex_index[name]] = value
            elif kind == "tube":
                if design_mode != "tube":
                    raise ValueError("tube-level defect weights are only valid in tube design mode")
                weights[:, :, :, tube_index[name]] = value
            else:
                raise ValueError(f"Unsupported defect weight kind: {kind}")
        except KeyError as exc:
            raise ValueError(f"Unknown entity for defect weight #{index}: {name}") from exc
        except Exception as exc:
            raise ValueError(f"Invalid defect weight #{index}: {exc}") from exc

    return weights


def parse_complex_lines(text, strand_map):
    complexes = []
    raw_lines = [line.strip() for line in text.splitlines() if line.strip()]
    if not raw_lines:
        raise ValueError("Complex mode requires at least one complex definition.")

    for idx, line in enumerate(raw_lines, start=1):
        tokens = [token.strip() for token in line.split("+") if token.strip()]
        if not tokens:
            raise ValueError(f"Invalid complex definition on line {idx}: {line}")
        missing = [token for token in tokens if token not in strand_map]
        if missing:
            raise ValueError(
                f"Unknown strand name(s) in complex line {idx}: {', '.join(missing)}"
            )
        strands = [strand_map[token] for token in tokens]
        complexes.append(Complex(strands, name=f"complex_{idx}"))

    return complexes


def stringify_complex(complex_obj):
    return " + ".join(strand.name for strand in complex_obj.strands)


def flatten_sequence(complex_obj):
    return "".join(str(strand) for strand in complex_obj.strands)


def sanitize_structure_for_rnaplot(structure):
    return str(structure).replace("+", "")


def safe_name(name):
    return re.sub(r"[^A-Za-z0-9_.-]+", "_", name).strip("_") or "structure"


def render_rnaplot_svg(complex_obj, structure, suffix):
    if not ENABLE_RNAPLOT:
        return None

    seq_name = safe_name(f"{stringify_complex(complex_obj)}_{suffix}")
    sequence = flatten_sequence(complex_obj)
    structure_text = sanitize_structure_for_rnaplot(structure)

    with tempfile.TemporaryDirectory(prefix="rnaplot-") as tmpdir:
        tmp_path = Path(tmpdir)
        input_path = tmp_path / f"{seq_name}.seq"
        output_path = tmp_path / f"{seq_name}_ss.svg"

        input_path.write_text(
            f">{seq_name}\n{sequence}\n{structure_text}\n",
            encoding="utf-8",
            newline="\n",
        )

        process = subprocess.run(
            [RNA_PLOT_CMD, "-f", "svg", input_path.name],
            cwd=tmpdir,
            capture_output=True,
            text=True,
            timeout=20,
        )

        if process.returncode != 0:
            raise RuntimeError(process.stderr.strip() or "RNAplot failed.")

        if not output_path.exists():
            raise RuntimeError("RNAplot finished without producing an SVG file.")

        return output_path.read_text(encoding="utf-8")


def render_structure_svg(complex_obj, structure, suffix):
    structure_text = str(structure)
    multistrand = len(complex_obj.strands) > 1 and "+" in structure_text
    should_try_split = STRUCTURE_PLOT_MODE == "split" or (
        STRUCTURE_PLOT_MODE != "rnaplot" and multistrand
    )
    split_error = None

    if should_try_split:
        try:
            return render_split_strands_svg(
                [str(strand) for strand in complex_obj.strands],
                structure_text,
                title=safe_name(f"{stringify_complex(complex_obj)}_{suffix}"),
            )
        except Exception as exc:
            split_error = exc
            if STRUCTURE_PLOT_MODE == "split":
                raise RuntimeError(f"Split-strand layout failed: {exc}") from exc

    if STRUCTURE_PLOT_MODE == "split":
        if split_error is not None:
            raise RuntimeError(f"Split-strand layout failed: {split_error}") from split_error
        return None

    try:
        return render_rnaplot_svg(complex_obj, structure, suffix)
    except Exception as exc:
        if split_error is not None:
            raise RuntimeError(
                f"Split-strand layout failed: {split_error}; RNAplot fallback failed: {exc}"
            ) from exc
        raise


def serialize_structures(complex_obj, items, plot_kind=None):
    output = []
    for index, item in enumerate(items or []):
        row = {
            "structure": str(item.structure),
            "energy": round(float(item.energy), 6),
            "stack_energy": round(float(item.stack_energy), 6),
        }
        if plot_kind and index == 0:
            try:
                row["rnaplot_svg"] = render_structure_svg(complex_obj, item.structure, plot_kind)
            except Exception as exc:
                row["rnaplot_error"] = str(exc)
        output.append(row)
    return output


def serialize_samples(items):
    return [str(item) for item in (items or [])]


def serialize_pairs(pair_matrix, preview_limit=24):
    if pair_matrix is None:
        return None

    dense = numpy.asarray(pair_matrix.to_array(), dtype=float)
    preview = dense[:preview_limit, :preview_limit]
    return {
        "shape": list(dense.shape),
        "preview": numpy.round(preview, 6).tolist(),
        "preview_size": int(preview_limit),
        "preview_truncated": bool(dense.shape[0] > preview_limit),
    }


def parse_structure_pairs(structure):
    stack = []
    pairs = []
    compact_index = -1
    index_map = []

    for char in str(structure):
        if char == "+":
            continue
        compact_index += 1
        index_map.append(compact_index)
        if char == "(":
            stack.append(compact_index)
        elif char == ")":
            if not stack:
                raise ValueError("Unbalanced structure: missing opening bracket.")
            left = stack.pop()
            pairs.append((left, compact_index))

    if stack:
        raise ValueError("Unbalanced structure: missing closing bracket.")

    pairs.sort()
    return pairs, compact_index + 1


def build_structure_probabilities(structure, pair_matrix):
    if pair_matrix is None:
        return None

    dense = numpy.asarray(pair_matrix.to_array(), dtype=float)
    mfe_pairs, structure_length = parse_structure_pairs(structure)

    if dense.shape[0] != dense.shape[1]:
        return None
    if structure_length != dense.shape[0]:
        return None

    pair_partner = {}
    for left, right in mfe_pairs:
        pair_partner[left] = right
        pair_partner[right] = left

    residue_probabilities = []
    for index in range(structure_length):
        partner = pair_partner.get(index)
        if partner is None:
            probability = float(dense[index, index])
        else:
            probability = float(dense[index, partner])
        residue_probabilities.append(round(probability, 6))

    pair_probabilities = [
        {
            "i": left + 1,
            "j": right + 1,
            "probability": round(float(dense[left, right]), 6),
        }
        for left, right in mfe_pairs
    ]

    return {
        "length": structure_length,
        "residue_probabilities": residue_probabilities,
        "pair_probabilities": pair_probabilities,
    }


def serialize_complex_result(complex_obj, data, pairs_preview_size=24):
    payload = {
        "name": getattr(complex_obj, "name", None),
        "display": stringify_complex(complex_obj),
        "strand_names": [strand.name for strand in complex_obj.strands],
        "sequence": flatten_sequence(complex_obj),
        "strand_lengths": [len(str(strand)) for strand in complex_obj.strands],
    }

    if data.pfunc is not None:
        payload["pfunc"] = format(data.pfunc, "g") if isinstance(data.pfunc, Decimal) else str(data.pfunc)
    if data.free_energy is not None:
        payload["free_energy_kcal_mol"] = round(float(data.free_energy), 6)
    if data.ensemble_size is not None:
        payload["ensemble_size"] = int(data.ensemble_size)
    if data.mfe_stack is not None:
        payload["mfe_stack_kcal_mol"] = round(float(data.mfe_stack), 6)
    if data.mfe is not None:
        payload["mfe"] = serialize_structures(complex_obj, data.mfe, plot_kind="mfe")
    if data.subopt is not None:
        payload["subopt"] = serialize_structures(complex_obj, data.subopt)
    if data.sample is not None:
        payload["sample"] = serialize_samples(data.sample)
    if data.pairs is not None:
        payload["pairs"] = serialize_pairs(data.pairs, preview_limit=pairs_preview_size)
        if payload.get("mfe"):
            payload["mfe"][0]["pair_probability_annotations"] = build_structure_probabilities(
                payload["mfe"][0]["structure"], data.pairs
            )

    return payload


def attach_pairs_result(payload, pair_data, pairs_preview_size=24):
    if pair_data is None or pair_data.pairs is None:
        return payload

    payload["pairs"] = serialize_pairs(pair_data.pairs, preview_limit=pairs_preview_size)
    if payload.get("mfe"):
        payload["mfe"][0]["pair_probability_annotations"] = build_structure_probabilities(
            payload["mfe"][0]["structure"], pair_data.pairs
        )
    return payload


def serialize_target_defect_rows(df, limit=None):
    if df is None:
        return []
    rows = []
    for record in df.to_dict("records"):
        row = {}
        for key, value in record.items():
            if key in {"tube", "complex"}:
                continue
            if isinstance(value, numpy.floating):
                row[key] = float(value)
            else:
                row[key] = value
        rows.append(row)
    if limit is not None:
        rows = rows[:limit]
    return rows


def serialize_design_result(
    design_result,
    design_mode,
    model_summary,
    design_options,
    target_rows,
    tube_rows,
    ordered_domains,
    ordered_strands,
):
    designed_domains = []
    designed_domain_map = getattr(design_result, "domains", {}) or {}
    for item in ordered_domains:
        designed_domain = designed_domain_map.get(item["object"])
        designed_domains.append(
            {
                "name": item["name"],
                "constraint": item["constraint"],
                "sequence": str(designed_domain) if designed_domain is not None else None,
                "length": len(str(designed_domain)) if designed_domain is not None else None,
            }
        )

    designed_strands = []
    for item in ordered_strands:
        target_strand = item["object"]
        analysis_strand = design_result.to_analysis[target_strand]
        designed_strands.append(
            {
                "name": item["name"],
                "constraint": item["constraint"],
                "constraint_kind": item.get("constraint_kind", "sequence_constraint"),
                "definition": item.get("definition", item["constraint"]),
                "sequence": str(analysis_strand),
                "length": len(str(analysis_strand)),
            }
        )

    target_complexes = []
    for target in target_rows:
        target_complex = target["object"]
        analysis_complex = design_result.to_analysis[target_complex]
        target_complexes.append(
            {
                "name": target["name"],
                "display": stringify_complex(analysis_complex),
                "strand_names": list(target["strands"]),
                "structure": target["structure"],
                "sequence": flatten_sequence(analysis_complex),
                "target_concentration_M": target.get("target_concentration_M"),
            }
        )

    concentration_rows = []
    concentration_table = getattr(getattr(design_result, "concentrations", None), "table", None)
    if design_mode == "tube" and concentration_table is not None:
        for record in concentration_table.to_dict("records"):
            concentration_rows.append(
                {
                    "complex_name": record["complex_name"],
                    "concentration_M": float(record["concentration"]),
                    "target_concentration_M": float(record["target_concentration"]),
                    "nucleotides": int(record["nucleotides"]),
                }
            )
        concentration_rows.sort(key=lambda row: row["concentration_M"], reverse=True)

    return {
        "workflow": "design",
        "mode": design_mode,
        "model": model_summary,
        "compute": ["design"],
        "options": design_options,
        "strands": designed_strands,
        "complexes": target_complexes,
        "total_complex_count": len(target_complexes),
        "displayed_complex_count": len(target_complexes),
        "tube": {
            "name": tube_rows[0]["name"] if tube_rows else "design_tube",
            "max_size": tube_rows[0]["max_size"] if tube_rows else design_options["off_target_max_size"],
            "complex_concentrations": concentration_rows[: design_options["result_limit"]],
            "total_complex_concentrations": len(concentration_rows),
        } if design_mode == "tube" else None,
        "design": {
            "ensemble_defect": float(design_result.ensemble_defect),
            "domains": designed_domains,
            "stats": {
                key: float(value) if isinstance(value, (int, float, numpy.floating)) else value
                for key, value in design_result.stats.items()
            },
            "job_options": {
                "f_stop": design_options["stop_condition"],
                "seed": design_options["seed"],
                "wobble_mutations": design_options["wobble_mutations"],
                "max_time_seconds": design_options["max_time_seconds"],
            },
            "target_tubes": [
                {
                    "name": row["name"],
                    "max_size": row["max_size"],
                    "on_targets": row["on_targets"],
                }
                for row in tube_rows
            ],
            "targets": target_complexes,
            "defects": {
                "tubes": serialize_target_defect_rows(getattr(design_result.defects, "tubes", None)),
                "complexes": serialize_target_defect_rows(getattr(design_result.defects, "complexes", None)),
                "tube_complexes": serialize_target_defect_rows(getattr(design_result.defects, "tube_complexes", None)),
            },
            "concentrations": concentration_rows[: design_options["result_limit"]],
        },
    }


def sort_concentrations(complex_concentrations):
    rows = []
    for complex_obj, value in complex_concentrations.items():
        rows.append(
            {
                "display": stringify_complex(complex_obj),
                "strand_names": [strand.name for strand in complex_obj.strands],
                "concentration_M": float(value),
            }
        )
    rows.sort(key=lambda row: row["concentration_M"], reverse=True)
    return rows


def run_job_payload(payload):
    workflow = payload.get("workflow", "analysis")
    mode = payload.get("mode", "tube")
    model_input = payload.get("model") or {}
    model = build_model(model_input)
    model_summary = build_model_summary(model_input)

    if workflow == "design":
        design_options = parse_design_options(payload)
        design_job_options = DesignOptions(
            f_stop=design_options["stop_condition"],
            seed=design_options["seed"],
            wobble_mutations=design_options["wobble_mutations"],
            max_time=design_options["max_time_seconds"],
        )
        design_domain_map, design_domains = build_design_domains(payload.get("design_domains") or [])
        target_strand_map, design_strands = build_design_strands(
            payload.get("strands") or [],
            domain_map=design_domain_map,
        )
        target_rows, target_complex_map = parse_design_complexes(payload, target_strand_map)
        hard_constraints = build_hard_constraints(payload, design_domain_map, target_strand_map)
        soft_constraints = build_soft_constraints(
            payload,
            design_domain_map,
            target_strand_map,
            target_complex_map,
        )
        tube_rows = []
        design_tubes = []
        if mode == "tube":
            tube_rows, design_tubes = parse_design_tubes(
                payload,
                target_rows,
                target_complex_map,
                design_options["off_target_max_size"],
            )
        validate_design_object_names(design_domains, design_strands, target_rows, tube_rows)
        defect_weights = build_defect_weights(
            payload,
            mode,
            design_domains,
            design_strands,
            target_rows,
            tube_rows,
        )

        if mode == "tube":
            tube_name = ((payload.get("tube") or {}).get("name") or "design_tube").strip() or "design_tube"
            design_job = tube_design(
                tubes=design_tubes,
                model=model,
                options=design_job_options,
                hard_constraints=hard_constraints,
                soft_constraints=soft_constraints,
                defect_weights=defect_weights,
            )
        elif mode == "complex":
            design_job = complex_design(
                complexes=[row["object"] for row in target_rows],
                model=model,
                options=design_job_options,
                hard_constraints=hard_constraints,
                soft_constraints=soft_constraints,
                defect_weights=defect_weights,
            )
        else:
            raise ValueError(f"Unsupported design mode: {mode}")

        results = design_job.run(trials=design_options["trials"])
        best_result = min(results, key=lambda item: float(item.ensemble_defect))
        return serialize_design_result(
            best_result,
            mode,
            model_summary,
            design_options,
            target_rows,
            tube_rows,
            design_domains,
            design_strands,
        )

    if workflow != "analysis":
        raise ValueError(f"Unsupported workflow: {workflow}")

    compute = parse_compute(payload)
    options = parse_options(payload)
    nupack_options = {
        key: value
        for key, value in options.items()
        if key not in {"result_limit", "pairs_preview_size"}
    }
    strand_map, strands = build_strands(payload.get("strands") or [])

    if mode == "tube":
        wants_pairs = "pairs" in compute
        tube_compute = [item for item in compute if item != "pairs"]
        if not tube_compute:
            tube_compute = ["pfunc"]
        tube_cfg = payload.get("tube") or {}
        max_size = int(tube_cfg.get("max_size", 2))
        if max_size < 1 or max_size > 8:
            raise ValueError("max_size must be between 1 and 8.")
        strand_concentrations = {
            item["object"]: unit_to_molar(item["concentration"], item["unit"])
            for item in strands
        }
        tube = Tube(
            strands=strand_concentrations,
            complexes=SetSpec(max_size=max_size),
            name=(tube_cfg.get("name") or "tube1").strip() or "tube1",
        )
        result = tube_analysis([tube], model=model, compute=tube_compute, options=nupack_options)
        concentration_rows = sort_concentrations(result[tube].complex_concentrations)
        total_complex_count = len(concentration_rows)
        result_limit = options["result_limit"]
        displayed_concentration_rows = concentration_rows[:result_limit]
        concentration_by_display = {
            row["display"]: row["concentration_M"] for row in displayed_concentration_rows
        }
        pair_data_by_display = {}
        if wants_pairs:
            displayed_complex_objects = [
                complex_obj
                for complex_obj in result.complexes
                if stringify_complex(complex_obj) in concentration_by_display
            ]
            if displayed_complex_objects:
                pair_result = complex_analysis(
                    displayed_complex_objects,
                    model=model,
                    compute=["pairs"],
                    options=nupack_options,
                )
                pair_data_by_display = {
                    stringify_complex(complex_obj): data
                    for complex_obj, data in pair_result.complexes.items()
                }
        complexes = []
        for complex_obj, data in result.complexes.items():
            display_name = stringify_complex(complex_obj)
            if display_name not in concentration_by_display:
                continue
            row = serialize_complex_result(
                complex_obj,
                data,
                pairs_preview_size=options["pairs_preview_size"],
            )
            if wants_pairs:
                attach_pairs_result(
                    row,
                    pair_data_by_display.get(display_name),
                    pairs_preview_size=options["pairs_preview_size"],
                )
            row["concentration_M"] = concentration_by_display[row["display"]]
            complexes.append(row)
        complexes.sort(key=lambda row: row["concentration_M"], reverse=True)
        return {
            "workflow": "analysis",
            "mode": "tube",
            "model": model_summary,
            "compute": compute,
            "options": options,
            "strands": [
                {
                    "name": item["name"],
                    "sequence": item["sequence"],
                    "concentration": item["concentration"],
                    "unit": item["unit"],
                    "concentration_M": unit_to_molar(item["concentration"], item["unit"]),
                }
                for item in strands
            ],
            "tube": {
                "name": tube.name,
                "max_size": max_size,
                "complex_concentrations": displayed_concentration_rows,
                "total_complex_concentrations": total_complex_count,
            },
            "complexes": complexes,
            "total_complex_count": total_complex_count,
            "displayed_complex_count": len(complexes),
        }

    if mode == "complex":
        complexes = parse_complex_lines(payload.get("complexes_text", ""), strand_map)
        result = complex_analysis(complexes, model=model, compute=compute, options=nupack_options)
        rows = [
            serialize_complex_result(
                complex_obj,
                data,
                pairs_preview_size=options["pairs_preview_size"],
            )
            for complex_obj, data in result.complexes.items()
        ]
        rows.sort(key=lambda row: row["display"])
        total_complex_count = len(rows)
        rows = rows[: options["result_limit"]]
        return {
            "workflow": "analysis",
            "mode": "complex",
            "model": model_summary,
            "compute": compute,
            "options": options,
            "strands": [
                {
                    "name": item["name"],
                    "sequence": item["sequence"],
                }
                for item in strands
            ],
            "complexes": rows,
            "total_complex_count": total_complex_count,
            "displayed_complex_count": len(rows),
        }

    raise ValueError(f"Unsupported analysis mode: {mode}")


def prune_jobs(now=None):
    now = now or time.time()
    expired_ids = []
    for job_id, job in JOB_STORE.items():
        updated_at = job.get("updated_at", job.get("created_at", now))
        if now - updated_at > JOB_TTL_SECONDS:
            expired_ids.append(job_id)
    for job_id in expired_ids:
        JOB_STORE.pop(job_id, None)

    if len(JOB_STORE) > JOB_MAX_COUNT:
        keep_ids = sorted(
            JOB_STORE,
            key=lambda job_id: JOB_STORE[job_id].get("created_at", 0),
            reverse=True,
        )[:JOB_MAX_COUNT]
        keep_set = set(keep_ids)
        for job_id in list(JOB_STORE):
            if job_id not in keep_set:
                JOB_STORE.pop(job_id, None)


def set_job_data(job):
    if redis_enabled():
        redis_client().setex(job_key(job["job_id"]), JOB_TTL_SECONDS, json.dumps(job, ensure_ascii=False))
        return
    with JOB_LOCK:
        prune_jobs(job.get("updated_at"))
        JOB_STORE[job["job_id"]] = dict(job)


def update_job_data(job_id, **updates):
    if redis_enabled():
        client = redis_client()
        current = get_job(job_id)
        if current is None:
            current = {"job_id": job_id, "created_at": time.time()}
        current.update(updates)
        current["updated_at"] = time.time()
        client.setex(job_key(job_id), JOB_TTL_SECONDS, json.dumps(current, ensure_ascii=False))
        return current

    with JOB_LOCK:
        current = JOB_STORE.get(job_id)
        if current is None:
            current = {"job_id": job_id, "created_at": time.time()}
        current.update(updates)
        current["updated_at"] = time.time()
        prune_jobs(current["updated_at"])
        JOB_STORE[job_id] = current
        return dict(current)


def create_job(payload):
    job_id = uuid4().hex
    now = time.time()
    job = {
        "job_id": job_id,
        "status": "queued",
        "error": None,
        "result": None,
        "created_at": now,
        "updated_at": now,
        "payload": payload,
    }
    set_job_data(job)
    log_event(f"accepted job_id={job_id} mode={payload.get('mode', 'tube')}")

    if redis_enabled():
        redis_client().lpush(JOB_QUEUE_KEY, job_id)
    else:
        thread = threading.Thread(target=_run_job, args=(job_id, payload), daemon=True)
        thread.start()
    return job_id


def _run_job(job_id, payload):
    started_at = time.time()
    log_event(f"running job_id={job_id}")
    update_job_data(job_id, status="running")
    if redis_enabled():
        redis_client().sadd(JOB_RUNNING_KEY, job_id)
    try:
        result = run_job_payload(payload)
        elapsed = round(time.time() - started_at, 3)
        update_job_data(job_id, status="success", result=result, payload=None, elapsed_seconds=elapsed)
        log_event(f"success job_id={job_id} elapsed={elapsed}s")
    except Exception as exc:
        elapsed = round(time.time() - started_at, 3)
        update_job_data(
            job_id,
            status="error",
            error={
                "message": str(exc),
                "traceback": traceback.format_exc(),
            },
            payload=None,
            elapsed_seconds=elapsed,
        )
        log_event(f"error job_id={job_id} elapsed={elapsed}s message={exc}")
    finally:
        if redis_enabled():
            redis_client().srem(JOB_RUNNING_KEY, job_id)


def get_job(job_id):
    if redis_enabled():
        raw = redis_client().get(job_key(job_id))
        if raw is None:
            return None
        job = json.loads(raw)
        job.pop("payload", None)
        return job
    with JOB_LOCK:
        prune_jobs()
        job = JOB_STORE.get(job_id)
        if job is None:
            return None
        output = dict(job)
        output.pop("payload", None)
        return output


def prune_shares(now=None):
    if SHARE_MAX_COUNT < 1:
        return
    if redis_enabled():
        client = redis_client()
        extra = client.zcard(SHARE_INDEX_KEY) - SHARE_MAX_COUNT
        if extra <= 0:
            return
        stale_ids = client.zrange(SHARE_INDEX_KEY, 0, extra - 1)
        if stale_ids:
            client.delete(*(share_key(share_id) for share_id in stale_ids))
            client.zrem(SHARE_INDEX_KEY, *stale_ids)
        return

    keep_ids = sorted(
        SHARE_STORE,
        key=lambda share_id: SHARE_STORE[share_id].get("created_at", 0),
        reverse=True,
    )[:SHARE_MAX_COUNT]
    keep_set = set(keep_ids)
    for share_id in list(SHARE_STORE):
        if share_id not in keep_set:
            SHARE_STORE.pop(share_id, None)


def create_share(record):
    payload = record.get("payload")
    result = record.get("result")
    if not isinstance(payload, dict):
        raise ValueError("Share payload must contain an input payload object.")
    if not isinstance(result, dict):
        raise ValueError("Share payload must contain a result object.")

    now = time.time()
    share_id = uuid4().hex[:16]
    item = {
        "id": share_id,
        "created_at": now,
        "created_at_iso": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime(now)),
        "payload": payload,
        "result": result,
        "result_summary": record.get("result_summary") or {},
    }

    if redis_enabled():
        client = redis_client()
        client.set(share_key(share_id), json.dumps(item, ensure_ascii=False))
        client.zadd(SHARE_INDEX_KEY, {share_id: now})
        prune_shares(now)
        return item

    with SHARE_LOCK:
        SHARE_STORE[share_id] = item
        prune_shares(now)
        return dict(item)


def get_share(share_id):
    share_id = str(share_id or "").strip()
    if not re.fullmatch(r"[0-9a-fA-F]{8,64}", share_id):
        return None
    if redis_enabled():
        raw = redis_client().get(share_key(share_id))
        return json.loads(raw) if raw else None
    with SHARE_LOCK:
        item = SHARE_STORE.get(share_id)
        return dict(item) if item else None


def run_worker_loop():
    if not redis_enabled():
        raise RuntimeError("Worker mode requires NP_REDIS_URL and the redis package.")

    client = redis_client()
    worker_count = max(1, WORKER_CONCURRENCY)
    log_event(
        f"Starting worker loop on Redis queue {JOB_QUEUE_KEY} "
        f"(concurrency={worker_count}, per_job_threads={max(1, PER_JOB_THREAD_LIMIT)})"
    )

    with ThreadPoolExecutor(max_workers=worker_count, thread_name_prefix="np-job") as executor:
        inflight = set()

        while True:
            finished = {future for future in inflight if future.done()}
            if finished:
                inflight.difference_update(finished)
                for future in finished:
                    future.result()

            if len(inflight) >= worker_count:
                done, _ = wait(inflight, return_when=FIRST_COMPLETED, timeout=1)
                inflight.difference_update(done)
                for future in done:
                    future.result()
                continue

            item = client.brpop(JOB_QUEUE_KEY, timeout=2)
            if not item:
                continue
            _, job_id = item
            raw = client.get(job_key(job_id))
            if raw is None:
                continue
            job = json.loads(raw)
            payload = job.get("payload")
            if payload is None:
                continue
            inflight.add(executor.submit(_run_job, job_id, payload))


EXAMPLE_PAYLOAD = {
    "workflow": "analysis",
    "mode": "tube",
    "model": {
        "material": "rna",
        "ensemble": "stacking",
        "celsius": 37,
        "sodium": 1.0,
        "magnesium": 0.0,
    },
    "compute": ["pfunc", "mfe", "pairs"],
    "options": {
        "num_sample": 20,
        "energy_gap": 1.0,
        "sparsity_fraction": 1.0,
        "sparsity_threshold": 0.0,
        "single_mfe": False,
        "result_limit": 25,
        "pairs_preview_size": 24,
    },
    "strands": [
        {"name": "A", "sequence": "AGUCUAGGAU", "concentration": 1.0, "unit": "uM"},
        {"name": "B", "sequence": "UUAACCCACG", "concentration": 2.0, "unit": "uM"},
    ],
    "tube": {"name": "tube1", "max_size": 2},
    "complexes_text": "A\nB\nA+B",
    "design": {
        "trials": 1,
        "result_limit": 25,
        "off_target_max_size": 2,
        "stop_condition": 0.02,
        "seed": 0,
        "wobble_mutations": False,
        "max_time_seconds": 0,
    },
    "design_domains": [
        {"name": "a", "sequence": "N10"},
        {"name": "b", "sequence": "N10"},
    ],
    "design_complexes": [
        {
            "name": "AB_target",
            "strands": "A+B",
            "structure": "(10+)10",
        }
    ],
    "design_tubes": [
        {
            "name": "Tube 1",
            "max_size": 2,
            "on_targets": [
                {"complex": "AB_target", "concentration": 1.0, "unit": "uM"},
            ],
        }
    ],
    "hard_constraints": [],
    "soft_constraints": [],
    "defect_weights": [],
}

# Use A=a and B=~a so the full-duplex target "(10+)10" is base-pair consistent.
DESIGN_TUBE_EXAMPLE_PAYLOAD = {
    "workflow": "design",
    "mode": "tube",
    "model": {
        "material": "rna",
        "ensemble": "stacking",
        "celsius": 37,
        "sodium": 1.0,
        "magnesium": 0.0,
    },
    "compute": ["pfunc", "mfe", "pairs"],
    "options": {
        "num_sample": 20,
        "energy_gap": 1.0,
        "sparsity_fraction": 1.0,
        "sparsity_threshold": 0.0,
        "single_mfe": False,
        "result_limit": 25,
        "pairs_preview_size": 24,
    },
    "strands": [
        {"name": "A", "sequence": "a", "concentration": 1.0, "unit": "uM"},
        {"name": "B", "sequence": "~a", "concentration": 1.0, "unit": "uM"},
    ],
    "tube": {"name": "Design Tube", "max_size": 2},
    "complexes_text": "A+B",
    "design": {
        "trials": 1,
        "result_limit": 25,
        "off_target_max_size": 2,
        "stop_condition": 0.05,
        "seed": 1,
        "wobble_mutations": False,
        "max_time_seconds": 0,
    },
    "design_domains": [
        {"name": "a", "sequence": "N10"},
    ],
    "design_complexes": [
        {
            "name": "AB_target",
            "strands": "A+B",
            "structure": "(10+)10",
        }
    ],
    "design_tubes": [
        {
            "name": "Design Tube",
            "max_size": 2,
            "on_targets": [
                {"complex": "AB_target", "concentration": 1.0, "unit": "uM"},
            ],
        }
    ],
    "hard_constraints": [],
    "soft_constraints": [],
    "defect_weights": [],
}

DESIGN_COMPLEX_EXAMPLE_PAYLOAD = {
    **DESIGN_TUBE_EXAMPLE_PAYLOAD,
    "mode": "complex",
    "tube": {"name": "Design Tube", "max_size": 1},
    "design": {
        **DESIGN_TUBE_EXAMPLE_PAYLOAD["design"],
        "off_target_max_size": 1,
    },
    "design_tubes": [],
}


def get_example_payload(query):
    params = parse_qs(query)
    workflow = (params.get("workflow") or ["analysis"])[0]
    mode = (params.get("mode") or ["tube"])[0]
    if workflow == "design":
        if mode == "complex":
            return DESIGN_COMPLEX_EXAMPLE_PAYLOAD
        return DESIGN_TUBE_EXAMPLE_PAYLOAD
    return EXAMPLE_PAYLOAD


class AppHandler(BaseHTTPRequestHandler):
    server_version = "NPReplica/0.1"

    def do_GET(self):
        parsed = urlparse(self.path)

        if parsed.path == "/":
            self._respond_file(INDEX_PATH, cache_control="public, max-age=60, stale-while-revalidate=86400")
            return

        if parsed.path in {"/favicon.svg", "/favicon.ico"}:
            self._respond_file(FAVICON_PATH, cache_control="public, max-age=86400")
            return

        if parsed.path == "/design-guide.html":
            self._respond_file(GUIDE_PATH, cache_control="public, max-age=3600")
            return

        if parsed.path == "/health":
            stats = job_stats()
            self._respond(
                *json_bytes(
                    {
                        "status": "ok",
                        "job_backend": "redis" if redis_enabled() else "memory",
                        "run_mode": RUN_MODE,
                        "queue_depth": queue_size(),
                        "jobs_running": stats["running"],
                        "jobs_queued": stats["queued"],
                        "job_ttl_seconds": JOB_TTL_SECONDS,
                        "worker_concurrency": max(1, WORKER_CONCURRENCY),
                        "per_job_thread_limit": max(1, PER_JOB_THREAD_LIMIT),
                        "nupack_threads": int(getattr(nupack_config, "threads", PER_JOB_THREAD_LIMIT)),
                        "nupack_cache_gb": float(getattr(nupack_config, "cache", NUPACK_CACHE_GB)),
                    }
                )
            )
            return

        if parsed.path == "/api/example":
            self._respond(*json_bytes(get_example_payload(parsed.query)))
            return

        if parsed.path.startswith("/api/jobs/"):
            job_id = parsed.path.rsplit("/", 1)[-1]
            job = get_job(job_id)
            if job is None:
                self._respond(*json_bytes({"error": "Job not found"}, status=HTTPStatus.NOT_FOUND))
                return
            self._respond(*json_bytes(job))
            return

        if parsed.path.startswith("/api/shares/"):
            share_id = parsed.path.rsplit("/", 1)[-1]
            share = get_share(share_id)
            if share is None:
                self._respond(*json_bytes({"error": "Share not found"}, status=HTTPStatus.NOT_FOUND))
                return
            self._respond(*json_bytes({"status": "success", "share": share}))
            return

        self._respond(*json_bytes({"error": "Not found"}, status=HTTPStatus.NOT_FOUND))

    def do_POST(self):
        parsed = urlparse(self.path)

        if parsed.path != "/api/analyze":
            if parsed.path == "/api/jobs":
                try:
                    length = int(self.headers.get("Content-Length", "0"))
                    raw = self.rfile.read(length)
                    payload = json.loads(raw.decode("utf-8"))
                    job_id = create_job(payload)
                    self._respond(*json_bytes({"status": "accepted", "job_id": job_id}, status=HTTPStatus.ACCEPTED))
                except Exception as exc:
                    self._respond(
                        *json_bytes(
                            {
                                "status": "error",
                                "error": str(exc),
                                "traceback": traceback.format_exc(),
                            },
                            status=HTTPStatus.BAD_REQUEST,
                        )
                    )
                return
            if parsed.path == "/api/shares":
                try:
                    length = int(self.headers.get("Content-Length", "0"))
                    if length > SHARE_MAX_BYTES:
                        raise ValueError(f"Share payload is too large. Limit is {SHARE_MAX_BYTES} bytes.")
                    raw = self.rfile.read(length)
                    record = json.loads(raw.decode("utf-8"))
                    share = create_share(record)
                    self._respond(
                        *json_bytes(
                            {
                                "status": "success",
                                "share_id": share["id"],
                                "url": f"/?share={share['id']}",
                                "max_count": SHARE_MAX_COUNT,
                            },
                            status=HTTPStatus.CREATED,
                        )
                    )
                except Exception as exc:
                    self._respond(
                        *json_bytes(
                            {
                                "status": "error",
                                "error": str(exc),
                                "traceback": traceback.format_exc(),
                            },
                            status=HTTPStatus.BAD_REQUEST,
                        )
                    )
                return
            self._respond(*json_bytes({"error": "Not found"}, status=HTTPStatus.NOT_FOUND))
            return

        try:
            length = int(self.headers.get("Content-Length", "0"))
            raw = self.rfile.read(length)
            payload = json.loads(raw.decode("utf-8"))
            result = run_job_payload(payload)
            self._respond(*json_bytes({"status": "success", "result": result}))
        except Exception as exc:
            self._respond(
                *json_bytes(
                    {
                        "status": "error",
                        "error": str(exc),
                        "traceback": traceback.format_exc(),
                    },
                    status=HTTPStatus.BAD_REQUEST,
                )
            )

    def log_message(self, format_, *args):
        print(f"{self.address_string()} - {format_ % args}")

    def _respond_file(self, path, cache_control):
        etag = file_etag(path)
        if self.headers.get("If-None-Match") == etag:
            self.send_response(HTTPStatus.NOT_MODIFIED)
            self.send_header("ETag", etag)
            self.send_header("Cache-Control", cache_control)
            self.end_headers()
            return
        status, content_type, body = html_bytes(path)
        self._respond(status, content_type, body, cache_control=cache_control, extra_headers={"ETag": etag})

    def _respond(self, status, content_type, body, cache_control=None, extra_headers=None):
        self.send_response(status)
        self.send_header("Content-Type", content_type)
        self.send_header("Content-Length", str(len(body)))
        self.send_header("Cache-Control", cache_control or "no-store, max-age=0")
        if cache_control is None:
            self.send_header("Pragma", "no-cache")
            self.send_header("Expires", "0")
        for key, value in (extra_headers or {}).items():
            self.send_header(key, value)
        self.end_headers()
        self.wfile.write(body)


def main():
    apply_thread_limits()
    if RUN_MODE == "worker":
        run_worker_loop()
        return
    server = ThreadingHTTPServer((HOST, PORT), AppHandler)
    log_event(
        f"Serving NP replica on http://{HOST}:{PORT} "
        f"(jobs via {'redis' if redis_enabled() else 'memory'}, "
        f"per_job_threads={max(1, PER_JOB_THREAD_LIMIT)})"
    )
    server.serve_forever()


if __name__ == "__main__":
    main()