[Indigo]<Cedar>PlotPackage>RealVecImpl.mesa!4

RealVecImpl.mesa, an implementation for dealing with vectors of real numbers

Last Modified On December 16, 1983 2:00 pm By Paul Rovner

DIRECTORY

PriorityQueue USING[Ref, Predict, Insert, Empty, Remove],

Real USING[FixI, RoundI],

RealVec USING[Handle, Object];

RealVecImpl: MONITOR -- protects sortee

IMPORTS PriorityQueue, Real

EXPORTS RealVec =

{ OPEN RealVec;

PUBLIC SIGNALS

LengthFault: PUBLIC ERROR[vec: Handle, expectedLength: NAT] = CODE;

PUBLIC PROCs

nElements <= h.length.

Result has averaged values as per some reasonable algorithm

Smooth: PUBLIC PROC[h: Handle, nElements: NAT] RETURNS[ans: Handle] =

{ IF nElements = h.length THEN RETURN[h];

a, b, c... are in the old space, i, j, k ... in the new

ans ← NEW[Object[nElements]];

FOR j: NAT IN [0..nElements)

x: REAL = j;

startT: REAL ← x / nElements;

endT: REAL ← (x + 1.0) / nElements;

a: NAT ← Real.FixI[startT * h.length];

b: NAT ← Real.RoundI[endT * h.length]; -- may be too small by 1

dtOld: REAL ← 1.0/h.length;

dtNew: REAL ← endT - startT;

IF b*dtOld < endT THEN b ← b + 1;

IF a = b THEN ERROR;

ans.elements[j] ← 0.0;

now add up the contributions of each old bucket that the new bucket

overlaps

FOR z: NAT IN [a..b) -- left edge index (from 0) identifies the bucket

value

DO --add in the value in the zth bucket times the fraction of the time

covered by the zth bucket in the new bucket

ans.elements[j] ← ans.elements[j]

+ h.elements[z] -- the value in the zth bucket

* (MIN[dtOld, MIN[endT, (z + 1)*dtOld]

- MAX[startT, z*dtOld]])/dtNew

ENDLOOP;

ENDLOOP};

convenient constructors

All: PUBLIC PROC[length: NAT, value: REAL ← 0.0] RETURNS[ans: Handle] =

{ans ← NEW[Object[length]];

FOR i: NAT IN [0..length) DO ans.elements[i] ← value ENDLOOP};

h[0] = 0, h[n] = n

IndexVec: PUBLIC PROC[length: NAT] RETURNS[ans: Handle] =

{ans ← NEW[Object[length]];

FOR i: NAT IN [0..length) DO ans.elements[i] ← i ENDLOOP};

convenient composition and extraction operations

SubVec: PUBLIC PROC[h: Handle, firstIndex, lastIndex: NAT]

RETURNS[ans: Handle] =

{length: NAT ← (IF h = NIL THEN 0 ELSE h.length);

IF lastIndex < firstIndex OR h = NIL OR lastIndex >= h.length

THEN ERROR LengthFault[h, length];

length ← lastIndex - firstIndex + 1;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h.elements[firstIndex + i] ENDLOOP};

ConCat: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF LOOPHOLE[LOOPHOLE[h1.length, CARDINAL]

+ LOOPHOLE[h2.length, CARDINAL], INTEGER] < 0

THEN ERROR LengthFault[NIL, 0]; -- overflow check

length ← h1.length + h2.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..h1.length)

DO ans.elements[i] ← h1.elements[i] ENDLOOP;

FOR i: NAT IN [0..h2.length)

DO ans.elements[h1.length + i] ← h2.elements[i] ENDLOOP};

h.elements[FixI[perms.elements[i]]] is put into ans.elements[i]

e.g., Permute[h, SortOrder[h]] produces a sorted version of h

Permute: PUBLIC PROC[h: Handle, perms: Handle] RETURNS[ans: Handle] =

{ length: NAT ← perms.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO x: NAT = Real.FixI[perms.elements[i]];

IF x >= h.length THEN ERROR LengthFault[h, x+1];

ans.elements[i] ← h.elements[x];

ENDLOOP;

};

sortee: Handle;

sortPred: SAFE PROC[x, y: REF ANY, data: REF ANY]

RETURNS[BOOLEAN] = TRUSTED

{i: NAT = LOOPHOLE--NARROW--[x, REF NAT]^;

j: NAT = LOOPHOLE--NARROW--[y, REF NAT]^;

RETURN[sortee.elements[i] < sortee.elements[j]]};

perms[i] is the index in h of the ith smallest element

i.e. (h.elements[perms[i]] <= h.elements[perms[i+1]])

SortOrder: PUBLIC ENTRY PROC[h: Handle] RETURNS[perms: Handle] =

{ENABLE UNWIND => NULL;

length: NAT ← h.length;

q: PriorityQueue.Ref;

sortee ← h;

q ← PriorityQueue.Predict[pred: sortPred, size: h.length];

perms ← NEW[Object[length]];

FOR i: NAT IN [0..h.length)

DO PriorityQueue.Insert[q, NEW[NAT ← i]] ENDLOOP;

FOR i: NAT IN [0..h.length)

DO perms.elements[i] ← NARROW[PriorityQueue.Remove[q], REF NAT]^

ENDLOOP;

IF NOT PriorityQueue.Empty[q] THEN ERROR;

};

element-by-element operations

Add: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];

length ← h1.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h1.elements[i] + h2.elements[i] ENDLOOP};

Subtract: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];

length ← h1.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h1.elements[i] - h2.elements[i] ENDLOOP};

Multiply: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];

length ← h1.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h1.elements[i] * h2.elements[i] ENDLOOP};

Divide: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: Handle] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];

length ← h1.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h1.elements[i] / h2.elements[i] ENDLOOP};

DotProduct: PUBLIC PROC[h1, h2: Handle] RETURNS[ans: REAL] =

{ length: NAT;

IF h1 = NIL OR h2 = NIL THEN ERROR LengthFault[NIL, 0];

IF h1.length # h2.length THEN ERROR LengthFault[h2, h1.length];

length ← h1.length;

ans ← 0.0;

FOR i: NAT IN [0..length)

DO ans ← ans + h1.elements[i] * h2.elements[i] ENDLOOP};

scalar element-by-element operations

ScalarAdd: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =

{ length: NAT;

IF h = NIL THEN RETURN[NIL];

length ← h.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h.elements[i] + s ENDLOOP};

ScalarSubtract: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =

{ length: NAT;

IF h = NIL THEN RETURN[NIL];

length ← h.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h.elements[i] - s ENDLOOP};

ScalarMultiply: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =

{ length: NAT;

IF h = NIL THEN RETURN[NIL];

length ← h.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h.elements[i] * s ENDLOOP};

ScalarDivide: PUBLIC PROC[h: Handle, s: REAL] RETURNS[ans: Handle] =

{ length: NAT;

IF h = NIL THEN RETURN[NIL];

length ← h.length;

ans ← NEW[Object[length]];

FOR i: NAT IN [0..length)

DO ans.elements[i] ← h.elements[i] / s ENDLOOP};

Equal: PUBLIC PROC[h1, h2: Handle] RETURNS[BOOLEAN] =

{ IF h1 = NIL OR h2 = NIL THEN RETURN[h1 = h2];

FOR i: NAT IN [0..h1.length)

DO IF NOT almost[h1.elements[i], h2.elements[i]]

THEN RETURN[FALSE]

ENDLOOP;

RETURN[TRUE]};

almost: PROC[p, q: REAL] RETURNS[BOOLEAN] = INLINE

{RETURN[MAX[p, q] = 0.0 OR ABS[p - q]/MAX[p, q] < 0.00001]};