================================================================================
                         PBITLANG LANGUAGE SPECIFICATION
                                  Version 1.0
================================================================================

1. INTRODUCTION
================================================================================

PbitLang is a domain-specific language for expressing Hamiltonians and energy
functions for thermodynamic computing systems. It compiles to PHAL Hamiltonian
objects for execution on P-bit hardware or simulators.

Design Goals:
  - Scientific rigor: Enforce physical correctness at compile time
  - Expressiveness: Support common physics models naturally
  - Safety: Prevent physically invalid configurations
  - Clarity: Make energy functions readable and maintainable


2. LEXICAL STRUCTURE
================================================================================

2.1 Character Set
-----------------
PbitLang source files are UTF-8 encoded. Identifiers use ASCII letters,
digits, and underscores.

2.2 Keywords
------------
Reserved words (case-sensitive):

  Types:        spin, coupling, field, lattice, graph, real, int, bool
  Spin types:   ising, binary, potts, clock, continuous
  Declarations: hamiltonian, system, param, const, let, import, from
  Constructs:   sum, product, over, where, for, in, if, else
  Lattice:      chain, square, triangular, honeycomb, kagome, cubic
  Boundary:     periodic, open, fixed
  Operators:    and, or, not, mod
  Physics:      neighbors, next_neighbors, all_pairs, distance
  Constraints:  require, assert, symmetric, positive, normalized

2.3 Operators
-------------
  Arithmetic:   + - * / ^ %
  Comparison:   == != < > <= >=
  Assignment:   = :=
  Range:        .. ...
  Subscript:    [ ]
  Access:       .
  Type:         :
  Annotation:   @

2.4 Literals
------------
  Integers:     42, -17, 0x2A, 0b101010
  Reals:        3.14, 2.5e-3, -1.0
  Booleans:     true, false
  Strings:      "hello", 'world'

2.5 Comments
------------
  Line:         // comment
                # comment
  Block:        /* multi-line
                   comment */
  Doc:          /// documentation comment


3. TYPE SYSTEM
================================================================================

3.1 Primitive Types
-------------------

  int           Signed integer
  real          Double-precision floating point
  bool          Boolean (true/false)
  complex       Complex number (for future quantum extension)

3.2 Spin Types
--------------

  ising         Classical Ising spin: s ∈ {-1, +1}
  binary        Binary variable: x ∈ {0, 1}
  potts(q)      q-state Potts: σ ∈ {0, 1, ..., q-1}
  clock(q)      q-state clock: θ = 2πk/q, k ∈ {0, ..., q-1}
  continuous    Continuous angle: θ ∈ [0, 2π)

Spin types encode their domain and valid operations:
  - ising supports multiplication s_i * s_j
  - binary supports x_i * x_j and (1 - x_i)
  - potts supports delta(σ_i, σ_j) (Kronecker delta)
  - clock supports cos(θ_i - θ_j)
  - continuous supports cos(θ_i - θ_j), sin(θ_i - θ_j)

3.3 Aggregate Types
-------------------

  array[T, n]       Fixed-size array of n elements of type T
  matrix[T, n, m]   n×m matrix of type T
  tensor[T, dims]   Multi-dimensional tensor

  lattice[S, L]     Spins of type S arranged on lattice L
  graph[S, G]       Spins of type S on graph G

3.4 Coupling Types
------------------

  coupling[S]       Coupling appropriate for spin type S
                    - ising/binary: real scalar (J_ij)
                    - potts: q×q matrix
                    - continuous: real (for cos interaction)

3.5 Physical Units (Optional)
-----------------------------

  Annotate values with physical units for documentation:
    @unit("meV")
    @unit("kelvin")
    @unit("tesla")

  Units are checked for consistency but not converted.


4. DECLARATIONS
================================================================================

4.1 Parameters
--------------

Parameters are values provided at compile time:

  param n: int                    // Required parameter
  param J: real = 1.0             // Default value
  param q: int where q >= 2       // Constrained parameter

4.2 Constants
-------------

Constants are computed at compile time:

  const pi: real = 3.14159265359
  const k_B: real = 8.617333262e-5  // eV/K

4.3 Lattice Declarations
------------------------

  lattice L = chain(n)                          // 1D chain
  lattice L = chain(n, boundary=periodic)       // Periodic boundary
  lattice L = square(nx, ny)                    // 2D square lattice
  lattice L = triangular(n)                     // Triangular lattice
  lattice L = honeycomb(n)                      // Honeycomb (graphene)
  lattice L = kagome(n)                         // Kagome lattice
  lattice L = cubic(nx, ny, nz)                 // 3D cubic
  lattice L = from_adjacency(A)                 // Custom from matrix
  lattice L = from_edges(edges)                 // Custom from edge list

4.4 Spin System Declarations
----------------------------

  system S: ising on chain(n)                   // Ising spins on chain
  system S: potts(3) on triangular(n)           // 3-state Potts
  system S: binary on graph(G)                  // Binary on custom graph


5. HAMILTONIAN DEFINITIONS
================================================================================

5.1 Basic Syntax
----------------

  hamiltonian Name(params) -> spin_type on lattice {
      // Energy terms
  }

5.2 Energy Terms
----------------

Sum over interactions:

  // Nearest-neighbor Ising
  sum((i,j) in neighbors(L)) { -J * s[i] * s[j] }

  // All pairs with distance weighting
  sum((i,j) in all_pairs(L) where i < j) { -J/distance(i,j)^2 * s[i]*s[j] }

  // External field
  sum(i in sites(L)) { -h[i] * s[i] }

5.3 Built-in Interaction Functions
----------------------------------

For lattices:
  neighbors(L)         Nearest-neighbor pairs
  next_neighbors(L)    Next-nearest-neighbor pairs
  all_pairs(L)         All pairs (i,j) where i < j
  sites(L)             All site indices
  distance(i, j)       Euclidean distance on lattice
  coordination(i)      Number of neighbors of site i

For spins:
  delta(σ_i, σ_j)      Kronecker delta (Potts)
  cos(θ_i - θ_j)       Cosine coupling (XY/clock)
  dot(S_i, S_j)        Dot product (Heisenberg)

5.4 Composition
---------------

Hamiltonians can be composed:

  hamiltonian Combined = H1 + H2           // Sum
  hamiltonian Scaled = 2.0 * H1            // Scalar multiplication
  hamiltonian Subsystem = H1 on subset(L)  // Restriction


6. STANDARD LIBRARY MODELS
================================================================================

6.1 Ising Models
----------------

  // Ferromagnetic Ising on chain
  @model
  hamiltonian IsingChain(n: int, J: real = 1.0, h: real = 0.0)
      -> ising on chain(n, periodic) {
      coupling: sum((i,j) in neighbors) { -J * s[i] * s[j] }
      field:    sum(i in sites) { -h * s[i] }
  }

  // 2D Ising on square lattice
  @model
  hamiltonian Ising2D(nx: int, ny: int, J: real = 1.0, h: real = 0.0)
      -> ising on square(nx, ny, periodic) {
      coupling: sum((i,j) in neighbors) { -J * s[i] * s[j] }
      field:    sum(i in sites) { -h * s[i] }
  }

  // Sherrington-Kirkpatrick spin glass
  @model
  hamiltonian SK(n: int, J_matrix: matrix[real, n, n])
      -> ising on complete(n)
      where symmetric(J_matrix) {
      coupling: sum((i,j) in all_pairs where i < j) { -J_matrix[i,j] * s[i]*s[j] }
  }

6.2 Potts Model
---------------

  @model
  hamiltonian Potts(n: int, q: int, J: real = 1.0)
      -> potts(q) on square(n, n, periodic)
      where q >= 2 {
      coupling: sum((i,j) in neighbors) { -J * delta(σ[i], σ[j]) }
  }

6.3 XY Model
------------

  @model
  hamiltonian XY(n: int, J: real = 1.0)
      -> clock(256) on square(n, n, periodic) {  // Discretized XY
      coupling: sum((i,j) in neighbors) { -J * cos(θ[i] - θ[j]) }
  }

6.4 Hopfield Network
--------------------

  @model
  hamiltonian Hopfield(patterns: matrix[ising, p, n])
      -> ising on complete(n) {
      // Hebbian weights: J_ij = (1/n) * sum_μ ξ_i^μ * ξ_j^μ
      let J[i,j] = (1.0/n) * sum(μ in 0..p) { patterns[μ,i] * patterns[μ,j] }
      coupling: sum((i,j) in all_pairs where i < j) { -J[i,j] * s[i] * s[j] }
  }

6.5 QUBO
--------

  @model
  hamiltonian QUBO(Q: matrix[real, n, n])
      -> binary on complete(n) {
      // E(x) = sum_ij Q_ij x_i x_j
      energy: sum((i,j) in all_pairs) { Q[i,j] * x[i] * x[j] }
             + sum(i in sites) { Q[i,i] * x[i] }
  }

6.6 MAX-SAT
-----------

  @model
  hamiltonian MaxSAT(clauses: list[clause])
      -> binary on variables(clauses) {
      // Each unsatisfied clause contributes +1 to energy
      energy: sum(c in clauses) { unsatisfied(c, x) }
  }

6.7 Restricted Boltzmann Machine
--------------------------------

  @model
  hamiltonian RBM(n_visible: int, n_hidden: int,
                  W: matrix[real, n_visible, n_hidden],
                  a: array[real, n_visible],
                  b: array[real, n_hidden])
      -> binary on bipartite(n_visible, n_hidden) {
      visible_bias: sum(i in visible) { -a[i] * v[i] }
      hidden_bias:  sum(j in hidden) { -b[j] * h[j] }
      interaction:  sum((i,j) in edges) { -W[i,j] * v[i] * h[j] }
  }


7. PHYSICS VALIDATION
================================================================================

7.1 Compile-Time Checks
-----------------------

The compiler enforces:

  1. Symmetry: Coupling matrices must be symmetric for undirected models
     - J[i,j] == J[j,i] for Ising
     - Error if asymmetric matrix provided without @directed annotation

  2. Diagonal: Self-coupling J[i,i] typically zero (flagged as warning)

  3. Bounds: Constrained parameters checked against declared bounds
     - param q: int where q >= 2  // Enforced

  4. Type Consistency:
     - Ising spins: multiplication produces real
     - Potts spins: must use delta function
     - Binary spins: quadratic terms only

  5. Lattice Compatibility:
     - neighbors() requires lattice structure
     - complete() graph for all-to-all

7.2 Runtime Checks
------------------

Generated code includes optional runtime validation:

  @validate(strict)   // All checks, slower
  @validate(fast)     // Minimal checks
  @validate(none)     // No runtime checks

Checks include:
  - Matrix dimension matching
  - Value range validation
  - NaN/Inf detection

7.3 Physical Warnings
---------------------

Compiler emits warnings for:

  - Frustrated systems (antiferro on triangular)
  - Near-critical temperature parameters
  - Large condition number in coupling matrix
  - Potential numerical instability


8. COMPILATION
================================================================================

8.1 Compilation Targets
-----------------------

  pbitlang compile model.pbit --target=phal    # PHAL Hamiltonian (default)
  pbitlang compile model.pbit --target=python  # Python code
  pbitlang compile model.pbit --target=json    # Serialized form
  pbitlang compile model.pbit --target=latex   # LaTeX equations

8.2 Compilation Process
-----------------------

  1. Lexical Analysis    Source → Tokens
  2. Parsing             Tokens → AST
  3. Type Checking       AST → Typed AST
  4. Physics Validation  Typed AST → Validated AST
  5. Optimization        Simplification, constant folding
  6. Code Generation     Validated AST → Target

8.3 Optimization Passes
-----------------------

  - Constant folding
  - Common subexpression elimination
  - Sparse matrix detection
  - Symmetry exploitation


9. EXAMPLES
================================================================================

9.1 Simple Ising Chain
----------------------

  // File: ising_chain.pbit

  hamiltonian IsingChain(n: int, J: real, h: real)
      -> ising on chain(n, periodic) {

      // Nearest-neighbor coupling
      coupling: sum((i,j) in neighbors) {
          -J * s[i] * s[j]
      }

      // External magnetic field
      field: sum(i in sites) {
          -h * s[i]
      }
  }

9.2 Spin Glass
--------------

  // File: spin_glass.pbit

  import random from stdlib

  hamiltonian EdwardsAnderson(n: int, disorder: real)
      -> ising on cubic(n, n, n, periodic) {

      // Random ±J couplings
      @quenched
      let J_bond = random.choice([-1.0, 1.0])

      coupling: sum((i,j) in neighbors) {
          -J_bond * s[i] * s[j]
      }
  }

9.3 Max-Cut Problem
-------------------

  // File: maxcut.pbit

  hamiltonian MaxCut(edges: list[(int, int)], weights: list[real])
      -> ising on graph(edges) {

      // Minimize: -sum_{(i,j) in E} w_ij * (1 - s_i * s_j) / 2
      // Equivalent to maximize cut weight

      energy: sum(k in 0..len(edges)) {
          let (i, j) = edges[k]
          -weights[k] * (1 - s[i] * s[j]) / 2
      }
  }


10. ERROR MESSAGES
================================================================================

PbitLang provides clear, educational error messages:

  Error [E0042]: Asymmetric coupling matrix
    --> model.pbit:15:5
     |
  15 |     let J = [[0, 1, 2], [1, 0, 1], [3, 1, 0]]
     |             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     |
     = note: Ising coupling matrices must be symmetric: J[i,j] = J[j,i]
     = note: Found J[0,2] = 2 but J[2,0] = 3
     = help: For directed interactions, use @directed annotation

  Warning [W0017]: Frustrated system detected
    --> model.pbit:8:1
     |
   8 | hamiltonian AFTriangular(n: int, J: real)
     | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
     |
     = note: Antiferromagnetic coupling on triangular lattice is geometrically
             frustrated. Ground state is highly degenerate.
     = help: This may cause slow thermalization. Consider using parallel
             tempering for sampling.


================================================================================
                              END OF SPECIFICATION
================================================================================
