local linear = require("linear")

-- Calculate an inner product
local x = linear.tolinear({ 1, 2, 3 })
local y = linear.tolinear({ 3, 2, 1 })
print("x^T y", linear.dot(x, y))  -- should be 10

-- Generate and analyze normally distributed random numbers with mean 10, std 5, i.e., N(10,5)
local x = linear.vector(100)
linear.normal(x)   -- N(0,1)
linear.scal(x, 5)  -- N(0,5)
linear.inc(x, 10)  -- N(10,5)
print("N(10,5)", x[1], x[2], x[3], "...")
print("mean, std", linear.mean(x), linear.std(x, 1))  -- should approximate 10, 5

-- Solve a system of linear equations: x + 2y = 7, 2x - y = 9
local A = linear.tolinear({ { 1, 2 }, { 2, -1 } })
local B = linear.matrix(2, 1)
local b = linear.tvector(B, 1)  -- column vector
b[1], b[2] = 7, 9
linear.gesv(A, B)
print("solutions", b[1], b[2])  -- should be 5, 1

-- Calculate pairwise Pearson correlation coefficients of 10k uniform random numbers
local A = linear.matrix(10000, 3)
linear.uniform(A)
local C = linear.matrix(3, 3)
linear.corr(A, C)
print("correlations", C[1][2], C[1][3], C[2][3])  -- should be close to zero

-- Approximate the normal distribution through percentiles, and compare with its CDF
local x = linear.vector(100000)
linear.normal(x)                        -- 100k random variables from N(0,1)
local percentiles = linear.vector(101)  -- use 101 sample percentiles
linear.ranks(100, percentiles, "zq")    -- set percentile ranks, i.e., [0.00, 0.01, ..., 0.99, 1.00]
linear.quantile(x, percentiles)         -- calculate the actual sample percentiles
local r = linear.vector(#x)
linear.copy(x, r)                       -- copy values of x to r
linear.rank(percentiles, r)             -- calculate the rank of each value of x
for i = 1, 3 do
	local rank = r[i]                          -- rank within percentiles should approximate CDF
	local cdf = linear.normalcdf(x[i])         -- actual CDF
	print("cdf", x[i], rank, cdf, rank - cdf)  -- difference should be close to 0
end

-- Make a cubic spline interpolant of the Runge function, and analyze its mean absolute error
local function runge (x)
	return 1 / (1 + 25 * x^2)
end
local x, y = { }, { }
for i = -5, 5 do
	table.insert(x, i / 5)
	table.insert(y, runge(i / 5))
end
x, y = linear.tolinear(x), linear.tolinear(y)
local sp = linear.spline(x, y, "not-a-knot")
local d = { }
for i = -100, 100 do
	table.insert(d, math.abs(sp(i / 100) - runge(i / 100)))
end
print("mae", "not-a-knot", linear.mean(linear.tolinear(d)))  -- should be close to 0

-- As we know the derivative, a clamped spline should be even better
local function rungep (x)
	return (-50 * x) / (1 + 25 * x^2)^2
end
sp = linear.spline(x, y, "clamped", nil, rungep(-1), rungep(1))
d = { }
for i = -100, 100 do
	table.insert(d, math.abs(sp(i / 100) - runge(i / 100)))
end
print("mae", "clamped", linear.mean(linear.tolinear(d)))  -- should be closer to 0