Adds 'symmetric' family of restrictions.

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mrf2d
 Type: Package
 Title: Markov Random Field Models for Image Analysis
-Version: 0.5.0.9999
+Version: 0.6.0.9999
 Authors@R: 
     person("Victor", "Freguglia", email = "victorfreguglia@gmail.com",
              role = c("aut", "cre"), 

NEWS.md

@@ -1,3 +1,7 @@
+# mrf2d 0.6
+
+  * Adds `symmetric` family for equality constraints of the potentials array.
+
 # mrf2d 0.5.1
 
   * Added `vec_description()` function to explicitly describe what relative positions and interactions are associated with each element of the summarized vector of parameters/sufficient statistics in each parameter restriction family.

R/array_vector_conversions.R

@@ -49,6 +49,10 @@ is_valid_array <- function(arr, family){
       m[1,1] == 0
     }))
 
+  } else if(family == "symmetric") {
+    all(apply(arr, MARGIN = 3, function(m) {
+      isSymmetric(m) & m[1,1] == 0
+    }))
   } else {
     stop("'", family, "' is not an implemented family.")
   }
@@ -76,7 +80,7 @@ vec_to_array <- function(vec, family, C, n_R){
     })))
 
   } else if(family == "dif") {
-    # Potential associated with zero difference (diagonal) is the sum of others.
+    # Potential associated with zero difference (diagonal) is 0.
     # The order is from -C to -1 then 1 to C.
     if(length(vec) != n_R*2*C) { stop("'vec' must have length n_R*2*C for family 'dif'.")}
     sanitize_theta(simplify2array(lapply(1:n_R, function(i){
@@ -98,6 +102,22 @@ vec_to_array <- function(vec, family, C, n_R){
              nrow = C+1)
     })))
 
+  } else if(family == "symmetric") {
+    if(length(vec) != (n_R*((C+2)*(C+1)/2 - 1))) {
+      stop("'vec' must have length n_R*((C+2)*(C+1)/2 -1) for family 'symmetric'.")}
+    # Parameters are filled by Columns.
+    sanitize_theta(simplify2array(lapply(1:n_R, function(i) {
+      m <- diag(C+1)*0
+      t <- 1
+      sub_vec <- c(0,vec[(1 + (i-1)*(((C+1)*(C+2)/2)-1)) : (i*((C+1)*(C+2)/2) - i)])
+      for(j in 1:(C+1)){
+        for(l in j:(C+1)){
+          m[j,l] <- m[l,j] <- sub_vec[t]
+          t <- t+1
+        }
+      }
+      return(m)
+    })))
   } else {
     stop("'", family, "' is not an implemented family.")
   }
@@ -137,6 +157,12 @@ array_to_vec <- function(arr, family){
       return(as.vector(m)[-1])
     }))
 
+  } else if(family == "symmetric") {
+    as.vector(apply(arr, MARGIN = 3, function(m){
+      out <- m[lower.tri(m, diag = TRUE)]
+      out[-1]
+    }))
+
   } else {
     stop("'", family, "' is not an implemented family.")
   }

R/family.R

@@ -53,4 +53,4 @@
 #' \url{https://arxiv.org/abs/2006.00383}
 NULL
 
-mrf2d_families <- c("onepar", "oneeach", "absdif", "dif", "free")
+mrf2d_families <- c("onepar", "oneeach", "absdif", "dif", "symmetric", "free")

R/mrfout.R

@@ -85,7 +85,13 @@ summary.mrfout <- function(object, ...){
     } else if (family == "free"){
       cat("Interaction for pairs of values:", "\n")
       cat("Position|",
-          sprintf("%-6s",paste0("(", rep(0:C, C+1), ",", rep(0:C, each = C+1), ")")[-1]),
+          sprintf("%-6s",paste0(rep(0:C, C+1), ",", rep(0:C, each = C+1))[-1]),
+          " Rel. Contribution")
+      cat("\n")
+    } else if (family == "symmetric"){
+      cat("Interaction for pairs of values:", "\n")
+      cat("Position|",
+          sprintf("%-6s",vec_description(mrfi, family, C)$interaction[1:((C+2)*(C+1)/2 - 1)]),
           " Rel. Contribution")
       cat("\n")
     }
@@ -113,7 +119,9 @@ plot.mrfout <- function(x, ...){
   } else if(x$family == "dif"){
     colnames(df) <- (-C:C)[-(C+1)]
   } else if(x$family == "free"){
-    colnames(df) <- paste0("(", rep(0:C, C+1), ",", rep(0:C, each = C+1), ")")[-1]
+    colnames(df) <- paste0(rep(0:C, C+1), ",", rep(0:C, each = C+1))[-1]
+  } else if(x$family == "symmetric"){
+    colnames(df) <- vec_description(mrfi(1), x$family, C)$interaction[1:((C+2)*(C+1)/2 - 1)]
   }
   df <- cbind(x$mrfi@Rmat, df)
   colnames(df)[1:2] <- c("rx", "ry")

R/suf_stat.R

@@ -32,6 +32,11 @@ suf_stat <- function(arr, family){
       return(as.vector(m)[-1])
     }))
 
+  } else if(family == "symmetric"){
+    as.vector(apply(arr, MARGIN = 3, function(m){
+      (t(m)+(m*lower.tri(m)))[lower.tri(m, diag = TRUE)][-1]
+    }))
+
   } else {
     stop("'", family, "' is not an implemented family.")
   }
@@ -164,20 +169,29 @@ vec_description <- function(mrfi, family, C){
 
     } else if(family == "absdif"){
         ints <- paste0("abs.dif. ",1:C)
-        res <- data.frame(position = as.factor(rep(pos, each = C)),
-                          interaction = as.factor(rep(ints, times = length(pos))))
+        res <- data.frame(position = rep(pos, each = C),
+                          interaction = rep(ints, times = length(pos)))
 
     } else if(family == "dif"){
         ints <- paste0("dif. ", c(-C:-1,1:C))
-        res <- data.frame(position = as.factor(rep(pos, each = 2*C)),
-                          interaction = as.factor(rep(ints, times = length(pos))))
+        res <- data.frame(position = rep(pos, each = 2*C),
+                          interaction = rep(ints, times = length(pos)))
+
 
     } else if(family == "free"){
         arr <- array(dim=c(C+1, C+1, length(pos)))
         ints <- paste0(rep(0:C, times = C+1), ",", rep(0:C, each = C+1))
         ints <- ints[ints != "0,0"]
-        res <- data.frame(position = as.factor(rep(pos, each = (C+1)^2 - 1)),
-                          interaction = as.factor(rep(ints, times = length(pos))))
+        res <- data.frame(position = rep(pos, each = (C+1)^2 - 1),
+                          interaction = rep(ints, times = length(pos)))
+
+    } else if(family == "symmetric"){
+        m <- matrix(nrow = C+1, ncol = C+1)
+        ints <- paste0(row(m)[lower.tri(m,TRUE)]-1, ",", col(m)[lower.tri(m, TRUE)]-1)
+        ints <- ints[ints != "0,0"]
+        res <- data.frame(position = rep(pos, each = (C+1)*(C+2)/2 - 1),
+                          interaction = rep(ints, times = length(pos)))
+
 
     } else {
         stop("'", family, "' is not an implemented family.")

tests/testthat/test-mrfout.R

@@ -29,4 +29,9 @@ test_that("mrfout methods work", {
   prnt <- capture.output(summary(a))
   expect_true(sum(grepl("Interaction", prnt)) > 0)
   expect_is(plot(a), "ggplot")
+
+  a <- fit_pl(Z_potts[1:30, 1:30], mrfi(1), family = "symmetric")
+  prnt <- capture.output(summary(a))
+  expect_true(sum(grepl("Interaction", prnt)) > 0)
+  expect_is(plot(a), "ggplot")
 })

tests/testthat/test-smr_functions.R

@@ -4,6 +4,7 @@ test_that("smr_array works", {
   expect_length(smr_array(theta_potts, "absdif"), 2*2)
   expect_length(smr_array(theta_potts, "dif"), 2*2*2)
   expect_length(smr_array(theta_potts, "free"), 2*(3^2 - 1))
+  expect_length(smr_array(theta_potts, "symmetric"), 2*(3*2 - 1))
 })
 
 test_that("smr_stat works", {
@@ -12,6 +13,7 @@ test_that("smr_stat works", {
   expect_length(smr_stat(Z_potts, mrfi(), "absdif"), 2*2)
   expect_length(smr_stat(Z_potts, mrfi(), "dif"), 2*2*2)
   expect_length(smr_stat(Z_potts, mrfi(), "free"), 2*(3^2 - 1))
+  expect_length(smr_stat(Z_potts, mrfi(), "symmetric"), 2*(3*2 - 1))
 })
 
 test_that("expand_array", {
@@ -22,5 +24,6 @@ test_that("expand_array", {
   expect_equivalent(1:4, smr_array(expand_array(1:4, "absdif", mrfi(1), 2), "absdif"))
   expect_equivalent(1:8, smr_array(expand_array(1:8, "dif", mrfi(1), 2), "dif"))
   expect_equivalent(1:16, smr_array(expand_array(1:16, "free", mrfi(1), 2), "free"))
+  expect_equivalent(1:10, smr_array(expand_array(1:10, "symmetric", mrfi(1), 2), "symmetric"))
 })
 

tests/testthat/test-suf_stat.R

@@ -5,6 +5,7 @@ test_that("sufficient statistics computing works", {
   expect_equal(length(suf_stat(pc, "absdif")), 2*2)
   expect_equal(length(suf_stat(pc, "dif")), 2*4)
   expect_equal(length(suf_stat(pc, "free")), 2*(3*3 -1))
+  expect_equal(length(suf_stat(pc, "symmetric")), 2*(3*2 -1))
 
   Z2 <- Z_potts
   Z2[10,] <- NA
@@ -14,6 +15,7 @@ test_that("sufficient statistics computing works", {
   expect_equal(length(suf_stat(pc, "absdif")), 2*2)
   expect_equal(length(suf_stat(pc, "dif")), 2*4)
   expect_equal(length(suf_stat(pc, "free")), 2*(3*3 -1))
+  expect_equal(length(suf_stat(pc, "symmetric")), 2*(3*2 -1))
 })
 
 test_that("summarized vector description works", {
@@ -23,9 +25,11 @@ test_that("summarized vector description works", {
   expect_is(vec_description(mrfi(1) + c(3,4), "oneeach", 1), "data.frame")
   expect_is(vec_description(mrfi(1) + c(3,4), "absdif", 1), "data.frame")
   expect_is(vec_description(mrfi(1) + c(3,4), "absdif", 2), "data.frame")
+  expect_is(vec_description(mrfi(1) + c(3,4), "dif", 1), "data.frame")
   expect_is(vec_description(mrfi(1) + c(3,4), "dif", 2), "data.frame")
-  expect_is(vec_description(mrfi(1) + c(3,4), "dif", 2), "data.frame")
-  expect_is(vec_description(mrfi(1) + c(3,4), "free", 2), "data.frame")
+  expect_is(vec_description(mrfi(1) + c(3,4), "free", 1), "data.frame")
   expect_is(vec_description(mrfi(1) + c(3,4), "free", 2), "data.frame")
+  expect_is(vec_description(mrfi(1) + c(3,4), "symmetric", 1), "data.frame")
+  expect_is(vec_description(mrfi(1) + c(3,4), "symmetric", 2), "data.frame")
   expect_error(vec_description(mrfi(1) + c(3,4), "invalidname", 2))
 })