fix: correct COI formulation to measure price erosion over time

The fundamental error was treating COI as instantaneous margin × alpha.
The corrected formulation is:

    COI = E[p_start] - p_transaction

This measures price erosion over time, capturing how agents using
multiple sessions gather information and drive prices down.

Key changes:
- Add coi.py with COIWindow, COITracker, and compute_multi_session_coi
- Add separability.py with KL-divergence behavioral classification
- Update simplified_env.py to track initial prices and compute windowed COI
- Add corrected COI metrics (coi_*_corrected) alongside legacy metrics

The new approach:
1. Tracks prices at episode start as E[p] (expected price)
2. Computes transaction prices as p (actual sale price)
3. Measures leak as the difference (price erosion)
4. Includes order statistic erosion (Theorem 1: N agents -> min price)

.gitignore

@@ -22,5 +22,3 @@ sim/rl/behavior_loader/*.png
 sim/rl/behavior_loader/*.svg
 sim/rl/behavior_loader/*.pdf
 tests/e2e/node_modules/**
-lab/case/thesis/runs*/
-sim/case/thesis_simplified/runs*/

engine/engine.py

@@ -1,66 +0,0 @@
-from sys import platform
-import numpy as np
-from .lib.demand import generate_demand, estimate_demand
-from .lib.behavior import sample_behavior
-from logging import INFO, getLogger
-logger = getLogger(__name__)
-logger.setLevel(INFO)
-
-
-
-class MarketEngine():
-    def __init__(self,
-                 alpha = 0.5,
-                 N = 100,
-                 demand_distribution = (50, 10),
-                 demand_sampling_function = np.random.normal):
-        self.Nagents = int(N*alpha)
-        self.Nhumans = int(N*(1-alpha))
-        self.demand = (demand_sampling_function, demand_distribution)
-
-    def act(self, prices):
-        demand = generate_demand(prices, *self.demand)
-        sample_n = lambda n, human: [sample_behavior(demand, human=human) for _ in range(n)]
-        human_t, agent_t = sample_n(100, True), sample_n(100, False)
-        trajectories = human_t + agent_t
-        demand_estimate = estimate_demand(trajectories)
-        return demand_estimate
-
-    def measure(self):
-        pass
-
-class PricingEngine():
-    def __init__(self,
-                 ) -> None:
-        pass
-
-    def act(self, demand):
-        return np.random.uniform(low=25, high=100, size=10)
-
-
-
-class Limbo():
-    def __init__(self,
-                 platform,
-                 market
-                 ) -> None:
-        self.platform_turn = True
-        self.platform = platform
-        self.market = market
-        self.output = None
-
-    def step(self):
-        # we could code golf this a little bit
-        if self.platform_turn:
-            self.output = self.platform.act(self.output)
-        else:
-            self.output = self.market.act(self.output)
-        print(self.output)
-        self.platform_turn = not self.platform_turn
-
-if __name__ == "__main__":
-    platform = PricingEngine()
-    market = MarketEngine()
-    limbo = Limbo(platform, market)
-    for _ in range(10):
-        limbo.step()

engine/lib/__init__.py

@@ -1,3 +0,0 @@
-from .demand import generate_demand, estimate_demand
-from .behavior import sample_behavior
-from .render import DashboardRenderer, style_axis

engine/lib/behavior.py

@@ -1,47 +0,0 @@
-from sim.rl.behavior_loader.models import BehaviorModel, AgentBehaviorModel, aggregate_event_transitions
-import pandas as pd
-import numpy as np
-from .demand import generate_demand
-
-base_dir = "/home/velocitatem/Documents/Projects/PHANTOM/experiments"
-human_dir, agent_dir = f"{base_dir}/collected_data/", f"{base_dir}/agents/collected_data/"
-
-_cache = {}  # lazy cache for models and base pivots
-
-def _get_base_pivot(human: bool):
-    key = 'human' if human else 'agent'
-    if key not in _cache:
-        model = BehaviorModel(human_dir) if human else AgentBehaviorModel(agent_dir)
-        mdp = model.build_MDP()
-        _cache[key] = pd.DataFrame(aggregate_event_transitions(mdp)).fillna(0.0)
-    return _cache[key]
-
-def adjust_behavior_to_condition(condition, transition_matrix):
-    # expand NxN transition matrix to (N*P)x(N*P) weighted by demand condition
-    cond_norm = condition / np.sum(condition)
-    n_products = len(condition)
-    base_vals = transition_matrix.values
-    base_cols, base_rows = transition_matrix.columns.tolist(), transition_matrix.index.tolist()
-
-    # expand via kronecker-like tiling: each cell becomes a P*P block weighted by outer product of cond_norm
-    expanded = np.kron(base_vals, np.outer(cond_norm, cond_norm))
-    new_cols = [f"{c}_product{p}" for c in base_cols for p in range(n_products)]
-    new_rows = [f"{r}_product{p}" for r in base_rows for p in range(n_products)]
-    return pd.DataFrame(expanded, index=new_rows, columns=new_cols)
-
-def sample_behavior(condition, human=True, max_len=40):
-    base_pivot = _get_base_pivot(human)
-    adjusted_transitions = adjust_behavior_to_condition(condition, base_pivot)
-
-    trajectory = [np.random.choice(adjusted_transitions.index)]
-    while len(trajectory) < max_len or 'checkout' in trajectory[-1]:
-        probs = adjusted_transitions.loc[trajectory[-1]].values
-        sample = np.random.choice(adjusted_transitions.columns, p=probs/np.sum(probs) if np.sum(probs) > 0 else None)
-        trajectory.append(sample)
-    return trajectory
-
-if __name__ == "__main__":
-    t=sample_behavior(generate_demand(np.array([10,20,30])), human=True)
-    print(t)
-    t=sample_behavior(generate_demand(np.array([10,20,30])), human=False)
-    print(t)

engine/lib/demand.py

@@ -1,45 +0,0 @@
-import logging
-import numpy as np
-from logging import getLogger
-logger = getLogger(__name__)
-
-def generate_demand(prices, distribution_method = np.random.normal, distribution_params = (50.0, 10.0)):
-    # assumption 1: each product has an intrinsic valuation drawn from a normal distribution centered at 50
-    product_valuations = distribution_method(*distribution_params, size=len(prices))
-    # assumption 2: demand decreases as price increases, following a simple linear model
-    demand = np.maximum(0, product_valuations - prices) # demand cannot be negative
-    total = np.sum(demand)
-    demand = demand / total * 100 if total > 0 else demand  # normalize to percentage, avoid div by zero
-    logger.info(f"Generated demand for prices {prices}: {demand} with valuations from distribution {distribution_params}")
-    return demand
-
-def estimate_demand(trajectories):
-    demand_estimate = {}
-    for traj in trajectories:
-        for event in traj:
-            if 'view_product' in event:
-                product_id = int(event.split('_')[-1].replace('product', ''))
-                demand_estimate[product_id] = demand_estimate.get(product_id, 0) + 1
-    total_views = sum(demand_estimate.values())
-    for product_id in demand_estimate:
-        demand_estimate[product_id] = (demand_estimate[product_id] / total_views) * 100  # normalize to percentage
-    return demand_estimate
-
-# Example usage
-if __name__ == "__main__":
-    np.random.seed(42)
-    prices = np.array([20.0, 35.0, 50.0, 65.0])
-    demand = generate_demand(prices)
-    print("Generated Demand:", demand)
-    from .behavior import sample_behavior
-    N, alphat =200, 0.1
-    trajectories = []
-    for _ in range(int(N*(1 - alphat))):
-        trajectories.append(sample_behavior(demand, human=True))
-    for _ in range(int(N*alphat)):
-        trajectories.append(sample_behavior(demand, human=False))
-    demand_estimate = estimate_demand(trajectories)
-    print("Estimated Demand from Behavior:", demand_estimate)
-    delta = {k: demand_estimate.get(k, 0) - demand[i] for i, k in enumerate(range(len(prices)))}
-    delta = np.mean([np.abs(v) for v in delta.values()])
-    print("Demand Delta:", delta)

engine/lib/render.py

@@ -1,126 +0,0 @@
-"""rendering logic for PHANTOM environment dashboard"""
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib.gridspec import GridSpec
-
-
-def style_axis(ax, title: str = None, xlabel: str = None, ylabel: str = None):
-    ax.spines['top'].set_visible(False)
-    ax.spines['right'].set_visible(False)
-    if title: ax.set_title(title, fontsize=11, fontweight='bold', pad=8)
-    if xlabel: ax.set_xlabel(xlabel, fontsize=9)
-    if ylabel: ax.set_ylabel(ylabel, fontsize=9)
-
-
-class DashboardRenderer:
-    """stateful renderer for PHANTOM market dynamics visualization"""
-
-    def __init__(self):
-        self.fig = None
-        self.gs = None
-
-    def render(self, env) -> None:
-        if self.fig is None:
-            plt.ion()
-            self.fig = plt.figure(figsize=(14, 10))
-            self.gs = GridSpec(3, 3, figure=self.fig, hspace=0.35, wspace=0.3,
-                               left=0.07, right=0.95, top=0.92, bottom=0.08)
-            plt.show(block=False)
-
-        self.fig.clear()
-        self.fig.suptitle(f'PHANTOM  Market Dynamics  [t={env._step_count}, a={env.alpha:.2f}]',
-                          fontsize=14, fontweight='bold')
-
-        demand_mat = np.array(env._demand_history).T
-        price_mat = np.array(env._price_history).T
-        elasticity = env._compute_elasticity()
-
-        self._render_scatter(env)
-        self._render_elasticity_bar(env, elasticity)
-        self._render_session_pie(env)
-        self._render_price_heatmap(price_mat)
-        self._render_demand_heatmap(demand_mat)
-        self._render_correlation(env.n_products, price_mat, demand_mat)
-        self._render_revenue(env)
-
-        self.fig.canvas.draw_idle()
-        self.fig.canvas.flush_events()
-
-    def _render_scatter(self, env):
-        ax = self.fig.add_subplot(self.gs[0, 0])
-        prices_flat = np.array(env._price_history).flatten()
-        demands_flat = np.array(env._demand_history).flatten()
-        product_ids = np.tile(np.arange(env.n_products), len(env._price_history))
-        ax.scatter(prices_flat, demands_flat, c=product_ids, cmap='plasma', alpha=0.6, s=15, edgecolors='none')
-        if len(prices_flat) > 1:
-            z = np.polyfit(prices_flat, demands_flat, 1)
-            p_line = np.linspace(prices_flat.min(), prices_flat.max(), 50)
-            ax.plot(p_line, np.polyval(z, p_line), '--', lw=1.5, alpha=0.8)
-        style_axis(ax, "Price-Demand Relationship", "Price ($)", "Demand")
-
-    def _render_elasticity_bar(self, env, elasticity):
-        ax = self.fig.add_subplot(self.gs[0, 1])
-        ax.barh(range(env.n_products), elasticity, alpha=0.8)
-        ax.axvline(0, lw=0.8, alpha=0.5)
-        ax.axvline(-1, lw=1, ls='--', alpha=0.5)
-        ax.set_yticks(range(env.n_products))
-        ax.set_yticklabels([f'P{i}' for i in range(env.n_products)], fontsize=7)
-        style_axis(ax, "Price Elasticity", "(dQ/dP)(P/Q)", None)
-
-    def _render_session_pie(self, env):
-        ax = self.fig.add_subplot(self.gs[0, 2])
-        n_h, n_a = env.market.Nhumans, env.market.Nagents
-        wedges, _ = ax.pie([n_h, n_a], startangle=90, wedgeprops={'linewidth': 2, 'edgecolor': 'white'})
-        ax.legend(wedges, [f'H ({n_h})', f'A ({n_a})'], loc='lower center', fontsize=8,
-                  frameon=False, bbox_to_anchor=(0.5, -0.05))
-        ax.set_title("Session Mix", fontsize=11, fontweight='bold')
-
-    def _render_price_heatmap(self, price_mat):
-        ax = self.fig.add_subplot(self.gs[1, :2])
-        im = ax.imshow(price_mat, aspect='auto', cmap='viridis', origin='lower')
-        style_axis(ax, "Price Heatmap P(product, t)", "Step", "Product")
-        cbar = self.fig.colorbar(im, ax=ax, fraction=0.03, pad=0.02)
-        cbar.set_label('$', fontsize=8)
-
-    def _render_demand_heatmap(self, demand_mat):
-        ax = self.fig.add_subplot(self.gs[1, 2])
-        im = ax.imshow(demand_mat, aspect='auto', cmap='Blues', origin='lower')
-        style_axis(ax, "Demand Q(product, t)", "Step", None)
-        self.fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02)
-
-    def _render_correlation(self, n_products, price_mat, demand_mat):
-        ax = self.fig.add_subplot(self.gs[2, 0])
-        if price_mat.shape[1] > 2:
-            corr = np.corrcoef(price_mat, demand_mat)[:n_products, n_products:]
-            im = ax.imshow(corr, cmap='RdBu', vmin=-1, vmax=1, aspect='auto')
-            ax.set_xticks(range(n_products))
-            ax.set_yticks(range(n_products))
-            ax.set_xticklabels([f'Q{i}' for i in range(n_products)], fontsize=6)
-            ax.set_yticklabels([f'P{i}' for i in range(n_products)], fontsize=6)
-            self.fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02)
-        style_axis(ax, "Price-Demand Correlation", None, None)
-
-    def _render_revenue(self, env):
-        ax = self.fig.add_subplot(self.gs[2, 1:])
-        n_steps = len(env._revenue_history)
-        demand_std = [np.std(d) for d in env._demand_history]
-        ax.fill_between(range(n_steps), env._revenue_history, alpha=0.3)
-        ax.plot(env._revenue_history, linewidth=2, label='Revenue')
-        ax.set_xlim(0, max(n_steps, 1))
-        ax.set_ylim(0, max(env._revenue_history) * 1.1 if env._revenue_history else 1)
-
-        ax2 = ax.twinx()
-        ax2.plot(range(n_steps), demand_std, linewidth=2, ls='-', alpha=0.9, label='sigma(Demand)')
-        d_min, d_max = min(demand_std), max(demand_std)
-        margin = (d_max - d_min) * 0.2 if d_max > d_min else 0.5
-        ax2.set_ylim(max(0, d_min - margin), d_max + margin)
-        ax2.set_ylabel('Demand sigma', fontsize=9)
-
-        style_axis(ax, "Revenue & Demand Dispersion", "Step", "Revenue ($)")
-        ax.legend(loc='upper left', fontsize=7, frameon=False)
-        ax2.legend(loc='upper right', fontsize=7, frameon=False)
-
-    def close(self):
-        if self.fig:
-            plt.close(self.fig)
-            self.fig = None

engine/studies/factors.py

@@ -1,34 +0,0 @@
-"""shared factor definitions for experimental designs"""
-import numpy as np
-from dataclasses import dataclass, field
-from typing import Callable, Any
-
-@dataclass
-class Factor:
-    name: str
-    levels: list
-    primary: bool = True  # full cross vs sampled
-
-# demand functions with compatible signatures
-def demand_linear(mu, sigma, size): return np.maximum(0, np.random.normal(mu, sigma, size))
-def demand_uniform(mu, sigma, size): return np.random.uniform(mu - sigma, mu + sigma, size)
-def demand_exponential(mu, sigma, size): return np.random.exponential(mu, size)
-def demand_logistic(mu, sigma, size): return np.random.logistic(mu, sigma, size)
-
-DEMAND_FUNCTIONS = {
-    "linear": demand_linear,
-    "uniform": demand_uniform,
-    "exponential": demand_exponential,
-    "logistic": demand_logistic,
-}
-
-FACTORS = [
-    Factor("demand_fn", list(DEMAND_FUNCTIONS.keys()), primary=True),
-    Factor("alpha", [0.1, 0.3, 0.5, 0.7], primary=True),
-    Factor("n_products", [5, 15, 30, 50], primary=True),
-    Factor("demand_mu", [30.0, 50.0, 70.0], primary=False),
-    Factor("demand_sigma", [5.0, 10.0, 20.0], primary=False),
-    Factor("N", [100, 500, 1000], primary=False),
-]
-
-SEEDS_PER_CONFIG = 5

engine/studies/full_factorial.py

@@ -1,89 +0,0 @@
-"""full factorial design - all factor combinations"""
-import sys
-sys.path.insert(0, "..")
-import logging
-from itertools import product
-import json
-import hashlib
-from pathlib import Path
-from concurrent.futures import ProcessPoolExecutor
-from .factors import FACTORS, DEMAND_FUNCTIONS, SEEDS_PER_CONFIG
-
-logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
-log = logging.getLogger(__name__)
-
-def generate_configs():
-    """generate all factor combinations with seeds"""
-    all_levels = [f.levels for f in FACTORS]
-    names = [f.name for f in FACTORS]
-
-    configs = []
-    for combo in product(*all_levels):
-        base = {names[i]: combo[i] for i in range(len(names))}
-        for seed in range(SEEDS_PER_CONFIG):
-            cfg = {**base, "seed": seed}
-            cfg["id"] = hashlib.md5(json.dumps(cfg, sort_keys=True).encode()).hexdigest()[:8]
-            configs.append(cfg)
-    return configs
-
-def run_single(cfg: dict) -> dict:
-    """execute one experiment config, return metrics"""
-    from engine.wrapper import PHANTOM
-    import numpy as np
-
-    np.random.seed(cfg["seed"])
-    demand_fn = DEMAND_FUNCTIONS[cfg["demand_fn"]]
-
-    env = PHANTOM(
-        n_products=cfg["n_products"],
-        alpha=cfg["alpha"],
-        N=cfg["N"],
-    )
-    env.market.demand = (demand_fn, (cfg["demand_mu"], cfg["demand_sigma"]))
-
-    obs, _ = env.reset()
-    total_reward, steps = 0.0, 0
-
-    for _ in range(100):
-        action = env.action_space.sample()
-        obs, reward, term, trunc, _ = env.step(action)
-        total_reward += reward
-        steps += 1
-        if term: break
-
-    env.close()
-    return {
-        "id": cfg["id"],
-        "config": cfg,
-        "total_reward": total_reward,
-        "avg_reward": total_reward / steps,
-        "steps": steps,
-    }
-
-def run_study(max_workers: int = None, output: str = "results_full.jsonl"):
-    configs = generate_configs()
-    log.info(f"full factorial: {len(configs)} configs ({len(configs)//SEEDS_PER_CONFIG} unique × {SEEDS_PER_CONFIG} seeds)")
-
-    results = []
-    with ProcessPoolExecutor(max_workers=max_workers) as ex:
-        for i, result in enumerate(ex.map(run_single, configs)):
-            results.append(result)
-            if (i+1) % 100 == 0: log.info(f"progress: {i+1}/{len(configs)}")
-
-    Path(output).write_text("\n".join(json.dumps(r) for r in results))
-    log.info(f"wrote {len(results)} results to {output}")
-    return results
-
-if __name__ == "__main__":
-    import argparse
-    p = argparse.ArgumentParser()
-    p.add_argument("--workers", type=int, default=None)
-    p.add_argument("--output", default="results_full.jsonl")
-    p.add_argument("--dry-run", action="store_true", help="only show design size")
-    args = p.parse_args()
-
-    configs = generate_configs()
-    log.info(f"design: {len(configs)} runs | factors: {[f.name for f in FACTORS]} | levels: {[len(f.levels) for f in FACTORS]}")
-
-    if not args.dry_run:
-        run_study(args.workers, args.output)

engine/studies/mixed_lh.py

@@ -1,106 +0,0 @@
-"""mixed design: full factorial on primary factors, latin hypercube on secondary"""
-import sys
-sys.path.insert(0, "..")
-import logging
-from itertools import product
-import json
-import hashlib
-from pathlib import Path
-from concurrent.futures import ProcessPoolExecutor
-import numpy as np
-from scipy.stats.qmc import LatinHypercube
-from factors import FACTORS, DEMAND_FUNCTIONS, SEEDS_PER_CONFIG
-
-logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
-log = logging.getLogger(__name__)
-
-LH_SAMPLES = 10
-
-def generate_configs(lh_samples: int = LH_SAMPLES):
-    primary = [f for f in FACTORS if f.primary]
-    secondary = [f for f in FACTORS if not f.primary]
-
-    primary_grid = list(product(*[f.levels for f in primary]))
-    lhs = LatinHypercube(d=len(secondary), seed=42)
-
-    configs = []
-    for p_combo in primary_grid:
-        samples = lhs.random(n=lh_samples)
-        for s in samples:
-            sec_vals = {
-                secondary[i].name: secondary[i].levels[int(s[i] * len(secondary[i].levels))]
-                for i in range(len(secondary))
-            }
-            base = {primary[i].name: p_combo[i] for i in range(len(primary))}
-            base.update(sec_vals)
-
-            for seed in range(SEEDS_PER_CONFIG):
-                cfg = {**base, "seed": seed}
-                cfg["id"] = hashlib.md5(json.dumps(cfg, sort_keys=True).encode()).hexdigest()[:8]
-                configs.append(cfg)
-    return configs
-
-def run_single(cfg: dict) -> dict:
-    from engine.wrapper import PHANTOM
-    import numpy as np
-
-    np.random.seed(cfg["seed"])
-    demand_fn = DEMAND_FUNCTIONS[cfg["demand_fn"]]
-
-    env = PHANTOM(
-        n_products=cfg["n_products"],
-        alpha=cfg["alpha"],
-        N=cfg["N"],
-    )
-    env.market.demand = (demand_fn, (cfg["demand_mu"], cfg["demand_sigma"]))
-
-    obs, _ = env.reset()
-    total_reward, steps = 0.0, 0
-
-    for _ in range(100):
-        action = env.action_space.sample()
-        obs, reward, term, trunc, _ = env.step(action)
-        total_reward += reward
-        steps += 1
-        if term: break
-
-    env.close()
-    return {
-        "id": cfg["id"],
-        "config": cfg,
-        "total_reward": total_reward,
-        "avg_reward": total_reward / steps,
-        "steps": steps,
-    }
-
-def run_study(max_workers: int = None, output: str = "results_mixed.jsonl", lh_samples: int = LH_SAMPLES):
-    configs = generate_configs(lh_samples)
-    n_primary_cells = int(np.prod([len(f.levels) for f in FACTORS if f.primary]))
-    log.info(f"mixed LH: {len(configs)} configs ({n_primary_cells} primary × {lh_samples} LH × {SEEDS_PER_CONFIG} seeds)")
-
-    results = []
-    with ProcessPoolExecutor(max_workers=max_workers) as ex:
-        for i, result in enumerate(ex.map(run_single, configs)):
-            results.append(result)
-            if (i+1) % 100 == 0: log.info(f"progress: {i+1}/{len(configs)}")
-
-    Path(output).write_text("\n".join(json.dumps(r) for r in results))
-    log.info(f"wrote {len(results)} results to {output}")
-    return results
-
-if __name__ == "__main__":
-    import argparse
-    p = argparse.ArgumentParser()
-    p.add_argument("--workers", type=int, default=None)
-    p.add_argument("--output", default="results_mixed.jsonl")
-    p.add_argument("--lh-samples", type=int, default=10)
-    p.add_argument("--dry-run", action="store_true", help="only show design size")
-    args = p.parse_args()
-
-    primary = [f for f in FACTORS if f.primary]
-    secondary = [f for f in FACTORS if not f.primary]
-    configs = generate_configs(args.lh_samples)
-    log.info(f"design: {len(configs)} runs | primary: {[f.name for f in primary]} | secondary (LH): {[f.name for f in secondary]}")
-
-    if not args.dry_run:
-        run_study(args.workers, args.output, args.lh_samples)

engine/train.py

@@ -1,45 +0,0 @@
-from stable_baselines3 import SAC
-from stable_baselines3.common.callbacks import EvalCallback, BaseCallback
-from .wrapper import PHANTOM
-
-
-class RenderCallback(BaseCallback):
-    """Renders environment on every step for live visualization."""
-    def __init__(self, env: PHANTOM):
-        super().__init__()
-        self.env = env
-
-    def _on_step(self) -> bool:
-        self.env.render()
-        return True
-
-
-env = PHANTOM(n_products=10, alpha=0.3, render_mode="human")
-eval_env = PHANTOM(n_products=10, alpha=0.3, render_mode=None)
-
-model = SAC(
-    "MultiInputPolicy",
-    env,
-    verbose=1,
-    learning_rate=3e-4,
-    buffer_size=50000,
-    batch_size=256,
-    tau=0.005,
-    gamma=0.99,
-)
-
-render_cb = RenderCallback(env)
-eval_cb = EvalCallback(eval_env, eval_freq=1000, n_eval_episodes=5, verbose=1)
-
-model.learn(total_timesteps=50000, callback=[render_cb, eval_cb])
-model.save("phantom_sac")
-
-# test trained policy
-env = PHANTOM(n_products=10, alpha=0.3, render_mode="human")
-obs, _ = env.reset()
-for _ in range(100):
-    action, _ = model.predict(obs, deterministic=True)
-    obs, reward, term, trunc, _ = env.step(action)
-    env.render()
-    if term or trunc: break
-env.close()

engine/wrapper.py

@@ -1,118 +0,0 @@
-import gymnasium as gym
-from gymnasium import spaces
-import numpy as np
-from .engine import Limbo, MarketEngine, PricingEngine
-from .lib.render import DashboardRenderer
-
-
-class PHANTOM(gym.Env):
-    """Gymnasium wrapper for the Limbo pricing-market simulation. Platform sets prices, market responds with demand."""
-    metadata = {"render_modes": ["human", "ansi"]}
-
-    def __init__(self,
-                 n_products: int = 10,
-                 alpha: float = 0.3,
-                 N: int = 100,
-                 price_bounds: tuple = (10.0, 150.0),
-                 lambda_coi: float = 0.1,
-                 render_mode: str = None):
-        super().__init__()
-        self.n_products = n_products
-        self.price_bounds = price_bounds
-        self.lambda_coi = lambda_coi
-        self.render_mode = render_mode
-        self.alpha = alpha
-        self.N = N
-
-        self.market = MarketEngine(alpha=alpha, N=N)
-        self._platform_stub = PricingEngine()
-        self._limbo = Limbo(self._platform_stub, self.market)
-
-        self.action_space = spaces.Box(
-            low=price_bounds[0], high=price_bounds[1],
-            shape=(n_products,), dtype=np.float32
-        )
-        self.observation_space = spaces.Dict({
-            "demand": spaces.Box(low=0.0, high=100.0, shape=(n_products,), dtype=np.float32),
-            "prices": spaces.Box(low=price_bounds[0], high=price_bounds[1], shape=(n_products,), dtype=np.float32),
-        })
-
-        self._prices = None
-        self._demand = None
-        self._step_count = 0
-        self._demand_history = []
-        self._price_history = []
-        self._revenue_history = []
-        self._renderer = None
-
-    def _get_obs(self) -> dict:
-        demand_arr = np.array([self._demand.get(i, 0.0) for i in range(self.n_products)], dtype=np.float32)
-        return {"demand": demand_arr, "prices": self._prices.astype(np.float32)}
-
-    def _compute_reward(self, prices: np.ndarray, demand: dict) -> float:
-        revenue = np.sum(prices * np.array([demand.get(i, 0.0) for i in range(self.n_products)]))
-        # TODO: implement supra-competitive price punishment
-        return float(revenue)
-
-    def _record_history(self):
-        demand_arr = np.array([self._demand.get(i, 0.0) for i in range(self.n_products)])
-        self._demand_history.append(demand_arr)
-        self._price_history.append(self._prices.copy())
-        self._revenue_history.append(np.sum(self._prices * demand_arr))
-
-    def reset(self, seed=None, options=None):
-        super().reset(seed=seed)
-        self._prices = np.random.uniform(*self.price_bounds, size=self.n_products)
-        self._demand = self.market.act(self._prices)
-        self._step_count = 0
-        self._demand_history, self._price_history, self._revenue_history = [], [], []
-        self._record_history()
-        return self._get_obs(), {}
-
-    def step(self, action: np.ndarray):
-        self._prices = np.clip(action, *self.price_bounds)
-        self._demand = self.market.act(self._prices)
-        self._step_count += 1
-        self._record_history()
-
-        reward = self._compute_reward(self._prices, self._demand)
-        terminated = self._step_count >= 100
-
-        return self._get_obs(), reward, terminated, False, {"step": self._step_count}
-
-    def _compute_elasticity(self) -> np.ndarray:
-        """point elasticity: e = (dQ/dP) * (P/Q) via finite differences, clipped to [-5, 5]"""
-        if len(self._price_history) < 2:
-            return np.zeros(self.n_products)
-        p, q = np.array(self._price_history), np.array(self._demand_history)
-        dp, dq = np.diff(p, axis=0), np.diff(q, axis=0)
-        valid = np.abs(dp) > 0.5
-        with np.errstate(divide='ignore', invalid='ignore'):
-            elasticity = np.where(valid, (dq / dp) * (p[:-1] / np.maximum(q[:-1], 1.0)), 0.0)
-            elasticity = np.nan_to_num(np.clip(elasticity, -5.0, 5.0), nan=0.0)
-        return np.mean(elasticity, axis=0) if len(elasticity) > 0 else np.zeros(self.n_products)
-
-    def render(self):
-        if self.render_mode == "human":
-            if self._renderer is None:
-                self._renderer = DashboardRenderer()
-            self._renderer.render(self)
-        elif self.render_mode == "ansi":
-            return f"step={self._step_count}, prices={self._prices}, demand={self._demand}"
-        return None
-
-    def close(self):
-        if self._renderer:
-            self._renderer.close()
-            self._renderer = None
-
-
-if __name__ == "__main__":
-    env = PHANTOM(n_products=15, alpha=0.3, N=100, render_mode="human")
-    obs, _ = env.reset()
-    for step in range(100):
-        action = env.action_space.sample()
-        obs, reward, term, trunc, info = env.step(action)
-        env.render()
-        if term: break
-    env.close()

lab/case/thesis/coi.py

@@ -0,0 +1,378 @@
+"""Cost of Information (COI) computation for thesis pricing simulation.
+
+Implements the corrected COI formulation:
+
+    COI = E[p] - p
+
+where:
+- E[p] = expected price BEFORE information revelation (window start price)
+- p = actual transaction price (price at which sales occur)
+
+The fundamental insight is that COI should measure PRICE EROSION over time,
+not instantaneous margin leakage. When agents explore across sessions:
+1. They reveal demand signals that drive platform price adjustments
+2. Coordinated agents can find the minimum price across their session pool
+3. The price path from window start to transaction captures information leakage
+
+Key components:
+- COIWindow: Windowed price erosion measurement over K steps
+- compute_coi_window: Per-episode COI from session-level transactions
+- coi_erosion: Order statistic erosion (Theorem 1: N agents -> min price)
+
+This fixes the fundamental error of treating COI as instantaneous margin × alpha.
+"""
+from __future__ import annotations
+from dataclasses import dataclass, field
+from typing import Dict, List, TYPE_CHECKING
+import numpy as np
+
+if TYPE_CHECKING:
+    from .simplified import Session
+
+EPS = 1e-10
+
+
+@dataclass
+class COIWindow:
+    """Windowed COI measurement capturing price erosion over time.
+
+    Attributes:
+        policy: Platform's intended COI (prices at window start - cost)
+        agent: Realized COI for agents (prices at transaction - cost)
+        leak: COI leakage = policy - agent (price erosion due to exploration)
+        survival_ratio: Fraction of intended COI that survives (agent/policy)
+        policy_by_product: Per-product policy COI
+        agent_by_product: Per-product agent COI
+        demand_weights: Demand weights used for aggregation
+    """
+    policy: float = 0.0          # E[p] - c at window start
+    agent: float = 0.0           # p_transaction - c
+    leak: float = 0.0            # policy - agent = price erosion
+    survival_ratio: float = 1.0  # agent / policy
+    policy_by_product: np.ndarray = field(default_factory=lambda: np.zeros(1))
+    agent_by_product: np.ndarray = field(default_factory=lambda: np.zeros(1))
+    demand_weights: np.ndarray = field(default_factory=lambda: np.zeros(1))
+
+    def to_dict(self) -> Dict[str, float]:
+        return {
+            'coi_policy': self.policy,
+            'coi_agent': self.agent,
+            'coi_leak': self.leak,
+            'coi_survival': self.survival_ratio,
+        }
+
+
+def compute_coi_window(
+    sessions: List["Session"],
+    costs: np.ndarray,
+    demand_mapping: Dict[str, float] = None,
+    window_prices: np.ndarray = None,
+) -> COIWindow:
+    """Compute COI from session data using the corrected formulation.
+
+    COI = E[p_start] - p_transaction
+
+    This measures how much the platform's pricing power eroded during the window.
+    Price at window start represents E[p] (what we expected to charge).
+    Transaction prices represent p (what we actually charged).
+
+    Args:
+        sessions: List of sessions with events containing price_seen and purchases
+        costs: Product costs array
+        demand_mapping: Optional session_id -> demand proxy mapping
+        window_prices: Optional explicit window start prices (otherwise use first seen)
+
+    Returns:
+        COIWindow with erosion metrics
+    """
+    if not sessions:
+        n = len(costs)
+        zeros = np.zeros(n)
+        return COIWindow(policy=0.0, agent=0.0, leak=0.0, survival_ratio=1.0,
+                        policy_by_product=zeros, agent_by_product=zeros, demand_weights=zeros)
+
+    n = len(costs)
+    demand_mapping = demand_mapping or {}
+
+    # Track prices seen at start (E[p]) and transaction prices (p)
+    first_prices = np.zeros(n)  # first price seen per product (window start proxy)
+    transaction_prices = np.zeros(n)  # prices at which purchases occurred
+    transaction_counts = np.zeros(n)
+    view_counts = np.zeros(n)
+    demand_weights = np.zeros(n)
+
+    for sess in sessions:
+        sid = sess.sid
+        sess_demand = demand_mapping.get(sid, 1.0)
+
+        for e in sess.events:
+            pidx = e.product_idx
+            if pidx < 0 or pidx >= n:
+                continue
+
+            price_seen = float(e.price_seen)
+
+            # Track first price seen (proxy for E[p] at window start)
+            if view_counts[pidx] == 0:
+                first_prices[pidx] = price_seen
+            view_counts[pidx] += 1
+
+            # Track transaction prices
+            if e.action == "purchase":
+                transaction_prices[pidx] += price_seen
+                transaction_counts[pidx] += 1
+                demand_weights[pidx] += sess_demand
+
+    # Compute per-product COI
+    # Policy COI: what we intended to charge (first seen price - cost)
+    policy_by_product = np.zeros(n)
+    agent_by_product = np.zeros(n)
+
+    for i in range(n):
+        if view_counts[i] > 0:
+            # Use explicit window prices if provided, else first seen
+            start_price = window_prices[i] if window_prices is not None else first_prices[i]
+            policy_by_product[i] = max(0, start_price - costs[i])
+
+        if transaction_counts[i] > 0:
+            avg_transaction = transaction_prices[i] / transaction_counts[i]
+            agent_by_product[i] = max(0, avg_transaction - costs[i])
+
+    # Aggregate with demand weighting
+    total_demand = np.sum(demand_weights) + EPS
+    weights = demand_weights / total_demand
+
+    # Only count products with transactions for fair comparison
+    active_mask = transaction_counts > 0
+    if np.any(active_mask):
+        policy = float(np.sum(policy_by_product[active_mask] * weights[active_mask]) /
+                      (np.sum(weights[active_mask]) + EPS))
+        agent = float(np.sum(agent_by_product[active_mask] * weights[active_mask]) /
+                     (np.sum(weights[active_mask]) + EPS))
+    else:
+        # No transactions - use view-weighted policy COI
+        view_weights = view_counts / (np.sum(view_counts) + EPS)
+        policy = float(np.sum(policy_by_product * view_weights))
+        agent = policy  # No erosion without transactions
+
+    # Leak = price erosion due to information revelation
+    leak = max(0, policy - agent)
+    survival = agent / (policy + EPS) if policy > EPS else 1.0
+
+    return COIWindow(
+        policy=policy,
+        agent=agent,
+        leak=leak,
+        survival_ratio=float(np.clip(survival, 0, 1)),
+        policy_by_product=policy_by_product,
+        agent_by_product=agent_by_product,
+        demand_weights=demand_weights,
+    )
+
+
+def coi_erosion(policy_coi: float, agent_coi: float) -> float:
+    """Compute COI erosion rate: (policy - agent) / policy.
+
+    Returns the fraction of intended COI that was lost to information leakage.
+    0 = no erosion, 1 = complete erosion.
+    """
+    if policy_coi < EPS:
+        return 0.0
+    return float(np.clip((policy_coi - agent_coi) / policy_coi, 0, 1))
+
+
+def order_statistic_erosion(n_agents: int, price_std: float, base_margin: float = 1.0) -> float:
+    """Compute COI erosion from order statistic effect (Theorem 1).
+
+    When N agents independently query prices:
+    - Each sees a price p_i ~ N(μ, σ²)
+    - They coordinate to buy at min(p_1, ..., p_N)
+    - Expected minimum: μ - σ * E[order_stat]
+
+    As N -> ∞, E[min] -> p_min, so COI -> 0.
+
+    This quantifies the price discovery benefit of multiple sessions.
+
+    Args:
+        n_agents: Number of independent agent sessions
+        price_std: Standard deviation of price distribution
+        base_margin: Expected margin (μ - cost)
+
+    Returns:
+        Erosion rate in [0, 1]
+    """
+    if n_agents <= 1 or price_std < EPS:
+        return 0.0
+
+    # For standard normal order statistics, E[min of N] ≈ -Φ^{-1}(1/(N+1))
+    # For large N, this grows like sqrt(2 * log(N))
+    log_n = np.log(n_agents)
+    if log_n < 0.1:
+        return 0.0
+
+    # Extreme value theory: expected min shift
+    shift = price_std * (np.sqrt(2 * log_n) -
+                        (np.log(log_n) + np.log(4 * np.pi)) / (2 * np.sqrt(2 * log_n) + EPS))
+
+    # Erosion = shift / base_margin, capped at 1
+    return float(np.clip(shift / (base_margin + EPS), 0, 1))
+
+
+@dataclass
+class COITracker:
+    """Track COI over multiple windows for temporal analysis.
+
+    This addresses the user's insight: compute COI over K episodes to see
+    how prices change from window start to end.
+
+    If at start of window price is A and by end it's B, the difference
+    A - B represents COI leakage from exploratory sessions.
+    """
+    window_size: int = 10  # K episodes per window
+    _price_history: List[np.ndarray] = field(default_factory=list)
+    _transaction_history: List[np.ndarray] = field(default_factory=list)
+    _coi_history: List[float] = field(default_factory=list)
+
+    def add_step(self, prices: np.ndarray, transactions: np.ndarray = None):
+        """Record price observation for current step."""
+        self._price_history.append(prices.copy())
+        if transactions is not None:
+            self._transaction_history.append(transactions.copy())
+
+    def compute_window_coi(self, costs: np.ndarray) -> float:
+        """Compute COI over the current window.
+
+        COI = E[p_start] - E[p_end] for the window.
+        This captures price erosion due to information revelation.
+        """
+        if len(self._price_history) < 2:
+            return 0.0
+
+        # Get prices at window boundaries
+        window_start = max(0, len(self._price_history) - self.window_size)
+        start_prices = self._price_history[window_start]
+        end_prices = self._price_history[-1]
+
+        # COI = (start_price - cost) - (end_price - cost) = start_price - end_price
+        start_margin = np.mean(start_prices - costs)
+        end_margin = np.mean(end_prices - costs)
+
+        coi = max(0, start_margin - end_margin)
+        self._coi_history.append(coi)
+        return coi
+
+    def get_cumulative_erosion(self, costs: np.ndarray) -> float:
+        """Compute total COI erosion from first observation to now."""
+        if len(self._price_history) < 2:
+            return 0.0
+
+        initial = np.mean(self._price_history[0] - costs)
+        current = np.mean(self._price_history[-1] - costs)
+        return max(0, initial - current)
+
+    def get_erosion_trend(self) -> float:
+        """Get average COI per window (erosion rate)."""
+        if not self._coi_history:
+            return 0.0
+        return float(np.mean(self._coi_history))
+
+    def reset(self):
+        """Reset tracker for new episode."""
+        self._price_history.clear()
+        self._transaction_history.clear()
+        self._coi_history.clear()
+
+
+def compute_multi_session_coi(
+    sessions: List["Session"],
+    costs: np.ndarray,
+    alpha: float,
+    initial_prices: np.ndarray,
+) -> Dict[str, float]:
+    """Compute COI accounting for multi-session agent behavior.
+
+    This is the key fix for the fundamental error:
+    - Agents use different sessions to gather information
+    - Each session reveals price information
+    - Coordinated agents find the minimum across their session pool
+
+    The COI is computed as:
+    1. What platform intended to charge: initial_prices - costs
+    2. What agents actually paid: min(prices seen across sessions) - costs
+    3. Leak = (1) - (2)
+
+    Args:
+        sessions: All sessions in the episode
+        costs: Product costs
+        alpha: Contamination level (fraction of agent sessions)
+        initial_prices: Prices at episode start (E[p])
+
+    Returns:
+        Dictionary with COI metrics
+    """
+    n = len(costs)
+
+    # Separate agent and human sessions by ground truth label
+    agent_sessions = [s for s in sessions if s.actor == "A"]
+    human_sessions = [s for s in sessions if s.actor == "H"]
+
+    # Track prices seen by agents per product (for min finding)
+    agent_prices_seen: Dict[int, List[float]] = {i: [] for i in range(n)}
+    human_prices_paid: Dict[int, List[float]] = {i: [] for i in range(n)}
+
+    for sess in agent_sessions:
+        for e in sess.events:
+            if 0 <= e.product_idx < n:
+                agent_prices_seen[e.product_idx].append(e.price_seen)
+
+    for sess in human_sessions:
+        for e in sess.events:
+            if 0 <= e.product_idx < n and e.action == "purchase":
+                human_prices_paid[e.product_idx].append(e.price_seen)
+
+    # Compute COI components
+    policy_coi = float(np.mean(initial_prices - costs))  # E[p] - c
+
+    # Agent COI: they find the minimum price via exploration
+    agent_coi_by_product = np.zeros(n)
+    for i in range(n):
+        if agent_prices_seen[i]:
+            min_price = min(agent_prices_seen[i])
+            agent_coi_by_product[i] = max(0, min_price - costs[i])
+        else:
+            agent_coi_by_product[i] = initial_prices[i] - costs[i]
+
+    agent_coi = float(np.mean(agent_coi_by_product))
+
+    # Human COI: they pay whatever price is offered
+    human_coi_by_product = np.zeros(n)
+    for i in range(n):
+        if human_prices_paid[i]:
+            avg_price = np.mean(human_prices_paid[i])
+            human_coi_by_product[i] = max(0, avg_price - costs[i])
+        else:
+            human_coi_by_product[i] = initial_prices[i] - costs[i]
+
+    human_coi = float(np.mean(human_coi_by_product))
+
+    # Total leak: weighted by contamination
+    # Agents erode COI, humans pay full price
+    realized_coi = (1 - alpha) * human_coi + alpha * agent_coi
+    leak = policy_coi - realized_coi
+
+    # Order statistic effect: more agents = more erosion
+    n_agents = len(agent_sessions)
+    price_std = float(np.std(initial_prices))
+    order_erosion = order_statistic_erosion(n_agents, price_std, policy_coi)
+
+    return {
+        'policy_coi': policy_coi,
+        'agent_coi': agent_coi,
+        'human_coi': human_coi,
+        'realized_coi': realized_coi,
+        'leak': leak,
+        'order_stat_erosion': order_erosion,
+        'n_agent_sessions': n_agents,
+        'n_human_sessions': len(human_sessions),
+        'survival_ratio': realized_coi / (policy_coi + EPS) if policy_coi > EPS else 1.0,
+    }

lab/case/thesis/separability.py

@@ -0,0 +1,104 @@
+"""Behavioral separability for thesis human/agent classification.
+
+Implements KL-divergence based separability scoring (Eq 20-21):
+- Δ_H = D_KL(T̂' || T̄_H): divergence from human reference kernel
+- Δ_A = D_KL(T̂' || T̄_A): divergence from agent reference kernel
+- α̂(τ') = σ(β(Δ_H - Δ_A)): per-session contamination estimate
+"""
+from __future__ import annotations
+from typing import Dict, List, TYPE_CHECKING
+import numpy as np
+
+if TYPE_CHECKING:
+    from .simplified import Session
+
+
+# Reference transition kernels T̄_H, T̄_A estimated from real data (Eq 19)
+TRANS_H = {
+    "start": {"view": 0.85, "end": 0.15},
+    "view": {"detail": 0.4, "add_to_cart": 0.3, "view": 0.2, "end": 0.1},
+    "detail": {"add_to_cart": 0.5, "view": 0.3, "end": 0.2},
+    "add_to_cart": {"purchase": 0.6, "view": 0.25, "end": 0.15},
+    "purchase": {"end": 1.0},
+    "checkout": {"purchase": 0.8, "end": 0.2},
+    "hover": {"view": 0.5, "detail": 0.3, "end": 0.2},
+}
+
+TRANS_A = {
+    "start": {"view": 0.95, "end": 0.05},
+    "view": {"detail": 0.6, "view": 0.25, "add_to_cart": 0.1, "end": 0.05},
+    "detail": {"view": 0.5, "add_to_cart": 0.15, "detail": 0.3, "end": 0.05},
+    "add_to_cart": {"view": 0.4, "purchase": 0.2, "end": 0.4},
+    "purchase": {"end": 1.0},
+    "checkout": {"purchase": 0.3, "end": 0.7},
+    "hover": {"view": 0.6, "detail": 0.35, "end": 0.05},
+}
+
+
+def kl_div(p: Dict[str, float], q: Dict[str, float], eps: float = 1e-10) -> float:
+    """Compute KL(p || q) with smoothing."""
+    if not p or not q:
+        return 0.0
+    all_keys = set(p.keys()) | set(q.keys())
+    total = 0.0
+    for k in all_keys:
+        pk = p.get(k, eps)
+        qk = q.get(k, eps)
+        if pk > eps:
+            total += pk * np.log(pk / max(qk, eps))
+    return max(0.0, total)
+
+
+def build_kernel(events: List) -> Dict[str, Dict[str, float]]:
+    """Build empirical transition kernel from event sequence."""
+    trans: Dict[str, Dict[str, int]] = {}
+    prev = "start"
+    for e in events:
+        curr = getattr(e, 'action', None) or e.get('action', 'end') if isinstance(e, dict) else 'end'
+        trans.setdefault(prev, {})
+        trans[prev][curr] = trans[prev].get(curr, 0) + 1
+        prev = curr
+    # add terminal transition
+    trans.setdefault(prev, {})
+    trans[prev]["end"] = trans[prev].get("end", 0) + 1
+
+    # normalize to probabilities
+    kernel = {}
+    for s, dests in trans.items():
+        total = sum(dests.values())
+        kernel[s] = {d: c / total for d, c in dests.items()} if total > 0 else {"end": 1.0}
+    return kernel
+
+
+def compute_divergence(kernel: Dict[str, Dict[str, float]], ref_h: Dict = None, ref_a: Dict = None) -> tuple[float, float]:
+    """Compute Δ_H, Δ_A divergence from reference kernels (Eq 20-21)."""
+    ref_h = ref_h or TRANS_H
+    ref_a = ref_a or TRANS_A
+    delta_h = sum(kl_div(kernel.get(s, {}), ref_h.get(s, {})) for s in kernel) / max(len(kernel), 1)
+    delta_a = sum(kl_div(kernel.get(s, {}), ref_a.get(s, {})) for s in kernel) / max(len(kernel), 1)
+    return delta_h, delta_a
+
+
+def estimate_alpha(session: "Session", beta: float = 2.0) -> float:
+    """Estimate per-session contamination α̂(τ') = σ(β(Δ_H - Δ_A)).
+
+    High Δ_H (far from human) and low Δ_A (close to agent) -> high α̂ (likely agent).
+    """
+    if not session.events:
+        return 0.5
+    kernel = build_kernel(session.events)
+    delta_h, delta_a = compute_divergence(kernel)
+
+    if delta_h + delta_a < 1e-6:
+        return 0.5
+
+    # sigmoid: high when trajectory is more divergent from human than agent
+    return 1.0 / (1.0 + np.exp(-beta * (delta_h - delta_a)))
+
+
+def batch_estimate_alpha(sessions: List["Session"]) -> tuple[float, List[float]]:
+    """Estimate aggregate and per-session contamination."""
+    if not sessions:
+        return 0.0, []
+    alphas = [estimate_alpha(s) for s in sessions]
+    return float(np.mean(alphas)), alphas

lab/case/thesis/simplified.py

@@ -8,6 +8,14 @@ Objects:
 - Demand proxy q_hat via weighted action aggregation
 - COI leakage penalty for agent reconnaissance
 - Limbo: alternating price/demand history for trajectory analysis
+
+COI Correction (Jan 2026):
+The fundamental COI formulation is:
+    COI = E[p_start] - p_transaction
+
+This measures price erosion over time, not instantaneous margin × alpha.
+Agents use multiple sessions to gather information and find minimum prices.
+The price path from episode start to transaction captures information leakage.
 """
 from __future__ import annotations
 from dataclasses import dataclass, field

lab/case/thesis/simplified_env.py

@@ -6,6 +6,14 @@ Supports multiple reward modes and contamination scenarios.
 Action: price multipliers [0.5, 1.5] applied to reference prices
 Observation: [prices, demand_agg, alpha_est, margins, position_proxy]
 Reward: configurable objective (revenue, profit, robust, coi-aware)
+
+COI Correction (Jan 2026):
+The fundamental COI formulation is now:
+    COI = E[p_start] - p_transaction
+
+This measures price erosion over time, not instantaneous margin × alpha.
+Agents using different sessions gather information and drive prices down.
+The COITracker now tracks prices over windows to capture this effect.
 """
 from __future__ import annotations
 from dataclasses import dataclass
@@ -20,7 +28,7 @@ except ImportError:
     HAS_GYM = False
 
 from .simplified import System, Session, Event, Limbo, put_prices_to_market, compute_demand, estimate_alpha
-from .coi import COIWindow, compute_coi_window, coi_erosion
+from .coi import COIWindow, compute_coi_window, coi_erosion, COITracker, compute_multi_session_coi
 
 
 @dataclass
@@ -73,6 +81,12 @@ class PricingEnv(gym.Env if HAS_GYM else object):
         self._episode_rewards: list[float] = []
         self._demand_agg = np.zeros(self.n)
 
+        # COI tracking: store initial prices for E[p] calculation
+        self._initial_prices: np.ndarray | None = None
+        self._coi_tracker = COITracker(window_size=10)
+        self._last_coi_metrics: Dict[str, float] = {}
+        self._last_window_coi: float = 0.0
+
         self.action_space = spaces.Box(low=0.5, high=1.5, shape=(self.n,), dtype=np.float32)
         obs_dim = self.n + self.n + 1 + 1 + self.n + 1  # prices + demand + alpha_hat + alpha + margins + t
         self.observation_space = spaces.Box(low=-np.inf, high=np.inf, shape=(obs_dim,), dtype=np.float32)
@@ -109,8 +123,29 @@ class PricingEnv(gym.Env if HAS_GYM else object):
         if self._last_prices is not None:
             vol_penalty = cfg.lambda_vol * float(np.mean(np.abs(prices - self._last_prices) / (sys.refs + 1e-6)))
 
+        # Track prices for windowed COI calculation
+        self._coi_tracker.add_step(prices)
+
+        # CORRECTED COI CALCULATION:
+        # COI = E[p_start] - p_transaction (price erosion over time)
+        # Use initial prices as E[p] and compute multi-session COI
+        coi_metrics = compute_multi_session_coi(
+            sessions=sys._last_sessions,
+            costs=sys.costs,
+            alpha=self._alpha,
+            initial_prices=self._initial_prices,
+        )
+        leak = float(coi_metrics['leak'])
+
+        # Also compute window-based COI for trend analysis
+        window_coi = self._coi_tracker.compute_window_coi(sys.costs)
+
+        # Store both for info dict
+        self._last_coi_metrics = coi_metrics
+        self._last_window_coi = window_coi
+
+        # For backward compatibility, also compute the old-style COI
         coi = compute_coi_window(sys._last_sessions, sys.costs, demand_mapping=demand)
-        leak = float(coi.leak)
 
         reward_fns = {
             "revenue": lambda: revenue,
@@ -127,6 +162,11 @@ class PricingEnv(gym.Env if HAS_GYM else object):
         self._t, self._alpha = 0, self.cfg.alpha_true
         self._last_prices, self._last_demand = None, None
         self._episode_rewards, self._demand_agg = [], np.zeros(self.n)
+
+        # COI tracking: store initial prices as E[p] for COI = E[p] - p calculation
+        self._initial_prices = self._sys.refs.copy()
+        self._coi_tracker.reset()
+
         return self._build_obs(), {"alpha_true": self._alpha, "alpha_est": self._sys.alpha,
                                    "costs": self._sys.costs.copy(), "refs": self._sys.refs.copy()}
 
@@ -150,6 +190,9 @@ class PricingEnv(gym.Env if HAS_GYM else object):
         n_agents = int(self._alpha * self.cfg.sessions_per_step)
         coi = compute_coi_window(self._sys._last_sessions, self._sys.costs, demand_mapping=demand)
 
+        # Corrected COI metrics (price erosion over time)
+        coi_m = self._last_coi_metrics
+
         info = {
             "alpha_true": self._alpha, "alpha_est": self._sys.alpha,
             "alpha_error": abs(self._alpha - self._sys.alpha),
@@ -157,9 +200,19 @@ class PricingEnv(gym.Env if HAS_GYM else object):
             "n_purchases": int(np.sum(purchases)),
             "avg_margin": float(np.mean((prices - self._sys.costs) / self._sys.costs)),
             "n_sessions": len(demand), "n_agents": n_agents, "price_std": float(np.std(prices)),
+            # Legacy COI metrics (for backward compatibility)
             "coi_erosion": coi_erosion(coi.policy, coi.agent),
             "coi_policy": float(coi.policy), "coi_agent": float(coi.agent),
             "coi_leakage": float(coi.leak), "coi_survival": float(coi.survival_ratio),
+            # CORRECTED COI metrics: E[p] - p (price erosion)
+            "coi_policy_corrected": float(coi_m.get('policy_coi', 0)),
+            "coi_agent_corrected": float(coi_m.get('agent_coi', 0)),
+            "coi_human_corrected": float(coi_m.get('human_coi', 0)),
+            "coi_realized": float(coi_m.get('realized_coi', 0)),
+            "coi_leak_corrected": float(coi_m.get('leak', 0)),
+            "coi_order_stat_erosion": float(coi_m.get('order_stat_erosion', 0)),
+            "coi_survival_corrected": float(coi_m.get('survival_ratio', 1.0)),
+            "coi_window": float(self._last_window_coi),
             "cumulative_reward": sum(self._episode_rewards), "step": self._t,
         }
         return self._build_obs(), reward, self._t >= self.cfg.max_steps, False, info

lab/case/thesis/train.py

@@ -65,7 +65,7 @@ class ExperimentConfig:
     n_envs: int = 4
     eval_freq: int = 5000
     n_eval_episodes: int = 10
-    log_dir: str = "sim/case/thesis_simplified/runs"
+    log_dir: str = "lab/case/thesis/runs"
     seed: int = 42
     n_products: int = 10
     max_steps: int = 200
@@ -312,7 +312,7 @@ def main():
     parser.add_argument("--n-products", type=int, default=10)
     parser.add_argument("--n-envs", type=int, default=4)
     parser.add_argument("--seed", type=int, default=42)
-    parser.add_argument("--log-dir", default="sim/case/thesis_simplified/runs")
+    parser.add_argument("--log-dir", default="lab/case/thesis/runs")
     parser.add_argument("--sweep", action="store_true", help="run contamination sweep")
     parser.add_argument("--compare", action="store_true", help="compare all baselines")
     parser.add_argument("--workers", type=int, default=None, help="max parallel workers for sweep (None=auto, 1=sequential)")

sim/case/__init__.py

@@ -1,2 +0,0 @@
-"""Case-specific simulations and experiments."""
-

sim/case/thesis_simplified/__init__.py

@@ -1,2 +0,0 @@
-"""Minimal thesis-aligned pricing simulation (self-contained)."""
-

sim/case/thesis_simplified/coi.py

@@ -1,125 +0,0 @@
-"""Cost of Information (COI) computation for thesis pricing system.
-
-Core KPI: COI = E[p_shown] - p_min measures pricing power from information asymmetry.
-Theorem 1 shows COI erodes as agent queries increase: as N->inf, p^(1)->p_min.
-"""
-from __future__ import annotations
-from dataclasses import dataclass
-from typing import Dict, List, TYPE_CHECKING
-import numpy as np
-
-if TYPE_CHECKING:
-    from .simplified import Session
-
-
-@dataclass(frozen=True)
-class COIWindow:
-    """Windowed COI metrics computed from realized price exposures.
-
-    policy: E[p_shown] - cost, the definition-level KPI
-    agent: E[p^(1)] - cost where p^(1) is min price under agent querying
-    leak: max(policy - agent, 0), observable gap from reconnaissance
-    survival_ratio: agent/policy, fraction of pricing power retained
-    """
-    policy: float
-    agent: float
-    leak: float
-    survival_ratio: float
-    policy_by_product: np.ndarray
-    agent_by_product: np.ndarray
-    demand_weights: np.ndarray
-
-
-def aggregate_prices(sessions: List["Session"], mode: str = "all") -> Dict[int, List[float] | float]:
-    """Unified price aggregation across sessions.
-
-    mode: "all" returns all prices per product, "min_per_session" returns min price per session per product,
-          "min_across" returns single min price per product
-    """
-    if mode == "min_across":
-        mins: Dict[int, float] = {}
-        for s in sessions:
-            for e in s.events:
-                pidx, price = int(e.product_idx), float(e.price_seen)
-                mins[pidx] = min(mins.get(pidx, price), price)
-        return mins
-    elif mode == "min_per_session":
-        result: Dict[int, List[float]] = {}
-        for s in sessions:
-            by_p: Dict[int, float] = {}
-            for e in s.events:
-                pidx, price = int(e.product_idx), float(e.price_seen)
-                by_p[pidx] = min(by_p.get(pidx, price), price)
-            for pidx, pmin in by_p.items():
-                result.setdefault(pidx, []).append(pmin)
-        return result
-    else:  # "all"
-        prices: Dict[int, List[float]] = {}
-        for s in sessions:
-            for e in s.events:
-                prices.setdefault(e.product_idx, []).append(float(e.price_seen))
-        return prices
-
-
-def demand_weights_by_product(sessions: List["Session"], demand_mapping: Dict[str, float], n_products: int) -> np.ndarray:
-    """Compute demand-weighted importance per product."""
-    w = np.zeros(n_products, dtype=float)
-    sessions_by_id = {s.sid: s for s in sessions}
-    for sid, q in demand_mapping.items():
-        sess = sessions_by_id.get(sid)
-        if sess and sess.events:
-            w[int(sess.events[0].product_idx)] += float(q)
-    total = float(np.sum(w))
-    return (w / total) if total > 0 else w
-
-
-def compute_coi_window(sessions: List["Session"], costs: np.ndarray, demand_mapping: Dict[str, float] | None = None) -> COIWindow:
-    """Compute COI metrics over session window.
-
-    Aggregates price exposures and computes policy-level vs agent-realized COI.
-    """
-    n = int(len(costs))
-    prices = aggregate_prices(sessions, mode="all")
-    agent_sessions = [s for s in sessions if s.actor == "A"]
-    agent_min = aggregate_prices(agent_sessions, mode="min_across") if agent_sessions else {}
-
-    policy_by = np.zeros(n, dtype=float)
-    agent_by = np.zeros(n, dtype=float)
-    seen = np.array([(i in prices) for i in range(n)], dtype=bool)
-    agent_seen = np.array([(i in agent_min) for i in range(n)], dtype=bool)
-
-    for pidx, ps in prices.items():
-        if 0 <= pidx < n and ps:
-            policy_by[pidx] = float(np.mean(ps) - float(costs[pidx]))
-    for pidx, pmin in agent_min.items():
-        if 0 <= pidx < n:
-            agent_by[pidx] = float(pmin - float(costs[pidx]))
-
-    agent_by[seen & ~agent_seen] = policy_by[seen & ~agent_seen]  # no erosion if no agent exposure
-
-    demand_w = demand_weights_by_product(sessions, demand_mapping, n) if demand_mapping else np.zeros(n, dtype=float)
-    has_weights = float(np.sum(demand_w)) > 0
-
-    if has_weights:
-        policy, agent = float(np.dot(demand_w, policy_by)), float(np.dot(demand_w, agent_by))
-    elif np.any(seen):
-        policy, agent = float(np.mean(policy_by[seen])), float(np.mean(agent_by[seen]))
-    else:
-        policy, agent = 0.0, 0.0
-
-    leak = float(max(policy - agent, 0.0))
-    survival = float(np.clip(agent / policy, 0.0, 1.0)) if policy > 0 else 0.0
-
-    return COIWindow(policy=policy, agent=agent, leak=leak, survival_ratio=survival,
-                     policy_by_product=policy_by, agent_by_product=agent_by, demand_weights=demand_w)
-
-
-def coi_erosion(coi_policy: float, coi_agent: float, eps: float = 1e-9) -> float:
-    """Thesis-consistent COI erosion: fraction of pricing power destroyed by agent queries.
-
-    erosion = 1 - (COI_agent / COI_policy)
-    When agents find low prices, COI_agent -> 0, erosion -> 1.
-    """
-    if coi_policy <= eps:
-        return 0.0
-    return float(np.clip(1.0 - (coi_agent / (coi_policy + eps)), 0.0, 1.0))

sim/case/thesis_simplified/experiments.py

@@ -1,325 +0,0 @@
-"""COI leakage experiments and policy comparisons.
-
-Demonstrates the core thesis contribution: COI erosion under agent contamination
-and recovery via robust pricing policies.
-
-Generates TensorBoard logs for:
-- COI erosion curves across contamination levels
-- Policy comparison (fixed vs adaptive vs RL)
-- Revenue/margin trade-offs
-"""
-from __future__ import annotations
-from dataclasses import dataclass
-from pathlib import Path
-from typing import Dict, List, Tuple
-import json
-import numpy as np
-
-try:
-    from torch.utils.tensorboard import SummaryWriter
-    HAS_TB = True
-except ImportError:
-    HAS_TB = False
-
-from .simplified_env import PricingEnv, EnvConfig, make_env
-from .simplified import System
-
-
-@dataclass
-class ExperimentResult:
-    """Container for experiment metrics."""
-    name: str
-    alpha: float
-    reward_mean: float
-    reward_std: float
-    coi_erosion: float
-    alpha_error: float
-    revenue: float
-    margin: float
-
-    def to_dict(self) -> dict:
-        return {k: getattr(self, k) for k in self.__dataclass_fields__}
-
-
-def theoretical_coi_erosion_curve(alphas: np.ndarray, n_sessions: int = 1000) -> np.ndarray:
-    """Theoretical COI erosion from Theorem 1 using order statistic model.
-
-    For N i.i.d. uniform queries on [p_min, p_max]:
-    E[p^(1)] = p_min + (p_max - p_min)/(N+1), so erosion = 1 - 2/(N+1)
-    """
-    erosions = []
-    for a in alphas:
-        n_agents = max(1, int(a * n_sessions))
-        erosions.append(1.0 - 2.0 / (n_agents + 1))
-    return np.array(erosions)
-
-
-def run_policy_episode(
-    env: PricingEnv,
-    policy_fn,
-    n_episodes: int = 10
-) -> Tuple[List[float], List[float], List[float], List[float]]:
-    """Run policy and collect per-step metrics."""
-    rewards, coi_erosions, alpha_errors, revenues = [], [], [], []
-
-    for _ in range(n_episodes):
-        obs, info = env.reset()
-        done = False
-        while not done:
-            action = policy_fn(obs, env.n)
-            obs, reward, terminated, truncated, info = env.step(action)
-            done = terminated or truncated
-            rewards.append(reward)
-            if 'coi_erosion' in info:
-                coi_erosions.append(info['coi_erosion'])
-            if 'alpha_true' in info and 'alpha_est' in info:
-                alpha_errors.append(abs(info['alpha_true'] - info['alpha_est']))
-            if 'revenue' in info:
-                revenues.append(info['revenue'])
-
-    return rewards, coi_erosions, alpha_errors, revenues
-
-
-class PolicyRegistry:
-    """Registry of baseline policies."""
-
-    @staticmethod
-    def fixed(obs: np.ndarray, n: int, margin: float = 0.15) -> np.ndarray:
-        return np.ones(n, dtype=np.float32) * (1.0 + margin)
-
-    @staticmethod
-    def random(obs: np.ndarray, n: int, rng: np.random.Generator = None) -> np.ndarray:
-        rng = rng or np.random.default_rng()
-        return rng.uniform(0.7, 1.3, n).astype(np.float32)
-
-    @staticmethod
-    def adaptive(obs: np.ndarray, n: int, base_margin: float = 0.15) -> np.ndarray:
-        """Reduce margins when alpha estimate is high."""
-        alpha_est = obs[2 * n] if len(obs) > 2 * n else 0.2
-        margin_scale = 1.0 - 0.4 * alpha_est
-        return np.ones(n, dtype=np.float32) * (1.0 + base_margin * margin_scale)
-
-    @staticmethod
-    def aggressive(obs: np.ndarray, n: int) -> np.ndarray:
-        """High margins, ignores contamination."""
-        return np.ones(n, dtype=np.float32) * 1.4
-
-    @staticmethod
-    def defensive(obs: np.ndarray, n: int) -> np.ndarray:
-        """Low margins, always cautious."""
-        return np.ones(n, dtype=np.float32) * 1.05
-
-    @staticmethod
-    def alpha_proportional(obs: np.ndarray, n: int, max_margin: float = 0.3) -> np.ndarray:
-        """Margin inversely proportional to estimated alpha."""
-        alpha_est = obs[2 * n] if len(obs) > 2 * n else 0.2
-        margin = max_margin * (1.0 - alpha_est)
-        return np.ones(n, dtype=np.float32) * (1.0 + margin)
-
-
-def run_contamination_sweep(
-    alphas: List[float],
-    policies: Dict[str, callable],
-    n_products: int = 10,
-    max_steps: int = 200,
-    n_episodes: int = 10,
-    seed: int = 42,
-    log_dir: str = None
-) -> Dict[str, List[ExperimentResult]]:
-    """Run policies across contamination levels."""
-
-    results = {name: [] for name in policies}
-    writer = SummaryWriter(Path(log_dir) / "sweep") if log_dir and HAS_TB else None
-
-    for alpha in alphas:
-        print(f"  alpha={alpha:.2f}", end=" ")
-        env_cfg = EnvConfig(
-            n_products=n_products, max_steps=max_steps,
-            alpha_true=alpha, reward_mode="robust", seed=seed)
-        env = make_env(env_cfg)
-
-        for name, policy_fn in policies.items():
-            rewards, coi_vals, alpha_errs, revenues = run_policy_episode(env, policy_fn, n_episodes)
-
-            result = ExperimentResult(
-                name=name, alpha=alpha,
-                reward_mean=float(np.mean(rewards)),
-                reward_std=float(np.std(rewards)),
-                coi_erosion=float(np.mean(coi_vals)) if coi_vals else 0.0,
-                alpha_error=float(np.mean(alpha_errs)) if alpha_errs else 0.0,
-                revenue=float(np.mean(revenues)) if revenues else 0.0,
-                margin=float(np.mean([policy_fn(np.zeros(3 * n_products + 3), n_products)]) - 1.0))
-
-            results[name].append(result)
-
-            if writer:
-                step = int(alpha * 100)
-                writer.add_scalar(f'{name}/reward', result.reward_mean, step)
-                writer.add_scalar(f'{name}/coi_erosion', result.coi_erosion, step)
-                writer.add_scalar(f'{name}/alpha_error', result.alpha_error, step)
-                writer.add_scalar(f'{name}/revenue', result.revenue, step)
-
-        print(f"done")
-
-    # add theoretical curve
-    if writer:
-        theo = theoretical_coi_erosion_curve(np.array(alphas))
-        for i, (a, e) in enumerate(zip(alphas, theo)):
-            writer.add_scalar('theoretical/coi_erosion', e, int(a * 100))
-        writer.close()
-
-    return results
-
-
-def run_coi_demonstration(log_dir: str = "sim/case/thesis_simplified/runs", seed: int = 42) -> Dict:
-    """Main COI demonstration experiment."""
-    print("=== COI Leakage Demonstration ===\n")
-
-    Path(log_dir).mkdir(parents=True, exist_ok=True)
-    writer = SummaryWriter(Path(log_dir) / "coi_demo") if HAS_TB else None
-
-    # theoretical erosion curve
-    print("1. Theoretical COI erosion (Theorem 1)")
-    alphas = np.linspace(0.0, 0.6, 13)
-    theo_erosion = theoretical_coi_erosion_curve(alphas, n_sessions=1000)
-
-    for a, e in zip(alphas, theo_erosion):
-        print(f"   alpha={a:.2f} -> erosion={e:.3f}")
-        if writer:
-            writer.add_scalar('theory/coi_erosion', e, int(a * 100))
-
-    # policy comparison
-    print("\n2. Policy comparison across contamination levels")
-    policies = {
-        'fixed': lambda obs, n: PolicyRegistry.fixed(obs, n),
-        'aggressive': PolicyRegistry.aggressive,
-        'defensive': PolicyRegistry.defensive,
-        'adaptive': PolicyRegistry.adaptive,
-        'alpha_proportional': PolicyRegistry.alpha_proportional,
-    }
-
-    sweep_alphas = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
-    results = run_contamination_sweep(
-        sweep_alphas, policies, n_products=10, max_steps=100,
-        n_episodes=5, seed=seed, log_dir=log_dir)
-
-    # summarize
-    print("\n3. Summary by policy")
-    for name, res_list in results.items():
-        avg_reward = np.mean([r.reward_mean for r in res_list])
-        avg_coi = np.mean([r.coi_erosion for r in res_list])
-        print(f"   {name:20s}: avg_reward={avg_reward:.2f}, avg_coi={avg_coi:.3f}")
-
-    # save results
-    output = {
-        'theoretical': {'alphas': alphas.tolist(), 'erosion': theo_erosion.tolist()},
-        'empirical': {name: [r.to_dict() for r in res_list] for name, res_list in results.items()}}
-
-    with open(Path(log_dir) / "coi_demo_results.json", 'w') as f:
-        json.dump(output, f, indent=2)
-
-    if writer:
-        writer.close()
-
-    print(f"\nResults saved to {log_dir}/coi_demo_results.json")
-    print(f"TensorBoard: tensorboard --logdir {log_dir}")
-
-    return output
-
-
-def run_reward_mode_comparison(log_dir: str = "sim/case/thesis_simplified/runs", seed: int = 42) -> Dict:
-    """Compare different reward modes."""
-    print("=== Reward Mode Comparison ===\n")
-
-    Path(log_dir).mkdir(parents=True, exist_ok=True)
-    writer = SummaryWriter(Path(log_dir) / "reward_modes") if HAS_TB else None
-
-    reward_modes = ["revenue", "profit", "robust", "coi_aware"]
-    alpha = 0.3  # moderate contamination
-
-    results = {}
-    for mode in reward_modes:
-        print(f"  mode={mode}", end=" ")
-        env_cfg = EnvConfig(
-            n_products=10, max_steps=200, alpha_true=alpha,
-            reward_mode=mode, seed=seed)
-        env = make_env(env_cfg)
-
-        rewards, coi_vals, _, revenues = run_policy_episode(
-            env, PolicyRegistry.adaptive, n_episodes=10)
-
-        results[mode] = {
-            'reward_mean': float(np.mean(rewards)),
-            'reward_std': float(np.std(rewards)),
-            'coi_erosion': float(np.mean(coi_vals)) if coi_vals else 0.0,
-            'revenue': float(np.mean(revenues)) if revenues else 0.0}
-
-        if writer:
-            for k, v in results[mode].items():
-                writer.add_scalar(f'{mode}/{k}', v, 0)
-
-        print(f"reward={results[mode]['reward_mean']:.2f}, coi={results[mode]['coi_erosion']:.3f}")
-
-    if writer:
-        writer.close()
-
-    with open(Path(log_dir) / "reward_mode_results.json", 'w') as f:
-        json.dump(results, f, indent=2)
-
-    return results
-
-
-def run_alpha_drift_experiment(log_dir: str = "sim/case/thesis_simplified/runs", seed: int = 42) -> Dict:
-    """Test policy robustness under non-stationary contamination."""
-    print("=== Alpha Drift Experiment ===\n")
-
-    Path(log_dir).mkdir(parents=True, exist_ok=True)
-    writer = SummaryWriter(Path(log_dir) / "alpha_drift") if HAS_TB else None
-
-    drift_rates = [0.0, 0.01, 0.02, 0.05]
-    results = {}
-
-    for drift in drift_rates:
-        print(f"  drift={drift:.2f}", end=" ")
-        env_cfg = EnvConfig(
-            n_products=10, max_steps=200, alpha_true=0.2,
-            alpha_drift=drift, reward_mode="robust", seed=seed)
-        env = make_env(env_cfg)
-
-        rewards, coi_vals, alpha_errs, _ = run_policy_episode(
-            env, PolicyRegistry.adaptive, n_episodes=10)
-
-        results[f'drift_{drift}'] = {
-            'reward_mean': float(np.mean(rewards)),
-            'coi_erosion': float(np.mean(coi_vals)) if coi_vals else 0.0,
-            'alpha_tracking_error': float(np.mean(alpha_errs)) if alpha_errs else 0.0}
-
-        if writer:
-            for k, v in results[f'drift_{drift}'].items():
-                writer.add_scalar(f'drift_{drift}/{k}', v, 0)
-
-        print(f"reward={results[f'drift_{drift}']['reward_mean']:.2f}, "
-              f"alpha_err={results[f'drift_{drift}']['alpha_tracking_error']:.3f}")
-
-    if writer:
-        writer.close()
-
-    return results
-
-
-if __name__ == "__main__":
-    import argparse
-    parser = argparse.ArgumentParser(description="Run COI experiments")
-    parser.add_argument("--exp", type=str, default="coi", choices=["coi", "reward", "drift", "all"])
-    parser.add_argument("--log-dir", type=str, default="sim/case/thesis_simplified/runs")
-    parser.add_argument("--seed", type=int, default=42)
-    args = parser.parse_args()
-
-    if args.exp == "coi" or args.exp == "all":
-        run_coi_demonstration(args.log_dir, args.seed)
-
-    if args.exp == "reward" or args.exp == "all":
-        run_reward_mode_comparison(args.log_dir, args.seed)
-
-    if args.exp == "drift" or args.exp == "all":
-        run_alpha_drift_experiment(args.log_dir, args.seed)

sim/case/thesis_simplified/separability.py

@@ -1,72 +0,0 @@
-"""Behavioral separability for human/agent detection.
-
-Computes divergence signals delta_H, delta_A from session trajectories using
-transition kernel estimation and KL divergence to prototype behavioral profiles.
-"""
-from __future__ import annotations
-from typing import Dict, List, Tuple, TYPE_CHECKING
-import numpy as np
-
-if TYPE_CHECKING:
-    from .simplified import Event, Session
-
-
-# prototype behavioral kernels for human vs agent sessions
-TRANS_H = {
-    "start": {"view": 0.85, "end": 0.15},
-    "view": {"detail": 0.4, "cart": 0.3, "view": 0.2, "end": 0.1},
-    "detail": {"cart": 0.5, "view": 0.3, "end": 0.2},
-    "cart": {"purchase": 0.6, "view": 0.25, "end": 0.15},
-    "purchase": {"end": 1.0},
-}
-
-TRANS_A = {
-    "start": {"view": 0.95, "end": 0.05},
-    "view": {"detail": 0.6, "view": 0.25, "cart": 0.1, "end": 0.05},
-    "detail": {"view": 0.5, "cart": 0.15, "detail": 0.3, "end": 0.05},
-    "cart": {"view": 0.4, "purchase": 0.2, "end": 0.4},
-    "purchase": {"end": 1.0},
-}
-
-
-def kl_div(p: Dict[str, float], q: Dict[str, float], eps: float = 1e-10) -> float:
-    """KL divergence D_KL(p || q) for discrete distributions."""
-    keys = set(p.keys()) | set(q.keys())
-    return sum(p.get(k, eps) * np.log((p.get(k, eps) + eps) / (q.get(k, eps) + eps)) for k in keys)
-
-
-def build_kernel(events: List["Event"]) -> Dict[str, Dict[str, float]]:
-    """Build empirical transition kernel T' from trajectory events."""
-    trans: Dict[str, Dict[str, int]] = {}
-    prev = "start"
-    for e in events:
-        curr = e.action
-        trans.setdefault(prev, {})
-        trans[prev][curr] = trans[prev].get(curr, 0) + 1
-        prev = curr
-    return {s: {d: c / sum(dsts.values()) for d, c in dsts.items()} for s, dsts in trans.items() if sum(dsts.values()) > 0}
-
-
-def compute_divergence(session: "Session") -> Tuple[float, float]:
-    """Compute divergence signals delta_H, delta_A for session.
-
-    delta_H = mean KL(T' || T_H) across states, measures distance to human prototype
-    delta_A = mean KL(T' || T_A) across states, measures distance to agent prototype
-    """
-    kernel = build_kernel(session.events)
-    if not kernel:
-        return 0.5, 0.5
-    delta_h = sum(kl_div(kernel.get(s, {}), TRANS_H.get(s, {})) for s in kernel) / len(kernel)
-    delta_a = sum(kl_div(kernel.get(s, {}), TRANS_A.get(s, {})) for s in kernel) / len(kernel)
-    return delta_h, delta_a
-
-
-def estimate_alpha(session: "Session", beta: float = 2.0) -> float:
-    """Per-session contamination estimate alpha_hat = sigma(beta*(delta_H - delta_A)).
-
-    Returns probability session is agent-generated based on behavioral divergence.
-    """
-    dh, da = compute_divergence(session)
-    if (dh + da) <= 0:
-        return 0.5
-    return 1.0 / (1.0 + np.exp(-beta * (dh - da)))

sim/case/thesis_simplified/summarize.py

@@ -1,168 +0,0 @@
-"""Summarize TensorBoard logs into comparison tables."""
-from __future__ import annotations
-import json
-import re
-from pathlib import Path
-from collections import defaultdict
-from dataclasses import dataclass
-import pandas as pd
-
-try:
-    from tensorboard.backend.event_processing.event_accumulator import EventAccumulator
-    HAS_TB = True
-except ImportError:
-    HAS_TB = False
-
-
-@dataclass
-class RunInfo:
-    algo: str
-    alpha: float
-    reward_mode: str
-    path: Path
-
-
-def parse_run_name(name: str) -> RunInfo | None:
-    """Extract algo, alpha, reward_mode from run directory name."""
-    # patterns: ppo_a0.20_robust, cmp_fixed_a0.20, sac_a0.90_robust
-    m = re.match(r'(cmp_)?(\w+)_a([\d.]+)_?(\w+)?', name)
-    if not m:
-        return None
-    prefix, algo, alpha, mode = m.groups()
-    return RunInfo(algo=algo, alpha=float(alpha), reward_mode=mode or 'robust', path=Path())
-
-
-def load_tb_scalars(log_dir: Path, tags: list[str], reduce: str = 'last') -> dict[str, float]:
-    """Load scalar values from TensorBoard event files."""
-    if not HAS_TB:
-        return {}
-    ea = EventAccumulator(str(log_dir))
-    ea.Reload()
-    results = {}
-    for tag in tags:
-        if tag in ea.Tags().get('scalars', []):
-            events = ea.Scalars(tag)
-            if not events:
-                continue
-            vals = [e.value for e in events]
-            if reduce == 'last':
-                results[tag] = vals[-1]
-            elif reduce == 'mean':
-                results[tag] = sum(vals) / len(vals)
-            elif reduce == 'max':
-                results[tag] = max(vals)
-            elif reduce == 'min':
-                results[tag] = min(vals)
-    return results
-
-
-def load_json_results(log_dir: Path) -> dict[str, float]:
-    """Load metrics from results.json if available."""
-    results_file = log_dir / 'results.json'
-    if results_file.exists():
-        with open(results_file) as f:
-            return json.load(f)
-    return {}
-
-
-def discover_runs(base_dir: Path) -> list[RunInfo]:
-    """Find all experiment runs in base directory."""
-    runs = []
-    for d in base_dir.iterdir():
-        if not d.is_dir():
-            continue
-        info = parse_run_name(d.name)
-        if info:
-            info.path = d
-            runs.append(info)
-    return runs
-
-
-def build_tables(runs: list[RunInfo], metrics: list[str], reduce: str = 'last') -> dict[str, dict[str, pd.DataFrame]]:
-    """Build pivot tables: reward_mode -> metric -> DataFrame[alpha x algo]."""
-    # collect data: {reward_mode: {metric: {(alpha, algo): value}}}
-    data = defaultdict(lambda: defaultdict(dict))
-
-    tb_tags = [f'economics/{m}' if m in ['revenue', 'profit', 'margin'] else f'coi/{m}' if m in ['erosion', 'leakage'] else f'alpha/{m}' for m in metrics]
-    tag_map = dict(zip(tb_tags, metrics))
-
-    for run in runs:
-        # try json first (final eval metrics)
-        jm = load_json_results(run.path)
-        tb = load_tb_scalars(run.path, tb_tags, reduce)
-
-        for tag, metric in tag_map.items():
-            val = None
-            json_key = f'{metric}_mean' if metric != 'reward' else 'reward_mean'
-            if json_key in jm:
-                val = jm[json_key]
-            elif tag in tb:
-                val = tb[tag]
-            if val is not None:
-                data[run.reward_mode][metric][(run.alpha, run.algo)] = val
-
-    # convert to DataFrames
-    tables = {}
-    for mode, metrics_data in data.items():
-        tables[mode] = {}
-        for metric, vals in metrics_data.items():
-            if not vals:
-                continue
-            alphas = sorted(set(a for a, _ in vals.keys()))
-            algos = sorted(set(al for _, al in vals.keys()))
-            df = pd.DataFrame(index=alphas, columns=algos, dtype=float)
-            for (a, al), v in vals.items():
-                df.loc[a, al] = v
-            df.index.name = 'alpha'
-            tables[mode][metric] = df
-    return tables
-
-
-def format_table(df: pd.DataFrame, fmt: str = '.3f') -> str:
-    """Format DataFrame as markdown table."""
-    return df.to_markdown(floatfmt=fmt)
-
-
-def summarize(base_dir: str = 'sim/case/thesis_simplified/runs',
-              metrics: list[str] | None = None,
-              reduce: str = 'last',
-              output: str | None = None) -> dict:
-    """Generate summary tables from experiment runs."""
-    base = Path(base_dir)
-    metrics = metrics or ['revenue', 'profit', 'margin', 'erosion', 'leakage']
-
-    runs = discover_runs(base)
-    if not runs:
-        print(f"No runs found in {base}")
-        return {}
-
-    print(f"Found {len(runs)} runs")
-    tables = build_tables(runs, metrics, reduce)
-
-    lines = []
-    for mode, metric_tables in sorted(tables.items()):
-        lines.append(f"\n# Reward Mode: {mode}\n")
-        for metric, df in sorted(metric_tables.items()):
-            lines.append(f"\n## {metric}\n")
-            lines.append(format_table(df))
-            lines.append("")
-
-    report = '\n'.join(lines)
-    print(report)
-
-    if output:
-        Path(output).write_text(report)
-        print(f"\nSaved to {output}")
-
-    return tables
-
-
-if __name__ == '__main__':
-    import argparse
-    p = argparse.ArgumentParser()
-    p.add_argument('--dir', default='sim/case/thesis_simplified/runs')
-    p.add_argument('--metrics', nargs='+', default=['revenue', 'profit', 'margin', 'erosion', 'leakage'])
-    p.add_argument('--reduce', default='last', choices=['last', 'mean', 'max', 'min'])
-    p.add_argument('--output', '-o', help='save markdown to file')
-    args = p.parse_args()
-    summarize(args.dir, args.metrics, args.reduce, args.output)

sim/rl/behavior_loader/models.py

@@ -226,7 +226,6 @@ if __name__ == "__main__":
 
     agent_model = AgentBehaviorModel(agent_dir)
     agent_mdp = agent_model.build_MDP()
-    print(agent_mdp)
     print(f"AGENT... Built MDP: {agent_mdp['num_states']} states, "
           f"{sum(len(t) for t in agent_mdp['transitions'].values())} transitions")
     if not agent_mdp['states']:
@@ -235,9 +234,6 @@ if __name__ == "__main__":
 
     human_evt = aggregate_event_transitions(human_mdp)
     agent_evt = aggregate_event_transitions(agent_mdp)
-    print(agent_evt)
-
-
     common = set(human_evt.keys()) & set(agent_evt.keys())
 
     if not common:

sim/rl/engine.py

@@ -76,7 +76,8 @@ class WildPricingEngine(BasePricingEngine):
 
     def compute_prices(self, current_prices: np.ndarray, observation: Dict[str, Any]) -> np.ndarray:
         self.step_count += 1
-        demand = _extract_demand(observation, self.c.product_catalogue_size)
+        # extract demand signal (from env observation) as proxy for sales
+        demand = observation.get('demand', np.zeros(self.c.product_catalogue_size, dtype=np.float32))
         return self._update_from_demand(current_prices, demand)
 
     def _update_from_demand(self, prices: np.ndarray, sold: np.ndarray) -> np.ndarray:
@@ -140,7 +141,7 @@ class SimpleDemandEngine(BasePricingEngine):
 
     def compute_prices(self, current_prices: np.ndarray, observation: Dict[str, Any]) -> np.ndarray:
         self.step_count += 1
-        demand = _extract_demand(observation, self.c.product_catalogue_size)
+        demand = observation.get('demand', np.zeros(self.c.product_catalogue_size, dtype=np.float32))
         if self.prev_demand is None:
             self.prev_demand = demand.copy()
             return current_prices.copy()
@@ -206,7 +207,7 @@ class ThompsonSamplingEngine(BasePricingEngine):
             lo = current_prices * 0.7
             hi = current_prices * 1.3
             self.price_grid = np.linspace(lo, hi, self.n_price_levels).T
-        demand = _extract_demand(observation, self.c.product_catalogue_size)
+        demand = observation.get('demand', np.zeros(self.c.product_catalogue_size, dtype=np.float32))
         # update beliefs based on last action
         if self.last_actions is not None:
             for i in range(self.c.product_catalogue_size):
@@ -225,14 +226,3 @@ class ThompsonSamplingEngine(BasePricingEngine):
             new_prices[i] = self.price_grid[i, actions[i]]
         self.last_actions = actions
         return np.clip(new_prices, self.c.system_min_price, self.c.system_max_price).astype(np.float32)
-
-
-def _extract_demand(observation: Dict[str, Any], n: int) -> np.ndarray:
-    if "elasticity" in observation and isinstance(observation["elasticity"], dict):
-        d = observation["elasticity"].get("demand")
-        if d is not None:
-            return np.asarray(d, dtype=np.float32)
-    d = observation.get("demand")
-    if d is not None:
-        return np.asarray(d, dtype=np.float32)
-    return np.zeros(n, dtype=np.float32)

sim/rl/environment.py

@@ -1,244 +1,682 @@
-from __future__ import annotations
-
-from dataclasses import dataclass
-from typing import Any, Dict, Optional, Tuple
-
+import gymnasium as gym
+from gymnasium import spaces
 import numpy as np
+from dataclasses import dataclass
+import pandas as pd
+from types import SimpleNamespace
+from typing import Optional, Dict, Any, List, Tuple
+
+from lib.separability import load_artifacts, score_session, estimate_alpha
+from sim.rl.behavior_loader.models import AgentBehaviorModel, BehaviorModel, aggregate_event_transitions
 
 try:
-    import gymnasium as gym
-    from gymnasium import spaces
-except ImportError as e:
-    raise ImportError("sim.rl.environment requires gymnasium") from e
+    import jax
+    from sim.rl.jax_core import JAX_AVAILABLE, compile_transitions, fallback_transitions, sample_sessions, compute_metrics
+    from sim.rl.jax_core import session_features, compute_session_transitions, compute_divergences, estimate_alpha_batch
+except ImportError:
+    JAX_AVAILABLE = False
 
-from sim.case.thesis_simplified.coi import COIWindow, coi_erosion, compute_coi_window
-from sim.case.thesis_simplified.separability import estimate_alpha as estimate_session_alpha
-from sim.case.thesis_simplified.simplified import Limbo, Session, put_prices_to_market
-from sim.rl.thesis_core import aggregate_demand_by_product, aggregate_purchases, constrain_prices
-
-
-@dataclass(frozen=True)
-class BusinessLogicConstraints:
-    product_catalogue_size: int = 100
-    max_steps: int = 2000
-    sessions_per_step: int = 250
+# "learner" agent learning to optimize pricing
+# "agent" part of environment creating demand signals that learner processes
 
+base_dir = "/home/velocitatem/Documents/Projects/PHANTOM/experiments"
+human_dir, agent_dir = f"{base_dir}/collected_data/", f"{base_dir}/agents/collected_data/"
+@dataclass
+class BusinessLogicConstraints():
+    max_price_adjustment: float = 0.30
     system_max_price: float = 500.0
     system_min_price: float = 1.0
-    max_price_adjustment: float = 0.30
-    min_margin_pct: float = 0.05
-
+    product_catalogue_size: int = 100
+    episode_length: int = 2000
+    sessions_per_step: int = 250
     agent_share: float = 0.2
-    alpha_drift: float = 0.0
-    alpha_bounds: tuple[float, float] = (0.0, 0.8)
-
+    agent_recon_multiplier: float = 6.0
+    agent_purchase_probability: float = 0.20
     coi_strength: float = 0.25
+    coi_threshold: float = 4.0
+    coi_sigmoid_temp: float = 1.25
+    base_human_demand: float = 0.08
+    base_agent_demand: float = 0.05
+    human_price_elasticity: float = -1.2 # assumptions here
+    agent_price_elasticity: float = -0.6
+    w_agent_loss: float = 1.0
     w_volatility: float = 5.0
     w_estimation_error: float = 0.25
-
     seed: int = 7
 
 
-def make_env(constraints: Optional[BusinessLogicConstraints] = None) -> "PHANTOMEnv":
-    return PHANTOMEnv(constraints=constraints or BusinessLogicConstraints())
+def _sigmoid(x: np.ndarray) -> np.ndarray:
+    return 1.0 / (1.0 + np.exp(-x))
+
+EVENT_PAGE_MAP = {
+    "session_start": "/",
+    "page_view": "/",
+    "view_item_page": "/products",
+    "learn_more_about_item": "/products/details",
+    "add_item_to_cart": "/cart",
+    "checkout_start": "/checkout",
+    "purchase_complete": "/checkout",
+    "session_end": "/checkout/success",
+}
+
+# map real collected event names to canonical simulation states
+EVENT_CANONICAL_MAP = {
+    "page_view": "session_start",
+    "hover_over_paragraph": "view_item_page",
+    "hover_over_title": "view_item_page",
+    "view_item_page": "view_item_page",
+    "learn_more_about_item": "learn_more_about_item",
+    "add_item_to_cart": "add_item_to_cart",
+    "checkout_start": "purchase_complete",
+    "remove_item": "view_item_page",
+}
+
+
+def _canonicalize_transitions(raw_trans: Dict[str, Dict[str, float]]) -> Dict[str, Dict[str, float]]:
+    """Map real event transition names to canonical simulation states."""
+    canonical: Dict[str, Dict[str, float]] = {}
+    for src, dsts in raw_trans.items():
+        src_canon = EVENT_CANONICAL_MAP.get(src, src)
+        if src_canon not in canonical:
+            canonical[src_canon] = {}
+        for dst, prob in dsts.items():
+            dst_canon = EVENT_CANONICAL_MAP.get(dst, dst)
+            canonical[src_canon][dst_canon] = canonical[src_canon].get(dst_canon, 0.0) + prob
+    # re-normalize after aggregation
+    for src in canonical:
+        total = sum(canonical[src].values())
+        if total > 0:
+            canonical[src] = {k: v / total for k, v in canonical[src].items()}
+    return canonical
+
+
+class BehavioralProfile:
+    """Synthetic Markov profile used to generate interaction sessions.
+    Uses aggregate_event_transitions from models.py to build transition kernels from real data."""
+
+    def __init__(self, actor: str, purchase_probs: np.ndarray):
+        self.actor = actor
+        self.purchase_probs = np.clip(purchase_probs, 0.0, 0.95)
+        self.states = [
+            "session_start",
+            "view_item_page",
+            "learn_more_about_item",
+            "add_item_to_cart",
+            "purchase_complete",
+            "session_end",
+        ]
+        model = AgentBehaviorModel(agent_dir) if actor == "agents" else BehaviorModel(human_dir)
+        mdp = model.build_MDP()
+        raw_trans = aggregate_event_transitions(mdp) if mdp.get("transitions") else {}
+        self.transitions = _canonicalize_transitions(raw_trans) if raw_trans else self._fallback_transitions()
+        self._ensure_terminal_states()
+        self.dwell_params = self._extract_dwell_params(mdp)
+
+    def _ensure_terminal_states(self):
+        # guarantee purchase_complete leads to session_end and session_start exists
+        if "purchase_complete" not in self.transitions:
+            self.transitions["purchase_complete"] = {"session_end": 1.0}
+        elif "session_end" not in self.transitions.get("purchase_complete", {}):
+            self.transitions["purchase_complete"]["session_end"] = 1.0
+            total = sum(self.transitions["purchase_complete"].values())
+            self.transitions["purchase_complete"] = {k: v/total for k, v in self.transitions["purchase_complete"].items()}
+        if "session_start" not in self.transitions:
+            self.transitions["session_start"] = {"view_item_page": 0.7, "learn_more_about_item": 0.2, "session_end": 0.1}
+
+    def _fallback_transitions(self) -> Dict[str, Dict[str, float]]:
+        return {
+            "session_start": {"view_item_page": 0.85, "session_end": 0.15},
+            "view_item_page": {"learn_more_about_item": 0.4, "add_item_to_cart": 0.3, "view_item_page": 0.2, "session_end": 0.1},
+            "learn_more_about_item": {"add_item_to_cart": 0.5, "view_item_page": 0.3, "session_end": 0.2},
+            "add_item_to_cart": {"purchase_complete": 0.6, "view_item_page": 0.25, "session_end": 0.15},
+            "purchase_complete": {"session_end": 1.0},
+        }
+
+    def _extract_dwell_params(self, mdp: Dict) -> Dict[str, Tuple[float, float]]:
+        state_vals = mdp.get("state_values", {})
+        params = {}
+        for state in self.states:
+            # try canonical and raw state names
+            val = state_vals.get(state, 0.5)
+            for raw, canon in EVENT_CANONICAL_MAP.items():
+                if canon == state and raw in state_vals:
+                    val = state_vals[raw]
+                    break
+            shape = 1.5 + val * 2.0
+            scale = 0.8 + (1.0 - val) * 1.2
+            params[state] = (shape, scale)
+        return params
+
+    def _transition_probs(self, state: str, product_idx: int) -> Dict[str, float]:
+        probs = dict(self.transitions.get(state, {"session_end": 1.0}))
+        if state == "add_item_to_cart":
+            base = probs.get("purchase_complete", 0.0)
+            demand_factor = float(self.purchase_probs[int(product_idx)])
+            if self.actor == "agents":
+                demand_factor *= 0.7
+            adjusted = np.clip(base * 0.5 + demand_factor * 0.5, 0.0, 0.95)
+            remainder = max(1e-6, 1.0 - adjusted)
+            other_total = sum(v for k, v in probs.items() if k != "purchase_complete")
+            scale = remainder / max(other_total, 1e-6)
+            for key in probs:
+                if key == "purchase_complete":
+                    probs[key] = adjusted
+                else:
+                    probs[key] = probs[key] * scale
+        total = sum(probs.values())
+        if total <= 0:
+            return {"session_end": 1.0}
+        return {state: val / total for state, val in probs.items()}
+
+    def sample_session(
+        self,
+        rng: np.random.Generator,
+        session_id: str,
+        prices: np.ndarray,
+        unit_cost: np.ndarray,
+    ) -> Tuple[List[Dict[str, Any]], List[SimpleNamespace]]:
+        """Generate a single session trajectory respecting business constraints."""
+        events: List[Dict[str, Any]] = []
+        feature_events: List[SimpleNamespace] = []
+        state = "session_start"
+        t = 0.0
+        product_idx = int(rng.integers(0, len(prices)))
+        product_id = f"product-{product_idx:04d}"
+
+
+        # enforce price >= cost constraint (lipschitz bound on pricing)
+        # This is a sort of last resort to not let an pricing learner go rogue
+        cost = float(unit_cost[product_idx])
+        constrained_price = max(float(prices[product_idx]), cost * 1.05)  # 5% min margin
+
+        while state != "session_end" and len(events) < 40:
+            if state != "session_start":
+                row = {
+                    "session_id": session_id,
+                    "actor": "agent" if self.actor == "agents" else "human",
+                    "eventName": state,
+                    "product_idx": product_idx,
+                    "productId": product_id,
+                    "price_offered": constrained_price,
+                    "price_paid": 0.0,
+                    "page": EVENT_PAGE_MAP.get(state, "/"),
+                    "ts": t,
+                    "unit_cost": cost,
+                    "base_price": float(prices[product_idx]),
+                }
+                if state == "purchase_complete":
+                    noise = float(rng.normal(0.0, 0.015))
+                    row["price_paid"] = max(constrained_price * (1.0 + noise), cost)
+                events.append(row)
+                feature_events.append(
+                    SimpleNamespace(
+                        eventName=row["eventName"],
+                        page=row["page"],
+                        productId=row["productId"],
+                        ts=row["ts"],
+                    )
+                )
+
+            transitions = self._transition_probs(state, product_idx)
+            next_state = rng.choice(list(transitions.keys()), p=list(transitions.values()))
+            shape, scale = self.dwell_params.get(state, (2.0, 1.0))
+            dwell = max(0.3, rng.gamma(shape=shape, scale=scale))
+            t += dwell
+            state = next_state
+
+        return events, feature_events
+
+
+def _load_behavioral_profile(actor: str, demand_forcing: np.ndarray) -> BehavioralProfile:
+    """returns a behavioral profile for generating synthetic sessions
+    actor: 'humans' or 'agents'
+    demand_forcing: per-product purchase probabilities used to weight interactions
+    """
+    return BehavioralProfile(actor, demand_forcing)
+
+
+class CommercePlatform:
+    """state management for the environment, simulates demand"""
+    def __init__(self, product_catalogue_size: int, max_price: float, min_price: float, constraints: BusinessLogicConstraints):
+        self.product_catalogue_size = product_catalogue_size
+        self.max_price = max_price
+        self.min_price = min_price
+        self.constraints = constraints
+        self.simulation_history: List[Dict[str, Any]] = []
+        self._rng = np.random.default_rng(constraints.seed)
+        self._last_interaction_df: pd.DataFrame = pd.DataFrame()
+        self.unit_cost = np.random.uniform(low=15.0, high=60.0, size=(self.product_catalogue_size,)).astype(np.float32)
+        self.base_price = np.random.uniform(low=60.0, high=140.0, size=(self.product_catalogue_size,)).astype(np.float32)
+        self.alpha_hat = constraints.agent_share
+        try:
+            self.separability_artifacts = load_artifacts()
+        except FileNotFoundError:
+            self.separability_artifacts = None
+
+    def setup_true_demand(self, prices: np.ndarray) -> Dict[str, np.ndarray]:
+        p = np.clip(prices, self.min_price, self.max_price)
+        cost = np.clip(self.unit_cost, self.min_price * 0.2, self.max_price)
+        margin = np.clip((p - cost) / np.maximum(cost, 1e-3), -0.9, 2.0)
+        # isoelastic demand approximation
+        human_prob = self.constraints.base_human_demand * np.exp(self.constraints.human_price_elasticity * margin)
+        agent_prob = self.constraints.base_agent_demand * np.exp(self.constraints.agent_price_elasticity * margin)
+        return {
+            "human_purchase_prob": np.clip(human_prob, 0.0, 0.95),
+            "agent_purchase_prob": np.clip(agent_prob, 0.0, 0.95),
+        }
+
+    def _simulate_sessions(self, prices: np.ndarray) -> Tuple[pd.DataFrame, Dict[str, Any]]:
+        demand = self.setup_true_demand(prices)
+        T = self.constraints.sessions_per_step
+        effective_share = float(np.clip(self.alpha_hat, 0.0, 0.95))
+        n_agent_sessions = max(1, int(round(T * effective_share)))
+        n_human_sessions = max(1, T - n_agent_sessions)
+
+        session_map = {
+            "humans": n_human_sessions,
+            "agents": n_agent_sessions,
+        }
+        pprob_map = {
+            "humans": demand["human_purchase_prob"],
+            "agents": demand["agent_purchase_prob"],
+        }
+
+        rows: List[Dict[str, Any]] = []
+        session_scores: List[Dict[str, float]] = []
+        demand_human = np.zeros_like(prices, dtype=np.float32)
+        demand_agent = np.zeros_like(prices, dtype=np.float32)
+
+        for actor, n_sessions in session_map.items():
+            profile = _load_behavioral_profile(actor, pprob_map[actor])
+            for idx in range(n_sessions):
+                session_id = f"{actor}_{idx:06d}"
+                session_rows, feature_events = profile.sample_session(
+                    self._rng, session_id, prices, self.unit_cost
+                )
+                rows.extend(session_rows)
+                if session_rows:
+                    df_session = pd.DataFrame(session_rows)
+                    purchases = df_session[df_session["eventName"] == "purchase_complete"]
+                    if not purchases.empty:
+                        counts = purchases.groupby("product_idx").size()
+                        if actor == "agents":
+                            demand_agent[counts.index.to_numpy(dtype=int)] += counts.to_numpy(dtype=np.float32)
+                        else:
+                            demand_human[counts.index.to_numpy(dtype=int)] += counts.to_numpy(dtype=np.float32)
+                if self.separability_artifacts and feature_events:
+                    score = score_session(feature_events, self.separability_artifacts)
+                    session_scores.append(score)
+
+        interactions_df = pd.DataFrame(rows)
+        diagnostics = {
+            "alpha_hat": float(self.alpha_hat),
+            "session_scores": session_scores,
+            "demand_human": demand_human,
+            "demand_agent": demand_agent,
+        }
+
+        if session_scores:
+            alphas = [
+                estimate_alpha(s["prob_agent"], s["delta_h"], s["delta_a"], temperature=2.0)
+                for s in session_scores
+            ]
+            mean_alpha = float(np.mean(alphas))
+            # exponential moving average for stability
+            self.alpha_hat = 0.7 * self.alpha_hat + 0.3 * mean_alpha
+            diagnostics.update(
+                {
+                    "alpha_hat": float(self.alpha_hat),
+                    "delta_h_mean": float(np.mean([s["delta_h"] for s in session_scores])),
+                    "delta_a_mean": float(np.mean([s["delta_a"] for s in session_scores])),
+                    "prob_agent_mean": float(np.mean([s["prob_agent"] for s in session_scores])),
+                }
+            )
+
+        self._last_interaction_df = interactions_df
+        return interactions_df, diagnostics
+
+    def compute_interaction_features(self, interaction_df: pd.DataFrame) -> Dict[str, float]:
+        if interaction_df.empty:
+            return {
+                "revenue_observed": 0.0,
+                "revenue_oracle": 0.0,
+                "agent_loss": 0.0,
+                "true_human_purchases": 0.0,
+                "true_agent_purchases": 0.0,
+                "mean_sale_price": 0.0,
+                "look_to_book": 0.0,
+                "coi": 0.0,
+                "expected_premium": 0.0,
+            }
+
+        purchases = interaction_df[interaction_df["eventName"] == "purchase_complete"]
+        human_purchases = purchases[purchases["actor"] == "human"]
+        agent_purchases = purchases[purchases["actor"] == "agent"]
+
+        revenue_observed = float(purchases["price_paid"].sum())
+        revenue_oracle = float(purchases["base_price"].sum())
+        agent_loss = float((agent_purchases["base_price"] - agent_purchases["price_paid"]).sum())
+
+        mean_sale_price = float(purchases["price_paid"].mean()) if not purchases.empty else 0.0
+        views = float((interaction_df["eventName"] == "view_item_page").sum())
+        look_to_book = float(views / (len(purchases) + 1e-6))
+        true_human = float(len(human_purchases))
+        true_agent = float(len(agent_purchases))
+
+        human_prices = human_purchases["price_offered"] if not human_purchases.empty else pd.Series(dtype=float)
+        human_costs = human_purchases["unit_cost"] if not human_purchases.empty else pd.Series(dtype=float)
+        human_base = human_purchases["base_price"] if not human_purchases.empty else pd.Series(dtype=float)
+        coi = 0.0
+        if not human_prices.empty and not human_costs.empty:
+            # COI = E[P] - p_min where p_min is cost, accounting for expected premium (base - realized)
+            margin = human_prices.mean() - human_costs.mean()
+            expected_premium = human_base.mean() - human_prices.mean() if not human_base.empty else 0.0
+            coi = float(np.maximum(0.0, margin - expected_premium * 0.5))
+
+        return {
+            "revenue_observed": revenue_observed,
+            "revenue_oracle": revenue_oracle,
+            "agent_loss": agent_loss,
+            "true_human_purchases": true_human,
+            "true_agent_purchases": true_agent,
+            "mean_sale_price": mean_sale_price,
+            "look_to_book": look_to_book,
+            "coi": coi,
+            "expected_premium": float(expected_premium) if not human_base.empty else 0.0,
+        }
+
+    def _session_feature_table(self, df: pd.DataFrame) -> pd.DataFrame:
+        """Extract per-session behavioral features for separability analysis."""
+        if df.empty:
+            return pd.DataFrame()
+        g = df.groupby("session_id", sort=False)
+        session_duration = g["ts"].max() - g["ts"].min()
+        total_interactions = g.size()
+        avg_time_between = g["ts"].apply(lambda x: float(np.diff(np.sort(x.to_numpy())).mean()) if len(x) > 1 else 0.0)
+        interaction_velocity = total_interactions / (session_duration + 1e-6)
+        views = g.apply(lambda x: int((x["eventName"] == "view_item_page").sum()), include_groups=False)
+        cart_adds = g.apply(lambda x: int((x["eventName"] == "add_item_to_cart").sum()), include_groups=False)
+        purchases = g.apply(lambda x: int((x["eventName"] == "purchase_complete").sum()), include_groups=False)
+        learn_more = g.apply(lambda x: int((x["eventName"] == "learn_more_about_item").sum()), include_groups=False)
+        conversion_rate = purchases / (views + 1e-6)
+        is_agent = g["actor"].apply(lambda s: bool((s == "agent").any()), include_groups=False)
+        # price sensitivity features
+        price_variance = g["price_offered"].var().fillna(0.0)
+        avg_price_seen = g["price_offered"].mean().fillna(0.0)
+        products_viewed = g["product_idx"].nunique()
+
+        return pd.DataFrame({
+            "session_duration_sec": session_duration.astype(float),
+            "avg_time_between_events": avg_time_between.astype(float),
+            "total_interactions": total_interactions.astype(int),
+            "interaction_velocity": interaction_velocity.astype(float),
+            "item_views": views.astype(int),
+            "cart_adds": cart_adds.astype(int),
+            "purchases": purchases.astype(int),
+            "learn_more_clicks": learn_more.astype(int),
+            "conversion_rate": conversion_rate.astype(float),
+            "price_variance": price_variance.astype(float),
+            "avg_price_seen": avg_price_seen.astype(float),
+            "products_viewed": products_viewed.astype(int),
+            "is_agent": is_agent.astype(bool),
+        }).reset_index()
+
+    def get_interaction_data(self) -> np.ndarray:
+        if self._last_interaction_df.empty:
+            return np.array([], dtype=object)
+        return self._last_interaction_df.to_dict(orient="records")
 
 
 class PHANTOMEnv(gym.Env):
-    metadata = {"render_modes": ["human", "ansi"]}
+    metadata = {"render_modes": []}
 
-    def __init__(self, constraints: Optional[BusinessLogicConstraints] = None):
+    def __init__(self, constraints: Optional[BusinessLogicConstraints] = None, use_jax: bool = True):
         super().__init__()
-        self.c = constraints or BusinessLogicConstraints()
-        self.n = int(self.c.product_catalogue_size)
+        self.constraints = constraints if isinstance(constraints, BusinessLogicConstraints) else BusinessLogicConstraints()
+        self.use_jax = use_jax and JAX_AVAILABLE
+        self.action_space = spaces.Box(low=-self.constraints.max_price_adjustment,
+                                       high=self.constraints.max_price_adjustment,
+                                       shape=(self.constraints.product_catalogue_size,), dtype=np.float32)
+        n_products = self.constraints.product_catalogue_size
+        self.observation_space = spaces.Dict({
+            "elasticity": spaces.Dict({
+                "price": spaces.Box(
+                    low=np.full((n_products,), self.constraints.system_min_price, dtype=np.float32),
+                    high=np.full((n_products,), self.constraints.system_max_price, dtype=np.float32),
+                    dtype=np.float32),
+                "demand": spaces.Box(
+                    low=np.zeros((n_products,), dtype=np.float32),
+                    high=np.full((n_products,), 1e6, dtype=np.float32),
+                    dtype=np.float32),
+            }),
+            "market": spaces.Dict({
+                "alpha_hat": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
+                "revenue_rate": spaces.Box(low=0.0, high=1e6, shape=(1,), dtype=np.float32),
+                "conversion_rate": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
+                "price_volatility": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
+            }),
+            "cost": spaces.Box(low=0.0, high=self.constraints.system_max_price, shape=(n_products,), dtype=np.float32),
+        })
+        self.commerce_platform = CommercePlatform(
+            product_catalogue_size=self.constraints.product_catalogue_size,
+            max_price=self.constraints.system_max_price,
+            min_price=self.constraints.system_min_price,
+            constraints=self.constraints)
+        self._rng = np.random.default_rng(self.constraints.seed)
+        self.t = 0
+        self._prev_prices: Optional[np.ndarray] = None
+        self.state: Dict[str, Any] = {}
+        self._jax_key = None
+        self._jax_trans = None
+        if self.use_jax:
+            self._jax_key = jax.random.PRNGKey(self.constraints.seed)
+            self._init_jax_transitions()
 
-        self._rng = np.random.default_rng(self.c.seed)
-        self._t = 0
-        self._alpha_true = float(self.c.agent_share)
-        self._alpha_hat = float(self.c.agent_share)
-        self._costs = np.zeros(self.n, dtype=np.float32)
-        self._refs = np.zeros(self.n, dtype=np.float32)
-        self._prices: Optional[np.ndarray] = None
-        self._last_sessions: list[Session] = []
-        self._last_coi: COIWindow | None = None
-        self._limbo = Limbo()
-
-        self.action_space = spaces.Box(
-            low=np.full((self.n,), self.c.system_min_price, dtype=np.float32),
-            high=np.full((self.n,), self.c.system_max_price, dtype=np.float32),
-            dtype=np.float32,
-        )
-        self.observation_space = spaces.Dict(
-            {
-                "elasticity": spaces.Dict(
-                    {
-                        "price": spaces.Box(
-                            low=np.full((self.n,), self.c.system_min_price, dtype=np.float32),
-                            high=np.full((self.n,), self.c.system_max_price, dtype=np.float32),
-                            dtype=np.float32,
-                        ),
-                        "demand": spaces.Box(
-                            low=np.zeros((self.n,), dtype=np.float32),
-                            high=np.full((self.n,), 1e9, dtype=np.float32),
-                            dtype=np.float32,
-                        ),
-                    }
-                ),
-                "market": spaces.Dict(
-                    {
-                        "alpha_hat": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
-                        "revenue_rate": spaces.Box(low=0.0, high=1e12, shape=(1,), dtype=np.float32),
-                        "conversion_rate": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
-                        "price_volatility": spaces.Box(low=0.0, high=1.0, shape=(1,), dtype=np.float32),
-                    }
-                ),
-                "cost": spaces.Box(
-                    low=np.zeros((self.n,), dtype=np.float32),
-                    high=np.full((self.n,), self.c.system_max_price, dtype=np.float32),
-                    dtype=np.float32,
-                ),
-            }
-        )
-
-    def _reset_catalogue(self) -> None:
-        self._costs = self._rng.uniform(15.0, 60.0, size=self.n).astype(np.float32)
-        margins = self._rng.uniform(0.2, 0.6, size=self.n).astype(np.float32)
-        self._refs = (self._costs * (1.0 + margins)).astype(np.float32)
-        self._prices = self._refs.copy()
-
-    def _observe_market(
-        self, prices: np.ndarray
-    ) -> tuple[list[Session], Dict[str, float], np.ndarray, np.ndarray, float, float, int]:
-        sessions, demand_map = put_prices_to_market(
-            prices,
-            costs=self._costs,
-            alpha=self._alpha_true,
-            n_sessions=int(self.c.sessions_per_step),
-            seed=int(self._rng.integers(0, 2**31 - 1)),
-        )
-        demand_by_product = aggregate_demand_by_product(sessions, demand_map, self.n)
-        purchases, revenue, cost, n_agents = aggregate_purchases(sessions, self._costs, self.n)
-        conversion = float(np.sum(purchases) / max(len(sessions), 1))
-        return sessions, demand_map, demand_by_product, purchases, revenue, cost, n_agents
-
-    def _update_alpha_hat(self, sessions: list[Session]) -> float:
-        scores = [estimate_session_alpha(s) for s in sessions if s.events]
-        if not scores:
-            return self._alpha_hat
-        alpha_step = float(np.mean(scores))
-        self._alpha_hat = 0.8 * self._alpha_hat + 0.2 * alpha_step
-        self._alpha_hat = float(np.clip(self._alpha_hat, 0.0, 1.0))
-        return self._alpha_hat
-
-    def _reward(self, prices: np.ndarray, revenue: float, cost: float, volatility: float) -> float:
-        profit = float(revenue - cost)
-        coi_leak = float(self._last_coi.leak) if self._last_coi else 0.0
-        alpha_err = abs(self._alpha_hat - self._alpha_true)
-        return profit - self.c.coi_strength * coi_leak - self.c.w_volatility * volatility - self.c.w_estimation_error * alpha_err
-
-    def _build_obs(
-        self,
-        prices: np.ndarray,
-        demand_by_product: np.ndarray,
-        revenue: float,
-        conversion: float,
-        volatility: float,
-    ) -> Dict[str, Any]:
-        return {
-            "elasticity": {"price": prices.astype(np.float32), "demand": demand_by_product.astype(np.float32)},
-            "market": {
-                "alpha_hat": np.array([self._alpha_hat], dtype=np.float32),
-                "revenue_rate": np.array([revenue], dtype=np.float32),
-                "conversion_rate": np.array([conversion], dtype=np.float32),
-                "price_volatility": np.array([volatility], dtype=np.float32),
-            },
-            "cost": self._costs.astype(np.float32),
-        }
+    def _init_jax_transitions(self):
+        try:
+            human_profile = _load_behavioral_profile("humans", np.ones(self.constraints.product_catalogue_size) * 0.1)
+            agent_profile = _load_behavioral_profile("agents", np.ones(self.constraints.product_catalogue_size) * 0.1)
+            self._jax_trans = compile_transitions(human_profile, agent_profile).to_jax()
+        except Exception:
+            self._jax_trans = fallback_transitions().to_jax()
 
     def reset(self, seed: Optional[int] = None, options: Optional[dict] = None):
         super().reset(seed=seed)
         if seed is not None:
             self._rng = np.random.default_rng(seed)
-        self._t = 0
-        self._alpha_true = float(np.clip(self.c.agent_share, *self.c.alpha_bounds))
-        self._alpha_hat = float(self.c.agent_share)
-        self._reset_catalogue()
-        self._limbo = Limbo()
-        self._last_sessions = []
-        self._last_coi = None
-
-        prices = self._prices if self._prices is not None else np.zeros(self.n, dtype=np.float32)
-        obs = self._build_obs(prices, np.zeros(self.n, dtype=np.float32), 0.0, 0.0, 0.0)
-        return obs, {"alpha_true": self._alpha_true}
-
-    def step(self, action: np.ndarray) -> Tuple[Dict[str, Any], float, bool, bool, Dict[str, Any]]:
-        if self._prices is None:
-            raise RuntimeError("reset() must be called before step()")
-
-        prev = self._prices
-        prices = constrain_prices(
-            prev,
-            np.asarray(action, dtype=np.float32),
-            costs=self._costs,
-            min_price=float(self.c.system_min_price),
-            max_price=float(self.c.system_max_price),
-            max_adjustment=float(self.c.max_price_adjustment),
-            min_margin_pct=float(self.c.min_margin_pct),
-        )
-        self._prices = prices
-        self._limbo.add_update("prices", prices)
-
-        sessions, demand_map, demand_by_product, purchases, revenue, cost, n_agents = self._observe_market(prices)
-        self._last_sessions = sessions
-        self._limbo.add_update("demand", demand_map)
-
-        self._update_alpha_hat(self._last_sessions)
-        self._last_coi = compute_coi_window(self._last_sessions, self._costs, demand_mapping=demand_map)
-
-        self._alpha_true = float(np.clip(self._alpha_true + self.c.alpha_drift, *self.c.alpha_bounds))
-        volatility = float(np.std((prices - prev) / (prev + 1e-6)))
-        reward = float(self._reward(prices, revenue, cost, volatility))
-        conversion = float(np.sum(purchases) / max(len(self._last_sessions), 1))
-
-        self._t += 1
-        terminated = self._t >= int(self.c.max_steps)
-
-        obs = self._build_obs(prices, demand_by_product, revenue, conversion, min(volatility, 1.0))
-        info = {
-            "step": self._t,
-            "reward": reward,
-            "revenue": float(revenue),
-            "profit": float(revenue - cost),
-            "n_sessions": int(self.c.sessions_per_step),
-            "n_agents": int(n_agents),
-            "alpha_true": float(self._alpha_true),
-            "alpha_hat": float(self._alpha_hat),
-            "alpha_error": float(abs(self._alpha_hat - self._alpha_true)),
-            "price_std": float(np.std(prices)),
-            "price_volatility": float(volatility),
+            self.commerce_platform._rng = np.random.default_rng(seed)
+            if self.use_jax:
+                self._jax_key = jax.random.PRNGKey(seed)
+        self.commerce_platform.alpha_hat = self.constraints.agent_share
+        self.t = 0
+        init_prices = self._rng.uniform(
+            low=60.0,
+            high=140.0,
+            size=(self.constraints.product_catalogue_size,),
+        ).astype(np.float32)
+        self.commerce_platform.unit_cost = self._rng.uniform(
+            low=15.0,
+            high=60.0,
+            size=(self.constraints.product_catalogue_size,),
+        ).astype(np.float32)
+        self.commerce_platform.base_price = init_prices.copy()
+        self._prev_prices = init_prices.copy()
+        self.state = {
+            "elasticity": {
+                "price": init_prices,
+                "demand": np.zeros((self.constraints.product_catalogue_size,), dtype=np.float32),
+            },
+            "market": {
+                "alpha_hat": np.array([self.constraints.agent_share], dtype=np.float32),
+                "revenue_rate": np.array([0.0], dtype=np.float32),
+                "conversion_rate": np.array([0.0], dtype=np.float32),
+                "price_volatility": np.array([0.0], dtype=np.float32),
+            },
+            "cost": self.commerce_platform.unit_cost.astype(np.float32),
         }
-        if self._last_coi is not None:
+        return self.state, {}
+
+    def _step_jax(self, new_prices: np.ndarray) -> Tuple[Dict, Dict]:
+        self._jax_key, subkey = jax.random.split(self._jax_key)
+        alpha = float(np.clip(self.commerce_platform.alpha_hat, 0.0, 0.95))
+        n_agent = max(1, int(self.constraints.sessions_per_step * alpha))
+        n_human = max(1, self.constraints.sessions_per_step - n_agent)
+        batch = sample_sessions(subkey, self._jax_trans, n_human, n_agent, len(new_prices))
+        sim = compute_metrics(batch, new_prices, self.commerce_platform.unit_cost, self.commerce_platform.base_price)
+        result = {"revenue_observed": sim.revenue, "revenue_oracle": sim.revenue_oracle,
+                  "agent_loss": sim.agent_loss, "coi": sim.coi, "look_to_book": sim.look_to_book,
+                  "mean_sale_price": sim.mean_sale_price, "true_human_purchases": sim.n_human_purchases,
+                  "true_agent_purchases": sim.n_agent_purchases}
+        diagnostics = {"demand_human": sim.demand_human, "demand_agent": sim.demand_agent, "alpha_hat": alpha}
+        return result, diagnostics
+
+    def step(self, action: np.ndarray):
+        self.t += 1
+        base_prices = self.state["elasticity"]["price"].astype(np.float32)
+        new_prices = np.clip(base_prices * (1.0 + action.astype(np.float32)),
+                           self.constraints.system_min_price,
+                           self.constraints.system_max_price).astype(np.float32)
+
+        self.state["elasticity"]["price"] = new_prices
+        if self.use_jax:
+            result, diagnostics = self._step_jax(new_prices)
+        else:
+            interactions_df, diagnostics = self.commerce_platform._simulate_sessions(new_prices)
+            result = self.commerce_platform.compute_interaction_features(interactions_df)
+        COI = float(result.get("coi", 0.0))
+
+        demand_vector = diagnostics.get("demand_human", np.zeros_like(new_prices)) + diagnostics.get(
+            "demand_agent", np.zeros_like(new_prices)
+        )
+        self.state["elasticity"]["demand"] = demand_vector.astype(np.float32)
+
+        volatility = 0.0 if self._prev_prices is None else \
+            float(np.mean(np.abs((new_prices - self._prev_prices) / (self._prev_prices + 1e-6))))
+        self._prev_prices = new_prices.copy()
+
+        # update market observation features
+        total_demand = float(np.sum(demand_vector))
+        total_purchases = float(result.get("true_human_purchases", 0.0) + result.get("true_agent_purchases", 0.0))
+        conv_rate = total_purchases / max(total_demand, 1.0)
+        self.state["market"] = {
+            "alpha_hat": np.array([float(diagnostics.get("alpha_hat", self.commerce_platform.alpha_hat))], dtype=np.float32),
+            "revenue_rate": np.array([float(result.get("revenue_observed", 0.0))], dtype=np.float32),
+            "conversion_rate": np.array([float(np.clip(conv_rate, 0.0, 1.0))], dtype=np.float32),
+            "price_volatility": np.array([float(volatility)], dtype=np.float32),
+        }
+        self.state["cost"] = self.commerce_platform.unit_cost.astype(np.float32)
+
+        # extract metrics with safe defaults for incomplete simulation
+        revenue_observed = float(result.get("revenue_observed", 0.0))
+        agent_loss = float(result.get("agent_loss", 0.0))
+
+        reward = (revenue_observed
+                  - COI
+                  - self.constraints.w_agent_loss * agent_loss
+                  - self.constraints.w_volatility * volatility
+                  - self.constraints.w_estimation_error)
+
+        terminated = self.t >= self.constraints.episode_length
+        info = {
+            "t": self.t,
+            "revenue_observed": revenue_observed,
+            "revenue_oracle": float(result.get("revenue_oracle", revenue_observed)),
+            "agent_loss": agent_loss,
+            "ux_volatility": volatility,
+            "look_to_book": float(result.get("look_to_book", 0.0)),
+            "mean_sale_price": float(result.get("mean_sale_price", 0.0)),
+            "true_human_purchases_total": float(result.get("true_human_purchases", 0.0)),
+            "true_agent_purchases_total": float(result.get("true_agent_purchases", 0.0)),
+            "coi": COI,
+            "alpha_hat": diagnostics.get("alpha_hat", self.commerce_platform.alpha_hat),
+            "mean_human_demand": float(np.mean(diagnostics.get("demand_human", np.zeros_like(new_prices)))),
+            "mean_agent_demand": float(np.mean(diagnostics.get("demand_agent", np.zeros_like(new_prices)))),
+        }
+        if "delta_h_mean" in diagnostics:
             info.update(
                 {
-                    "coi_policy": float(self._last_coi.policy),
-                    "coi_agent": float(self._last_coi.agent),
-                    "coi_leakage": float(self._last_coi.leak),
-                    "coi_survival": float(self._last_coi.survival_ratio),
-                    "coi_erosion": float(coi_erosion(self._last_coi.policy, self._last_coi.agent)),
+                    "delta_h_mean": diagnostics["delta_h_mean"],
+                    "delta_a_mean": diagnostics["delta_a_mean"],
+                    "prob_agent_mean": diagnostics["prob_agent_mean"],
                 }
             )
-        return obs, reward, terminated, False, info
+        return self.state, float(reward), terminated, False, info
 
-    def render(self, mode: str = "human") -> str | None:
-        if self._prices is None:
-            return None
-        out = (
-            f"t={self._t}/{self.c.max_steps} "
-            f"alpha_true={self._alpha_true:.3f} alpha_hat={self._alpha_hat:.3f} "
-            f"price_std={float(np.std(self._prices)):.2f}"
-        )
-        if mode == "human":
-            print(out)
-        return out
 
-    def close(self) -> None:
-        return
+if __name__ == "__main__":
+    import matplotlib.pyplot as plt
+    from collections import defaultdict
+
+    env = PHANTOMEnv(constraints=BusinessLogicConstraints())
+    obs, _ = env.reset(seed=42)
+    metrics = defaultdict(list)
+    total_reward = 0.0
+    done = False
+
+    while not done:
+        action = env.action_space.sample()
+        obs, reward, done, _, info = env.step(action)
+        total_reward += reward
+        p_mean = float(np.mean(obs["elasticity"]["price"]))
+        q_mean = float(np.mean(obs["elasticity"]["demand"]))
+        p_std = float(np.std(obs["elasticity"]["price"]))
+
+        metrics['t'].append(info['t'])
+        metrics['price_mean'].append(p_mean)
+        metrics['price_std'].append(p_std)
+        metrics['demand_mean'].append(q_mean)
+        metrics['revenue_observed'].append(info['revenue_observed'])
+        metrics['revenue_oracle'].append(info['revenue_oracle'])
+        metrics['agent_loss'].append(info['agent_loss'])
+        metrics['ux_volatility'].append(info['ux_volatility'])
+        metrics['look_to_book'].append(info['look_to_book'])
+        metrics['reward'].append(reward)
+        metrics['human_purchases'].append(info['true_human_purchases_total'])
+        metrics['agent_purchases'].append(info['true_agent_purchases_total'])
+        metrics['coi'].append(info.get('coi', 0.0))
+        metrics['alpha_hat'].append(info.get('alpha_hat', env.commerce_platform.alpha_hat))
+        metrics['mean_human_demand'].append(info.get('mean_human_demand', 0.0))
+        metrics['mean_agent_demand'].append(info.get('mean_agent_demand', 0.0))
+        metrics['delta_h_mean'].append(info.get('delta_h_mean', 0.0))
+        metrics['delta_a_mean'].append(info.get('delta_a_mean', 0.0))
+        metrics['prob_agent_mean'].append(info.get('prob_agent_mean', 0.0))
+
+        if info['t'] % 20 == 0 or done:
+            print(f"t={info['t']:03d} p={p_mean:6.2f}±{p_std:4.2f} q={q_mean:6.2f} "
+                  f"rev={info['revenue_observed']:7.2f} oracle={info['revenue_oracle']:7.2f} "
+                  f"loss={info['agent_loss']:6.2f} ux={info['ux_volatility']:.3f} "
+                  f"coi={info.get('coi', 0.0):6.2f} alpha={info.get('alpha_hat', 0.0):4.2f} "
+                  f"ltb={info['look_to_book']:5.2f} r={reward:7.2f}")
+
+    print(f"total_reward={total_reward:.2f}")
+
+    fig, axes = plt.subplots(3, 4, figsize=(18, 12))
+    fig.suptitle('PHANTOM Environment Run', fontsize=14, fontweight='bold')
+
+    plot_configs = [
+        ('price_mean', 'Mean Price', 'Price'),
+        ('demand_mean', 'Mean Demand (All)', 'Demand'),
+        ('mean_human_demand', 'Mean Human Demand', 'Count'),
+        ('mean_agent_demand', 'Mean Agent Demand', 'Count'),
+        ('revenue_observed', 'Revenue (Observed)', 'Revenue'),
+        ('agent_loss', 'Agent Loss (Oracle - Observed)', 'Loss'),
+        ('coi', 'Cost of Information', 'COI'),
+        ('alpha_hat', 'Estimated α̂', 'alpha'),
+        ('ux_volatility', 'UX Volatility (Price Change)', 'Volatility'),
+        ('look_to_book', 'Look-to-Book Ratio', 'Ratio'),
+        ('reward', 'Step Reward', 'Reward'),
+        ('prob_agent_mean', 'Avg Agent Probability', 'Probability'),
+    ]
+
+    for idx, (key, title, ylabel) in enumerate(plot_configs):
+        ax = axes[idx // 4, idx % 4]
+        ax.plot(metrics['t'], metrics[key], color='blue', alpha=0.7, linewidth=1.5)
+        ax.set_xlabel('Step')
+        ax.set_ylabel(ylabel)
+        ax.set_title(title, fontsize=10, fontweight='bold')
+        ax.grid(True, alpha=0.3)
+
+    plt.tight_layout()
+    plt.savefig('phantom_env_comparison.png', dpi=150, bbox_inches='tight')
+    print("Plot saved to phantom_env_comparison.png")
+    plt.show()