#!/usr/bin/octave -q
# http://people.ofset.org/~ckhung/b/ma/lademo.m

1;

function rv = pivot(A, r, c, varargin)
# Pivots matrix A around the element A(r,c), but only operates
# on the rows specified by range. If range is an illegal index
# (such as inf, nan, or 0), then operates on the entire matrix.
    s = size(A);
    # see http://wiki.octave.org/wiki.pl?VariableLengthArgumentLists
    if (nargin > 3)
	range = varargin{1};
    else
	range = [1:r-1,r+1:s(1)];
    endif
    A(r,:) = A(r,:) / A(r,c);
    for i = range
	if (i == r); error("pivot row must not occur in range"); endif
	A(i,:) = A(i,:) - A(r,:) * A(i,c);
    endfor;
    rv = A;
endfunction

global epsilon;
epsilon = 1e-5;

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

    s = size(A);
    r = 1;
    for c = 1:s(2)
	p = 0;
	for i = r:s(1)
	    if (abs(A(i,c)) > epsilon); p = i; break; endif
	endfor
	if (p <= 0); continue; endif
	A([r,p],:) = A([p,r],:);
	A = pivot(A, r, c, r+1:s(1));
	r = r + 1;
    endfor
    for r = r-1:-1:2
	p = 0;
	for i = 1:s(2)
	    if (abs(A(r,i)) > epsilon); p = i; break; endif
	endfor
	if (p <= 0); continue; endif
	A = pivot(A, r, p, 1:r-1);
    endfor
    rv = A;
endfunction

function rv = rand_permute(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)
    s = size(v);
    n = s(1) * s(2)
    if (n == 1); n = v; v = 1:v; endif
    for k = n:-1:2
	t = floor(rand() * k) + 1;
	v([t,k]) = v([k,t]);
    endfor
    rv = v;
endfunction

function rv = rand_mat(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"); endif
    A = zeros(m,n);
    p = rand_permute(n);
    d = [sort(p(1:r)), inf];
    i = 1;
    for j = 1:n
	if (j == d(i))
	    A(i,j) = 1;
	    ++i;
	else
	    A(1:i-1,j) = rand(i-1,1)*2 - 1;
	endif
    endfor
    rv = (rand(m,m)*2-1) * A;
endfunction

function rv = simplex(obj, A)
    global epsilon;
    T = [A; -obj];
    s = size(T);
    T = [T(:,1:s(2)-1), eye(s(1),s(1)), T(:,s(2))];
    s = size(T);
    nr = s(1) - 1;
    nc = s(2) - 1;
    while (1)
	bad = mini(T(1:nr,nc+1));
	if (T(bad,nc+1) >= 0); break; endif
	pc = mini(T(bad,1:nc));
	if (T(bad,pc) >= 0)
	    T
	    error("no solution");
	endif
	ratio = T(1:nr,nc+1) ./ T(1:nr,pc);
	pr = mini(ratio + (ratio <= 0) / epsilon);
	T = pivot(T, pr, pc);
    endwhile
    while (1)
	pc = mini(T(nr+1,1:nc));
	if (T(nr+1,pc) >= 0); break; endif
	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; endif
	T = pivot(T, pr, pc);
    endwhile
    rv = T;
endfunction

function rv = mini(c)
    rv = 1;
    for i = 2:prod(size(c))
	if (c(i) < c(rv))
	    rv = i;
	endif
    endfor
endfunction

# test
rand("seed", 7);
T = rand_mat(5,6,4) * 10
pivot(T,3,4,[2,5])
A = r_r_echelon(T)

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])