Search code examples
mathmatrixjuliasymbolic-math

Rewriting on Infinite Matrices


I've taken a project named with "Symbolic Linear Algebra" which is about doing basic operations on infinite matrices like addition, multiplication, accessing specific element etc. I will be implementing those on Julia.

For specifying those infinite matrices we'll have some mathematical cases like:

Cases

So the visual representation of matrix will be like:

visual_rep_matrix

For example let's say we want to find A + A' for this example. Here our cases change so we need to rewrite those cases to get desired output right ? I know Mathematica does this but how can I implement this? Yes, this was too general so let me ask some questions;

  • Let's start with taking cases as input. There can be many cases with different rules like if i % 2 == 0 or i == j like in this example how can I provide a generic input ?
  • Let's say that I'm done with input and I want to make those simple operations. How can I combine those cases in a programming language like Julia ?

I've wrote some non-generic dumb code to see how things will go so I will provide my code to apply minimum reproducible example but don't take it seriously, I think I'm just looking for a clue or a roadmap to get rid of the question marks in my head.

 using Parameters



struct inf_matrix
   mod_of :: Integer
   mod_value :: Integer
   i_coefficient :: Integer
   j_coefficient :: Integer
   value :: Integer
end


function single_demo(_mod_of :: Integer, _mod_value :: Integer, _i_coefficient :: Integer, _j_coefficient :: Integer, _value :: Integer)
   test_matrix = inf_matrix(_mod_of, _mod_value, _i_coefficient, _j_coefficient, _value)
   return test_matrix
end

function get_elem(st::inf_matrix ,i :: Integer, j :: Integer)
   #This function is not completed yet  
   if (i % st.mod_of == st.mod_value) && (2 * st.i_coefficient == j)
         return st.value;
   else
      return -1
   end
   

end

demo_1 = single_demo(2, 0 ,1, 2, 1)

println(get_elem(demo_1, 1, 0))

Any help would be appreciated.


Solution

  • Here is how you could do this

    import Base: getindex, +, *
    
    abstract type InfiniteMatrix end
    
    struct InfiniteIdentity <: InfiniteMatrix end
    
    getindex(::InfiniteIdentity, i, j) = i .== j'
    
    struct InfiniteConstant <: InfiniteMatrix 
        c
    end
    
    getindex(m::InfiniteConstant, i::Integer, j::Integer) = m.c
    getindex(m::InfiniteConstant, i, j) = fill(m.c, size(i)..., size(j)...)
    
    struct InfiniteMatrixFilter <: InfiniteMatrix
        condition::Function
        iftrue::InfiniteMatrix
        iffalse::InfiniteMatrix
    end
    
    getindex(m::InfiniteMatrixFilter, i, j) = ifelse.(m.condition.(i,j'), m.iftrue[i,j], m.iffalse[i,j]) 
    
    
    struct InfiniteMatrixFunction <: InfiniteMatrix
        f::Function
        args
    end
    
    getindex(m::InfiniteMatrixFunction, i, j) = m.f(getindex.(m.args, Ref(i), Ref(j))...)
    
    +(m1::InfiniteMatrix, m2::InfiniteMatrix) = InfiniteMatrixFunction(+, (m1, m2))
    *(n::Number, m::InfiniteMatrix) = InfiniteMatrixFunction(x -> n*x, (m,))
    
    
    julia> i = InfiniteIdentity()
    InfiniteIdentity()
    
    julia> c1 = InfiniteConstant(1)
    InfiniteConstant(1)
    
    julia> (2i+3c1)[1:5, 1:5]
    5×5 Array{Int64,2}:
     5  3  3  3  3
     3  5  3  3  3
     3  3  5  3  3
     3  3  3  5  3
     3  3  3  3  5
    
    julia> m = InfiniteMatrixFilter((i,j) -> i%2 == 0, c1, 0c1)
    InfiniteMatrixFilter(var"#43#44"(), InfiniteConstant(1), InfiniteMatrixFunction(var"#41#42"{Int64}(0), (InfiniteConstant(1),)))
    
    julia> m[1:5, 1:5]
    5×5 Array{Int64,2}:
     0  0  0  0  0
     1  1  1  1  1
     0  0  0  0  0
     1  1  1  1  1
     0  0  0  0  0
    

    (this is only a proof of concept and it's not optimized or bugfree)