Kiwis - von Hamstern und Früchten

Die Welt scheint keine Fourier-Motzkin-Elimination zu brauchen. Zumindest habe ich keine brauchbare Implementierung in den Weiten des Internets, die google abdeckt, gefunden.

Da ich aber als Hausaufgabe die Fourier-Motzkin-Elimination in Matlab programmieren soll, stelle ich sie allen, die sonst so danach suchen, zur Verfügung:

%% Fourier-Motzkin-Elimination %%
function[D, d] = fourmotz(A,b,j)

% Input: c = e_j: die Projektionsrichtung
% P(A,b): ein Polyeder
% H = {x in K^n | x_j = 0} = {x in K^n | c^T x = 0}
% Output: P(D, d), so dass {x in K^n | Dx <= d, x_j = 0}
% die orthogonale Projektion von P(A, b) auf H ist

[m,n] = size(A);
e = zeros(n,1);
e(j) = 1;

% (1) Partitioniere die Zeilenindexmenge M = {1, 2, . . . ,m} von A wie
% folgt:
% N = {i in M | A_ij < 0},
% Z = {i in M | A_ij = 0},
% P = {i in M | A_ij > 0}.

vneg = (A(:,j) < 0); % negative Einträge in der j.ten Spalte
vzer = (A(:,j) == 0); % null Einträge in der j.ten Spalte
vpos = (A(:,j) > 0); % positive Einträge in der j.ten Spalte

x = 1;
y = 1;

for i = 1:m
for j = 1:m
if (vneg(i) == 1 && vpos(j) == 1)
negpos(x,1) = i;
negpos(x,2) = j;
x = x+1;
end
end

if vzer(i) == 1
zero(y) = i;
y = y+1;
end
end

neg = sum(vneg);
zer = sum(vzer);
pos = sum(vpos);

% (2) Setze r = |Z v (N × P)| und sei p : {1, 2, . . . , r} -> Z v (N × P)
% eine Bijektion.

r = zer + (neg * pos);

%(3) Für i = 1, 2, . . . , r führe aus:

for i = 1:zer
% (a) Ist p(i) in Z, dann setze D_i· = A_p(i)· , d_i = b_p(i)
D(i,:) = A(zero(i),:);
d(i) = b(zero(i));
end

for i = (zer+1):r,
% (b) Ist p(i) = (s, t) in N × P, dann setze
% Di· = (At· c)As· - (As· c)At·,
% di = (At· c)bs - (As· c)bt.
s = negpos(i-zer,1)
t = negpos(i-zer,2)

D(i,:) = ((A(t,:) * e) * A(s,:)) - ((A(s,:) * e) * A(t,:));
d(i) = ((A(t,:) * e) * b(s)) - ((A(s,:) * e) * b(t));

end