Rob writes:
>Any thought on the matter of Object!s potentially being self modifying???
This happens all the time. I don't think it's controversial, is it?
One nice thing that objects do is allow you to use a lot of persistent local
variables.
Here's an implementation I did of an algorithm from Numerical Recipes in C.
Don't ask me how the algorithm works - it's too complicated for me!
There are two functions in this object, DECOMP and BACKSUB. The purpose is
to find solutions for various values of B for the equation A * X = B, where
A is a square matrix, and X and B are vectors. DECOMP has to be run once
for every value of A, and then BACKSUB needs to be run to solve for X given B.
lu: make object! [
indx: ; holds info on permutations
lu: ; holds LU decomposed matrix produced by DECOMP
x: ; holds solution to A * X = B given by BACKSUB
d: ; holds permutation parity
n: ; holds size of square matrix
none
decomp: func [
{do an LU decomposition of square matrix A - from NR in C p. 46-7}
a [block!]
/local i imax j k big dum sum vv
][
lu: copy/deep a
n: length? a
vv: head insert/dup copy [] none n
indx: copy vv
d: 1 ; ends up -1 for odd number of row changes
i: 1 while [ i <= n ] [ ; loop over rows to get implicit scaling
big: 0
j: 1 while [ j <= n ] [
big: max abs lu/:i/:j big
j: j + 1
]
if zero? big [make error! "singular matrix in LU-DECOMP"]
poke vv i 1.0 / big ; save scaling
i: i + 1
]
j: 1 while [ j <= n ] [ ; go column by column
i: 1 while [ i < j ] [ ; do rows above the diagonal
sum: lu/:i/:j
k: 1 while [ k < i ] [
sum: sum - (lu/:i/:k * lu/:k/:j)
k: k + 1
]
poke lu/:i j sum
i: i + 1
]
big: 0
i: j while [ i <= n ] [ ; find largest pivot on or below diagonal
sum: lu/:i/:j
k: 1 while [ k < j ] [
sum: sum - (lu/:i/:k * lu/:k/:j)
k: k + 1
]
poke lu/:i j sum
if big <= dum: vv/:i * abs sum [
big: dum ; set big to largest pivot relative to its row
imax: i
]
i: i + 1
]
if j <> imax [
swap-items lu imax j ; swap rows
d: (- d)
poke vv imax vv/:j ; swap scaling factors
]
poke indx j imax ; save swap info for BACKSUB
if zero? lu/:j/:j [ poke lu/:j j 1e-20 ]
if j <> n [
dum: 1.0 / lu/:j/:j ; divide by pivot element
i: j + 1 while [ i <= n ] [
poke lu/:i j lu/:i/:j * dum
i: i + 1
]
]
j: j + 1
]
lu
]
backsub: func [
{perform back substution on LU to solve for X given B}
b [block!] "vector to solve"
/local a i ii ip j sum
][
if any [not lu not indx] [make error! "no LU yet to solve with"]
x: copy b
ii: 0
i: 1 while [ i <= n ] [
ip: indx/:i
sum: x/:ip
poke x ip x/:i
either zero? ii [
if not zero? sum [ii: i]
][
j: ii while [ j <= (i - 1) ] [
sum: sum - (lu/:i/:j * x/:j)
j: j + 1
]
]
poke x i sum
i: i + 1
]
i: n while [ i >= 1 ] [
sum: x/:i
j: i + 1 while [ j <= n ] [
sum: sum - (lu/:i/:j * x/:j)
j: j + 1
]
poke x i sum / lu/:i/:i
i: i - 1
]
x
]
]
So, the object LU starts out life like this:
>> source lu
lu:
make object! [
indx: none
lu: none
x: none
d: none
n: none
.... ; function code snipped
]
Now let's:
>> a: [[1 2 3] [4 5 6] [3 6 6]]
== [[1 2 3] [4 5 6] [3 6 6]]
>> lu/decomp a
== [[4 5 6] [0.75 2.25 1.5] [0.25 0.333333333333333 1]]
and take a look at LU again:
>> source lu
lu:
make object! [
indx: [2 3 3]
lu: [[4 5 6] [0.75 2.25 1.5] [0.25 0.333333333333333 1]]
x: none
d: 1
n: 3
....
]
In a sense LU has modified itself. LU/LU, LU/INDX and LU/N are all used by
LU/BACKSUB in its calculations. Now let's:
>> lu/backsub [7 6 3]
== [-1.66666666666667 -4.66666666666667 6]
Hmmm, did that work?
>> lu/x
== [-1.66666666666667 -4.66666666666667 6]
>> multiply-matrices a lu/x
== [[7] [6] [3]]
Guess so. All this is still really rough around the edges, though.
Anyway, what I'm trying to illustrate here, is that objects can be used as a
kind of workbench when you've got a problem that is convenient to do in
stages, and you need to hang on to some intermediate values.
LU/DECOMP can't tell the difference between INDX and LU as words local to its
object and global words. In fact LU/DECOMP doesn't even realize it's part of
an object. But from our bird's-eye perspective we can view the environment
that LU/DECOMP functions in, and watch how LU/DECOMP modifies that
environment. I think it's pretty damn cool to be able to view that whole
local environment with SOURCE.
After implementing an algorithm like this, it's a good idea to do
>> huh * data? ; see %huh.r on www.rebol.org
to make sure you don't have any words "leaking out" into the global context.
Don't know if this is anything like what you were asking about, but it
seemed like a good excuse to make a couple of comments on objects that I
had stored up!
See you,
Eric
