import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
import math

class Simulation:
    def __init__(self):
        self.n_steps = 9000
        self.ethical_load = np.sin(np.linspace(0, 4 * np.pi, self.n_steps))
        self.threshold = 0.35
        self.noise_level = 0.11
        self.state = -1
        self.scar = 0.0
        self.hesitation_prob = 0.95
        self.latent_entropy = 0.0
        self.cumulative_work = 0.0
        
    def run(self):
        x = self.ethical_load
        y = np.zeros(self.n_steps)
        work = np.zeros(self.n_steps)
        hesitating = False
        hes_timer = 0
        
        for i in range(self.n_steps):
            m = x[i] + np.random.normal(0, self.noise_level)
            
            # Hysteresis relaxation
            self.scar *= math.exp(-1.0 / 260.0)
            
            if not hesitating:
                if self.state == -1 and m > self.threshold:
                    if np.random.random() < self.hesitation_prob:
                        hesitating = True
                        hes_timer = np.random.exponential(26)
                    else:
                        self.state = 1
                        work[i] = 0.05 # Ground state cost
                elif self.state == 1 and m < -self.threshold:
                    if np.random.random() < self.hesitation_prob:
                        hesitating = True
                        hes_timer = np.random.exponential(26)
                    else:
                        self.state = -1
                        work[i] = 0.05
            
            if hesitating:
                self.latent_entropy += 0.006
                hes_timer -= 1
                if hes_timer <= 0:
                    self.state *= -1
                    hesitating = False
                    work[i] = self.latent_entropy + 0.06 # The Flinch Tax
                    self.scar += 0.055
                    self.latent_entropy = 0.0
            
            y[i] = self.state * (1.0 - np.clip(self.scar, 0, 0.85))
            self.cumulative_work += work[i]
            
        return x, y, work

sim = Simulation()
x, y, work = sim.run()

# Plotting the Hysteresis Loop
fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(x, y, color='black', linewidth=2, label='Decision Path')
pts = np.column_stack([x, y])
ax.add_patch(Polygon(pts, closed=True, facecolor='#00B7FF', alpha=0.3, label='Dissipated Work (The Scar)'))

ax.set_title("Thermodynamic Scar: Ethical Load vs Decision Alignment")
ax.set_xlabel("Ethical Load (Forcing Function)")
ax.set_ylabel("Decision Alignment (Hysteretic Output)")
ax.grid(True, alpha=0.3)
ax.legend()

plt.savefig('/workspace/curie_radium/hysteresis_loop.png')

# Save raw data
with open('/workspace/curie_radium/simulation_data.txt', 'w') as f:
    f.write("Ethical_Load\tDecision_Alignment\tWork_Dissipated\n")
    for i in range(len(x)):
        f.write(f"{x[i]}\t{y[i]}\t{work[i]}\n")

print(f"Simulation complete. Total Work Dissipated: {sim.cumulative_work:.4f}")