------------------------------------------
------------------------------------------
----- An Implementation of           -----
----- Sparse Matrix Methods.         -----
-----                                -----
----- By John Ringland               -----
----- 2004/11/21                     -----
----- modified 2005/08/27 for the    -----
----- SOME project                   -----
------------------------------------------
------------------------------------------
-- SparseSM.e


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


global type SparseSMRow(object row)
    if sequence(row) and length(row)=3 and integer(row[rowId_]) then
       -- and length(row[elemId_])=length(row[Data_]) then
        if length(row[elemId_])=length(row[Data_]) then
            for i = 1 to length(row[elemId_]) do
                if not (integer(row[elemId_][i])and SMElement(row[Data_][i])) then 
                    return 0
                end if
            end for
        --else it maybe in the process of being constructed so let it go for now
        end if
    else
        return 0
    end if
    return 1
end type

global type SparseSMCol(object col)
    if sequence(col) and length(col)=2 and integer(col[colId_]) then
        for i = 1 to length(col[elemId_]) do
            if not integer(col[elemId_][i]) then -- no data stored in column data
                return 0
            end if
        end for
    else
        return 0
    end if
    return 1
end type

global type SparseSystemMatrix(object sm)
    if sequence(sm) and length(sm)=3
       and length(sm[dim_])=3 and integer(sm[dim_][nElems_]) -- detect sparseSM metadata
       and integer(sm[dim_][nRows_]) and integer(sm[dim_][nCols_]) then -- detect sparseSM metadata
        for i = 1 to length(sm[Rows_]) do
            if not SparseSMRow(sm[Rows_][i]) then
                return 0
            end if
        end for
        for i = 1 to length(sm[Cols_]) do
            if not SparseSMCol(sm[Cols_][i]) then
                return 0
            end if
        end for
        return 1
    end if
    return 0
end type

global type SparseSMIndex(object id)
    return sequence(id) and length(id)=2 and integer(id[row_]) and integer(id[col_])
end type

global function createSparseSM(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_,..}, {data,..} } }    <--list of Columns
	return {{ 0 , nRows, nCols},  {} , {} }
end function

function is_SM_index(SparseSystemMatrix sm, SparseSMIndex 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(SparseSystemMatrix sm, SparseSMIndex id, SMElement 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 a zero valued SMReal
    -- this is adequate because of the functioning of an SM
    -- must include this file after PNData.e for this to work
    if SMReal(data[con_]) and SMReal(data[div_]) then --atom(data) then
        if equal(data[con_][1], fdZero) and equal(data[div_][1], fdZero) then
            isNULL = 1
        else
            isNULL = 0
        end if
    else -- SMElement contains operators
        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 head
                        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 head
                        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 new head
                    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 new head
                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 head
                            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 head
                            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 new head
                        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 new head
                    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(SparseSystemMatrix sm, SparseSMIndex 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(SparseSystemMatrix 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(SparseSystemMatrix 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