circulant implementation to obtain structure of CGMRF

Author: ErickChacon

README.md

@@ -14,3 +14,7 @@ Implementation of GMRF for spatial analysis.
 ## Graphs
 
 ![graphs](figures/03-graph.png)
+
+## Todo
+
+- [ ] Check structure_base to return always `Integers`

src/GMRFs.jl

@@ -1,7 +1,7 @@
 module GMRFs
 
 using Meshes
-import SparseArrays: sparse, sparsevec, spdiagm, spzeros
+import SparseArrays: sparse, sparsevec, spdiagm, spzeros, findnz, SparseVector, SparseMatrixCSC
 import SuiteSparse # this is only for fixes
 import FFTW
 import Distributions: InverseGamma, Gamma, Distributions

src/gmrf/cgmrf.jl

@@ -31,7 +31,16 @@ CGMRF(domain::CartesianGrid, order::Integer, κ::Real, δ::Real) =
 Base.length(d::CGMRF) = length(d.base)
 scale(d::CGMRF) = d.κ
 structure_base(d::CGMRF) = d.base
-# structure(d::CGMRF) = structure(d.g; δ = d.δ, order = d.order, circular = true)
+
+function structure(d::CGMRF)
+    # get base and transpose for 2d
+    base = structure_base(d)
+    if base isa SparseMatrixCSC
+        base = sparse(base')
+    end
+    # convert to circulant
+    spcirculant(base)
+end
 
 ## Random generator
 

src/utils.jl

@@ -1,4 +1,4 @@
-for filename in ["adjacency.jl", "difference.jl", "structure.jl"]
+for filename in ["adjacency.jl", "difference.jl", "structure.jl", "circulant.jl"]
     include(joinpath("utils", filename))
 end
 

src/utils/circulant.jl

@@ -0,0 +1,40 @@
+"""
+    spcirculant(x)
+
+Return a sparse circulant or block-circulant matrix based on the input sparse vector or
+matrix `x`.
+
+The size of the resulting matrix is `n×n`, where `n` represents the total number of
+elements in `x`. If `x` is a `SparseVector`, the function generates a `n×n` sparse
+circulant matrix. Conversely, if `x` is a `SparseMatrixCSC`, the function creates a `n×n`
+sparse block circulant matrix.
+"""
+
+function spcirculant(x::SparseVector)
+    # get indices and values
+    n = length(x)
+    (is, vals) = findnz(x)
+
+    # compute the indices for the full structure matrix
+    Is = repeat([i for i in 1:n], length(vals))
+    Js = [mod1(i -is[k]+1, n) for i in 1:n, k in 1:length(vals)]
+    Vs = repeat(vals, inner = n)
+
+    # spare structure matrix
+    sparse(Is, vec(Js), Vs, n, n)
+end
+
+function spcirculant(x::SparseMatrixCSC)
+    # get base and values
+    (n1,n2) = size(x)
+    (is, js, vals) = findnz(x)
+
+    # compute the indices for the full structure matrix
+    Is = repeat(vec([i + n1 * (b-1) for i in 1:n1, b in 1:n2]), length(vals))
+    Js = [mod1(n1*(js[k]-1) + n1 * (b-1) + mod1(i -is[k]+1, n1), n1*n2) for i in 1:n1, b in 1:n2, k in 1:length(vals)]
+    Vs = repeat(vals, inner = n1 * n2)
+
+    # spare structure matrix
+    sparse(Is, vec(Js), Vs, n1 * n2, n1 * n2)
+end
+

test/cgmrf.jl

@@ -14,6 +14,7 @@
         @test length(X) == n
         @test GMRFs.scale(X) == kappa[i]
         @test GMRFs.structure_base(X) == GMRFs.structure_base(grid, order = order[i], δ = δ[i])
+        @test GMRFs.structure(X) == GMRFs.structure(grid, δ = δ[i], order = order[i], circular = true)
         # rand and logpdf: single
         x = rand(X)
         @test length(x) == n
@@ -41,6 +42,7 @@ end
         @test length(X) == n1 * n2
         @test GMRFs.scale(X) == kappa[i]
         @test GMRFs.structure_base(X) == GMRFs.structure_base(grid, order = order[i], δ = δ[i])
+        @test GMRFs.structure(X) == GMRFs.structure(grid, δ = δ[i], order = order[i], circular = true)
         # rand and logpdf: single
         x = rand(X)
         @test length(x) == n1 * n2

test/utils.jl

@@ -112,12 +112,12 @@ end
     # circular: structure
     S = GMRFs.structure(grid, order = 1, circular = true)
     @test S == GMRFs.structure(graphc, order = 1)
-    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 1))
+    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 1)')
     S = GMRFs.structure(grid, order = 2, circular = true)
     @test S == GMRFs.structure(graphc, order = 2)
-    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 2))
+    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 2)')
     S = GMRFs.structure(grid, order = 3, circular = true)
     @test S == GMRFs.structure(graphc, order = 3)
-    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 3))
+    @test S[1, :] == vec(GMRFs.structure_base(grid, order = 3)')
 
 end