#!/usr/bin/env python3
"""Minimal 2D Ising Monte Carlo model for classroom demonstration."""

from __future__ import annotations

import argparse
import numpy as np


def initialize(n: int, rng: np.random.Generator) -> np.ndarray:
    return rng.choice(np.array([-1, 1], dtype=np.int8), size=(n, n))


def metropolis_sweep(spins: np.ndarray, temperature: float, field: float, rng: np.random.Generator) -> None:
    n = spins.shape[0]
    for _ in range(n * n):
        i = rng.integers(0, n)
        j = rng.integers(0, n)
        s = spins[i, j]
        neighbors = (
            spins[(i + 1) % n, j]
            + spins[(i - 1) % n, j]
            + spins[i, (j + 1) % n]
            + spins[i, (j - 1) % n]
        )
        delta_e = 2.0 * s * (neighbors + field)
        if delta_e <= 0.0 or rng.random() < np.exp(-delta_e / temperature):
            spins[i, j] = -s


def magnetization(spins: np.ndarray) -> float:
    return float(spins.mean())


def main() -> None:
    parser = argparse.ArgumentParser()
    parser.add_argument("--n", type=int, default=48)
    parser.add_argument("--temperature", type=float, default=1.8)
    parser.add_argument("--field", type=float, default=0.0)
    parser.add_argument("--sweeps", type=int, default=200)
    parser.add_argument("--seed", type=int, default=2)
    args = parser.parse_args()

    rng = np.random.default_rng(args.seed)
    spins = initialize(args.n, rng)
    for step in range(args.sweeps):
      metropolis_sweep(spins, args.temperature, args.field, rng)
      if step % max(1, args.sweeps // 10) == 0:
          print(f"sweep {step:4d}: M = {magnetization(spins): .3f}")
    print(f"final: M = {magnetization(spins): .3f}")


if __name__ == "__main__":
    main()