Examples

All examples files are under examples/ directory. Run any of them with:

physika examples/<name>.phyk

Add --print-ast or --print-code to inspect the generated AST or PyTorch source.

---

Data Types

Physika supports three core numeric types ℝ for Real number, ℤ for integer and ℂ for complex number which can be used to declare scalar, arrays, tensors, etc.

Numeric: ℝ and ℤ

# -----------------------
# Integer type (ℤ)
# -----------------------
a : ℤ = 10
b : ℤ = 3
z_add = a + b
z_array: Z[2] = [1, 2]

a
b
z_add
z_array


# -----------------------
# Real type (ℝ)
# -----------------------
x : ℝ = 3.14
y : ℝ = 2
r_mul : ℝ = x * y

x
y
r_mul


# ---------------------------
# working with both ℤ and ℝ
# ---------------------------
z_number: Z = 1
r_number: R = 2.0
result = z_number * r_number
result


# -----------------------
# Negative values
# -----------------------
neg_int : ℤ = -7
neg_float : ℝ = -3.14
neg_array: R[3] = [-1, -2.0, -3]

neg_int
neg_float
neg_array

Complex: ℂ

# -------------------------------------------
# Basic complex declarations
# -------------------------------------------
x : ℂ = 3 + 1j
y : ℂ = 5 + 3j
x
y

# -------------------------------------------
# Addition and Subtraction
# -------------------------------------------
complex_add: ℂ = x + y
complex_add

x - y



# -------------------------------------------
# Multiplication
# -------------------------------------------
# (3 + 1j) * (5 + 3j)
# 3 * (5 + 3j) + 1j * (5 + 3j)
# 15 + 9j + 5j + 3j**2
# 15 + 14j + 3j**2
# since, j**2 = -1
# 15 + 14j - 3
# 12 + 14j
# -------------------------------------------

complex_mul: ℂ = x * y
complex_mul


# -------------------------------------------
# complex in function (calculating magnitude)
# -------------------------------------------
# x = 3 + 1j 
# formula -> square_root(a**2 + b**2)
# a = 3
# b = 1
# square_root(9 + 1)
# square_root(10)
# -> 3.162
# -------------------------------------------
def magnitude(z: ℂ): ℂ:
    return abs(z)

magnitude(x)




# -------------------------------------------
# Array of complex numbers
# -------------------------------------------

array_imag: ℂ[3] = [1j, 2j, 3j]

array_complex: ℂ[3] = [1 + 9j, 7 + 2j, 3 + 5j]

nested_complex: ℂ[2, 2] = [
    [1 + 2j, 3 + 4j],
    [4 + 9j, 7 + 2j]
]


array_imag
array_complex
nested_complex

Arrays

x : \mathbb{R}[6] = [1, 2, 3, 5, 6, 7]
y : ℝ[3] = x[0:3] + x[0:3]
z : \R[3] = y + [1, 3, 4]

x
y
z

u0: R[4, 4] = [
    [0.00, 0.00, 0.00, 0.00],
    [0.00, 0.75, 0.75, 0.00],
    [0.00, 0.75, 0.75, 10.00],
    [0.00, 0.00, 0.00, 0.00]
]

u00 : ℝ = u0[0, 0]
u01 : ℝ = u0[0, 1]
u02 : ℝ = u0[0, 2]
u03 : ℝ = u0[0, 3]

u10 : ℝ = u0[1, 0]
u11 : ℝ = u0[1, 1]
u12 : ℝ = u0[1, 2]
u13 : ℝ = u0[1, 3]

u20 : ℝ = u0[2, 0]
u21 : ℝ = u0[2, 1]
u22 : ℝ = u0[2, 2]
u23 : ℝ = u0[2, 3]

u30 : ℝ = u0[3, 0]
u31 : ℝ = u0[3, 1]
u32 : ℝ = u0[3, 2]
u33 : ℝ = u0[3, 3]

A : ℝ[3, 3] = [
    [1, 0, 0],
    [0, 2, 0],
    [0, 0, 3]
]

A00 : ℝ = A[0, 0]
A11 : ℝ = A[1, 1]
A22 : ℝ = A[2, 2]

T : ℝ[2, 3, 4] = [
    [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]],
    [[13, 14, 15, 16], [17, 18, 19, 20], [21, 22, 23, 24]]
]

T0 : ℝ[3, 4] = T[0]
T12 : ℝ[4] = T[1, 2]
T000 : ℝ = T[0, 0, 0]
T123 : ℝ = T[1, 2, 3]
T012 : ℝ = T[0, 1, 2]


# -------------------------------------------------------
# Examples of updating arrays through indexing
# -------------------------------------------------------

# -----------------------
# Program level
# -----------------------

prog_1d : ℝ[2] = [1.0, 1.0]
prog_2d : ℝ[2, 2] = [[1.0, 1.0], [1.0, 1.0]]
prog_3d : ℝ[2, 2, 2] = [[[1.0, 1.0], [1.0, 1.0]], [[1.0, 1.0], [1.0, 1.0]]]

prog_1d[1] = 2
prog_2d[1, 1] = 2
prog_3d[1, 1, 1] = 2

prog_1d
prog_2d
prog_3d


# -----------------------
# Function level
# -----------------------

func_1d : ℝ[2] = [1.0, 1.0]
func_2d : ℝ[2, 2] = [[1.0, 1.0], [1.0, 1.0]]
func_3d : ℝ[2, 2, 2] = [[[1.0, 1.0], [1.0, 1.0]], [[1.0, 1.0], [1.0, 1.0]]]

def update_1d_array(x: R[m]): R[m]:
    x[1] = 3
    return x

def update_3d_array(x: R[m, n, o]): R[m, n, o]:
    x[1, 1, 1] = 3
    return x

def update_2d_array(x: R[m, n]): R[m, n]:
    x[1, 1] = 3
    return x

update_1d_array(func_1d)
update_2d_array(func_2d)
update_3d_array(func_3d)

Matrices

# 2x3 matrix
A: R[2,3] = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
A

# Scalar multiplication (broadcasting)
B: R[2,3] = 2.0 * A
B

# Element-wise addition
C: R[2,3] = A + B
C

scale: R = 2.0

def transform(M: R[2,2]): R[2,2]:
    return M * scale + 1.0

X: R[2,2] = [[1.0, 2.0], [3.0, 4.0]]
result: R[2,2] = transform(X)
result

# Dot product of two vectors
a: R[3] = [1.0, 2.0, 3.0]
b: R[3] = [4.0, 5.0, 6.0]

result: R = a @ b
result

# Matrix-vector multiplication
A: R[2,3] = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
x: R[3, 1] = [[1.0], [2.0], [3.0]]

result: R[2, 1] = A @ x
result

Tensors (contravariant / covariant)

# Contravariant vector v^α (like a column vector in standard basis)
v: R[+3] = [1.0, 2.0, 3.0]
v

# Covariant vector ωα (like a row vector, dual space)
omega: R[-3] = [1.0, 2.0, 3.0]
omega

# Linear map: T^α_β applied to vector v^β
T: R[+3,-3] = [[2.0, 0.0, 0.0],
                [0.0, 3.0, 0.0],
                [0.0, 0.0, 4.0]]

v: R[+3, +1] = [[1.0], [2.0], [3.0]]

# Contraction: T^α_β v^β = w^α
w: R[+3, +1] = T @ v
w

# Works - single line
#A: R[2,2] = [[1.0, 2.0], [3.0, 4.0]]

# Also works - multi-line
#B: R[2,2] = [
#    [1.0, 2.0],
#    [3.0, 4.0]
#]

# Also works - with tensor variance
#T: R[+3,-3] = [
#    [1.0, 0.0, 0.0],
#    [0.0, 2.0, 0.0],
#    [0.0, 0.0, 3.0]
#]

#T

# Linear map: T^α_β applied to vector v^β
T: R[+3,-3] = [[2.0, 0.0, 0.0], # Add support for ten
               [0.0, 3.0, 0.0],
               [0.0, 0.0, 4.0]]

v: R[+3, +1] = [[1.0], [2.0], [3.0]]

# Contraction: T^α_β v^β = w^α
w: R[+3, +1] = T @ v
w

T

---

Control Flow Operators & Differentiation

Conditionals

# Examples of if-else in different contexts

# 1. Inside a function 
# if x > 0:
#   f(x) = x**2
# else:
#   f(x) = -x
def f(x: ℝ): ℝ:
    if x > 0:
        return x * x
    else:
        return - x

a : ℝ = 2.0
b : ℝ = -1.5
f(a)
grad(f(a), a)
f(b)
grad(f(b), b)

# 2. if-else with assignments inside a function
def clamp(x: ℝ): ℝ:
    if x > 1:
        y = 1
    else:
        y = x
    if y < - 1:
        y = - 1
    return y

clamp(3)
clamp(0.5)
clamp(-5)

# 3. Nested if-else inside a function
# ``classify`` returns -1, 0, 1, or 2 based on four regions
def classify(x: ℝ): ℝ:
    if x > 0.5:
        if x > 1.5:
            return 2
        else:
            return 1
    else:
        if x < - 0.5:
            return - 1
        else:
            return 0

classify(2)
classify(0)
classify(-1)


# 4. if-else inside a class method body 
# PiecewiseNet: y = x**2 if x > 0 else -x
class PiecewiseNet(threshold: ℝ):
    def λ(x: ℝ) -> ℝ:
        if x > threshold:
            y = x * x
        else:
            y = - x
        return y
    def loss(pred: ℝ, target: ℝ) -> ℝ:
        return (pred - target)**2

net = PiecewiseNet(0)
a : ℝ = 2.0
net(a)
grad(net(a), a)
b : ℝ = -1.5
net(b)
grad(net(b), b)


# 5. At top-level
x : ℝ = 0.3
if x > 0.5:
    y = 3 * (x - 0.75)**2
else:
    y = x**2 + 2

y

Differentiable conditional

# Differentiable conditional

def f(x: ℝ): ℝ:
    if x > 0.0:
        return x * x
    else:
        return - x

# positive bracnh
a : ℝ = 3.0
f(a)
grad(f(a), a)

# negative branch
b : ℝ = - 2.0
f(b)
grad(f(b), b)

Differentiable sin / cos

# -------------------------------------------------
# single arg functions
# -------------------------------------------------

def torch_funcs_with_scalar_R(x: ℝ): ℝ[m]:
    result_sin: ℝ = sin(x)
    result_cos: ℝ = cos(x)
    result_exp: ℝ = exp(x)
    result_sqrt: ℝ = sqrt(x)
    result_log: ℝ = log(x)
    result_abs: ℝ = abs(x)
    return [result_sin, result_cos, result_exp, result_sqrt, result_log, result_abs]

def check_diff_torch_funcs(x: ℝ): ℝ:
    results = grad(torch_funcs_with_scalar_R, x)
    return results

x: ℝ = 1.0
torch_funcs_with_scalar_R(x)
check_diff_torch_funcs(5.0)


# -------------------------------------------------
# multiple args functions
# -------------------------------------------------

x_matrix: ℝ[3, 3] = [
    [1, -1, 0],
    [-1, 0, 0],
    [0, 0, 0]
]
rolled_neg: ℝ[3, 3] = roll(x_matrix, -1)
rolled_pos: ℝ[3, 3] = roll(x_matrix, 1)

rolled_neg
rolled_pos



# -------------------------------------------------
# Differentiable sin/cos in conditional statement
# if x > 0:
#     f(x) = cos(x) -> f'(x) = -sin(x)
# else:
#     f(x) = sin(x) -> f'(x) = -cos(x)
# -------------------------------------------------

def f(x: ℝ): ℝ:
    if x > 0:
        return cos(x)
    else:
        return sin(x)

x0 : ℝ = - 1.5
f(x0)
grad(f(x0), x0)

x1 : ℝ = - 0.5
f(x1)
grad(f(x1), x1)

x2 : ℝ = 0.5
f(x2)
grad(f(x2), x2)

x3 : ℝ = 1.5
f(x3)
grad(f(x3), x3)

x4 : ℝ = 3.14
f(x4)
grad(f(x4), x4)

Threshold optimisation (gradient descent)

# Threshold Optimisation:
# A piecewise loss with two regimes separated at t = 0.5:
#  L(t) =  if t > 0.5:  
#               3·(t - 0.75)² + 0.1   -> HIGH regime, global min at t = 0.75
#          else:
#               t**2  + 2             -> LOW regime,  local  min at t= 0.0

def L(t: ℝ): ℝ:
    if t > 0.5:
        return 3 * (t - 0.75) **2 + 0.1
    else:
        return t**2 + 2

# Trajectory A: High regime (t > 0.5)
# Converges to the global minimum t* = 0.75
t0 : ℝ = 0.9
L(t0)
grad(L(t0), t0)

t1 : ℝ = t0 - 0.2 * grad(L(t0), t0)
L(t1)
grad(L(t1), t1)

t2 : ℝ = t1 - 0.2 * grad(L(t1), t1)
L(t2)
grad(L(t2), t2)

t3 : ℝ = t2 - 0.2 * grad(L(t2), t2)
L(t3)
grad(L(t3), t3)
t3

# Trajectory B: LOW regime (t <= 0.5)
s0 : ℝ = 0.3
L(s0)
grad(L(s0), s0)

s1 : ℝ = s0 - 0.2 * grad(L(s0), s0)
L(s1)
grad(L(s1), s1)

s2 : ℝ = s1 - 0.2 * grad(L(s1), s1)
L(s2)
grad(L(s2), s2)

s3 : ℝ = s2 - 0.2 * grad(L(s2), s2)
L(s3)
grad(L(s3), s3)

s4 : ℝ = s3 - 0.2 * grad(L(s3), s3)
L(s4)
grad(L(s4), s4)
s4

For loops and for-expressions

# Physika for-loops

# 1. Top-level imperative loops
# Implicit range
# Size inferred from array indexing in body

# 1.1 Sum of array elements
arr: ℝ[5] = [1, 2, 3, 4, 5]
total: ℝ = 0
for i:
    total += arr[i]
total

# 1.2 Sum of squares
X: ℝ[4] = [1, 2, 3, 4]
sum_sq: ℝ = 0
for i:
    sum_sq += X[i] ** 2
sum_sq

# 1.3 MSE
y: ℝ[4] = [2, 4, 6, 8]
mse: ℝ = 0
for i:
    mse += (X[i] - y[i]) ** 2
mse

# 1.4 Index assignment
src : ℝ[5] = [1, 2, 3, 4, 5]
dst : ℝ[5] = [0, 0, 0, 0, 0]
for i:
    dst[i] = src[i] * src[i]
dst

# Explicit range
# for i: ℕ(end) or for i: ℕ(start, end)

# 1.5 ℕ(start, end)
start : ℝ = 10
end : ℝ = 20
total: ℝ = 0
for i: ℕ(start, end):
    total += i
total

# 1.6 Nested range
# outer loop variable used as inner bound
n: ℝ = 10
for i: ℕ(n):
    for j: ℕ(i, 10):
        total += i
        total += j
total

# 2. Compact for-expressions
# Explicit size
# 2.1 1D array construction
a : ℝ[5] = for i : ℕ(5) -> i * 1
a

# 2.2 Mathematic operations
cos_wave : ℝ[6] = for i : ℕ(6) → cos(i * 0.5)
cos_wave

# 2.3 Nested for-loops to produce multi-dimensional arrays
# 2D matrix
add : ℝ[3, 4] = for i : ℕ(3) → for j : ℕ(4) → i + j
add

# 2.4 3D tensor
t : ℝ[2, 3, 4] = for i : ℕ(2) → for j : ℕ(3) → for k : ℕ(4) → i + j + k
t

# For-loop range inferred from indexed array
# 2.5 Single array transform
arr : ℝ[5] = [1.0, 2.0, 3.0, 4.0, 5.0]
doubled : ℝ[5] = for i → arr[i] * 2.0
doubled

# 2.6 Element-wise product of two arrays
u : ℝ[4] = [1.0, 2.0, 3.0, 4.0]
v : ℝ[4] = [4.0, 3.0, 2.0, 1.0]
dot_elems : ℝ[4] = for i → u[i] * v[i]
dot_elems

# Inside functions

# 2.7 Outer product: C[i,j] = u[i] * v[j]
def outer_product(u : ℝ[n], v : ℝ[m]): ℝ[n, m]:
    return for i → for j → u[i] * v[j]

p : ℝ[3] = [1, 2, 3]
q : ℝ[4] = [10, 20, 30, 40]
outer_product(p, q)

# 2.8 Dot product as reduction with sum()
x : ℝ[4] = [1, 0, 0, 0]
y : ℝ[4] = [0, 1, 0, 0]
dot : ℝ = sum(for i → x[i] * y[i])
dot

# 3. Function body with implicit range (for i:)
# 3.1 Loop assignment
def get_last(arr : ℝ[n]): ℝ:
    cur : ℝ = 0
    for i:
        cur = arr[i]
    return cur
vals : ℝ[5] = [3, 1, 4, 1, 5]
get_last(vals)

# 4. Function body with explicit range (for i: ℕ(n):)
# 4.1 Single level iteration
def iter_prod(n : ℝ): ℝ:
    total : ℝ = 0
    for i: ℕ(n):
        total += i * 1
    return total
iter_prod(10)

# 4.2 Start and end from parameters
def partial_sum(arr : ℝ[n], low : ℝ, high : ℝ): ℝ:
    total : ℝ = 0
    for i: ℕ(low, high):
        total += arr[i]
    return total
data : ℝ[6] = [1, 2, 3, 4, 5, 6]
partial_sum(data, 2, 5)

# 5. Function body with accumulation loops
# 5.1 2 iter variables
def outer_accum(u : ℝ[n], v : ℝ[m]): ℝ[n, m]:
    C : ℝ[n, m]
    for i j:
        C[i, j] += u[i] * v[j]
    return C
outer_accum(p, q)

# 5.2 3 iter variables
def matmul_physika(A : ℝ[n, m], B : ℝ[m, o]): ℝ[n, o]:
    C : ℝ[n, o]
    for i j k:
        C[i, j] += A[i, k] * B[k, j]
    return C

A : ℝ[2, 4] = [[1, 2, 3, 4], [0, 1, 1, 2]]
B : ℝ[4, 1] = [[1], [0], [0], [2]]
matmul_physika(A, B)

# 5.3 Sequential loops
# Computes (A@B)@A
A2 : ℝ[2, 2] = [[1, 2], [0, 1]]
B2 : ℝ[2, 2] = [[1, 0], [0, 2]]
def chain_mm(A : ℝ[m, n], B : ℝ[n, o]): ℝ[m, o]:
    C : ℝ[2, 2]
    for i j k:
        C[i, j] += A[i, k] * B[k, j]
    D : ℝ[2, 2]
    for i j k:
        D[i, j] += C[i, k] * A[k, j]
    return D
chain_mm(A2, B2)



# 5.5 Four iteration variables and one reduction dimension
# 3D tensor contraction T[i,j,l] = Σ_k A[i,k]*B[k,j]*C[k,l]
C_mat : ℝ[2, 3] = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
def tensor_contraction(A : ℝ[m, n], B : ℝ[n, o], C : ℝ[n, r]): ℝ[m, o, r]:
    T : ℝ[m, o, r]
    for i j l k:
        T[i, j, l] += A[i, k] * B[k, j] * C[k, l]
    return T
tensor_contraction(A2, B2, C_mat)

# 6. For loops with conditionals

# 6.1 For loop inside if-else branch
# Branch selects which computation to run
def sum_or_sum_sq(arr : ℝ[n], sq : ℝ): ℝ:
    if sq > 0.0:
        return sum(for i → arr[i] ** 2)
    else:
        return sum(for i → arr[i])
w : ℝ[4] = [1.0, 2.0, 3.0, 4.0]
sum_or_sum_sq(w, 1.0)
sum_or_sum_sq(w, 0.0 - 1.0)

# 6.2 body_for inside if/else — each branch runs a different loop
def abs_sum(arr : ℝ[n]): ℝ:
    total : ℝ = 0.0
    for i:
        total += arr[i]
    if total > 0.0:
        return total
    else:
        return 0.0 - total
pos : ℝ[4] = [1.0, 2.0, 3.0, 4.0]
neg : ℝ[4] = [0.0 - 1.0, 0.0 - 2.0, 0.0 - 3.0, 0.0 - 4.0]
abs_sum(pos)
abs_sum(neg)

# 6.3 If (one-sided) inside for loop body
# Accumulate only elements above zero
def sum_positive(arr : ℝ[n]): ℝ:
    total : ℝ = 0.0
    for i:
        if arr[i] > 0.0:
            total += arr[i]
    return total
mixed : ℝ[5] = [1.0, 0.0 - 2.0, 3.0, 0.0 - 4.0, 5.0]
sum_positive(mixed)

# 6.4 If-else inside for loop body
# Sum of absolute values
def sum_abs(arr : ℝ[n]): ℝ:
    total : ℝ = 0.0
    for i:
        if arr[i] > 0.0:
            total += arr[i]
        else:
            total += 0.0 - arr[i]
    return total
sum_abs(mixed)

# 6.5 If-else inside body_for
# Count elements that exceed a threshold
def count_above(arr : ℝ[n], thresh : ℝ): ℝ:
    count : ℝ = 0.0
    for i:
        if arr[i] > thresh:
            count += 1
    return count
data : ℝ[6] = [1, 5, 2, 8, 3, 7]
count_above(data, 3)

# 6.6 If-else inside body_for_range
# Count elements in a slice that exceed a threshold
def count_above_range(arr : ℝ[n], lo : ℝ, hi : ℝ, thresh : ℝ): ℝ:
    count : ℝ = 0.0
    for i: ℕ(lo, hi):
        if arr[i] > thresh:
            count += 1
    return count
data : ℝ[6] = [1, 5, 2, 8, 3, 7]
count_above_range(data, 1, 5, 3)

# 6.7 Deeply nested for/if — 10 levels
# Alternates body_if_else → body_for → loop_if_else → loop_for_range → ...
def deep_nest(arr : ℝ[n]): ℝ:
    n : R = 3
    a : ℝ = 0
    if a > - 1:
        for i: N(2):
            if arr[i] > 0:
                for j: N(2):
                    if arr[j] < 100:
                        for k: N(2):
                            if arr[k] != 0:
                                for l: N(2):
                                    if a < 10:
                                        a += 1
                                    else:
                                        a += 2
                            else:
                                a = a + 3
                    else:
                        a = a + 4
            else:
                a = a + 5
    else:
        a = - 1
    return a
arr2 : ℝ[3] = [1, 2, 4]
deep_nest(arr2)

# 7. Top-level nested for/if/else
# Same structures as section 6 but written as a top-level program.

# 7.1 If inside explicit range for loop
arr3 : ℝ[5] = [1, - 2, 3, - 4, 5]
pos_sum : ℝ = 0
for i: ℕ(5):
    if arr3[i] > 0:
        pos_sum += arr3[i]
pos_sum

# 7.2 If-else inside explicit range for loop
abs_total : ℝ = 0.0
for i: ℕ(5):
    if arr3[i] > 0.0:
        abs_total += arr3[i]
    else:
        abs_total += 0.0 - arr3[i]
abs_total

# 7.3 If inside implicit-range for
pos_sum2 : ℝ = 0.0
for i:
    if arr3[i] > 0.0:
        pos_sum2 += arr3[i]
pos_sum2

# 7.4 Nested implicit for inside explicit conditional for
mat : ℝ[3] = [1.0, 2.0, 3.0]
vec : ℝ[3] = [4.0, 5.0, 6.0]
dot2 : ℝ = 0.0
for i: ℕ(3):
    if mat[i] > 0.0:
        dot2 += mat[i] * vec[i]
dot2

# 7.5 Deeply nested for/if/else at top level (6 levels)
arr4 : ℝ[2] = [1.0, 4.0]
a : ℝ = 0.0
for i: ℕ(2):
    if arr4[i] > 0.0:
        for j: ℕ(2):
            if arr4[j] < 100.0:
                for k: ℕ(2):
                    if arr4[k] != 0.0:
                        if a < 10.0:
                            a += 1
                        else:
                            a += 2
                    else:
                        a += 3
            else:
                a += 4
    else:
        a += 5
a

# 8. List-comprehension style for-expressions with conditionals (top-level)

# 8.1 Explicit single var
# Select computation via if/else, produce ℝ[5]
vals : ℝ[5] = [1, -2, 3, 4, 5]
flag : ℝ = 1
res : ℝ[5] = for i → vals[i]
if flag > 0:
    res = for i : ℕ(5) -> vals[i] ** 3
else:
    res = for i : ℕ(5) -> vals[i]
res
# Sum of cubes when positive, sum values otherwise
total_c : ℝ = sum(for i → res[i])
total_c


# 8.3 Nested explicit two iter var for loop
# Build a 3×4 weight matrix conditioned on scale
# scale > 1 → W[i,j] = (i + j) * scale; else W[i,j] = i + j
scale : ℝ = 2
W : ℝ[3, 4] = for i : ℕ(3) -> for j : ℕ(4) -> i + j
if scale > 1:
    W = for i : ℕ(3) → for j : ℕ(4) → (i + j) * scale
else:
    W = for i : ℕ(3) → for j : ℕ(4) → i + j
W

# 8.4 Implicit two-var
# Element wise product of two arrays
# for i:
#   row_sums[i] = Σ_j u[i]*v[j]
u2 : ℝ[3] = [1, 2, 3]
v2 : ℝ[3] = [4, 5, 6]
row_sums : ℝ[3] = for i → sum(for j → u2[i] * v2[j])
row_sums

# 8.5 Implicit multiple loop variables (function level)
def outer_product(u : ℝ[n], v : ℝ[m]): ℝ[n, m]:
    results : ℝ[n, m]
    for i j:
        results[i, j] = u[i] * v[j]
    return results

x: R[4] = [1, 2, 3, 4]
y: R[4] = [0, 5, 6, 7]
outer_product(x, y)

# 8.6 Mixed implicit/explicit
data2 : ℝ[4] = [10, 20, 30, 40]
norm_flag : ℝ = 1
normed : ℝ[4] = for i → data2[i]
if norm_flag > 0:
    normed = for i : ℕ(4) → data2[i] * (1 / (i + 1))
else:
    normed = for i : ℕ(4) → data2[i]
normed

# 9. Indexing array (1d, 2d) inside for loop

# helper functions
def get_array_length(x: ℝ[m]): ℝ:
    total: ℝ = 0
    for i:
        curr = x[i]
        total += 1
    return total

def get_2d_array_num_rows(x: ℝ[m, n]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

# 9.1 Top-level 1d array indexing
sample_1d_array: ℝ[3] = [1, 2, 3]
length_array: ℝ = get_array_length(sample_1d_array)
for i:ℕ(length_array):
    sample_1d_array[i] = i*2
sample_1d_array

# 9.2 Top-level 2d array indexing
sample_2d_array: ℝ[2, 2] = [
    [1, 1],
    [1, 1]
]
rows: R = get_2d_array_num_rows(sample_2d_array)
cols: R = get_array_length(sample_2d_array[0])
for i:ℕ(rows):
    for j:ℕ(cols):
        sample_2d_array[i, j] = j*2
sample_2d_array

# 9.3 Top-level 3d array indexing
sample_3d_array: ℝ[2, 2, 2] = [
    [
        [1, 2],
        [1, 2]
    ],
    [
        [1, 2],
        [1, 2]
    ]
]
for i:ℕ(2):
    for j:ℕ(2):
        for k:ℕ(2):
            sample_3d_array[i, j, k] = j*2 + k
sample_3d_array

# 9.4 1d array indexing inside function body
def manipulate_1d_array(x: ℝ[m]): ℝ[n]:
    m = get_array_length(x)
    for i:ℕ(m):
        x[i] = i*2
    return x
arr1d: ℝ[3] = [1, 2, 3]
manipulate_1d_array(arr1d)

# 9.5 2d array indexing inside function body
def manipulate_2d_array(x: ℝ[m, n]): ℝ[m, n]:
    m = get_2d_array_num_rows(x)
    n = get_array_length(x[0])
    for i:ℕ(m):
        for j:ℕ(n):
            x[i, j] = j*2
    return x
arr2d: ℝ[2, 2] = [
    [1, 1],
    [1, 1]
]
manipulate_2d_array(arr2d)

def manipulate_3d_array(x: ℝ[m, n, o]): ℝ[m, n, o]:
    for i:ℕ(2):
        for j:ℕ(2):
            for k:ℕ(2):
                x[i, j, k] = i*2 + j + k
    return x
arr3d: ℝ[2, 2, 2] = [
    [
        [3, 2],
        [1, 1]
    ],
    [
        [1, 4],
        [1, 2]
    ]
]
manipulate_3d_array(arr3d)

Differentiable for loops

All loop forms are differentiable. grad() returns a gradient vector for scalar-output functions and a full Jacobian matrix for vector-output functions.

# Examples for differentiable for-loops in Physika
# Gradient of for-expr
# f(s) = sum(for i : ℕ(4) → s * i) = s(0+1+2+3) = 6s
# f'(s) = 6
def sum_for_expr(s : ℝ): ℝ:
    return sum(for i : ℕ(4) → s * i)
s0 : ℝ = 2.0
sum_for_expr(s0)
grad(sum_for_expr(s0), s0)

# Gradient of implicit for-loop
# g(s) = sum(s * arr[i] for i) = sum([1s, 2s, 3s, 4s])= s * (1+2+3+4) = 10*s
# g'(s) = 10
def dot_with_arr(s : ℝ): ℝ:
    a3 : ℝ[4] = [1.0, 2.0, 3.0, 4.0]
    result : ℝ = 0.0
    for i:
        result += s * a3[i]
    return result
s1 : ℝ = 1.0
dot_with_arr(s1)
grad(dot_with_arr(s1), s1)

# Gradient of multiple iter variable loops
# h(s) = sum(s * A @ I) = s * sum(A) = s(1+2+3+4) = 10s
# h'(s) = 10
def matmul_scale(s : ℝ): ℝ:
    A3 : ℝ[2, 2] = [[1.0, 2.0], [3.0, 4.0]]
    I : ℝ[2, 2] = [[1.0, 0.0], [0.0, 1.0]]
    C3 : ℝ[2.0, 2.0]
    for i j k:
        C3[i, j] += s * A3[i, k] * I[k, j]
    return sum(C3)
s2 : ℝ = 1.0
matmul_scale(s2)
grad(matmul_scale(s2), s2)

# Gradient of nested for i: ℕ(n) / for j: ℕ(i, 10) loop
# k(s) = s * Σ {i=0}^{9} Σ{j=i}^{9} (i + j)
# k'(s) = Σ {i=0}^{9} Σ_{j=i}^{9} (i + j) = 495
def nested_sum(s : ℝ): ℝ:
    result : ℝ = 0
    for i: ℕ(10):
        for j: ℕ(i, 10):
            result += s * i * 1.0 + s * j * 1.0
    return result
s3 : ℝ = 1.0
nested_sum(s3)
grad(nested_sum(s3), s3)


# Gradients for tensors using for loops
# scale_vec: f(x) = [x, 2x, 3x]
# Jac(f)(x) = [1, 2, 3]
def scale_vec(x : ℝ): ℝ[3]:
    return for i : ℕ(3) → x * (i + 1)
s : ℝ = 2.0
scale_vec(s)
grad(scale_vec(s), s)

# sq_vec: f(x) = [x², 2x², 3x², 4x²]
# Jac(f)(x) = [2x, 4x, 6x, 8x]
def sq_vec(x : ℝ): ℝ[4]:
    return for i : ℕ(4) → x ** 2 * (i + 1)
sv : ℝ = 3.0
sq_vec(sv)
grad(sq_vec(sv), sv)

# cos_freqs: f(x) = [cos(x), cos(2x), cos(3x), cos(4x)]
# Jac(f)(x) = [-sin(x), -2sin(2x), -3sin(3x), -4sin(4x)]
def cos_freqs(x : ℝ): ℝ[4]:
    return for i : ℕ(4) → cos(x * (i + 1))
x : ℝ = 0.5
cos_freqs(x)
grad(cos_freqs(x), x)

# elementwise_sq: f(x)[i] = x[i]²
# f(x) ∈ ℝ[n]
# Jac(f)(x) = diag(2·x) ∈ ℝ[n,n]
def elementwise_sq(x : ℝ[n]): ℝ[n]:
    return for i → x[i] ** 2
ev : ℝ[3] = [1.0, 2.0, 3.0]
elementwise_sq(ev)
grad(elementwise_sq(ev), ev)

Factorial (recursive)

# Factorial Example

# Recursive Function with if/else
def fact(n: ℝ): ℝ:
    if n == 0.0:
        return 1.0
    else:
        return n * fact(n - 1.0)

fact(0.0)   # 1
fact(1.0)   # 1
fact(2.0)   # 2
fact(3.0)   # 6
fact(4.0)   # 24
fact(5.0)   # 120
fact(10.0)  # 3628800

---

Classes

Physika classes

# Physika classes

class Vec:
    x : ℝ
    y : ℝ
    def dot(other : Vec) : ℝ:
        return this.x * other.x + this.y * other.y
    def scale(s : ℝ) : Vec:
        return Vec(this.x * s, this.y * s)
    def norm_sq() : ℝ:
        return this.x * this.x + this.y * this.y

class Particle:
    pos : ℝ[2]
    vel : ℝ[2]
    mass : ℝ
    def kinetic_energy() : ℝ:
        return 0.5 * this.mass * sum(this.vel * this.vel)
    def step(force : ℝ[2], dt : ℝ) : Particle:
        acc : ℝ[2] = force * (1.0 / this.mass)
        new_vel : ℝ[2] = this.vel + acc * dt
        new_pos : ℝ[2] = this.pos + this.vel * dt
        return Particle(new_pos, new_vel, this.mass)

# Vec
a = Vec(3.0, 4.0)
b = Vec(1.0, 0.0)
a.x
a.y
dot_ab : ℝ = a.dot(b)
dot_ab
c = a.scale(4)
c.x
c.y


# Particle object
pos0 : ℝ[2] = [0.0, 10.0]
vel0 : ℝ[2] = [1.0, 0.0]
gravity : ℝ[2] = [0.0, -9.81]
p = Particle(pos0, vel0, 9.0)
ke0 : ℝ = p.kinetic_energy()
ke0
p1 = p.step(gravity, 0.5)
p1.pos
p2 = p1.step(gravity, 0.5)
p2.pos

# Computing grads through Physika class methods
# grad(Particle.kinetic_enery, v)
def ke_wrt_vel(vel : ℝ[2]) : ℝ:
    particle : Particle = Particle(pos0, vel, 1.0)
    return particle.kinetic_energy()

v : ℝ[2] = [2.0, 3.4]
ke0_v : ℝ = ke_wrt_vel(v)
ke0_v
dKE_dv : ℝ[2] = grad(ke_wrt_vel(v), v)
dKE_dv

# grad of kinetic energy w.r.t. a velocity component:
def ke_vy(vy : ℝ) : ℝ:
    # v is composed of [vx, vy]
    # Lets differentiate kinetic energy wrt vy component
    p : Particle = Particle(pos0, [1.0, vy], 2.0)
    return p.kinetic_energy()

vy0 : ℝ = 3.0
ke_vy(vy0)
grad(ke_vy(vy0), vy0)    # d(0.5*2*(1**2 + vy**2))/d(vy) = 2*vy = 6

# grad(Vec.norm_sq, x)
def norm_sq_wrt_x(x : ℝ) : ℝ:
    vec : Vec = Vec(x, 4.0)
    return vec.norm_sq()

x0 : ℝ = 3.0
norm_sq_wrt_x(x0)
grad(norm_sq_wrt_x(x0), x0)

# Directly computing grad of a method w.r.t. the input:
x1 : ℝ = 5.0
vec : Vec = Vec(x1, 4.0)
grad(vec.norm_sq(), x1)    # d(x**2 + 4**2)/dx = 2x = 10

x1 : ℝ = 5.0
vec : Vec = Vec(x1, 4.0)
grad(vec.x, x1)    # d(x)/dx = 1

Neural network classes

# ============================================================
# Fully Connected Network in Physika
# ============================================================

# ------------------------------------------------------------
# Helper functions 
# ------------------------------------------------------------

def get_1d_array_length(x: ℝ[m]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

def get_2d_array_num_rows(x: ℝ[m, n]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

def zero_2d_array(rows: ℝ, cols: ℝ): ℝ[m, n]:
    results: ℝ[rows, cols] = for i:ℕ(rows) -> for j:N(cols) -> j*0
    return results


# ------------------------------------------------------------
# Activation: σ(x) = 1/(1 + e^(-x))
# ------------------------------------------------------------
def sigma(x: ℝ[a,b]): ℝ[a,b]:
    rows: ℝ = get_2d_array_num_rows(x)
    cols: ℝ = get_1d_array_length(x[0])
    results: ℝ[cols, rows] = zero_2d_array(rows, cols)
    for i:ℕ(rows):
        for j:ℕ(cols):
            results[i,j] = 1.0 / (1.0 + exp(0.0 - x[i,j]))
    return results



# ============================================================
# Example 1: One-Layer Network (simple version)
# ============================================================
class OneLayerNet:
    W0: ℝ[2,3]
    c0: ℝ[2, 1] 
    w1: ℝ[1, 2]
    b1: ℝ
    def λ(x: ℝ[3, 1]) → ℝ:
        z = sigma(W0 @ x + c0)
        results = w1 @ z + b1
        return results[0,0]
    def loss(y: ℝ, target: ℝ) → ℝ:
        return (y - target) ** 2.0

W0: ℝ[2,3] = [[0.1, 0.2, 0.3], [0.4, 0.5, 0.6]]
c0: ℝ[2,1] = [[0.1],[0.2]]
w1: ℝ[1,2] = [[0.7,0.8]]
b1: ℝ = 0.3
net1 = OneLayerNet(W0, c0, w1, b1)

net1([1.0, 2.0, 3.0])

# ============================================================
# Example 2: Generalized n-Layer Network
# ============================================================

class FullyConnectedNetwork:
    W: ℝ[2,3,3]
    B: ℝ[2,3,1]
    w: ℝ[1,3]
    b: ℝ 
    n: ℕ
    def λ(x: ℝ[3, 1]) → ℝ[m, n]:
        for k:
            x = sigma(W[k] @ x + B[k])
        return w @ x + b
    def loss(y: ℝ, target: ℝ) → ℝ:
        return (y - target) ** 2.0

# 2-layer network: ℝ[3] -> ℝ[3] -> ℝ[3] -> ℝ
W: ℝ[2,3,3] = [[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]],
               [[0.2, 0.3, 0.4], [0.5, 0.6, 0.7], [0.8, 0.9, 0.1]]]
B: ℝ[2,3,1] = [
    [[0.1],[0.2],[0.3]],
    [[0.1],[0.2],[0.3]]
]
w: ℝ[1,3] = [[0.5,0.5,0.5]]
b: ℝ = 0.1

net2 = FullyConnectedNetwork(W, B, w, b, 2)

net2([[1.0], [2.0], [3.0]])
net2([[0.0], [0.0], [0.0]])
net2([[1.0], [1.0], [1.0]])

Training a fully connected network

# ============================================================
# Training a Fully Connected Network in Physika
# ============================================================


# ------------------------------------------------------------
# Helper functions 
# ------------------------------------------------------------

def get_1d_array_length(x: ℝ[m]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

def get_2d_array_num_rows(x: ℝ[m, n]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

def zero_2d_array(rows: ℝ, cols: ℝ): ℝ[m, n]:
    results: ℝ[rows, cols] = for i:ℕ(rows) -> for j:N(cols) -> j*0
    return results

# ------------------------------------------------------------
# Activation: σ(x) = 1/(1 + e^(-x))
# ------------------------------------------------------------

def sigma(x: ℝ[a,b]): ℝ[a,b]:
    rows: ℝ = get_2d_array_num_rows(x)
    cols: ℝ = get_1d_array_length(x[0])
    results: ℝ[cols, rows] = zero_2d_array(rows, cols)
    for i:ℕ(rows):
        for j:ℕ(cols):
            results[i,j] = 1.0 / (1.0 + exp(0.0 - x[i,j]))
    return results


# ------------------------------------------------------------
# Network Definition: 2-layer FCN with custom MSE loss
# R[3] -> R[3] -> R[3] -> R
# ------------------------------------------------------------
class FullyConnectedNetwork:
    W: ℝ[2,3,3]
    B: ℝ[2,3,1]
    w: ℝ[1,3]
    b: ℝ
    n: ℕ
    def λ(x: ℝ[3, 1]) → ℝ:
        for k:
            x = sigma(W[k] @ x + B[k])
        results: ℝ[m, n] = w @ x + b
        return results[0, 0]
    def loss(y: ℝ, target: ℝ) → ℝ:
        return (y - target) ** 2.0


# ------------------------------------------------------------
# Training Data
# ------------------------------------------------------------
X: ℝ[4,3,1] = [
    [[1.0],[0.0],[0.0]],
    [[0.0],[1.0],[0.0]],
    [[0.0],[0.0],[1.0]],
    [[1.0],[1.0],[1.0]]
]
y: ℝ[4] = [0.2, 0.4, 0.6, 0.9]

# ------------------------------------------------------------
# Create Network (2 hidden layers)
# ------------------------------------------------------------
W: ℝ[2,3,3] = [[[0.1, 0.2, 0.3], [0.4, 0.5, 0.6], [0.7, 0.8, 0.9]],
               [[0.2, 0.3, 0.4], [0.5, 0.6, 0.7], [0.8, 0.9, 0.1]]]
B: ℝ[2,3,1] = [
    [[0.1],[0.2],[0.3]],
    [[0.1],[0.2],[0.3]]
]
w: ℝ[1,3] = [[0.5,0.5,0.5]]
b: ℝ = 0.1

net = FullyConnectedNetwork(W, B, w, b, 2)

# ------------------------------------------------------------
# Evaluate Before Training (uses class-defined loss)
# ------------------------------------------------------------
loss_before = evaluate(net, X, y)
loss_before

# ------------------------------------------------------------
# Train Network (uses class-defined loss)
# ------------------------------------------------------------
epochs : ℕ = 1000
lr : ℝ = 0.1 
net_trained = train(net, X, y, epochs, lr)

# ------------------------------------------------------------
# Evaluate After Training
# ------------------------------------------------------------
loss_after = evaluate(net_trained, X, y)
loss_after

# ------------------------------------------------------------
# Test Predictions (targets: 0.2, 0.4, 0.6, 0.9)
# ------------------------------------------------------------
net_trained([1.0, 0.0, 0.0])
net_trained([0.0, 1.0, 0.0])
net_trained([0.0, 0.0, 1.0])
net_trained([1.0, 1.0, 1.0])

---

Physics Simulations

Note

The following three examples open a GUI window (matplotlib / PyVista).

Harmonic oscillator

# ============================================================
# Harmonic Oscillator - Differentiable Time Evolution
# ============================================================
# Equation: 
#   m * d²x/dt² = -k * x
# 
# Angular frequence:
#   ω = √(k / m)
#
# General solution (sine-cosine form):
#   x(t) = A * cos(ω t) + B * sin(ω t)
# 
# Velocity:
#   v(t) = dx/dt = -Aω sin(ω t) + Bω cos(ω t)
# 
# Initial conditions:
#   x(0) = x0
#   v(0) = v0
#
# Apply initial conditions:
#   x(0) = A = x0
#   v(0) = Bω = v0 -> B = v0 / ω
# 
# Final solution:
#   x(t) = x0 * cos(ω t) + (v0 / ω) * sin(ω t)
# ============================================================


# ------------------------------------------------------------
# Helper functions
# ------------------------------------------------------------

def get_1d_array_length(x: ℝ[m]): ℝ:
    total: ℝ = 0
    temp: ℝ = 0
    for i:
        temp = x[i]
        total += 1
    return total

def zero_1d_array(len: ℝ): ℝ[m]:
    results: ℝ[len] = for i:ℕ(len) -> i*0
    return results

def zero_2d_array(rows: ℝ, cols: ℝ): R[m, n]:
    results: ℝ[rows, cols] = for i:ℕ(rows) -> for j:ℕ(cols) -> j*0
    return results

# ------------------------------------------------------------
# Gaussian elimination method 
# ------------------------------------------------------------

def solve(A: ℝ[m, m], b: ℝ[m]) : ℝ[m]:
    n: ℝ = get_1d_array_length(b)
    Aug: ℝ[m, m] = zero_2d_array(n, n+1)
    result: ℝ[m] = zero_1d_array(n)
    # create augmented matrix
    for i:ℕ(n):
        for j:ℕ(n):
            Aug[i, j] = A[i, j]
        Aug[i, n] = b[i]
    # convert into upper triangular matrix
    for pivot:ℕ(n):
        for i:ℕ(pivot+1, n):
            factor = Aug[i, pivot] / Aug[pivot, pivot]
            for j:ℕ(pivot, n+1):
                Aug[i, j] = Aug[i, j] - factor * Aug[pivot, j]
    # back substitution
    for i:ℕ(n):
        ri = n - 1 - i
        result[ri] = Aug[ri, n]
        for j:ℕ(ri+1, n):
            result[ri] = result[ri] - Aug[ri, j] * result[j]
        result[ri] = result[ri] / Aug[ri, ri]
    return result


# ------------------------------------------------------------
# Time Evolution U(k, m, t, x0, v0)
# ------------------------------------------------------------
def U(k: ℝ, m: ℝ, t: ℝ, x0: ℝ, v0: ℝ): ℝ:
    omega: ℝ = (k / m) ** 0.5
    A: ℝ[2, 2] = [[1.0, 0.0], [0.0, omega]]
    B: ℝ[2] = [x0, v0]
    coeffs: ℝ[2] = solve(A, B)
    a: ℝ = coeffs[0]
    b: ℝ = coeffs[1]
    return a * cos(omega * t) + b * sin(omega * t)

# ------------------------------------------------------------
# Physical Constants
# ------------------------------------------------------------
k: ℝ = 1.0 # TODO: Support units
m: ℝ = 1.0

# Initial Conditions
x0: ℝ = 1.0
v0: ℝ = 0.0

# ------------------------------------------------------------
# Time Evolution: U(k, m, t, x0, v0)
# For x0=1, v0=0: x(t) = cos(ωt)
# ------------------------------------------------------------
# TODO: Add support for π symbol and value
U(k, m, 0.0 , x0, v0) # t = 0 * π
U(k, m, 1.5708, x0, v0) # t = π / 2
U(k, m, 3.1416, x0, v0) # t = π
U(k, m, 4.7124, x0, v0) # t = 3 * π / 2
U(k, m, 6.2832, x0, v0) # t = 2 * π

# ------------------------------------------------------------
# Different ICs: x0=0, v0=1
# Solution: x(t) = sin(t) for ω=1
# ------------------------------------------------------------
U(k, m, 0.0, 0.0, 1.0)
U(k, m, 1.5708, 0.0, 1.0)
U(k, m, 3.1416, 0.0, 1.0)


# ------------------------------------------------------------
# Animate the oscillator using U(t) and velocity is computed by
# differentiation
# ------------------------------------------------------------
animate(U, k, m, x0, v0, 0.0, 31.1416) # t = 10π


# ------------------------------------------------------------
# Velocity via grad(U, t) = dU/dt
# For x0=1, v0=0: x(t) = cos(t)
#                  v(t) = -sin(t)
# Expected: v(0) = 0,
#            v(π/2) = -1,
#            v(π) = 0
# ------------------------------------------------------------
t0: ℝ = 0.0
grad(U(k, m, t0, x0, v0), t0)

t1: ℝ = 1.5708
grad(U(k, m, t1, x0, v0), t1)

t2: ℝ = 3.1416
grad(U(k, m, t2, x0, v0), t2)

Simple pendulum (RK4)

# ======================================================
# RK4 ODE Solver in Physika Applied to Simple Pendulum
# ======================================================
# ----------------------------------
# Simple Pendulum System:
# dθ/dt = ω,  dω/dt = -(g/L)·sin(θ)
# State vector: x = [θ, ω]
# ----------------------------------
def pendulum(x: ℝ[2]): ℝ[2]:
    return [x[1], 0.0 - (9.81 / 1.0) * sin(x[0])]

# ----------------
# RK4 ODE Solver
# ----------------
class RK4:
    dt: ℝ
    def λ(x: ℝ[2]) → ℝ[2]:
        k1 = pendulum(x)
        k2 = pendulum(x + 0.5 * dt * k1)
        k3 = pendulum(x + 0.5 * dt * k2)
        k4 = pendulum(x + dt * k3)
        x = x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
        return x

# ------------
# Parameters
# ------------
dt: ℝ = 0.01

# ------------------------------------------------------------
# Solve: θ₀ = 0.5 rad, ω₀ = 0.0 rad/s (1000 steps = 10s)
# ------------------------------------------------------------
solver = RK4(dt)
solver([0.5, 0.0])

# Different initial conditions
solver([1.0, 0.0])

# ------------------------------------------------------------
# Visualize oscillations
# ------------------------------------------------------------
step = RK4(dt)
simulate(step, [0.5, 0.0], 1000, dt)

Spring pendulum (RK4)

# ============================================================
# Spring Pendulum with RK4
# ============================================================
# State: x = [r, θ, dr/dt, dθ/dt]
# Equations of motion (Lagrangian mechanics):
#   dr/dt   = x[2]
#   dθ/dt   = x[3]
#   d²r/dt² = r*dθ² - g*cos(θ) - k*(r - L₀)
#   d²θ/dt² = (-2*dr*dθ - g*sin(θ)) / r
# Parameters: g = 9.81, k = 50.0, L₀ = 1.0

# ------------------------------------------------------------
# Spring pendulum dynamics
# ------------------------------------------------------------
def spring_pendulum(x: ℝ[4]): ℝ[4]:
    return [x[2], x[3], x[0] * x[3] * x[3] - 9.81 * cos(x[1]) - 50.0 * (x[0] - 1.0), (0.0 - 2.0 * x[2] * x[3] - 9.81 * sin(x[1])) / x[0]]

# ------------------------------------------------------------
# RK4 Integrator
# ------------------------------------------------------------
class RK4:
    dt: ℝ
    def λ(x: ℝ[4]) → ℝ[4]:
        k1 = spring_pendulum(x)
        k2 = spring_pendulum(x + 0.5 * dt * k1)
        k3 = spring_pendulum(x + 0.5 * dt * k2)
        k4 = spring_pendulum(x + dt * k3)
        x = x + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
        return x

# --------
# Solve
# -------
dt = 0.01
solver = RK4(dt)

# r₀ = 1.2, θ₀ = 0.3 rad, dr₀ = 0, dθ₀ = 0
solver([1.2, 0.3, 0.0, 0.0])

# r₀ = 1.3, θ₀ = 0.8 rad
solver([1.3, 0.8, 0.0, 0.0])

# ----------
# Visualize
# ----------
step = RK4(dt)
simulate(step, [1.2, 0.3, 0.0, 0.0], 500, dt)

Hamiltonian Neural Network (HNN)

# ============================================================
# Hamiltonian Neural Network (HNN) in Physika
# ============================================================
# Learning the dynamics of a Simple Harmonic Oscillator
# True Hamiltonian: H(q, p) = (q² + p²) / 2
# Hamilton's equations: dq/dt = ∂H/∂p = p, dp/dt = -∂H/∂q = -q

# ------------------------------------------------------------
# Activation Function
# ------------------------------------------------------------
def tanh(x: ℝ): ℝ:
    return (exp(x) - exp(0.0 - x)) / (exp(x) + exp(0.0 - x))

# ------------------------------------------------------------
# Hamiltonian Network with Physics-Informed Loss
# Input: x = [q, p] ∈ R[2], Output: H ∈ R
# Loss enforces Hamilton's equations in a single expression
# ------------------------------------------------------------


class HamiltonianNet:
    W1: ℝ[M,N]
    b1: ℝ[M,1]
    w2: ℝ[1,M]
    b2: ℝ
    def λ(x: ℝ[N, 1]) → ℝ[2, 1]:
        h: ℝ[1, 1] = w2 @ tanh(W1 @ x + b1) + b2
        dh_grad = grad(h, x)
        dh_dp: R = dh_grad[1, 0]
        dh_dq: R = -dh_grad[0, 0]
        return [[dh_dp], [dh_dq]]
     def loss(symplectic_gradient: ℝ[2], target: ℝ[N]) → ℝ:
        lo : ℝ = (symplectic_gradient[0] - target[0]) ** 2.0 + (symplectic_gradient[1] - target[1]) ** 2.0
        return lo

# ------------------------------------------------------------
# Training Data: Simple Harmonic Oscillator
# x = [q, p], target = [dq/dt, dp/dt] = [p, -q]
# ------------------------------------------------------------
X: ℝ[8,2,1] = [
    [[0.0],[1.0]],
    [[1.0],[0.0]],
    [[0.0],[-1.0]],
    [[-1.0],[0.0]],
    [[0.5],[0.5]],
    [[-0.5],[-0.5]],
    [[0.7],[-0.7]],
    [[-0.7],[0.7]]
]

y: ℝ[8,2,1] = [
    [[1.0],[0.0]],
    [[0.0],[-1.0]],
    [[-1.0],[0.0]],
    [[0.0],[1.0]],
    [[0.5],[-0.5]],
    [[-0.5],[0.5]],
    [[-0.7],[-0.7]],
    [[0.7],[0.7]]
]

# ------------------------------------------------------------
# Initialize Hamiltonian Network (weights and bias can be hanfled in backend)
# ------------------------------------------------------------
W1: ℝ[16,2] = [[0.5, 0.1], [0.1, 0.5], [0.3, 0.3], [0.4, 0.2],
               [0.2, 0.4], [0.1, 0.1], [0.3, 0.1], [0.1, 0.3],
               [0.2, 0.2], [0.4, 0.4], [0.5, 0.3], [0.3, 0.5],
               [0.2, 0.1], [0.1, 0.2], [0.4, 0.1], [0.1, 0.4]]
b1: ℝ[16,1] = [
    [0.0],[0.0],[0.0],[0.0],
    [0.0],[0.0],[0.0],[0.0],
    [0.0],[0.0],[0.0],[0.0],
    [0.0],[0.0],[0.0],[0.0]
]

w2: ℝ[1,16] = [[
 0.1,0.1,0.1,0.1,
 0.1,0.1,0.1,0.1,
 0.1,0.1,0.1,0.1,
 0.1,0.1,0.1,0.1
]]
b2: ℝ = 0.0

H_net = HamiltonianNet(W1, b1, w2, b2)

# ------------------------------------------------------------
# Evaluate before training
# ------------------------------------------------------------
loss_before = evaluate(H_net, X, y)
loss_before

# ------------------------------------------------------------
# Train HNN using the defined loss in class
# ------------------------------------------------------------

epochs : ℕ = 500
lr : ℝ = 0.01 
H_trained = train(H_net, X, y, epochs, lr)

# ------------------------------------------------------------
# Evaluate after training
# ------------------------------------------------------------
loss_after = evaluate(H_trained, X, y)
loss_after

Symbolic support

Physika supports symbolic math via sympy, letting you declare Symbol, Function, and Equation types and use built-ins like diff, subs, lambdify and symbolic_solve to derive and solve expressions analytically.

Note

Physika uses SymPy for symbolic computation.

Core concepts

• Symbol — symbolic variables used in expressions

• Function — symbolic functions of one or more variables

• Equation — mathematical relationships between expressions

Operations

• diff(expr, var) — compute derivative w.r.t. a variable

• subs(expr, old, new) — substitute part of an expression

• lambdify(expr, vars) — convert to a numerical function

• symbolic_solve(eq, var) — solve equations analytically

#-----------------------------------------------
# Basic use case of symbolic maths with physika
#-----------------------------------------------

# declaration
x, y : Symbol
u : Function
x

# create expression
f = x**2 + y**2
f

# function call with symbols
expr = u(x, y)
expr

# substitution
subs(f, x, 3.0, y, 4.0)

# diff
f = x**3 + 2*x**2 + x
diff(f, x)

# lambdify example

expr = x**2 + y**2
f = lambdify([x, y], expr)  
f(3.0, 4.0)

# solve a linear equation
eq: Equation := 2.0*x + 3.0 = 7.0
symbolic_solve(eq, x)

Greek Letters and Scientifc Notation

Physika treats Greek letters as valid symbols and variables, allowing you to write physics and mathematics in a natural, notation-friendly style. Variables like α, β, μ, σ, λ and all standard Greek letters are supported.

Scientific notation is also supported natively — values like 1e5, 2.5e-3, or 6.674e-11 are valid numeric literals.

Note

Δ is reserved for the Laplacian operator and cannot be used as an symbol/variable.

# basic declarations
α: R = 1.0
β: R = 2.0
x: R = 1e5
y: R = 3e5


# numeric operations
results = α + β
results
z = x + y
z

# create array
greek_letters_array = [α, β]
greek_letters_array


# check gradients correctness
def f(x: R[m]): R[m]:
    return x**2 + 1

μ: R[1] = [2]
grad(f, μ) # should return 4

Import statements

Physika supports importing functions, classes from other Physika modules.

# ----------------------------------
# example with single import
# ----------------------------------

from factorial import fact

x: ℝ = 1.0
fact_results: ℝ = fact(x)
fact_results

# ----------------------------------
# example with multiple imports
# ----------------------------------

from diff_functions import f, torch_funcs_with_scalar_R

torch_funcs_results: ℝ[6] = torch_funcs_with_scalar_R(x)
torch_funcs_results

f_results: ℝ = f(x)
f_results