scientific-domains

The scientific-domains worked example.

Run it from sema/:

sema check examples/scientific-domains
SEMA_STRICT=1 sema run examples/scientific-domains
sema assure examples/scientific-domains --grade silver

Source

src/main.sema

"""Public scalar scientific domains with explicit, checked semantics.

These values are native Sema domains rather than opaque Python handles. Their
operations therefore have the same typed failures, tagged interchange, and
tree-walker/VM behavior as the rest of the language.
"""

import math
import latex

assure silver

equation kinetic_energy(m: any, v: any) -> any:
    return 1 / 2 * m * v^2

equation bounded_number_theory() -> any:
    return (factorint(360), next_prime(100), prev_prime(100), divisor_count(360))

equation bounded_prime_sequence() -> any:
    return (prime_nth(25), prime(1), prime_count(100), primepi(541))

equation bounded_category_type() -> any:
    category := [2, [0, 1, 0], [0, 1, 1], [0, 1], [[0, -1, -1], [-1, 1, 2], [2, -1, -1]]]
    return (category_validate(category), category_compose(category, 1, 2), category_hom(category, 0, 1), category_isomorphisms(category), functor_validate(category, category, [0, 1], [0, 1, 2]), natural_transformation_validate(category, category, [0, 1], [0, 1, 2], [0, 1], [0, 1, 2], [0, 1]), stlc_type(["lam", "Nat", ["var", 0]], []), stlc_normalize(["app", ["lam", "Nat", ["var", 0]], ["nat", 7]], 1), stlc_definitional_equal(["app", ["lam", "Nat", ["var", 0]], ["nat", 7]], ["nat", 7], 1), category_validate([1, [0], [0], [0], [[-1]]]))

equation bounded_morphism_composition() -> any:
    # `∘` composes functions right-to-left: (f ∘ g)(x) = f(g(x)), so the right
    # operand runs first. It builds a value rather than applying anything, so a
    # composite can be named and reused like any other function.
    double(x) := 2 * x
    shift(x) := x + 3
    forward := double ∘ shift
    backward := shift ∘ double
    # `∘` is left-associative and composition is associative, so both groupings
    # of one chain agree. Either operand may be a native unary kernel: `sqrt`
    # composes exactly like a user-defined equation.
    left_grouped := (sqrt ∘ double) ∘ shift
    right_grouped := sqrt ∘ (double ∘ shift)
    return (forward(5), backward(5), left_grouped(5), right_grouped(5))

equation bounded_constraints() -> any:
    domains := [[0, 1, 2], [0, 1, 2]]
    relations := [[[0, 1], [[0, 2], [1, 1], [2, 0]]]]
    return (csp(domains, relations, 8), solve_constraints(domains, relations, 8), sat(domains, relations), sat([[0, 1]], [[[0], []]]))

def bounded_units():
    energy = quantity(2.0, "N") * quantity(3.0, "m")
    bounds = quantity_uncertainty(quantity(10.0, "m"), quantity(0.1, "m"))
    return (dimension(energy), dimension_check(energy, "J"), equivalent_unit("J", "N*m"), simplify_unit(physical_constant("c")), unit_prefix("k"), convert_unit(physical_constant("c"), "m*s^-1"), convert_unit(bounds[0], "m"), convert_unit(bounds[1], "m"), len(base_units()), len(derived_units()))

equation bounded_interpolation() -> any:
    return (interpolate([0, 1, 2], [1, 3, 7], 1 / 2), polynomial_interpolate([0, 1, 2], [1, 3, 7]), polynomial_interpolate([0, 1, 2], [1, 2, 3]))

equation bounded_interpolation_real(xs: any, ys: any, x: any) -> any:
    return (interpolate(xs, ys, x), polynomial_interpolate(xs, ys))

equation bounded_statistics() -> any:
    return (mean([1, 2, 3]), variance([1, 2, 3]), correlation([1, 2, 3], [2, 4, 6]), entropy([0.25, 0.75]), js_divergence([0.25, 0.75], [0.5, 0.5]))

equation bounded_convolution() -> any:
    return convolve([1, 2, 3], [4, 5])

equation bounded_cross_correlation() -> any:
    return correlate([1, 2, 3], [4, 5])

equation bounded_normal() -> any:
    return (normal_pdf(0, 0, 1), normal_logpdf(0, 0, 1), normal_cdf(0, 0, 1), normal_sf(0, 0, 1), normal_logcdf(0, 0, 1), normal_logsf(0, 0, 1))

equation bounded_normal_tails() -> any:
    return (normal_cdf(9, 0, 1), normal_sf(9, 0, 1), normal_logcdf(-40, 0, 1), normal_logsf(40, 0, 1))

equation bounded_normal_quantiles() -> any:
    return (normal_ppf(0.5, 0, 1), normal_ppf(0.975, 1, 2), normal_logppf(-0.6931471805599453, 0, 1), normal_logppf(-800, 0, 1))

equation bounded_rounding(values: any, matrix: any) -> any:
    return (floor(7 / 2), ceil(0 - 7 / 2), round([3 / 2, 5 / 2, 0 - 5 / 2]), trunc((0 - 7 / 2, 7 / 2)), fract([7 / 2, 0 - 7 / 2]), floor(values), round(matrix), floor("x"), subst(fract("x"), "x", 0 - 7 / 2))

equation bounded_svd() -> any:
    return (svd([[3.0, 0.0], [0.0, 2.0], [0.0, 0.0]]), rank([[1e-200, 0.0], [0.0, 5e-201]]))

def bounded_sparse():
    equation:
        s := sparse(2, 2, [1, 0], [1, 0], [4.0, 2.0])
        applied := sparse.matmul(s, [1.0, 1.0])
        solved := sparse.solve(s, [2.0, 8.0])
        transposed := sparse.transpose(s)
        product := sparse.matmul(s, transposed)
        factors := sparse.lu(s)
        fitted := sparse.lstsq(s, [2.0, 8.0])
        eigen := sparse.eigs(s)
    return (s, applied, solved, transposed, product, factors, fitted, eigen)

def bounded_complex_linalg():
    a = tensor([[complex(2.0, 1.0), complex(0.0, 0.0)], [complex(0.0, 0.0), complex(1.0, -1.0)]])
    x = tensor([complex(1.0, 0.0), complex(0.0, 1.0)])
    b = tensor([complex(2.0, 2.0), complex(4.0, -4.0)])
    real = tensor([[2.0, 0.0], [0.0, 4.0]])
    equation:
        product := matvec(a, x)
        solved := solve(a, product)
        square := matmul(a, a)
        promoted := solve(real, b)
    return (product, solved, square, promoted, dtype(square))

test "complex arithmetic and elementary dispatch":
    z = complex(3.0, 4.0)
    check abs(z) == 5.0
    check (z * z.conj).re == 25.0
    check math.exp(complex(0.0, math.pi)).re < -0.999999999999

test "certified intervals enclose every represented result":
    x = interval(1.0, 2.0)
    y = interval(3.0, 4.0)
    product = x * y
    check product.lo <= 3.0 and product.hi >= 8.0
    check 1.5 in x
    check interval(1.25, 1.75) in x

test "quaternion geometry preserves vector norm":
    quarter_turn = quaternion(0.7071067811865476, 0.0, 0.0, 0.7071067811865476)
    rotated = rotate(quarter_turn, [1.0, 0.0, 0.0])
    check abs(rotated[0]) < 0.000000000000001
    check abs(rotated[1] - 1.0) < 0.000000000000001
    check abs(rotated[2]) < 0.000000000000001

test "modular division and negative powers require an inverse":
    a = modint(3, 7)
    check (a ** -1).value == 5
    check (a / modint(2, 7)).value == 5
    check (a ** 6).value == 1

test "decimal arithmetic carries an explicit context":
    left = decimal("1.005", precision=3, rounding="half_even")
    right = decimal("0.005", precision=3, rounding="half_even")
    total = left + right
    check total == decimal("1.01", precision=3, rounding="half_even")
    check total.precision == 3
    check total.rounding == "half_even"

test "bounded prime navigation and divisor evidence are exact":
    evidence = bounded_number_theory()
    check evidence[0] == [(2, 3), (3, 2), (5, 1)]
    check evidence[1] == 101
    check evidence[2] == 97
    check evidence[3] == 24

test "bounded prime indexing and counting are exact and one based":
    evidence = bounded_prime_sequence()
    check evidence[0] == 97
    check evidence[1] == 2
    check evidence[2] == 25
    check evidence[3] == 100

test "bounded categories and STLC preserve exact laws and normal forms":
    evidence = bounded_category_type()
    check evidence[0] and evidence[4] and evidence[5]
    check evidence[1] == 2 and evidence[2] == [2]
    check evidence[3] == [[0, 0], [1, 1]]
    check evidence[6] == ["->", "Nat", "Nat"]
    check evidence[7] == ["nat", 7]
    check evidence[8] and not evidence[9]

test "function composition applies right to left and associates":
    evidence = bounded_morphism_composition()
    # shift(5) = 8 runs first, then double(8) = 16.
    check evidence[0] == 16
    # The flipped order is a different function: double(5) = 10, shift(10) = 13.
    check evidence[1] == 13
    # sqrt(double(shift(5))) = sqrt(16), under either grouping of the chain.
    check evidence[2] == 4.0
    check evidence[2] == evidence[3]

test "bounded pure table constraints enumerate exact models":
    evidence = bounded_constraints()
    check evidence[0] == [[0, 2], [1, 1], [2, 0]]
    check evidence[0] == evidence[1]
    check evidence[2] and not evidence[3]

test "bounded SI units expose dimensions conversion and explicit uncertainty":
    evidence = bounded_units()
    check evidence[0] == [2, 1, -2, 0, 0, 0, 0]
    check evidence[1] and evidence[2]
    check evidence[3] == "m*s^-1"
    check evidence[4] == 1000.0 and evidence[5] == 299792458.0
    check abs(evidence[6] - 9.9) < 0.000000000000001
    check abs(evidence[7] - 10.1) < 0.000000000000001
    check evidence[8] == 7 and evidence[9] == 9

test "bounded interpolation keeps exact lanes exact":
    evidence = bounded_interpolation()
    check evidence[0] == QQ(7, 4)
    check evidence[1] == [1, 1, 1]
    check evidence[2] == [1, 1, 0]

test "float interpolation inputs select the strict finite-real lane":
    evidence = bounded_interpolation_real([0.0, 1.0], [1.0, 3.0], 0.5)
    check evidence[0] == 2.0
    check evidence[1] == [1.0, 2.0]

test "bounded statistics declare population and information semantics":
    evidence = bounded_statistics()
    check evidence[0] == 2.0
    check abs(evidence[1] - 0.6666666666666666) < 0.000000000000001
    check abs(evidence[2] - 1.0) < 0.000000000000001
    check evidence[3] > 0.0
    check evidence[4] > 0.0

test "bounded convolution is full and rank one":
    output = bounded_convolution()
    check output[0] == 4.0
    check output[1] == 13.0
    check output[2] == 22.0
    check output[3] == 15.0

test "bounded cross-correlation uses ascending full-lag order":
    output = bounded_cross_correlation()
    check output[0] == 5.0
    check output[1] == 14.0
    check output[2] == 23.0
    check output[3] == 12.0

test "bounded Normal evaluation is scalar and tail-aware":
    output = bounded_normal()
    tails = bounded_normal_tails()
    quantiles = bounded_normal_quantiles()
    check abs(output[0] - 0.3989422804014327) < 0.000000000000001
    check abs(output[1] + 0.9189385332046727) < 0.000000000000001
    check output[2] == 0.5
    check output[3] == 0.5
    check abs(output[4] + 0.6931471805599453) < 0.000000000000001
    check abs(output[5] + 0.6931471805599453) < 0.000000000000001
    check tails[0] == 1.0
    check tails[1] > 0.0
    check tails[2] < -800.0
    check tails[3] < -800.0
    check abs(tails[2] - tails[3]) < 0.000000000000001
    check quantiles[0] == 0.0
    check abs(quantiles[1] - 4.919927969080108) < 0.000000000001
    check quantiles[2] == 0.0
    check abs(quantiles[3] + 39.88469483825668) < 0.000000000001

test "equation rounding is exact symbolic and shape preserving":
    output = bounded_rounding(tensor([1.9, -1.1]), tensor([[1.5, 2.5], [-1.5, -2.5]]))
    check output[0] == 3
    check output[1] == -3
    check output[2] == [2, 2, -2]
    check output[3] == (-3, 3)
    check output[4] == [QQ(1, 2), QQ(-1, 2)]
    check output[5][0] == 1.0 and output[5][1] == -2.0
    check output[6][0][0] == 2.0 and output[6][1][1] == -2.0
    check output[7] == "floor(x)"
    check output[8] == QQ(-1, 2)

test "reduced SVD is reconstructable and matrix rank is scale relative":
    output = bounded_svd()
    factors = output[0]
    check factors[0][0][0] == 1.0 and factors[0][1][1] == 1.0
    check factors[1][0] == 3.0 and factors[1][1] == 2.0
    check factors[2][0][0] == 1.0 and factors[2][1][1] == 1.0
    check output[1] == 2

test "bounded sparse CSR composes, factors, fits, and solves":
    evidence = bounded_sparse()
    s = evidence[0]
    check s.schema == "sema.sparse-matrix/v1"
    check s.format == "csr"
    check s.rows == 2 and s.cols == 2 and s.nnz == 2
    check s.indptr == [0, 1, 2]
    check s.indices == [0, 1]
    check s.values == [2.0, 4.0]
    check evidence[1][0] == 2.0 and evidence[1][1] == 4.0
    check evidence[2][0] == 1.0 and evidence[2][1] == 2.0
    check evidence[3].values == [2.0, 4.0]
    check evidence[4].values == [4.0, 16.0]
    check evidence[5][0].format == "csr" and evidence[5][1].format == "csr"
    check evidence[6][0] == 1.0 and evidence[6][1] == 2.0
    check evidence[7][4] and evidence[7][0] == 4.0

test "bounded complex dense products solve and promote exactly":
    evidence = bounded_complex_linalg()
    check dtype(evidence[0]) == "complex" and shape(evidence[0]) == [2]
    check sum(evidence[0]) == complex(3.0, 2.0)
    check sum(evidence[1]) == complex(1.0, 1.0)
    check sum(evidence[2]) == complex(3.0, 2.0)
    check sum(evidence[3]) == complex(2.0, 0.0)
    check evidence[4] == "complex"

test "latex renders formulas and serializes exact and symbolic values":
    block = latex.render("\\frac{1}{2}")
    check len(block) > 0
    check block == " 1 \n───\n 2 \n"
    check latex.of(QQ(1, 2)) == "\\frac{1}{2}"
    check latex.of(kinetic_energy) == "\\frac{1}{2} m v^{2}"
    pretty = latex.render(latex.of(kinetic_energy))
    check len(pretty) > 0
    check pretty == " 1    \n───mv²\n 2    \n"

def main() !{observe.record}:
    z = math.sqrt(complex(-4.0, 0.0))
    bounds = math.exp(interval(0.0, 1.0))
    orientation = quaternion(1.0, 0.0, 0.0, 0.0)
    residue = modint(17, 5)
    amount = decimal("12.345", precision=4, rounding="half_up")
    sparse_evidence = bounded_sparse()
    complex_evidence = bounded_complex_linalg()
    report = {
        "complex": z,
        "interval": bounds,
        "quaternion": orientation,
        "modint": residue,
        "decimal": amount,
        "number_theory": bounded_number_theory(),
        "prime_sequence": bounded_prime_sequence(),
        "category_type": bounded_category_type(),
        "constraints": bounded_constraints(),
        "units": bounded_units(),
        "interpolation": bounded_interpolation(),
        "statistics": bounded_statistics(),
        "convolution": bounded_convolution(),
        "cross_correlation": bounded_cross_correlation(),
        "normal": bounded_normal(),
        "rounding": bounded_rounding(tensor([1.9, -1.1]), tensor([[1.5, 2.5], [-1.5, -2.5]])),
        "svd": bounded_svd(),
        "sparse": sparse_evidence[0],
        "complex_linalg": (sum(complex_evidence[1]), complex_evidence[4]),
    }
    log.info("scientific domains", report=report)
    return report

Reflected API

main

Public scalar scientific domains with explicit, checked semantics.

These values are native Sema domains rather than opaque Python handles. Their operations therefore have the same typed failures, tagged interchange, and tree-walker/VM behavior as the rest of the language.

def bounded_units

def bounded_units()

def bounded_sparse

def bounded_sparse()

def bounded_complex_linalg

def bounded_complex_linalg()

def main

def main() !{observe.record}

Effects !{observe.record}