main.py

import json
import math
import random
from copy import deepcopy
from pathlib import Path

import ollama


# ============================================================
# Configuration
# ============================================================

LLM_MODEL = "gpt-oss:20b"

MAX_ROUNDS = 8
STOP_PROBABILITY = 0.95
NOISE_STD = 0.15

MODEL_MISMATCH_THRESHOLD = 2.5

random.seed(42)

OUTPUT_DIR = Path("runs_v04")
OUTPUT_DIR.mkdir(exist_ok=True)


# ============================================================
# Candidate boiling mass-transfer models
# ============================================================

MODEL_REGISTRY = {
    "M0": {
        "description": "linear interfacial mass-transfer closure",
        "base": "linear",
        "a": 0.80,
        "corrections": [],
    },

    "M1": {
        "description": "linear + quadratic mass-transfer closure",
        "base": "linear",
        "a": 0.80,
        "corrections": [
            {
                "type": "quadratic",
                "coefficient": 0.004,
            }
        ],
    },

    "M2": {
        "description": "saturating rational mass-transfer closure",
        "base": "rational",
        "a": 0.80,
        "b": 0.01,
        "corrections": [],
    },
}


def model_to_string(model_name):

    spec = MODEL_REGISTRY[model_name]

    if spec["base"] == "linear":
        expression = f"{spec['a']:.6g} * ΔT"

    elif spec["base"] == "rational":
        expression = (
            f"{spec['a']:.6g} * ΔT "
            f"/ (1 + {spec['b']:.6g} * ΔT)"
        )

    else:
        raise ValueError(
            f"Unknown base model: {spec['base']}"
        )

    for correction in spec["corrections"]:

        if correction["type"] == "quadratic":
            c = correction["coefficient"]
            expression += f" + ({c:.6g}) * ΔT^2"

        elif correction["type"] == "linear":
            c = correction["coefficient"]
            expression += f" + ({c:.6g}) * ΔT"

        elif correction["type"] == "constant":
            c = correction["coefficient"]
            expression += f" + ({c:.6g})"

    return expression


def physics_model(model_name, delta_T):

    spec = MODEL_REGISTRY[model_name]

    if spec["base"] == "linear":

        y = spec["a"] * delta_T

    elif spec["base"] == "rational":

        y = (
            spec["a"] * delta_T
            / (1.0 + spec["b"] * delta_T)
        )

    else:

        raise ValueError(
            f"Unknown base model: {spec['base']}"
        )

    for correction in spec["corrections"]:

        correction_type = correction["type"]
        coefficient = correction["coefficient"]

        if correction_type == "quadratic":

            y += coefficient * delta_T**2

        elif correction_type == "linear":

            y += coefficient * delta_T

        elif correction_type == "constant":

            y += coefficient

        else:

            raise ValueError(
                f"Unknown correction: {correction_type}"
            )

    return y


# ============================================================
# Hidden boiling physical world
# ============================================================

def hidden_physics(delta_T: float) -> float:
    """
    Synthetic boiling interfacial mass-transfer world.

    The intelligence system never sees this equation.

    The hidden world contains a nonlinear mass-transfer
    contribution that is not represented exactly by the
    initial candidate model class.
    """

    return (
        0.80 * delta_T
        + 0.002 * delta_T**2
    )


def query_hidden_world(delta_T: float) -> dict:

    clean = hidden_physics(delta_T)

    noise = random.gauss(
        0.0,
        NOISE_STD,
    )

    return {
        "delta_T": delta_T,
        "observed_mass_transfer": clean + noise,
        "noise_std": NOISE_STD,
    }


# ============================================================
# Scientific state
# ============================================================

state = {

    "scientific_question": (
        "Determine an adequate constitutive closure for "
        "boiling interfacial mass transfer as a function "
        "of interfacial thermal driving ΔT. "
        "Detect failure of the initial closure class and "
        "construct a revised executable closure if required."
    ),

    "physical_quantity": (
        "normalized interfacial mass-transfer response"
    ),

    "candidate_models": {},

    "allowed_delta_T": [
        2,
        5,
        8,
        12,
        16,
        20,
        24,
        28,
    ],

    "evidence": [],

    "posterior": {},

    "round": 0,

    "revision_history": [],
}


def synchronize_state_models():

    state["candidate_models"] = {
        model_name: model_to_string(model_name)
        for model_name in MODEL_REGISTRY
    }


def reset_posterior():

    n_models = len(MODEL_REGISTRY)

    state["posterior"] = {
        model_name: 1.0 / n_models
        for model_name in MODEL_REGISTRY
    }


synchronize_state_models()
reset_posterior()


# ============================================================
# Deterministic physics tools
# ============================================================

def prediction_table(state):

    table = {}

    for delta_T in state["allowed_delta_T"]:

        table[delta_T] = {}

        for model in state["candidate_models"]:

            table[delta_T][model] = physics_model(
                model,
                delta_T,
            )

    return table


def discrimination_scores(state):
    """
    Rank unused thermal conditions according to the minimum
    pairwise separation between executable mass-transfer
    closures, normalized by observational noise.
    """

    used = {
        obs["delta_T"]
        for obs in state["evidence"]
    }

    scores = {}

    for delta_T in state["allowed_delta_T"]:

        if delta_T in used:
            continue

        predictions = [
            physics_model(
                model,
                delta_T,
            )
            for model in state["candidate_models"]
        ]

        if len(predictions) < 2:
            scores[delta_T] = 0.0
            continue

        pairwise = []

        for i in range(len(predictions)):

            for j in range(
                i + 1,
                len(predictions),
            ):

                separation = abs(
                    predictions[i]
                    - predictions[j]
                )

                pairwise.append(
                    separation / NOISE_STD
                )

        scores[delta_T] = min(pairwise)

    return scores


# ============================================================
# Bayesian evidence update
# ============================================================

def gaussian_log_likelihood(
    observed,
    predicted,
    sigma,
):

    z = (
        observed - predicted
    ) / sigma

    return -0.5 * z**2


def update_posterior(
    state,
    observation,
):

    old = state["posterior"]

    log_weights = {}

    for model in state["candidate_models"]:

        prediction = physics_model(
            model,
            observation["delta_T"],
        )

        log_likelihood = gaussian_log_likelihood(
            observation["observed_mass_transfer"],
            prediction,
            observation["noise_std"],
        )

        prior = max(
            old.get(model, 1e-300),
            1e-300,
        )

        log_weights[model] = (
            math.log(prior)
            + log_likelihood
        )

    max_log_weight = max(
        log_weights.values()
    )

    weights = {
        model: math.exp(
            value - max_log_weight
        )
        for model, value
        in log_weights.items()
    }

    normalizer = sum(
        weights.values()
    )

    return {
        model: value / normalizer
        for model, value
        in weights.items()
    }


def recompute_posterior_from_all_evidence():

    reset_posterior()

    for observation in state["evidence"]:

        state["posterior"] = (
            update_posterior(
                state,
                observation,
            )
        )


# ============================================================
# Model-class adequacy
# ============================================================

def model_mismatch_scores(state):

    if len(state["evidence"]) < 3:
        return {}

    scores = {}

    for model in state["candidate_models"]:

        residuals = []

        for obs in state["evidence"]:

            predicted = physics_model(
                model,
                obs["delta_T"],
            )

            residual = (
                obs["observed_mass_transfer"]
                - predicted
            )

            residuals.append(
                residual
            )

        rmse = math.sqrt(
            sum(
                r**2
                for r in residuals
            )
            / len(residuals)
        )

        scores[model] = (
            rmse / NOISE_STD
        )

    return scores


def best_model_by_mismatch(state):

    scores = model_mismatch_scores(
        state
    )

    if not scores:
        return None, None

    best_model = min(
        scores,
        key=scores.get,
    )

    return (
        best_model,
        scores[best_model],
    )


# ============================================================
# Residual analysis
# ============================================================

def residual_table(
    state,
    baseline_model,
):

    table = []

    for obs in state["evidence"]:

        prediction = physics_model(
            baseline_model,
            obs["delta_T"],
        )

        residual = (
            obs["observed_mass_transfer"]
            - prediction
        )

        table.append({
            "delta_T":
                obs["delta_T"],

            "observed_mass_transfer":
                obs["observed_mass_transfer"],

            "baseline_prediction":
                prediction,

            "residual":
                residual,
        })

    return table


def fit_constant_correction(
    state,
    baseline_model,
):

    residuals = []

    for obs in state["evidence"]:

        x = obs["delta_T"]

        residuals.append(
            obs["observed_mass_transfer"]
            - physics_model(
                baseline_model,
                x,
            )
        )

    return (
        sum(residuals)
        / len(residuals)
    )


def fit_linear_correction(
    state,
    baseline_model,
):

    numerator = 0.0
    denominator = 0.0

    for obs in state["evidence"]:

        x = obs["delta_T"]

        residual = (
            obs["observed_mass_transfer"]
            - physics_model(
                baseline_model,
                x,
            )
        )

        numerator += residual * x
        denominator += x**2

    if denominator == 0:
        return 0.0

    return numerator / denominator


def fit_quadratic_correction(
    state,
    baseline_model,
):

    numerator = 0.0
    denominator = 0.0

    for obs in state["evidence"]:

        x = obs["delta_T"]

        residual = (
            obs["observed_mass_transfer"]
            - physics_model(
                baseline_model,
                x,
            )
        )

        numerator += (
            residual * x**2
        )

        denominator += x**4

    if denominator == 0:
        return 0.0

    return (
        numerator / denominator
    )


def correction_fit_rmse(
    state,
    baseline_model,
    correction_type,
    coefficient,
):

    errors = []

    for obs in state["evidence"]:

        x = obs["delta_T"]

        baseline = physics_model(
            baseline_model,
            x,
        )

        if correction_type == "constant":

            correction = coefficient

        elif correction_type == "linear":

            correction = (
                coefficient * x
            )

        elif correction_type == "quadratic":

            correction = (
                coefficient * x**2
            )

        else:

            raise ValueError(
                correction_type
            )

        prediction = (
            baseline + correction
        )

        errors.append(
            obs["observed_mass_transfer"]
            - prediction
        )

    return math.sqrt(
        sum(
            error**2
            for error in errors
        )
        / len(errors)
    )


def deterministic_revision_search(
    state,
    baseline_model,
):
    """
    Search a deliberately small mathematical correction space.

    The LLM may interpret the residual structure, but the
    executable revised closure is selected and fitted by
    deterministic numerical tools.
    """

    candidates = {}

    constant_c = fit_constant_correction(
        state,
        baseline_model,
    )

    candidates["constant"] = {
        "coefficient": constant_c,
        "rmse": correction_fit_rmse(
            state,
            baseline_model,
            "constant",
            constant_c,
        ),
    }

    linear_c = fit_linear_correction(
        state,
        baseline_model,
    )

    candidates["linear"] = {
        "coefficient": linear_c,
        "rmse": correction_fit_rmse(
            state,
            baseline_model,
            "linear",
            linear_c,
        ),
    }

    quadratic_c = fit_quadratic_correction(
        state,
        baseline_model,
    )

    candidates["quadratic"] = {
        "coefficient": quadratic_c,
        "rmse": correction_fit_rmse(
            state,
            baseline_model,
            "quadratic",
            quadratic_c,
        ),
    }

    best_type = min(
        candidates,
        key=lambda name:
            candidates[name]["rmse"],
    )

    return (
        best_type,
        candidates[best_type]["coefficient"],
        candidates,
    )


# ============================================================
# Executable theory revision
# ============================================================

def register_revised_model(
    parent_model,
    correction_type,
    coefficient,
):

    new_index = len(
        MODEL_REGISTRY
    )

    new_model_name = (
        f"M{new_index}"
    )

    new_spec = deepcopy(
        MODEL_REGISTRY[parent_model]
    )

    new_spec["description"] = (
        f"revised boiling mass-transfer closure "
        f"derived from {parent_model}"
    )

    new_spec["corrections"].append({
        "type":
            correction_type,

        "coefficient":
            coefficient,
    })

    MODEL_REGISTRY[
        new_model_name
    ] = new_spec

    synchronize_state_models()

    return new_model_name


# ============================================================
# LLM interface
# ============================================================

def ask_agent(
    system_prompt,
    user_prompt,
):

    response = ollama.chat(
        model=LLM_MODEL,
        messages=[
            {
                "role": "system",
                "content": system_prompt,
            },
            {
                "role": "user",
                "content": user_prompt,
            },
        ],
    )

    return (
        response["message"]["content"]
    )


# ============================================================
# Scientific agents
# ============================================================

def proposer(
    state,
    table,
):

    prompt = f"""
SCIENTIFIC STATE

{json.dumps(state, indent=2)}

EXECUTABLE BOILING MASS-TRANSFER MODEL PREDICTIONS

{json.dumps(table, indent=2)}

You are the scientific hypothesis proposer.

The problem is constitutive modeling of interfacial mass
transfer in boiling.

The numerical predictions were computed by an external
physics tool and are authoritative.

Do NOT perform new arithmetic.
Do NOT invent additional models.

Using the current evidence:

1. identify which mass-transfer closures remain plausible,
2. explain their mathematical differences,
3. state what uncertainty remains.

Do not assign a microscopic boiling mechanism to a
mathematical term unless the evidence identifies it.

Be concise.
"""

    return ask_agent(
        "You are a rigorous boiling-physics hypothesis agent.",
        prompt,
    )


def critic(
    state,
    proposal,
    scores,
):

    prompt = f"""
SCIENTIFIC STATE

{json.dumps(state, indent=2)}

PROPOSER

{proposal}

COMPUTED TEST DISCRIMINATION SCORES

{json.dumps(scores, indent=2)}

You are an independent scientific critic evaluating
candidate boiling interfacial mass-transfer closures.

The numerical values were computed externally.
Do not recompute them.

Assess:

1. overclaiming,
2. evidence sufficiency,
3. surviving alternatives,
4. whether another thermal condition is required,
5. whether preference among existing closures could hide
   model-class inadequacy.

Do not invent microscopic boiling physics.
"""

    return ask_agent(
        "You are a skeptical boiling-physics reviewer.",
        prompt,
    )


def select_next_test(
    state,
    proposal,
    critique,
    scores,
):

    ranked = sorted(
        scores.items(),
        key=lambda x: x[1],
        reverse=True,
    )

    prompt = f"""
CURRENT SCIENTIFIC STATE

{json.dumps(state, indent=2)}

PROPOSER

{proposal}

CRITIC

{critique}

AVAILABLE THERMAL CONDITIONS RANKED BY
COMPUTED MODEL DISCRIMINATION

{json.dumps(ranked, indent=2)}

Select ONE available Delta T condition for the next
synthetic boiling mass-transfer observation.

Return ONLY the numerical Delta T value.
"""

    answer = ask_agent(
        "You select informative falsification tests.",
        prompt,
    )

    allowed = list(
        scores.keys()
    )

    for value in sorted(
        allowed,
        reverse=True,
    ):

        if str(value) in answer:
            return value

    return ranked[0][0]


def theory_revision_agent(
    state,
    best_model,
    mismatch_score,
    residuals,
    revision_candidates,
):

    prompt = f"""
SCIENTIFIC STATE

{json.dumps(state, indent=2)}

BEST CURRENT EXECUTABLE MASS-TRANSFER CLOSURE

{best_model}

NORMALIZED MODEL-MISMATCH SCORE

{mismatch_score}

COMPUTED RESIDUALS

{json.dumps(residuals, indent=2)}

DETERMINISTIC REVISION SEARCH

{json.dumps(revision_candidates, indent=2)}

The current model class is inadequate.

You are the theory-revision component of an autonomous
boiling-physics discovery system.

The numerical fitting was performed by external tools.
Do NOT recompute coefficients.

Interpret the evidence.

Answer concisely:

MATHEMATICAL INFERENCE:
What residual structure is supported?

MODEL REVISION:
What minimal constitutive correction is justified?

PHYSICAL HYPOTHESIS:
What, if anything, can be inferred about missing
interfacial mass-transfer physics?

FALSIFICATION:
What observation would most strongly challenge the
revised closure?

Do not claim a specific microscopic boiling mechanism
unless the evidence actually identifies one.
"""

    return ask_agent(
        (
            "You revise inadequate boiling mass-transfer "
            "closures using falsifiable numerical evidence."
        ),
        prompt,
    )


def final_evaluator(state):

    mismatch_scores = (
        model_mismatch_scores(state)
    )

    prompt = f"""
FINAL SCIENTIFIC STATE

{json.dumps(state, indent=2)}

FINAL NORMALIZED MODEL-MISMATCH SCORES

{json.dumps(mismatch_scores, indent=2)}

You are the final scientific evaluator.

This is a synthetic benchmark for autonomous discovery
of a boiling interfacial mass-transfer closure.

Use only the supplied numerical evidence.

Answer:

1. Which executable closure is best supported?
2. Is it adequate within the observational uncertainty?
3. Was the original model class falsified?
4. Was a revised executable closure generated?
5. What should be tested next?

Do not invent microscopic physics.
"""

    return ask_agent(
        "You are an evidence-based scientific judge.",
        prompt,
    )


# ============================================================
# Persistent scientific memory
# ============================================================

trace = []


def save_state():

    with open(
        OUTPUT_DIR / "scientific_state.json",
        "w",
    ) as f:

        json.dump(
            state,
            f,
            indent=2,
        )

    with open(
        OUTPUT_DIR / "trace.json",
        "w",
    ) as f:

        json.dump(
            trace,
            f,
            indent=2,
        )

    with open(
        OUTPUT_DIR / "model_registry.json",
        "w",
    ) as f:

        json.dump(
            MODEL_REGISTRY,
            f,
            indent=2,
        )


# ============================================================
# Closed self-revising discovery loop
# ============================================================

print("\n" + "=" * 72)
print("BOILING INTELLIGENCE v0.4")
print("Executable Mass-Transfer Closure Discovery")
print("=" * 72)


original_model_class_failed = False
revision_generated = False
revision_count = 0

MAX_REVISIONS = 1


for round_id in range(
    1,
    MAX_ROUNDS + 1,
):

    state["round"] = round_id

    print("\n" + "=" * 72)
    print(f"ROUND {round_id}")
    print("=" * 72)

    table = prediction_table(
        state
    )

    scores = discrimination_scores(
        state
    )

    if not scores:

        print(
            "\nNo unused thermal conditions remain."
        )

        break


    # --------------------------------------------------------
    # 1. Hypothesis proposer
    # --------------------------------------------------------

    print("\n[1] HYPOTHESIS PROPOSER")

    proposal = proposer(
        state,
        table,
    )

    print(proposal)


    # --------------------------------------------------------
    # 2. Critic
    # --------------------------------------------------------

    print("\n[2] SCIENTIFIC CRITIC")

    critique = critic(
        state,
        proposal,
        scores,
    )

    print(critique)


    # --------------------------------------------------------
    # 3. Falsification-test designer
    # --------------------------------------------------------

    print("\n[3] TEST DESIGNER")

    delta_T = select_next_test(
        state,
        proposal,
        critique,
        scores,
    )

    print(
        f"Selected ΔT = {delta_T}"
    )


    # --------------------------------------------------------
    # 4. Synthetic boiling physical world
    # --------------------------------------------------------

    print(
        "\n[4] EXECUTABLE BOILING WORLD"
    )

    observation = query_hidden_world(
        delta_T
    )

    print(
        "Observed normalized mass transfer = "
        f"{observation['observed_mass_transfer']:.6f}"
    )


    # --------------------------------------------------------
    # 5. Evidence update
    # --------------------------------------------------------

    print("\n[5] EVIDENCE UPDATE")

    state["evidence"].append(
        observation
    )

    state["posterior"] = (
        update_posterior(
            state,
            observation,
        )
    )

    for model, probability in sorted(
        state["posterior"].items(),
        key=lambda x: x[1],
        reverse=True,
    ):

        print(
            f"{model}: "
            f"P = {probability:.4f}"
        )


    best_posterior_model = max(
        state["posterior"],
        key=state["posterior"].get,
    )

    best_probability = (
        state["posterior"][
            best_posterior_model
        ]
    )


    # --------------------------------------------------------
    # 6. Model adequacy
    # --------------------------------------------------------

    best_model, mismatch = (
        best_model_by_mismatch(
            state
        )
    )

    if mismatch is not None:

        print(
            "\nBest closure by adequacy: "
            f"{best_model}"
        )

        print(
            "Normalized mismatch = "
            f"{mismatch:.3f} sigma"
        )


    round_record = {

        "round":
            round_id,

        "proposal":
            proposal,

        "critique":
            critique,

        "test_scores":
            scores,

        "selected_delta_T":
            delta_T,

        "observation":
            observation,

        "posterior":
            state["posterior"].copy(),

        "best_model_by_adequacy":
            best_model,

        "mismatch_score":
            mismatch,
    }

    trace.append(
        round_record
    )

    save_state()


    print(
        "\nCurrent Bayesian best candidate:",
        best_posterior_model,
        f"(P={best_probability:.4f})",
    )


    # --------------------------------------------------------
    # 7. Open-set model-class failure
    # --------------------------------------------------------

    if (
        mismatch is not None
        and mismatch >= MODEL_MISMATCH_THRESHOLD
        and revision_count < MAX_REVISIONS
    ):

        print("\n" + "=" * 72)
        print("MODEL-CLASS FAILURE DETECTED")
        print("=" * 72)

        original_model_class_failed = True

        residuals = residual_table(
            state,
            best_model,
        )

        (
            correction_type,
            coefficient,
            revision_candidates,
        ) = deterministic_revision_search(
            state,
            best_model,
        )

        print(
            "\nDeterministic residual search:"
        )

        for name, result in (
            revision_candidates.items()
        ):

            print(
                f"{name:10s} "
                f"c={result['coefficient']:.6g} "
                f"RMSE={result['rmse']:.6g}"
            )


        print(
            "\n[6] THEORY REVISION AGENT"
        )

        revision_text = (
            theory_revision_agent(
                state,
                best_model,
                mismatch,
                residuals,
                revision_candidates,
            )
        )

        print(revision_text)


        # ----------------------------------------------------
        # 8. Compile revision into executable closure
        # ----------------------------------------------------

        print(
            "\n[7] EXECUTABLE MODEL REVISION"
        )

        new_model = register_revised_model(
            best_model,
            correction_type,
            coefficient,
        )

        revision_count += 1
        revision_generated = True

        revision_record = {

            "parent_model":
                best_model,

            "new_model":
                new_model,

            "correction_type":
                correction_type,

            "coefficient":
                coefficient,

            "equation":
                model_to_string(
                    new_model
                ),

            "agent_interpretation":
                revision_text,
        }

        state[
            "revision_history"
        ].append(
            revision_record
        )

        trace.append({
            "event":
                "executable_model_revision",

            **revision_record,
        })


        print(
            f"Created {new_model}"
        )

        print(
            "Executable closure:"
        )

        print(
            f"{new_model}: "
            f"{model_to_string(new_model)}"
        )


        # ----------------------------------------------------
        # Re-evaluate all accumulated evidence with M3
        # ----------------------------------------------------

        recompute_posterior_from_all_evidence()

        save_state()


        print(
            "\nRe-entering revised closure "
            "into falsification loop."
        )

        continue


    # --------------------------------------------------------
    # Stop only when revised model has survived testing
    # --------------------------------------------------------

    if (
        revision_generated
        and mismatch is not None
        and mismatch < MODEL_MISMATCH_THRESHOLD
        and round_id >= 4
    ):

        print("\n" + "=" * 72)
        print("REVISED CLOSURE SURVIVES FALSIFICATION")
        print("=" * 72)

        print(
            f"{best_model}: "
            f"{model_to_string(best_model)}"
        )

        break


# ============================================================
# Final assessment
# ============================================================

print("\n" + "=" * 72)
print("FINAL EVALUATION")
print("=" * 72)

assessment = final_evaluator(
    state
)

print(
    assessment
)


with open(
    OUTPUT_DIR / "final_assessment.txt",
    "w",
) as f:

    f.write(
        assessment
    )


with open(
    OUTPUT_DIR / "discovered_models.txt",
    "w",
) as f:

    for model_name in MODEL_REGISTRY:

        f.write(
            f"{model_name}: "
            f"{model_to_string(model_name)}\n"
        )


print("\n" + "=" * 72)

if revision_generated:

    print(
        "SELF-REVISING BOILING DISCOVERY LOOP COMPLETE"
    )

else:

    print(
        "BOILING DISCOVERY LOOP COMPLETE"
    )

print("=" * 72)