```------------------------------------------
------------------------------------------
----- An Implementation of           -----
----- Sparse Matrix Methods.         -----
-----                                -----
----- By John Ringland               -----
----- 2004/11/21                     -----
------------------------------------------
------------------------------------------
-- SparseSM.e

-- some constants to make the data accessing more human readable
--
global constant dim_           = 1
global constant nElems_    = 1
global constant nRows_     = 2
global constant nCols_     = 3
-- matrix row components
global constant Rows_          = 2
--constant rowId_     = 1
global constant elemId_    = 2
global constant Data_      = 3
-- matrix column components
global constant Cols_          = 3
--constant colId_     = 1
--constant elemId_    = 2
--constant Data_      = 3
-- matrix indexing components
global constant row_           = 1
global constant col_           = 2
-- BSearch result components
global constant index_         = 1
global constant found_         = 2

-- anything below this value is treated as if it was zero
constant threshold      = 1e-10

global function sparse_SM(integer nRows, integer nCols)
-- creates an empty sparse system matrix with the following structure
-- sm = {{nElems_, nRows_,nCols_},                  <--metadata
--       { {rowId_,{ {colId_, data},.. }},.. },     <--list of Rows
--       { {colId_,{rowId_,..}},.. }}               <--list of Columns
return {{  0  , nRows, nCols},  {} , {} }
end function

function is_SM_index(sequence sm, sequence id)
-- verifies that the index id={row,col} is valid for the matrix sm
if length(id)=2 then
if integer(id[row_]) and integer(id[col_]) then
if id[row_] > 0 and id[col_] > 0 then
if id[row_] <= sm[dim_][nRows_] and id[col_] <= sm[dim_][nCols_] then
return 1
end if
end if
end if
end if
return 0
end function

global function BSearch(sequence olist, integer target)
-- This function performs a binary search looking for the integer target,
-- olist must either be an ordered list of integers or a list of sequences
-- where the first element of each sequence is an integer and the sequences are ordered.
-- It returns {index, found} where if found=0 then index indicates the element prior to
-- where the target should have been, this allows for insertion if desired. If found=1
-- then index indicates the location of the target.

integer left, right, mid, excluded
left = 1
right = length(olist)
excluded = 0
mid = left + floor((right-left) / 2)
-- if olist is a list of integers
if integer(olist[1]) then
while mid != left do
if olist[mid] >= target then
-- search left half
right = mid
elsif olist[mid+1] <= target then
-- serch right half
left = mid+1
else -- target not in list
excluded = 1
exit -- mid is element prior to insertion point
end if
mid = left + floor((right-left) / 2)
end while
-- check left then right if not yet excluded
if not excluded then
if olist[left] = target then
return {left,1}
elsif olist[right] = target then
return {right,1}
else -- not in list
if olist[left] < target then
if olist[right] > target then
return {left, 0} -- insert between left and right
else
return {right, 0} -- insert after right, e.g as new tail
end if
else
return {left-1,0} -- insert to the left of left, e.g as new head
end if
end if
else
return {mid, 0}
end if
else -- if olist is a list of sequences
while mid != left do
if olist[mid][1] >= target then
-- search left half
right = mid
elsif olist[mid+1][1] <= target then
-- serch right half
left = mid+1
else -- target not in list
excluded = 1
exit -- mid is element prior to insertion point
end if
mid = left + floor((right-left) / 2)
end while
-- check left then right if not yet excluded
if not excluded then
if olist[left][1] = target then
return {left,1}
elsif olist[right][1] = target then
return {right,1}
else -- not in list
if olist[left][1] < target then
if olist[right][1] > target then
return {left, 0} -- insert between left and right
else
return {right, 0} -- insert after right, e.g as new tail
end if
else
return {left-1,0} -- insert to the left of left, e.g as new head
end if
end if
else
return {mid, 0}
end if
end if
end function

global function set_SM(sequence sm, sequence id, object data)
-- sets the sm matrix element specified by id={row,col} to data
-- to remove a matrix element set data=0 only, don't use {0,0} etc
-- This particular implementation stores data twice by both row and column.
-- This is assuming that data is a function Id which is just an atom
-- and it also assumes that systems will require efficient access to these id's
-- by both row and column, row when acquiring input, and column when changing
-- output transforms to create non-ergodic system networks.
integer rId,rcId, cId, crId, isNULL
sequence rSearch, rcSearch, cSearch, crSearch

-- verify that the index is valid
if not is_SM_index(sm,id) then
puts(1,"The SM index was invalid: id = ")
print(1, id)
puts(1,"\n")
--abort(2)
return sm
end if

-- this function assumes that the only NULL data is an integer zero
-- this is adequate because of the functioning of an SM
if atom(data) then
if data > threshold then
isNULL = 0
else
isNULL = 1
end if
else
isNULL = 0
end if
-- find where to put the data
if length(sm[Rows_]) then
-- search for row
rSearch = BSearch(sm[Rows_], id[row_])
rId = rSearch[index_]
rcSearch = {0,0} -- for a later 'if' test when we begin the column search
-- if row exists then
if rSearch[found_] then
--  search for element
rcSearch = BSearch(sm[Rows_][rId][elemId_], id[col_])
rcId = rcSearch[index_]
-- if element exists then modify or remove it
if rcSearch[found_] then
if isNULL then -- remove element
sm[dim_][nElems_] -= 1
-- if not last element then remove from sequence
if length(sm[Rows_][rId][elemId_]) > 1 then
if rcId = 1 then
sm[Rows_][rId][elemId_] =
sm[Rows_][rId][elemId_][2..length(sm[Rows_][rId][elemId_])]
sm[Rows_][rId][Data_] =
sm[Rows_][rId][Data_][2..length(sm[Rows_][rId][Data_])]
-- if tail
elsif rcId = length(sm[Rows_][rId][elemId_]) then
sm[Rows_][rId][elemId_] =
sm[Rows_][rId][elemId_][1..length(sm[Rows_][rId][elemId_])-1]
sm[Rows_][rId][Data_] =
sm[Rows_][rId][Data_][1..length(sm[Rows_][rId][Data_])-1]
-- if body
else
sm[Rows_][rId][elemId_] =
sm[Rows_][rId][elemId_][1..rcId-1]
& sm[Rows_][rId][elemId_][rcId+1..length(sm[Rows_][rId][elemId_])]
sm[Rows_][rId][Data_] =
sm[Rows_][rId][Data_][1..rcId-1]
& sm[Rows_][rId][Data_][rcId+1..length(sm[Rows_][rId][Data_])]
end if
else -- it is last element so remove whole sequence
-- if last row then set sm[Rows_] to empty
if rId = 1 then
sm[Rows_] = sm[Rows_][2..length(sm[Rows_])]
-- if tail
elsif rId = length(sm[Rows_]) then
sm[Rows_] = sm[Rows_][1..length(sm[Rows_])-1]
-- if body
else
sm[Rows_] = sm[Rows_][1..rId-1] & sm[Rows_][rId+1..length(sm[Rows_])]
end if
end if
else -- modify element
sm[Rows_][rId][Data_][rcId] = data
end if
else -- no element so create element if necessary
if not isNULL then -- create element
sm[dim_][nElems_] += 1
if rcId = 0 then
sm[Rows_][rId][elemId_] = {id[col_]} & sm[Rows_][rId][elemId_]
sm[Rows_][rId][Data_] = {data} & sm[Rows_][rId][Data_]
-- if tail
elsif rcId = length(sm[Rows_][rId][elemId_]) then
sm[Rows_][rId][elemId_] = sm[Rows_][rId][elemId_] & {id[col_]}
sm[Rows_][rId][Data_] = sm[Rows_][rId][Data_] & {data}
-- if body
else
sm[Rows_][rId][elemId_] =
sm[Rows_][rId][elemId_][1..rcId]
& {id[col_]}
& sm[Rows_][rId][elemId_][rcId+1..length(sm[Rows_][rId][elemId_])]
sm[Rows_][rId][Data_] =
sm[Rows_][rId][Data_][1..rcId]
& {data}
& sm[Rows_][rId][Data_][rcId+1..length(sm[Rows_][rId][Data_])]
end if
end if
end if
else -- no row so create row if necessary
if not isNULL then
sm[dim_][nElems_] += 1
if rId = 0 then
sm[Rows_] = {{id[row_],{id[col_]}, {data}}} & sm[Rows_]
-- if tail
elsif rId = length(sm[Rows_]) then
sm[Rows_] = sm[Rows_] & {{id[row_],{id[col_]}, {data}}}
-- if body
else
sm[Rows_] =
sm[Rows_][1..rId]
& {{id[row_],{id[col_]}, {data}}}
& sm[Rows_][rId+1..length(sm[Rows_])]
end if
end if
end if

-- only if you created an entry or you removed an entry then
-- you chould check the columns
if not rcSearch[found_] xor isNULL then
-- search for column
cSearch = BSearch(sm[Cols_], id[col_])
cId = cSearch[index_]
-- if column exists then
if cSearch[found_] then
--  search for element
crSearch = BSearch(sm[Cols_][cId][elemId_], id[row_])
crId = crSearch[index_]
-- if element exists then modify or remove it
if crSearch[found_] then
if isNULL then -- remove element
-- if not last element then remove from sequence
if length(sm[Cols_][cId][elemId_]) > 1 then
if crId = 1 then
sm[Cols_][cId][elemId_] =
sm[Cols_][cId][elemId_][2..length(sm[Cols_][cId][elemId_])]
-- if tail
elsif crId = length(sm[Cols_][cId][elemId_]) then
sm[Cols_][cId][elemId_] =
sm[Cols_][cId][elemId_][1..length(sm[Cols_][cId][elemId_])-1]
-- if body
else
sm[Cols_][cId][elemId_] =
sm[Cols_][cId][elemId_][1..crId-1]
& sm[Cols_][cId][elemId_][crId+1..length(sm[Cols_][cId][elemId_])]
end if
else -- it is last element so remove whole sequence
-- if last column then set sm[Cols_] to empty
if cId = 1 then
sm[Cols_] = sm[Cols_][2..length(sm[Cols_])]
-- if tail
elsif cId = length(sm[Cols_]) then
sm[Cols_] = sm[Cols_][1..length(sm[Cols_])-1]
-- if body
else
sm[Cols_] = sm[Cols_][1..cId-1] & sm[Cols_][cId+1..length(sm[Cols_])]
end if
end if
end if
else -- no element so create element if necessary
if not isNULL then -- create element
if crId = 0 then
sm[Cols_][cId][elemId_] = {id[row_]} & sm[Cols_][cId][elemId_]
-- if tail
elsif crId = length(sm[Cols_][cId][elemId_]) then
sm[Cols_][cId][elemId_] = sm[Cols_][cId][elemId_] & {id[row_]}
-- if body
else
sm[Cols_][cId][elemId_] =
sm[Cols_][cId][elemId_][1..crId]
& {id[row_]}
& sm[Cols_][cId][elemId_][crId+1..length(sm[Cols_][cId][elemId_])]
end if
end if
end if
else -- no column so create column if necessary
if not isNULL then
if cId = 0 then
sm[Cols_] = {{id[col_],{id[row_]}}} & sm[Cols_]
-- if tail
elsif cId = length(sm[Cols_]) then
sm[Cols_] = sm[Cols_] & {{id[col_],{id[row_]}}}
-- if body
else
sm[Cols_] =
sm[Cols_][1..cId]
& {{id[col_],{id[row_]}}}
& sm[Cols_][cId+1..length(sm[Cols_])]
end if
end if
end if
end if
else -- matrix is empty so simply insert the data
if not isNULL then
sm[dim_][nElems_] += 1
-- sm[Rows_][rowId_] = id[row_]
-- sm[Rows_][elemId_] = {id[col_]}
-- sm[Rows_][Data_] = {data}
sm[Rows_] = {{id[row_],{id[col_]}, {data}}}
-- sm[Cols_][colId_] = id[col_]
-- sm[Cols_][elemId_] = {id[row_]}
-- sm[Cols_][Data_] = {data}
sm[Cols_] = {{id[col_],{id[row_]}}}
end if
-- else if isNULL then do nothing
end if
return sm
end function

global function get_SM(sequence sm, sequence id)
-- retrieves the element of sm specified by the id={row,col}
-- assumes the index is valid to save time during runtime processing
-- but returns zero if index is not valid.
integer rId
sequence rSearch, rcSearch, smRowsrId

-- search for row
rSearch = BSearch(sm[Rows_], id[row_])
rId = rSearch[index_]
-- if row exists then
if rSearch[found_] then
smRowsrId = sm[Rows_][rId]
--  search for element
rcSearch = BSearch(smRowsrId[elemId_], id[col_])
-- if element exists then return its data
if rcSearch[found_] then
return smRowsrId[Data_][rcSearch[index_]]
end if
end if
return 0
end function

global function get_col_SM(sequence sm, integer col)
-- retrieves the column of sm specified by col
-- assumes the index is valid to save time during runtime processing
-- but returns zero if index is not valid.
sequence cSearch, smCols

-- search for the column
smCols = sm[Cols_]
cSearch = BSearch(smCols, col)
-- if column exists then
if cSearch[found_] then
return smCols[cSearch[index_]]
end if
return 0
end function

global function get_row_SM(sequence sm, integer row)
-- retrieves the row of sm specified by row
-- assumes the index is valid to save time during runtime processing
-- but returns zero if index is not valid.
sequence rSearch, smRows

-- search for row
smRows = sm[Rows_]
rSearch = BSearch(smRows, row)
-- if row exists then
if rSearch[found_] then
return smRows[rSearch[index_]]
end if
return 0
end function
```

www.Anandavala.info