------------------------------------------
------------------------------------------
----- A System Matrix Implementation -----
----- using Sparse Matrix Methods,   -----
----- and Energy Flow Processing     -----
-----                                -----
----- John Ringland                  -----
----- 2004/11/21                     -----
------------------------------------------
------------------------------------------

-- this file uses sparse matrix energy flow methods
-- to implement an ergodic pure metric sysWrapped SMN framework
-- using System Transforms.
-- Sparse matrix elements are either scalar multipliers or they are function Id's
-- and these are stored in both row and column entries of the sparse matrix.

with trace
include misc.e -- for the pretty_print function
include SparseSM.e -- for sparse system matrices

--without type_check
--with profile_time
--with profile

-- some constants to make the data accessing more human readable
--
-- System Model components
global constant model_         = 1
    global constant ST_        = 1 -- operator list
    global constant SM_        = 2 -- system matrix
    global constant SV_        = 3 -- state vector
    global constant CU_        = 4 -- update list
    --global constant MD_      = 4
-- matrix element components
global constant con_    = 1 -- converging
global constant div_    = 2 -- diverging


function union_SET(sequence set1, sequence set2)
    -- set1 and set2 are ordered lists of atoms with no duplicates
    -- this function performs the operation of the Union of the two sets by adding the 
    -- elements of the smallest set to to the largest set if they are not already present.
    -- It takes advantage of the ordering of the sets to optimise the procedure.
    sequence result, set, data
    integer rId, slen1, slen2, slen, dlen
    rId = 1
    slen1 = length(set1)
    slen2 = length(set2)
    if slen1 < slen2 then
        -- so that we iterate over the smallest set for efficiency
        data = set1
        set = set2
        slen = slen2
        dlen = slen1
    else
        data = set2
        set = set1
        slen = slen1
        dlen = slen2
    end if
    for i = 1 to dlen do
        -- uses rId so we search only that part of the set that the datum may be in
        -- this is possible because of the ordering of both sets.
        if rId > slen then
            set = set & data[i..dlen]
            exit
        end if
        result = BSearch(set[rId..slen],data[i])
        rId += result[index_]-1
        if not result[found_] then -- insert data[i]
            -- if new head
            if rId = 0 then
                set = {data[i]} & set
                slen += 1
                rId += 2
            -- if tail
            elsif rId = slen then
                -- because data is ordered the remainder can be appended
                set = set & data[i..dlen]
                exit
            -- if body
            else
                set = set[1..rId] & {data[i]} & set[rId+1..slen]
                slen += 1
                rId += 2
            end if
        else -- data[i] already in set
            rId += 1
        end if
    end for
	return set
end function

global function step_Model(sequence mdl, integer nSteps)
    -- steps the model forwards by nSteps
    sequence st, sm, sv, rUpdate, cUpdate, row, rowData, new_sv
    object acc, col, data, rUr, interfaces, output, input
    
    st = mdl[ST_]
    sm = mdl[SM_]
    sv = mdl[SV_]
    new_sv = sv
    cUpdate = mdl[CU_]
    -- for each moment in time
    for n = 1 to nSteps do
        rUpdate = {}
        -- cUpdate is an ordered list of colId’s which require updating
        -- for each element of cUpdate
        for c = 1 to length(cUpdate) do
            -- find the matrix column
            col = get_col_SM(sm,cUpdate[c])
            -- if the column existed in the matrix; for global outputs col=0
            if sequence(col) then
                col = col[elemId_] -- to keep only the rowIds
                -- form the union of each column with rUpdate
                rUpdate = union_SET(rUpdate,col)
            end if
        end for
        cUpdate = {} -- empty cUpdate now that it has been used
        -- rUpdate is now an ordered list of rowId’s which require updating
        -- for each element of rUpdate
        for r = 1 to length(rUpdate) do
            -- retrieve the Row from the sm
            rUr = rUpdate[r]
            row = get_row_SM(sm, rUr)
            acc = {}    
            -- for each matrix element in that row
            rowData= row[Data_]
            for e = 1 to length(rowData) do
                interfaces = rowData[e]
                if sequence(interfaces) then
                    -- process the interface function pipelines
                    -- separate input from output pipelines
                    output = interfaces[div_]
                    input = interfaces[con_]
                    if sequence(output) then
                        data = sv[row[elemId_][e]]
                        for f = length(output) to 1 by -1 do
                            data = call_func(output[f], {data})
                        end for
                    else
                        data = output * sv[row[elemId_][e]]
                    end if
                    if sequence(input) then
                        for f = length(input) to 1 by -1 do
                            data = call_func(input[f], {data})
                        end for
                    else
                        data = input * data
                    end if
                else
                    data = interfaces * sv[row[elemId_][e]]
                end if
                acc = acc & {data}
            end for 
            -- here we apply the system transform to the merged input
            if st[rUr] >= 0 then
                acc = call_func(st[rUr], {acc})
            elsif length(acc)=1 then
            -- prevent the accumulation of brackets around single datums
                acc = acc[1]
            end if
            -- compare the accumulated result with the sv element
            -- if different then
            if not equal(acc, sv[rUr]) then
                -- store the accumulated result in the corresponding sv element
                new_sv[rUr] = acc
                -- insert the rowId into cUpdate to indicate that the sv element has changed
                cUpdate = cUpdate & rUr
            end if
        end for
        sv = new_sv
    end for
    return {{st,sm,sv,cUpdate},nSteps}
end function

global procedure print_sv(sequence sv, sequence sv_names)
    for i = 1 to length(sv) do
        puts(1,"\n# ")
        puts(1, sv_names[i])
        puts(1, " #####\n")
        pretty_print(1,sv[i],{})
    end for
end procedure