File : lass.adb
with Messages, Reps.Ops, Lin.Ops, Ada.Unchecked_Deallocation, Store;
use Messages, Reps.Ops, Lin.Ops;
package body Lass is
pragma Optimize (Time);
package LassStorage is new Store(Index => L_Set);
use LassStorage;
Small: constant Rep := 1.0E-6;
Tiny: constant Rep := 1.0E-10;
NegSmall: constant Rep := -Small;
NegTiny: constant Rep := -Tiny;
MaxTrans: constant Rep := 1.0E+3; -- shouldn't have to translate more than that
Table_Size: Int renames Prim;
Volumes_Stored: Int renames Filled;
Volumes_Retrieved: Int2 := 0;
procedure Free is new Ada.Unchecked_Deallocation(SSV,SSVPtr);
procedure Free is new Ada.Unchecked_Deallocation(LSV,LSVPtr);
procedure Free is new Ada.Unchecked_Deallocation(Qivot,QivotPtr);
procedure Free is new Ada.Unchecked_Deallocation(Pivots,PivotsPtr);
procedure Free is new Ada.Unchecked_Deallocation(Polts,PoltsPtr);
procedure Free is new Ada.Unchecked_Deallocation(Goo,GooPtr);
procedure AllocGStuff is
begin
if GPol = null then GPol := new Polts; end if;
if GQiv = null then GQiv := new Qivot; end if;
if GPiv = null then GPiv := new Pivots; end if;
if GNum = null then GNum := new Pivots; end if;
if GCol = null then GCol := new SSV; end if;
if GRow = null then GRow := new LSV; end if;
if GMat = null then GMat := new Goo; end if;
end AllocGStuff;
procedure FreeGStuff is
begin
if GQiv /= null then Free(GQiv); end if;
if GPiv /= null then Free(GPiv); end if;
if GNum /= null then Free(GNum); end if;
if GPol /= null then Free(GPol); end if;
if GCol /= null then Free(GCol); end if;
if GRow /= null then Free(GRow); end if;
if GMat /= null then Free(GMat); end if;
end FreeGStuff;
procedure StorageInfo is
L: constant Integer := Integer(Rep(100)*Rep(Volumes_Stored)/Rep(2*Table_Size));
begin
Show("storage level: " & L'Img & " %");
Show("volumes stored: " & Volumes_Stored'Img);
Show("volumes retrieved: " & Volumes_Retrieved'Img);
end StorageInfo;
function NearPos(R: Rep) return Boolean is
begin
return R>Small;
end NearPos;
function NearZero(R: Rep) return Boolean is
begin
return R <= Small and then R >= NegSmall;
end NearZero;
function Cleanup(R: Rep) return Rep is
begin
if R >= Tiny or else R <= NegTiny then return R;
else return Zero; end if;
end Cleanup;
pragma Inline (NearPos,NearZero,Cleanup);
function FillPivot(M,D: Integer; P: PivotPtr; A: PoltPtr) return Boolean is
Jp: Integer;
R: Rep;
begin
for I in 1 .. M loop
declare
Ai: constant FacePtr := A(I)'Access;
begin
Jp := 0;
R := Zero;
for J in 1 .. D loop
if abs(Ai(J))>R then
Jp := J;
R := abs(Ai(J));
end if;
end loop;
if R=Zero then
if Ai(0)<Zero then return False; end if; -- no solution
P(I) := 0;
else
P(I) := Jp;
R := (One-Small)*R;
for N in 1 .. I-1 loop -- compare with earlier planes
Jp := P(N);
if Jp>0 and then abs(Ai(Jp))>R then
declare
An: constant FacePtr := A(N)'Access;
Rn: constant Rep := An(Jp);
Ri: constant Rep := Ai(Jp);
begin
for J in 1 .. D loop
if not NearZero(Ai(J)/Ri-An(J)/Rn) then goto OK; end if;
end loop;
if Ri>Zero then
if Rn>Zero then
if (Ai(0)/Ri)<(An(0)/Rn) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(Ai(0)/Ri-An(0)/Rn) then
return False;
end if;
elsif Rn<Zero then
if (Ai(0)/Ri)>(An(0)/Rn) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(An(0)/Rn-Ai(0)/Ri) then
return False;
end if;
<<OK>> null;
end;
end if;
end loop;
end if;
end;
end loop;
return True;
end FillPivot;
function Ortho(F,T: in FacePtr; H: in PoltPtr; M,D: in Integer) return Integer is
R: Rep;
Hn: constant FacePtr := H(M+1)'Access;
begin
Hn(1 .. D) := F(1 .. D);
for I in 1 .. M loop -- find orthogonal complement of H(M+1) to H(1 .. M)
declare
Hi: constant FacePtr := H(I)'Access;
begin
R := Zero;
for J in 1 .. D loop -- dot product
R := R+Hn(J)*Hi(J);
end loop;
R := R/Hi(0);
for J in 1 .. D loop
Hn(J) := Cleanup(Hn(J)-R*Hi(J));
end loop;
end;
end loop;
R := Zero;
for J in 1 .. D loop -- norm^2
R := R+Sqr(Hn(J));
end loop;
if Sqrt(R) <= abs(F(0))/MaxTrans then
return 0; -- H(M+1) close to span H(1) .. H(M)
end if;
Hn(0) := R;
if F(0)=Zero then -- no new translation needed
if M=0 then
for J in 1 .. D loop T(J) := Zero; end loop;
end if;
return 1;
end if;
R := F(0)/R;
if M=0 then
for J in 1 .. D loop
T(J) := R*Hn(J);
end loop;
else
for J in 1 .. D loop
T(J) := T(J)+R*Hn(J);
end loop;
end if;
return 1;
end Ortho;
procedure Translate(M,D: in Integer; P: in PivotPtr; A: in PoltPtr) is
N: Integer := 0;
R: Rep;
H: constant PoltPtr := GPol(1)'Access;
T: constant FacePtr := H(NumEqu)'Access; -- accumulates translations
begin
for I in 1 .. M loop
if P(I) /= 0 and then A(I)(0)=Zero then
N := N+1;
P(I) := -1;
end if;
end loop;
if N>0 then
if N >= D then return; end if;
N := 0;
for I in 1 .. M loop
if P(I)<0 then N := N+Ortho(A(I)'Access,T,H,N,D); end if;
end loop;
end if;
for I in 1 .. M loop
if P(I)>0 then
N := N+Ortho(A(I)'Access,T,H,N,D);
exit when N=D;
end if;
end loop;
for I in 1 .. M loop
if P(I)<0 then
A(I)(0) := Zero;
elsif P(I)>0 then
declare
Ai: constant FacePtr := A(I)'Access;
begin
R := Zero;
for J in 1 .. D loop
R := R+Ai(J)*T(J);
end loop;
Ai(0) := Cleanup(Ai(0)-R);
end;
end if;
end loop;
end Translate;
procedure Reduce(M,D,I: in Integer; P: in PivotPtr; A,B: in PoltPtr; Ib: out Integer) is
Save: constant Boolean := D>MinVolSave;
Ai: constant FacePtr := A(I)'Access;
Jp: constant Integer := P(I);
begin
Ib := 0;
for Ia in 1 .. M loop
if P(Ia) /= 0 and then Ia /= I then
Ib := Ib+1;
if Save then
GNum(D-1)(Ib) := GNum(D)(Ia); -- copy halfspace labels
end if;
declare
Aia: constant FacePtr := A(Ia)'Access;
Bib: constant FacePtr := B(Ib)'Access;
C: constant Rep := Aia(Jp)/Ai(Jp);
begin
for J in 0 .. Jp-1 loop
Bib(J) := Cleanup(Aia(J)-C*Ai(J));
end loop;
for J in Jp+1 .. D loop
Bib(J-1) := Cleanup(Aia(J)-C*Ai(J));
end loop;
end;
end if;
end loop;
end Reduce;
function Vol1(M: Integer) return Rep is
A: constant PoltPtr := GPol(1)'Access;
R: Rep := Rep'Last;
L: Rep := -R;
Q: Rep;
begin
for I in 1 .. M loop
Q := A(I)(1);
if Q>Tiny then
Q := A(I)(0)/Q;
if Q<R then
R := Q;
if R <= L+Tiny then return Zero; end if;
end if;
elsif Q<NegTiny then
Q := A(I)(0)/Q;
if Q>L then
L := Q;
if R <= L+Tiny then return Zero; end if;
end if;
elsif A(I)(0)<NegSmall then
return Zero;
end if;
end loop;
return R-L;
end Vol1;
function Vol2(M: Integer) return Rep is
A: constant PoltPtr := GPol(2)'Access;
P: constant PivotPtr := GPiv(2)'Access;
Pi1,Pi2: Boolean;
V: Rep := Zero;
Q,L,R: Rep;
begin
for I in 1 .. M loop -- fill pivot and check parallel planes
declare
Ai1: constant Rep := A(I)(1);
Ai2: constant Rep := A(I)(2);
begin
if Ai1=Zero and then Ai2=Zero then
if A(I)(0)<Zero then return Zero; end if; -- no solution
P(I) := 0;
else
if abs(Ai2)<abs(Ai1) then
P(I) := 1;
Pi1 := True;
Pi2 := abs(Ai2/Ai1)>(One-Small);
else
P(I) := 2;
Pi2 := True;
Pi1 := abs(Ai1/Ai2)>(One-Small);
end if;
for N in 1 .. I-1 loop -- compare with earlier planes
if P(N)=1 and then Pi1 then
declare
An1: constant Rep := A(N)(1);
begin
if NearZero(Ai2/Ai1-A(N)(2)/An1) then
if Ai1>Zero then
if An1>Zero then
if (A(I)(0)/Ai1)<(A(N)(0)/An1) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(A(I)(0)/Ai1-A(N)(0)/An1) then
return Zero;
end if;
elsif An1<Zero then
if (A(I)(0)/Ai1)>(A(N)(0)/An1) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(A(N)(0)/An1-A(I)(0)/Ai1) then
return Zero;
end if;
end if;
end;
elsif P(N)=2 and then Pi2 then
declare
An2: constant Rep := A(N)(2);
begin
if NearZero(Ai1/Ai2-A(N)(1)/An2) then
if Ai2>Zero then
if An2>Zero then
if (A(I)(0)/Ai2)<(A(N)(0)/An2) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(A(I)(0)/Ai2-A(N)(0)/An2) then
return Zero;
end if;
elsif An2<Zero then
if (A(I)(0)/Ai2)>(A(N)(0)/An2) then P(N) := 0; else P(I) := 0; end if;
elsif not NearPos(A(N)(0)/An2-A(I)(0)/Ai2) then
return Zero;
end if;
end if;
end;
end if;
end loop;
end if;
end;
end loop;
for I in 1 .. M loop -- loop over edges
declare
Ai: constant FacePtr := A(I)'Access;
Ai0: constant Rep := Ai(0);
begin
if Ai0 /= Zero then
if P(I)=1 then
R := Rep'Last;
L := -Rep'Last;
for Ia in 1 .. M loop -- Vol1
if P(Ia)>0 and then Ia /= I then
declare
Aia: constant FacePtr := A(Ia)'Access;
C: constant Rep := Aia(1)/Ai(1);
begin
Q := Aia(2)-C*Ai(2);
if Q>Tiny then
Q := (Aia(0)-C*Ai0)/Q;
if Q<R then
R := Q;
if R <= L+Tiny then goto ZeroVol; end if;
end if;
elsif Q<NegTiny then
Q := (Aia(0)-C*Ai0)/Q;
if Q>L then
L := Q;
if R <= L+Tiny then goto ZeroVol; end if;
end if;
elsif (Aia(0)-C*Ai0) < NegSmall then
goto ZeroVol;
end if;
end;
end if;
end loop;
V := V+(R-L)*Ai0/abs(Ai(1));
elsif P(I)=2 then
R := Rep'Last;
L := -Rep'Last;
for Ia in 1 .. M loop -- Vol1
if P(Ia)>0 and then Ia /= I then
declare
Aia: constant FacePtr := A(Ia)'Access;
C: constant Rep := Aia(2)/Ai(2);
begin
Q := Aia(1)-C*Ai(1);
if Q>Tiny then
Q := (Aia(0)-C*Ai0)/Q;
if Q<R then
R := Q;
if R <= L+Tiny then goto ZeroVol; end if;
end if;
elsif Q<NegTiny then
Q := (Aia(0)-C*Ai0)/Q;
if Q>L then
L := Q;
if R <= L+Tiny then goto ZeroVol; end if;
end if;
elsif (Aia(0)-C*Ai(0)) < NegSmall then
goto ZeroVol;
end if;
end;
end if;
end loop;
V := V+(R-L)*Ai0/abs(Ai(2));
end if;
end if;
<<ZeroVol>> null;
end;
end loop;
return Cleanup(V/Two);
end Vol2;
function VolFactor(VarSet: Int) return Rep is
Nw,J: Integer;
R,F: Rep;
VarNum: IVec(0 .. NumVar);
begin
UnPack(VarSet,VarNum,Nw);
R := One;
for Iw in 1 .. Nw loop
declare
D: constant Integer := NumVar-Iw+1;
I: constant Integer := GQiv(D);
Jp: constant Integer := GPiv(D)(I);
Ai: constant FacePtr := GPol(D)(I)'Access;
begin
R := R*Ai(Jp);
for Jw in 1 .. Nw loop
J := VarNum(Jw);
if J=0 then
GMat(Iw,Jw) := Zero;
else
GMat(Iw,Jw) := Ai(J);
if J>Jp then
VarNum(Jw) := J-1;
elsif J=Jp then
VarNum(Jw) := 0;
end if;
end if;
end loop;
end;
end loop;
Det(Nw,GMat.all,F);
return abs(R/F);
end VolFactor;
procedure FindVol(D: in Integer; V: out Rep; Found: out Boolean) is
I: Int;
begin
FindVol(GRow(D),I,V,Found);
if not Found then return; end if;
if V=Zero or else I=GCol(D) then
Volumes_Retrieved := Volumes_Retrieved+1;
elsif D<MinVolConv then
Found := False;
else
V := V*VolFactor(I);
Volumes_Retrieved := Volumes_Retrieved+1;
end if;
end FindVol;
function Vol(Ma,Da: in Integer) return Rep is
Found: Boolean;
Va: Rep;
begin
if Da=2 then return Vol2(Ma); end if;
if Da >= MinVolSave and then Da <= MaxVolSave then
FindVol(Da,Va,Found);
if Found then return Va; end if;
end if;
declare -- dimension d>2
A: constant PoltPtr := GPol(Da)'Access;
P: constant PivotPtr := GPiv(Da)'Access;
begin
if not FillPivot(Ma,Da,P,A) then return Zero; end if;
if Da >= MinDoTrans then Translate(Ma,Da,P,A); end if;
declare
Db: constant Integer := Da-1;
B: constant PoltPtr := GPol(Db)'Access;
Mb: Integer;
Vb: Rep;
begin
Va := Zero;
for I in 1 .. Ma loop
if A(I)(0) /= Zero and then P(I)>0 then
Reduce(Ma,Da,I,P,A,B,Mb);
if Da >= MinVolSave then
GCol(Db) := GCol(Da);
FillGap(P(I),GCol(Db));
GRow(Db) := GRow(Da);
AddElement(GNum(Da)(I),GRow(Db));
GQiv(Da) := I;
end if;
Vb := Vol(Mb,Db);
Va := Va+Vb*A(I)(0)/abs(A(I)(P(I)));
if Da=NumVar and then Verbose>0 then
Show("I =" & I'Img & " V =",Va/Rep(Da));
if Verbose>1 then StorageInfo; end if;
end if;
end if;
end loop;
end;
Va := Cleanup(Va/Rep(Da));
if Da >= MinVolSave and then Da <= MaxVolSave then
SaveVol(GRow(Da),GCol(Da),Va);
end if;
return Va;
end;
end Vol;
procedure Vol(A: in out RMatPtr; V: out Rep) is
-- constraints are a(i,1)*x_1 + ... + a(i,d)*x_d <= a(i,0)
begin
AllocGStuff;
for I in 1 .. NumEqu loop
GNum(NumVar)(I) := I; -- halfspace labels
for J in 0 .. NumVar loop
GPol(NumVar)(I)(J) := Cleanup(A(I,J));
end loop;
end loop;
Free(A);
GCol(NumVar) := Empty;
AddElement(0,GCol(NumVar));
GRow(NumVar) := Empty;
AddElement(0,GRow(NumVar));
Make_Empty_Table;
V := Vol(NumEqu,NumVar);
Free_Table;
FreeGStuff;
end Vol;
end Lass;