# vim: syntax=matlab
# http://people.ofset.org/~ckhung/b/ma/lademo.r

pivot = function(A, r, c, range) {
# Pivots matrix A around the element A[r;c], but only operates
# on the rows specified by range. If range is not specified,
# then operates on the entire matrix.
    if (! exist(range)) { range = [1:r-1, r+1:A.nr]; }
    A[r;] = A[r;] / A[r;c];
    for (i in range) {
	if (i == r) { error("pivot row must not occur in range"); }
	A[i;] = A[i;] - A[r;] * A[i;c];
    }
    return A;
}

epsilon = 1e-8;

r_r_echelon = function(A) {
# Computes reduced row echelon form of A
    global(epsilon);

    r = 1;
    for (c in 1:A.nc) {
	p = 0;
	for (i in r:A.nr) { if (abs(A[i;c]) > epsilon) { p = i; break } }
	if (p <= 0) { continue }
	A[r,p;] = A[p,r;];
	A = pivot(A, r, c, r+1:A.nr);
	r = r + 1;
    }
    for (r in r-1:2:-1) {
	p = 0;
	for (i in 1:A.nc) { if (abs(A[r;i]) > epsilon) { p = i; break } }
	if (p <= 0) { continue }
	A = pivot(A, r, p, 1:r-1);
    }
    return A;
}

rand_permute = function(v) {
# Returns a random permutation of the vector v (or of [1,2,...v] if
# v is a scalar) Can be useful too, if you need random combinations.
# (Just pick the first k elements)
    if (v.nr * v.nc == 1) { v = 1:v; }
    n = v.nr * v.nc;
    for (k in n:2:-1) {
	t = floor(rand() * k) + 1;
	v[t,k] = v[k,t];
    }
    return v;
}

rand_mat = function(m, n, r) {
# Returns a random matrix of size m*n and of rank r
    if (m > n) { error("matrix width must be >= matrix height"); }
    A = zeros(m,n);
    p = rand_permute(n);
    d = [sort(p[1:r]).val, inf()];
    i = 1;
    for (j in 1:n) {
	if (j == d[i]) {
	    A[i;j] = 1;
	    i++;
	else
	    A[1:i-1;j] = rand(i-1,1)*2 - 1;
	}
    }
    return (rand(m,m)*2-1) * A;
}

simplex = function(obj, A) {
    global (epsilon);
    T = [A; -obj];
    T = [T[;1:T.nc-1], eye(T.nr,T.nr), T[;T.nc]];
    nr = T.nr - 1;
    nc = T.nc - 1;
    while (1) {
	bad = mini(T[1:nr;nc+1]);
	if (T[bad;nc+1] >= 0) { break; }
	pc = mini(T[bad;1:nc]);
	if (T[bad;pc] >= 0) {
	    writem("stderr", T);
	    error("no solution");
	}
	ratio = T[1:nr;nc+1] ./ T[1:nr;pc];
	pr = mini(ratio + (ratio <= 0) / epsilon);
	T = pivot(T, pr, pc);
    }
    while (1) {
	pc = mini(T[nr+1;1:nc]);
	if (T[nr+1;pc] >= 0) { break; }
	ratio = T[1:nr;nc+1] ./ T[1:nr;pc];
	bad = T[1:nr;pc]<=0;
	pr = mini(ratio + bad / epsilon);
	if (bad[pr]) { break; }
	T = pivot(T, pr, pc);
    }
    return T;
}

smooth = function(A) {
# Finds all entries in A that are close to 0 and makes them 0
# (in order to get a cleaner output)
    global(epsilon);

    return A .* (abs(A) > epsilon) + 0;
}

# test
# srand(7)
# T = rand_mat(5,6,4) * 10
# pivot(T,3,4,[2,5])
# A = r_r_echelon(T)
# smooth(A)

# simplex([-2,-3,-375], [1,0,30; 0,1,25; -1,-1,-15;1,1,45])
# simplex([13,4,0], [-2,-1,-11; 1,1,10; 1/3,1,6; -1/4,-1,-4])