# bivology_mars_proto.py
# Bivology: Mars-inspired digital "organics" → life-like thresholds (prototype)
# Run: python bivology_mars_proto.py
# Optional: tweak PARAMS below to explore regimes

import random
import math
from collections import Counter, defaultdict
from dataclasses import dataclass, field

# =========================
# ====== PARAMETERS =======
# =========================
PARAMS = {
    "GRID_W": 80,
    "GRID_H": 48,
    "SEED": 42, # set to None for random
    "STEPS": 5000, # total ticks to simulate
    "INIT_MONOMERS": 3500, # starting 0/1 units
    "DIFFUSION_P": 0.60, # chance a particle moves each tick
    "POLYMERIZE_P": 0.10, # monomers join a chain when adjacent
    "CHAIN_BREAK_P": 0.001, # chains can break (abiotic decay)
    "TEMPLATE_REPL_P": 0.05, # chance a chain templates a neighbor
    "MUTATION_P": 0.01, # bit flip during replication
    "MIN_REPL_LEN": 6, # minimum length for templating
    "ENERGY_MAX": 20, # cell energy capacity
    "ENERGY_DIFFUSE": 0.25, # energy diffusion per tick
    "ENERGY_RECHARGE": 0.08, # recharge from global “geochemistry”
    "ENERGY_HARVEST_PER_SITE": 2, # what a chain can take per step (if available)
    "CHAIN_ENERGY_USE": 1, # upkeep energy per step per chain
    "CHAIN_LEN_METABOLIC": 5, # length at/above which metabolism is active
    "FRESH_SEED_RATE": 6, # monomers injected per tick (meteoritic input)
    "REPORT_EVERY": 200, # print metrics every N ticks
}

# =========================
# ====== DATA MODELS ======
# =========================
@dataclass
class Chain:
    bits: list # list of 0/1
    x: int
    y: int
    energy: int = 0

    def signature(self) -> str:
        # canonical text identity of a chain
        return "".join("1" if b else "0" for b in self.bits)

    def length(self) -> int:
        return len(self.bits)

@dataclass
class Cell:
    monomers: list = field(default_factory=list) # list of 0/1
    chains: list = field(default_factory=list) # list[Chain indices] (indices refer to global pool)
    energy: float = 0.0

# =========================
# ====== SIMULATION =======
# =========================
class World:
    def __init__(self, p):
        if p["SEED"] is not None:
            random.seed(p["SEED"])
        self.p = p
        self.grid = [[Cell() for _ in range(p["GRID_W"])] for _ in range(p["GRID_H"])]
        self.chains = [] # global chain store
        self._init_energy_field()
        self._seed_monomers()

    def _init_energy_field(self):
        # “Jezero-like” gradient: energy peaks along a meandering band (river analog)
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                # sinusoidal band across x with some y variance
                band = 0.5 + 0.5*math.sin((x/8.0) + math.sin(y/7.0))
                base = band * self.p["ENERGY_MAX"]
                self.grid[y][x].energy = min(self.p["ENERGY_MAX"], max(0.0, base * 0.6))

    def _seed_monomers(self):
        for _ in range(self.p["INIT_MONOMERS"]):
            x = random.randrange(self.p["GRID_W"])
            y = random.randrange(self.p["GRID_H"])
            self.grid[y][x].monomers.append(random.choice([0,1]))

    def step(self, t):
        self._inject_monomers()
        self._energy_recharge_and_diffuse()
        self._diffuse_particles()
        self._abiotic_polymerization()
        self._chains_metabolism_and_decay()
        self._templated_replication()

    # ----- processes -----
    def _inject_monomers(self):
        for _ in range(self.p["FRESH_SEED_RATE"]):
            x = random.randrange(self.p["GRID_W"])
            y = random.randrange(self.p["GRID_H"])
            self.grid[y][x].monomers.append(random.choice([0,1]))

    def _energy_recharge_and_diffuse(self):
        # recharge: global “geochemistry” adds a little to every cell
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                c = self.grid[y][x]
                c.energy = min(self.p["ENERGY_MAX"], c.energy + self.p["ENERGY_RECHARGE"])

        # simple diffusion: each cell shares a fraction with neighbors
        out = [[0.0 for _ in range(self.p["GRID_W"])] for _ in range(self.p["GRID_H"])]
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                share = self.grid[y][x].energy * self.p["ENERGY_DIFFUSE"]
                if share <= 0: 
                    continue
                nbrs = self._neighbors4(x,y)
                if not nbrs: 
                    continue
                give = share / len(nbrs)
                self.grid[y][x].energy -= share
                for nx, ny in nbrs:
                    out[ny][nx] += give
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                self.grid[y][x].energy = min(self.p["ENERGY_MAX"], self.grid[y][x].energy + out[y][x])

    def _diffuse_particles(self):
        # move monomers randomly; move chains as units
        moves = [(-1,0),(1,0),(0,-1),(0,1)]
        # monomers
        new_mons = [[[] for _ in range(self.p["GRID_W"])] for _ in range(self.p["GRID_H"])]
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                for m in self.grid[y][x].monomers:
                    if random.random() < self.p["DIFFUSION_P"]:
                        dx,dy = random.choice(moves)
                        nx,ny = (x+dx) % self.p["GRID_W"], (y+dy) % self.p["GRID_H"]
                        new_mons[ny][nx].append(m)
                    else:
                        new_mons[y][x].append(m)
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                self.grid[y][x].monomers = new_mons[y][x]

        # chains: move as a single object with small probability (slower)
        for idx, ch in enumerate(self.chains):
            if ch is None: 
                continue
            if random.random() < (self.p["DIFFUSION_P"] * 0.25):
                dx,dy = random.choice(moves)
                ch.x = (ch.x + dx) % self.p["GRID_W"]
                ch.y = (ch.y + dy) % self.p["GRID_H"]

    def _abiotic_polymerization(self):
        # if a cell has adjacent monomers, they can join to existing chains
        # or start a new chain abiotically
        for y in range(self.p["GRID_H"]):
            for x in range(self.p["GRID_W"]):
                cell = self.grid[y][x]
                if len(cell.monomers) < 2: 
                    continue
                if random.random() > self.p["POLYMERIZE_P"]:
                    continue

                # try to extend an existing chain in this cell (if any)
                chain_indices = [i for i,c in enumerate(self.chains) if (c is not None and c.x==x and c.y==y)]
                if chain_indices:
                    ci = random.choice(chain_indices)
                    ch = self.chains[ci]
                    # consume one monomer to extend, pick randomly 0/1 from cell
                    bit = cell.monomers.pop(random.randrange(len(cell.monomers)))
                    # append (simple rule; could also insert at start randomly)
                    ch.bits.append(bit)
                else:
                    # start a new chain from two monomers
                    m1 = cell.monomers.pop(random.randrange(len(cell.monomers)))
                    m2 = cell.monomers.pop(random.randrange(len(cell.monomers)))
                    new = Chain(bits=[m1,m2], x=x, y=y, energy=0)
                    self.chains.append(new)

    def _chains_metabolism_and_decay(self):
        # metabolism: chains of sufficient length harvest energy; all chains pay upkeep; chains can break
        for i,ch in enumerate(self.chains):
            if ch is None: 
                continue
            # upkeep
            ch.energy -= self.p["CHAIN_ENERGY_USE"]
            if ch.energy < -5:
                # chain disintegrates into monomers (death)
                cell = self.grid[ch.y][ch.x]
                cell.monomers.extend(ch.bits)
                self.chains[i] = None
                continue

            # harvest if long enough
            if ch.length() >= self.p["CHAIN_LEN_METABOLIC"]:
                cell = self.grid[ch.y][ch.x]
                take = min(self.p["ENERGY_HARVEST_PER_SITE"], int(cell.energy))
                if take > 0:
                    cell.energy -= take
                    ch.energy += take

            # abiotic decay: break at random site into two shorter chains (or drop monomers)
            if random.random() < self.p["CHAIN_BREAK_P"] and ch.length() >= 3:
                cut = random.randrange(1, ch.length()-1)
                left_bits = ch.bits[:cut]
                right_bits = ch.bits[cut:]
                # replace original:
                self.chains[i] = Chain(bits=left_bits, x=ch.x, y=ch.y, energy=max(0, ch.energy//2))
                # second piece becomes new chain (or monomers if too short)
                if len(right_bits) >= 2:
                    self.chains.append(Chain(bits=right_bits, x=ch.x, y=ch.y, energy=max(0, ch.energy//2)))
                else:
                    self.grid[ch.y][ch.x].monomers.extend(right_bits)

    def _templated_replication(self):
        # If two chains are in same cell and one is long enough, it can template the other (copy pattern)
        # Copy consumes local monomers matching needed bits; mutation may flip bits
        cell_chain_map = defaultdict(list)
        for idx,ch in enumerate(self.chains):
            if ch is None: 
                continue
            cell_chain_map[(ch.x,ch.y)].append(idx)

        for (x,y), idxs in cell_chain_map.items():
            if len(idxs) < 1:
                continue
            # pick a chain that might act as a template
            t_idx = random.choice(idxs)
            templ = self.chains[t_idx]
            if templ is None or templ.length() < self.p["MIN_REPL_LEN"]:
                continue
            if random.random() > self.p["TEMPLATE_REPL_P"]:
                continue

            cell = self.grid[y][x]
            # try to build a copy from available monomers in the cell
            copy_bits = []
            temp_bits = templ.bits
            # naive “need one monomer per bit” model
            local = cell.monomers
            if len(local) < len(temp_bits):
                continue

            # assemble a candidate using available monomers
            # we don't require matching monomer type here (abiotic looseness); mutation handles differences
            for k in range(len(temp_bits)):
                bit = local.pop(random.randrange(len(local)))
                # mutate with small probability
                if random.random() < self.p["MUTATION_P"]:
                    bit = 1 - bit
                copy_bits.append(bit)

            # birth: new chain starts with low energy
            newborn = Chain(bits=copy_bits, x=x, y=y, energy=0)
            self.chains.append(newborn)

    # ----- utilities -----
    def _neighbors4(self, x, y):
        w,h = self.p["GRID_W"], self.p["GRID_H"]
        return [((x-1)%w,y), ((x+1)%w,y), (x,(y-1)%h), (x,(y+1)%h)]

# =========================
# ====== METRICS ==========
# =========================
def life_detection_scale(stats):
    """
    LDS (0-6):
      0: Only random monomers, no chains
      1: Chains exist but are short and fleeting (median len < 4, half-life short)
      2: Stable chains (median len ≥ 4) persist across reports
      3: Metabolic coupling (≥20% chains length ≥ metabolic threshold)
      4: Replication signatures (growth in most-common signatures, newborn rate > 0)
      5: Selection/adaptation (one/few signatures dominate over time with rising median length)
      6: Ecosystem-like (diverse signatures co-exist, trophic-like: evidence of long-lived variety)
    """
    # Heuristics from current snapshot
    chains = stats["chain_count"]
    median_len = stats["median_len"]
    frac_metabolic = stats["frac_metabolic"]
    newborns = stats["newborns_recent"]
    sig_dom_frac = stats["sig_dom_frac"]
    unique_sigs = stats["unique_sigs"]

    if chains == 0:
        return 0
    if median_len < 4:
        return 1
    if median_len >= 4 and frac_metabolic < 0.2:
        return 2
    if frac_metabolic >= 0.2 and newborns == 0:
        return 3
    if newborns > 0 and sig_dom_frac < 0.6:
        return 4
    if sig_dom_frac >= 0.6 and unique_sigs <= 3:
        return 5
    return 6

def summarize(world, history_births):
    lengths = [ch.length() for ch in world.chains if ch is not None]
    chain_count = len(lengths)
    median_len = 0 if chain_count == 0 else sorted(lengths)[chain_count//2]
    metabolic = [L for L in lengths if L >= world.p["CHAIN_LEN_METABOLIC"]]
    frac_metabolic = 0.0 if chain_count == 0 else len(metabolic)/chain_count

    sigs = Counter(ch.signature() for ch in world.chains if ch is not None)
    unique_sigs = len(sigs)
    sig_dom_frac = 0.0 if chain_count == 0 else (sigs.most_common(1)[0][1] / chain_count)

    # newborns in last interval:
    newborns_recent = history_births[-1] if history_births else 0

    return {
        "chain_count": chain_count,
        "median_len": median_len,
        "frac_metabolic": frac_metabolic,
        "unique_sigs": unique_sigs,
        "sig_dom_frac": sig_dom_frac,
        "newborns_recent": newborns_recent
    }

# =========================
# ====== RUN LOOP =========
# =========================
def main():
    p = PARAMS
    w = World(p)
    births_tracker = []
    prev_chain_total = 0

    for t in range(1, p["STEPS"]+1):
        before = sum(1 for c in w.chains if c is not None)
        w.step(t)
        after = sum(1 for c in w.chains if c is not None)
        births = max(0, after - before)
        births_tracker.append(births)

        if t % p["REPORT_EVERY"] == 0:
            stats = summarize(w, [sum(births_tracker[-p["REPORT_EVERY"] :])])
            lds = life_detection_scale(stats)
            print(
                f"[t={t}] chains={stats['chain_count']:5d} "
                f"median_len={stats['median_len']:2d} "
                f"metab%={stats['frac_metabolic']*100:4.1f} "
                f"unique={stats['unique_sigs']:4d} "
                f"dom_frac={stats['sig_dom_frac']:.2f} "
                f"newborns={stats['newborns_recent']:4d} "
                f"LDS={lds}"
            )

if __name__ == "__main__":
    main()