From 8e2d21b71f062ee80188194a22f1b187dd57e7bb Mon Sep 17 00:00:00 2001 From: Freguglia Date: Mon, 9 Aug 2021 14:47:25 -0300 Subject: [PATCH] Adds 'symmetric' family of restrictions. --- DESCRIPTION | 2 +- NEWS.md | 4 ++++ R/array_vector_conversions.R | 28 +++++++++++++++++++++++++++- R/family.R | 2 +- R/mrfout.R | 12 ++++++++++-- R/suf_stat.R | 24 ++++++++++++++++++------ tests/testthat/test-mrfout.R | 5 +++++ tests/testthat/test-smr_functions.R | 3 +++ tests/testthat/test-suf_stat.R | 8 ++++++-- 9 files changed, 75 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5266c80..4b5bdcf 100644 --- a/DESCRIPTION +++ b/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"), diff --git a/NEWS.md b/NEWS.md index 78a25fb..9bfad12 100644 --- a/NEWS.md +++ b/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. diff --git a/R/array_vector_conversions.R b/R/array_vector_conversions.R index 1cc4673..091936a 100644 --- a/R/array_vector_conversions.R +++ b/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.") } diff --git a/R/family.R b/R/family.R index d2aeec5..608b815 100644 --- a/R/family.R +++ b/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") diff --git a/R/mrfout.R b/R/mrfout.R index de283ae..64b4bd4 100644 --- a/R/mrfout.R +++ b/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") diff --git a/R/suf_stat.R b/R/suf_stat.R index 50beed0..cdd7714 100644 --- a/R/suf_stat.R +++ b/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.") } @@ -144,8 +149,8 @@ expand_array <- function(theta_vec, family, mrfi, C){ } #' @rdname smr_stat -#' -#' @inheritParams expand_array +#' +#' @inheritParams expand_array #' @return A `data.frame` describing the relative position #' and interaction associated with each potential in the vector #' form in each row, in the same order. @@ -162,21 +167,28 @@ vec_description <- function(mrfi, family, C){ } else if(family == "absdif"){ ints <- paste0("abs.dif. ",1:C) - res <- data.frame(position = rep(pos, each = C), + 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 = rep(pos, each = 2*C), + 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 = rep(pos, each = (C+1)^2 - 1), + 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.") } diff --git a/tests/testthat/test-mrfout.R b/tests/testthat/test-mrfout.R index 8859334..2dcef6f 100644 --- a/tests/testthat/test-mrfout.R +++ b/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") }) diff --git a/tests/testthat/test-smr_functions.R b/tests/testthat/test-smr_functions.R index 6cf2bbf..2e502a8 100644 --- a/tests/testthat/test-smr_functions.R +++ b/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")) }) diff --git a/tests/testthat/test-suf_stat.R b/tests/testthat/test-suf_stat.R index d4cc01b..7b08b0a 100644 --- a/tests/testthat/test-suf_stat.R +++ b/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)) })