cobranum

Overview

cobranum is a Go plugin that adds a dense, typed N-dimensional array (1-D and 2-D) to Cobra, backed by a contiguous float64 buffer. Operations are vectorized and multi-core; on macOS they use Apple Accelerate (BLAS + LAPACK) and hand-written ARM NEON kernels. Arrays are opaque NDARRAY handles — inspect them with to_list, str or shape.

Download

Prebuilt plugins. Unzip to get cobranum.so, put it next to your program (or on the import path), and import "cobranum.so" — using the matching cobra binary.

platform	plugin	needs binary
macOS (Apple Silicon)	cobranum-0.11.1-darwin-arm64.zip	standard cobra-0.11.1-darwin-arm64
Linux (x86-64)	cobranum-0.11.0-linux-amd64.zip	cobra-0.11.0-linux-amd64-plugin (glibc ≥ 2.34)

Checksums are in MD5SUMS. The macOS/arm64 release binary is already plugin-capable; on Linux the default binary is fully static (no plugins), so the separate dynamically-linked plugin binary above is provided for plugin users.

Build from source (any platform)

Go plugins can't be cross-compiled and aren't supported on Windows, so on other platforms (Intel Mac, Linux arm64) build against your own cobra — one command, no extra tooling:

$ make build                              # a plugin-capable cobra binary
$ make plugin PLUGIN=./plugins/cobranum   # -> cobranum.so
$ ./dist/cobra your_program.cobra         # with cobranum.so on the path

Then import it and (optionally) give it a short alias:

import "cobranum.so"
let np = cobranum                      # a short alias, by convention

let a = np.array([1, 2, 3, 4])
print(np.to_list(np.mul(a, 10)))       # => [10, 20, 30, 40]

Creating arrays

Build arrays from (nested) lists, or generate them directly — generators never materialize a big Cobra list, so they are cheap even for millions of elements.

let v = np.array([1, 2, 3])                 # 1-D, shape [3]
let M = np.array([[1, 2, 3], [4, 5, 6]])    # 2-D, shape [2, 3]

np.zeros([2, 3])          # all zeros, shape [2, 3]
np.ones([4])             # [1, 1, 1, 1]
np.full([2, 2], 7.0)     # [[7, 7], [7, 7]]
np.eye(3)                # 3x3 identity
np.arange(0, 10, 2)      # [0, 2, 4, 6, 8]
np.linspace(0.0, 1.0, 5) # [0, 0.25, 0.5, 0.75, 1]

Inspecting & reshaping

let M = np.array([[1, 2, 3], [4, 5, 6]])
print(np.shape(M))                 # => [2, 3]
print(np.ndim(M), np.size(M))      # => 2  6
print(np.to_list(np.reshape(M, [3, 2])))   # => [[1, 2], [3, 4], [5, 6]]
print(np.str(np.eye(2)))           # ndarray(shape=[2, 2], data=[[1, 0], [0, 1]])
print(np.get(M, 1, 2))             # => 6   (element [1][2])
np.set(M, 0, 0, 99)                # in-place element write

Elementwise math & broadcasting

Arithmetic works elementwise on equal-shaped arrays, with a scalar, or with NumPy-style broadcasting (a row or column vector stretched across a matrix).

let a = np.array([1, 2, 3])
let b = np.array([10, 20, 30])

np.add(a, b)        # [11, 22, 33]
np.sub(b, a)        # [9, 18, 27]
np.mul(a, 2)        # scalar broadcast -> [2, 4, 6]
np.div(b, a)        # [10, 10, 10]
np.pow(a, 2)        # [1, 4, 9]
np.neg(a)           # [-1, -2, -3]

# unary math, elementwise
np.sqrt(np.array([1, 4, 9]))    # [1, 2, 3]
np.exp(np.array([0, 1]))        # [1, 2.718...]
np.log(np.array([1, 2.718281828]))
np.sin(a)   np.cos(a)   np.abs(np.array([-1, 2, -3]))

Broadcasting

let X = np.array([[1, 2, 3], [4, 5, 6]])     # shape [2, 3]
let row = np.array([10, 20, 30])             # shape [3]  -> stretched across rows
print(np.to_list(np.add(X, row)))            # [[11, 22, 33], [14, 25, 36]]

let col = np.array([[100], [200]])           # shape [2, 1] -> stretched across columns
print(np.to_list(np.add(X, col)))            # [[101, 102, 103], [204, 205, 206]]

Fused & in-place (performance)

Fused kernels do the whole computation in one pass (no temporaries); in-place variants (*_into) reuse a preallocated buffer, so hot loops never allocate.

np.hypot(a, b)            # sqrt(a*a + b*b), one fused NEON pass
np.fma(a, b, c)           # a*b + c
np.axpy(2.0, a, b)        # 2*a + b

# in-place: write the result into a preallocated dst (dst may alias an input)
let dst = np.zeros([3])
np.mul_into(dst, a, a)    # dst = a*a, no allocation
np.add_into(dst, dst, b)  # dst = dst + b
np.sqrt_into(dst, dst)    # dst = sqrt(dst)

Reductions & statistics

let s = np.array([4, 8, 15, 16, 23, 42])
np.sum(s)      # 108
np.mean(s)     # 18
np.min(s)      np.max(s)      np.prod(s)
np.var(s)      # 151.66...
np.std(s)      # 12.31...
np.norm(s)     # sqrt(sum of squares)
np.argmin(s)   np.argmax(s)   # indices of the extremes
np.to_list(np.cumsum(s))      # [4, 12, 27, 43, 66, 108]
np.to_list(np.clip(s, 10, 20))# [10, 10, 15, 16, 20, 20]

# axis reductions on a matrix (axis 0 = down columns, 1 = across rows)
let M = np.array([[1, 2, 3], [4, 5, 6]])
np.to_list(np.sum(M, 0))   # [5, 7, 9]
np.to_list(np.sum(M, 1))   # [6, 15]

Linear algebra

Matrix products, transpose, inverse, determinant and linear solves go through Accelerate (BLAS/LAPACK) on macOS.

let A = np.array([[1, 2], [3, 4]])
let B = np.array([[5, 6], [7, 8]])

np.to_list(np.matmul(A, B))     # [[19, 22], [43, 50]]
np.dot(np.array([1, 2, 3]), np.array([4, 5, 6]))   # 32  (1-D dot product)
np.to_list(np.transpose(A))     # [[1, 3], [2, 4]]
np.det(A)                       # -2.0
np.to_list(np.inv(A))           # [[-2, 1], [1.5, -0.5]]
np.trace(A)                     # 5
np.to_list(np.diag(A))          # [1, 4]
np.to_list(np.outer(np.array([1, 2]), np.array([3, 4])))  # [[3, 4], [6, 8]]

# solve A x = b
np.to_list(np.solve(np.array([[3, 2], [1, 2]]), np.array([5, 5])))  # [0, 2.5]

Scientific examples

Complete, runnable programs. Every snippet below was executed against cobranum; the comments show the real output. A few results, visualized:

1 · Descriptive statistics & z-score normalization

Standardize each column of a data matrix to zero mean and unit variance using broadcasting.

import "cobranum.so"
let np = cobranum

let X = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
let n = np.shape(X)[0]

let mu = np.mul(np.sum(X, 0), 1.0 / n)            # column means (row vector)
let dev = np.sub(X, mu)                            # broadcast subtract
let sd  = np.sqrt(np.mul(np.sum(np.mul(dev, dev), 0), 1.0 / n))
let Z   = np.div(np.sub(X, mu), sd)                # broadcast divide

print(np.to_list(Z))
# => [[-1.2247..., -1.2247..., -1.2247...],
#     [0, 0, 0],
#     [1.2247...,  1.2247...,  1.2247...]]

2 · Linear regression — normal equation

Fit y = Xw in closed form: w = (XᵀX)⁻¹ Xᵀy. The design matrix's first column is the bias term.

let X = np.array([[1, 1], [1, 2], [1, 3], [1, 4]])   # bias + feature
let y = np.array([[2], [3], [4], [5]])               # column vector

let Xt = np.transpose(X)
let w  = np.matmul(np.matmul(np.inv(np.matmul(Xt, X)), Xt), y)

print(np.to_list(w))   # => [[1.0], [1.0]]   (intercept 1, slope 1)

3 · Linear regression — gradient descent

The iterative form, useful when the matrix is large or ill-conditioned. Each step is two matrix products and a scaled subtraction.

let X = np.array([[1, 1], [1, 2], [1, 3], [1, 4]])
let y = np.array([[2], [3], [4], [5]])
let Xt = np.transpose(X)
let n  = np.shape(X)[0]
let lr = 0.05

let w = np.zeros([2, 1])
for step in range(2000)
  let pred = np.matmul(X, w)
  let err  = np.sub(pred, y)
  let grad = np.mul(np.matmul(Xt, err), 1.0 / n)
  w = np.sub(w, np.mul(grad, lr))
end
print(np.to_list(w))   # => [[~1.0], [~1.0]]

4 · Solving linear systems

Solve A x = b with a direct LAPACK solve (LU factorization) — faster and more stable than forming the inverse.

let A = np.array([[3, 2, -1], [2, -2, 4], [-1, 0.5, -1]])
let b = np.array([1, -2, 0])
print(np.to_list(np.solve(A, b)))   # => [1, -2, -2]

5 · Dominant eigenvalue by power iteration

Repeatedly multiply by the matrix and renormalize; the vector converges to the dominant eigenvector and the Rayleigh quotient to its eigenvalue.

let M = np.array([[2, 1], [1, 3]])
let v = np.array([[1], [1]])

for i in range(60)
  let Mv = np.matmul(M, v)
  v = np.div(Mv, np.norm(Mv))        # renormalize (scalar broadcast)
end

let lambda = np.matmul(np.transpose(v), np.matmul(M, v))   # 1x1
print(np.to_list(lambda))   # => [[3.6180...]]   ((5 + sqrt 5) / 2)

6 · Numerical integration (trapezoid rule)

Estimate π from ∫₀¹ 4/(1+x²) dx = π on a fine grid — fully array-native, no loop.

let n = 100000
let x  = np.linspace(0.0, 1.0, n)
let fx = np.div(4.0, np.add(1.0, np.mul(x, x)))      # 4 / (1 + x^2)
let h  = 1.0 / (n - 1)

let area = h * (np.sum(fx) - (np.get(fx, 0) + np.get(fx, n - 1)) / 2.0)
print(area)   # => 3.14159265...

7 · Solving an ODE (Euler integration)

Integrate the harmonic oscillator y' = K y with K = [[0,1],[-1,0]]. After a quarter period the state rotates from (1, 0) to about (0, -1).

let K  = np.array([[0, 1], [-1, 0]])
let dt = 0.01
let state = np.array([[1], [0]])

for i in range(157)                                  # ~ pi/2 in time
  state = np.add(state, np.mul(np.matmul(K, state), dt))
end
print(np.to_list(state))   # => [[~0.0], [~-1.0]]

8 · Markov chain stationary distribution

Iterate a probability row vector through the transition matrix until it stops changing.

let P = np.array([[0.9, 0.1], [0.5, 0.5]])
let v = np.array([[1.0, 0.0]])                       # start in state 0

for i in range(100)
  v = np.matmul(v, P)
end
print(np.to_list(v))   # => [[0.8333..., 0.1666...]]

9 · Polynomial least-squares fit

Fit a quadratic to noisy-looking data by least squares on a Vandermonde matrix. The points lie exactly on y = x² + 1, so the recovered coefficients are [1, 0, 1].

let xs = [0, 1, 2, 3, 4]
let rows = []
for xi in xs
  push(rows, [1.0, xi * 1.0, xi * xi * 1.0])         # [1, x, x^2]
end
let V  = np.array(rows)
let y  = np.reshape(np.array([1, 2, 5, 10, 17]), [5, 1])
let Vt = np.transpose(V)
let coef = np.matmul(np.matmul(np.inv(np.matmul(Vt, V)), Vt), y)

print(np.to_list(coef))   # => [[1.0], [~0.0], [1.0]]   (1 + 0*x + 1*x^2)

10 · Covariance matrix

Center the data and form Cᵀ C / n — the basis of PCA.

let D = np.array([[1, 2], [3, 6], [5, 10], [7, 14]])
let n = np.shape(D)[0]
let mu = np.mul(np.sum(D, 0), 1.0 / n)
let Dc = np.sub(D, mu)                                # center
let cov = np.mul(np.matmul(np.transpose(Dc), Dc), 1.0 / n)
print(np.to_list(cov))   # => [[5, 10], [10, 20]]

11 · Logistic sigmoid & a forward pass

Vectorized σ(z) = 1 / (1 + e⁻ᶻ) over a whole layer, then a one-layer logistic prediction σ(Xw).

let z = np.array([-2, -1, 0, 1, 2])
print(np.to_list(np.div(1.0, np.add(1.0, np.exp(np.neg(z))))))
# => [0.1192..., 0.2689..., 0.5, 0.7310..., 0.8807...]

let X = np.array([[1, 0.5], [1, -1.0], [1, 2.0]])    # samples (with bias)
let w = np.array([[0.2], [0.8]])                     # weights
let probs = np.div(1.0, np.add(1.0, np.exp(np.neg(np.matmul(X, w)))))
print(np.to_list(probs))                             # P(class = 1) per sample

12 · Pairwise vectors & physics — projectile range

Compute the range of a projectile for many launch angles at once: R = v²·sin(2θ)/g, no explicit loop.

let g = 9.81
let v = 20.0
let deg = np.linspace(10.0, 80.0, 8)                 # launch angles
let rad = np.mul(deg, 3.141592653589793 / 180.0)
let R = np.mul(np.sin(np.mul(rad, 2.0)), v * v / g)
print(np.to_list(R))                                 # range per angle (m)
print("max range at index", np.argmax(R))            # ~ 45 degrees

API reference

Every function takes and returns NDARRAY values (or numbers). Errors are catchable Cobra errors.

group	functions
Creation	array(list) · zeros(shape) · ones(shape) · full(shape, v) · eye(n) · arange(start, stop, step?) · linspace(a, b, count)
Inspect	shape(a) · ndim(a) · size(a) · reshape(a, shape) · to_list(a) · str(a) · get(a, i[, j]) · set(a, i[, j], v)
Elementwise	add · sub · mul · div · pow · neg (arrays, scalars, or broadcast)
Unary math	sqrt · exp · log · sin · cos · abs
Fused / in-place	hypot(a,b) · fma(a,b,c) · axpy(α,x,y) · add_into · sub_into · mul_into · div_into · sqrt_into · exp_into · log_into · sin_into · cos_into · abs_into · neg_into
Reductions	sum · mean · min · max · prod (+ axis for 2-D) · var · std · norm · argmin · argmax · cumsum · clip(a, lo, hi)
Linear algebra	matmul / dot · transpose · inv · det · solve · trace · diag · outer · concat

See the downloads for the plugin and the cobra binary, and the project's bench/cobranum/ for the NumPy benchmark suite.