--
-- ########################
-- # SMN core implementation
-- ########################
-- John Ringland 2005/08/27

include PNData.e -- includes FDData.e and SOME_Utilities.e

-- {a,b} => [a] column vector
--          [b]
global type StateVector(object sv)
    return sequence(sv) -- it may contain any objects, no need to test them
end type

global type Operator(object op)
    if integer(op) and (op >= -1) then
        return 1
    else
        return 0
    end if
end type

-- matrix element components
global constant con_    = 1 -- converging
global constant div_    = 2 -- diverging

global type SMElement(object sme)
    if sequence(sme) and length(sme)=2 
           and ((PNData(sme[con_]) and PNData(sme[div_]))
           or (Operator(sme[con_]) and Operator(sme[div_]))) then
        return 1
    elsif atom(sme) then
        return 1
    end if
    return 0
end type

-- include here so it gets SMElement definition
include SparseSM.e
include EnergyFlowUtils.e

-- {a,b} => [a b] row vector
global type SMRow(object smr)
    if SparseSMRow(smr) then
        for i = 1 to length(smr[Data_]) do
            if not SMElement(smr[Data_][i]) then
                return 0
            end if
        end for
    elsif sequence(smr) then -- should be an atom or a sequence containing Operators or PNData
        for i = 1 to length(smr) do
            if not SMElement(smr[i]) then
                return 0
            end if
        end for
    else
        return 0
    end if
    return 1
end type

-- {a,b} => [a] column vector
--          [b]
global type SMColumn(object smc)
    if SparseSMCol(smc) then
        return 1
    elsif sequence(smc) then -- should be an atom or a sequence containing Operators or PNData
        for i = 1 to length(smc) do
            if not SMElement(smc[i]) then
                return 0
            end if
        end for
    else
        return 0
    end if
    return 1
end type

-- {{a,b},{c,d}} => [a b] matrix represented as a column of rows
--                  [c d]
global type SystemMatrix(object sm)
    integer len
    if SparseSystemMatrix(sm) then -- detect sparseSM
        return 1
    end if
    
    if not sequence(sm) then 
        return 0 
    end if
    len = length(sm)
    for i = 1 to len do
        if not (SMRow(sm[i]) and length(sm[i])=len) then
            return 0
        end if
        -- the data may be any kind of object, no need to test them
    end for
    return 1
end type

-- smModel = {SM , SV}
global type smModel(object mdl)
    if sequence(mdl) and length(mdl) = 2 then
        return SystemMatrix(mdl[1]) and StateVector(mdl[2])
    end if
    return 0
end type

global function createSM(integer nRows, integer nCols, object data)
    if PNData(data) or Operator(data) then
        return repeat(repeat({data,data},nCols),nRows)
    end if
    return {} -- empty matrix
end function

global function createSV(integer nRows, object data)
    return repeat(data,nRows)
end function

global type OpSequence(object os)
    if sequence(os) then
        for i = 1 to length(os) do
            if not Operator(os[i]) then
                return 0
            end if
        end for
        return 1
    end if
    return 0
end type

function smPairwise(object out, SMElement d1, object d2)
    if atom(d1) and atom(d2) then
        return d1*d2
    elsif Operator(d1[con_]) and Operator(d1[div_]) then
        return call_func(d1[con_],{call_func(d1[div_],{d2})})
    else --if PNData(out) and PNData(d1) and PNData(d2) then
         -- type checking will ensure PNData'hood
        return smMultiply(out,d1[con_],smMultiply(out,d1[div_],d2))
    end if
end function
    
-- instead of 
-- M[i] use smGetRow(M,i)
global function smGetRow(SystemMatrix M, integer rowId)
    if SparseSystemMatrix(M) then
        return get_row_SM(M,rowId)
    else -- a full sequence row
        return M[rowId]
    end if
end function

-- equivalent to extracting the jth element out of each row
-- to form the jth column
global function smGetCol(SystemMatrix M, integer colId)
    SMColumn res res = {}
    
    if SparseSystemMatrix(M) then
        return get_col_SM(M,colId)
    else -- a full sequence row
        for i = 1 to length(M) do
            res = res & {M[i][colId]}
        end for
        return res
    end if
end function

-- instead of 
-- row[i] use smGetElemFromRow(row,i)
function smGetElemFromRow(SMRow row, integer colId)
    sequence rSearch
    if SparseSMRow(row) then
        -- search for element in row
        rSearch = BSearch(row[elemId_], colId)
        -- if element exists then
        if rSearch[found_] then
            return row[Data_][rSearch[index_]]
        else
            return {{fdZero},{fdZero}}
        end if
    else -- a full sequence row
        return row[colId]
    end if
end function

-- a general grouping function
function smGroup(SMRow row, StateVector sv)
    sequence out
    
    out = {}
    for i = 1 to length(sv) do
        out = out & {smPairwise(sv[i],smGetElemFromRow(row,i),sv[i])}
    end for
    return out
end function

-- a basic aggregate function
function smSummation_(sequence in) -- (PNData out, sequence of PNData (pairwise results from smGroup))
                                   -- or {atom out, sequence of atoms}
    object out
    if atom(in[1]) then
        out = 0
        for i = 1 to length(in[2]) do
            out += in[2][i]
        end for
    elsif PNData(in[1]) then
        out = eAssign(in[1],{fdZero})
        for i = 1 to length(in[2]) do
            out = smAdd(out,out, in[2][i])
        end for
    end if
    return out
end function
constant smSummation = routine_id ("smSummation_")

global type UpdateList(object cu)
    if sequence(cu) then
        for i = 1 to length(cu) do
            if not integer(cu[i]) then
                return 0
            end if
        end for
    else
        return 0
    end if
    return 1
end type


-- F for full processing
-- EF for energy flow processing
function smFGlobal_(sequence in) --(OpSequence A, SystemMatrix M, StateVector V, integer steps)
    OpSequence A
    SystemMatrix M
    integer len, steps
    StateVector new, curr

    A = in[1]
    M = in[2]
    curr = in[3]
    steps = in[4]
    
    len = length(curr)
    new = curr
    if length(A) = 1 then
        A = repeat(A[1], len)
    end if
    for t = 1 to steps do
        for i = 1 to len do
            new[i] = call_func(A[i],{{new[i],smGroup(smGetRow(M,i), curr)}})
        end for
        curr = new
    end for
    return new
end function
constant smFGlobal = routine_id ("smFGlobal_")


-- smEnergyFlowGlobal
-- implements the energy flow method on arbitrary matrices
-- smEFGlobal(OpSequence A, SystemMatrix M, StateVector V, integer steps, UpdateList cUpdate)
global function smEFGlobal_(sequence in) 
    OpSequence A
    SystemMatrix M
    UpdateList cUpdate, rUpdate
    SMColumn col
    StateVector new, curr
    integer len, rUr, steps
    
    A = in[1]
    M = in[2]
    curr = in[3]
    steps = in[4]
    cUpdate = in[5]

    len = length(curr)
    new = curr
    if length(A) = 1 then
        A = repeat(A[1], len)
    end if
    -- for each moment in time
    for n = 1 to steps 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 = smGetCol(M,cUpdate[c]) --get_col_SM(M,cUpdate[c])
            -- if the column existed in the matrix; for global outputs col=0
            if SMColumn(col) then
                -- to keep only the rowIds
                -- form the union of each column with rUpdate
                if SparseSystemMatrix(M) then
                    rUpdate = union_SET(rUpdate,col[elemId_])
                else -- full SystemMatrix
                    for i = 1 to length(col) do
                        if SMElement(col[i]) and ( ( atom(col[i]) and not equal(0,col[i]) ) or (  sequence(col[i]) and not equal(fdZero,col[i][div_]) ) ) then
                            rUpdate = union_SET(rUpdate,{i})
                        end if
                    end for
                end if
            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
            rUr = rUpdate[r]
            new[rUr] = call_func(A[rUr],{{new[rUr],smGroup(smGetRow(M,rUr), curr)}})
            -- compare the accumulated result with the sv element
            -- if different then
            if not equal(new[rUr], curr[rUr]) then
                -- insert the rowId into cUpdate to indicate that the sv element has changed
                cUpdate = cUpdate & rUr
            end if
        end for
        curr = new
    end for
    return new
end function
constant smEFGlobal = routine_id ("smEFGlobal_")


--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=-
-- in each of these two general functions one may use:

-- with this function one may implement any type of SMN
-- full (F) or energy flow (EF) processing depending on G
-- linear or nonlinear depending on A
-- full or sparse matrices depending on M
global function smGeneral(Operator G, OpSequence A, SystemMatrix M, StateVector V, integer steps, UpdateList cU)
    return call_func(G, {{A,M,V,steps,cU}})
end function

-- linear SMN
global function smLinear(Operator G, SystemMatrix M, StateVector V, integer steps, UpdateList cU)
    return smGeneral(G,{smSummation},M,V,steps,cU)
end function

-- full processing linear SMN
global function smFLinear(SystemMatrix M, StateVector V, integer steps)
    return smLinear(smFGlobal,M,V,steps,{})
end function

-- full processing general SMN
global function smFGeneral(OpSequence A, SystemMatrix M, StateVector V, integer steps)
    return smGeneral(smFGlobal,A,M,V,steps,{})
end function

-- energy flow processing linear SMN
global function smEFLinear(SystemMatrix M, StateVector V, integer steps, UpdateList cU)
    return smLinear(smEFGlobal,M,V,steps,cU)
end function

-- energy flow processing general SMN
global function smEFGeneral(OpSequence A, SystemMatrix M, StateVector V, integer steps, UpdateList cU)
    return smGeneral(smEFGlobal,A,M,V,steps,cU)
end function


--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=-

-- function used later in the test code
-- in[1] is output type
-- in[2] is the data, containing all three pair results
-- in[2][2] and in[2][3] are the desired datums
function NAND(sequence in)
    -- return a*b = 0
    if equal(smMultiply(in[1],in[2][2],in[2][3]), repeat(fdZero,length(in[1]))) then
        return {fdOne}
    else
        return {fdZero}
    end if
end function
constant NAND_id = routine_id("NAND")

procedure testSMN()
---------------------------------------------
-- A summary of the test procedure:
-- test linear SMN processing
    -- test non-sparse matrices
        -- test full processing with non-sparse linear Matrices
        -- test energy flow processing with non-sparse linear Matrices
    -- test sparse matrices
        -- test full processing with sparse linear Matrices
        -- test energy flow processing with sparse linear Matrices
-- test non-linear SMN processing
    -- test non-sparse matrices
        -- test full processing with non-sparse non-linear Matrices
        -- test energy flow processing with non-sparse non-linear Matrices
    -- test sparse matrices
        -- test full processing with sparse non-linear Matrices
        -- test energy flow processing with sparse non-linear Matrices
---------------------------------------------
-- let a = isLinear = 0 or 1   l means linear
--     b = isSparse = 0 or 1   1 means sparse
--     c = isEFlow  = 0 or 1   1 means uses energy flow processing
--
-- the permutation table is:
--------------------------------------------
--  a | b | c |    type     |   test order
--------------------------------------------
--  0 | 0 | 0 | nl, ns, nef |   5
--  0 | 0 | 1 | nl, ns,  ef |   6
--  0 | 1 | 0 | nl,  s, nef |   7
--  0 | 1 | 1 | nl,  s,  ef |   8   working
--  1 | 0 | 0 |  l, ns, nef |   1   working
--  1 | 0 | 1 |  l, ns,  ef |   2   working
--  1 | 1 | 0 |  l,  s, nef |   3
--  1 | 1 | 1 |  l,  s,  ef |   4
------------------------------------------
-- There are 2^3 = 8 basic types of SMN
-- So in this process there are 8 different ways in which
-- the kings of Edom can unite into a single king of Israel
------------------------------------
-- Test Code: ----------------------
------------------------------------
    SMElement r0, r1
    SystemMatrix M1, M2
    StateVector Vi, Vnew

    OpSequence st
    SystemMatrix sm
    StateVector sv
    UpdateList cUpdate
    object data
-- test linear SMN processing

    -- test non-sparse matrices
        
        r0 = {{fdZero},{fdZero}}
        r1 = {{fdOne},{fdOne}}
        M1 = {{r1,r0},{r0,r1}} -- identity
        M2 = {{r0,r1},{r1,r0}} -- flip
        Vi = {{defaultFDV(9)},{defaultFDV(15)}}
        -- test full processing with non-sparse linear Matrices
            puts(1,"\nTest:    Linear,    nonSparse,      nonEF   SMN\n")
            trace (1)
            ?{fdGet(first(Vi[1])),fdGet(first(Vi[2]))}
            Vnew = smFLinear(M1,Vi,1) -- identity
            puts(1,"After one iteration:\n")
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
            Vnew = smFLinear(M1,Vi,1000) -- identity
            puts(1,"After 1000 iterations, notice the drift due to FDData:\n")
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
            puts(1,"There's a steady drift upwards of the decimal digits.\n")
            Vnew = smFLinear(M2,Vi,1) -- flip
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
        -- test energy flow processing with non-sparse linear Matrices
            puts(1,"\nTest:    Linear,    nonSparse,         EF   SMN\n")
            trace(1)
            Vnew = smEFLinear(M1,Vi,1,{1,2}) -- identity
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
            Vnew = smEFLinear(M2,Vi,1,{1,2}) -- flip
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
    -- test sparse matrices
            M1 = createSparseSM(2,2)
            M1 = set_SM(M1, {1,1}, r1)
            M1 = set_SM(M1, {2,2}, r1)
            M2 = createSparseSM(2,2)
            M2 = set_SM(M2, {1,2}, r1)
            M2 = set_SM(M2, {2,1}, r1)
            Vi = {{defaultFDV(9)},{defaultFDV(15)}}
        
        -- test full processing with sparse linear Matrices
            puts(1,"\nTest:    Linear,       Sparse,      nonEF   SMN\n")
            trace(1)
            ?{fdGet(first(Vi[1])),fdGet(first(Vi[2]))}
            Vnew = smFLinear(M1,Vi,1)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
            Vnew = smFLinear(M2,Vi,1)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}

        -- test energy flow processing with sparse linear Matrices
            puts(1,"\nTest:    Linear,       Sparse,         EF   SMN\n")
            trace(1)
            cUpdate = {1,2} -- flag initial inputs for processing
            Vnew = smEFLinear(M1,Vi,1,cUpdate)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
            Vnew = smEFLinear(M2,Vi,1,cUpdate)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2]))}
             
-- test non-linear SMN processing

    -- test non-sparse matrices
        -- test full processing with non-sparse non-linear Matrices
            puts(1,"\nTest: nonLinear,    nonSparse,      nonEF   SMN\n")
        -- test energy flow processing with non-sparse non-linear Matrices
            puts(1,"\nTest: nonLinear,    nonSparse,         EF   SMN\n")
    -- test sparse matrices
        -- test full processing with sparse non-linear Matrices
            puts(1,"\nTest: nonLinear,       Sparse,      nonEF   SMN\n")
        -- test energy flow processing with sparse non-linear Matrices
            data = r1 -- input selector
            st = {NAND_id,-1,-1} -- OpSequence undefined for rows 2 and 3
            -- construct the Metric System Matrix element by element
            sm = createSparseSM(3,3)
            sm = set_SM(sm, {1,2}, data)
            sm = set_SM(sm, {1,3}, data)

            -- construct state vector, sv = {0,1,1}
            -- 1 NAND 1 = 0
            sv = {{fdZero},{fdOne},{fdOne}}
            puts(1,"\nTest: nonLinear,       Sparse,         EF   SMN\n")
            trace(1)
            ?{fdGet(first(sv[1])),fdGet(first(sv[2])),fdGet(first(sv[3]))}
            cUpdate = {2,3} -- flag initial inputs for processing
            Vnew = smEFGeneral(st,sm,sv,1,cUpdate)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2])),fdGet(first(Vnew[3]))}

            -- construct state vector, sv = {0,0,1}
            -- 0 NAND 1 = 1
            sv = {{fdZero},{fdZero},{fdOne}}
            ?{fdGet(first(sv[1])),fdGet(first(sv[2])),fdGet(first(sv[3]))}
            Vnew = smEFGeneral(st,sm,sv,1,cUpdate)
            ?{fdGet(first(Vnew[1])),fdGet(first(Vnew[2])),fdGet(first(Vnew[3]))}

end procedure

--testSMN()
-- It works!

--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=--=#=-

-- General SMN tools for construction, manipulation and visualistion of SM's and SV's

-- smLinear tools:
-- ###############
--
-- The SM and SV contains PNData