I found out that the problem
is in the so called refinement procedure of the algorithm which consists of
10 different functions, returning an adjacency matrix. I¹m calling the
refinement procedure with
M <- refine1(M, A, B, p_A, p_B, FAIL) (note that all parameters in call
refine1 have been defined previosly)
Then the following steps look like this
#####################################
# Refinement process #
#####################################
refine1 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 1")
# elim marks if there was eliminated a 1 (and changed to 0)
lst <- vector(mode = "numeric")
elim <- 0
i <- 1
refine2(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
refine2 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 2")
k <- 1
h <- 1
sc <- vector(mode = "numeric", length = p_A)
for (l in 1:p_A){
sc[l] <- 0
}
sc[1] <- 1
refine3(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
refine3 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 3")
# check if there is a 1 in the current row in adjacency matrix of graph1
which is on the same position of the 1 in sc
# sc is a binary string whith only one 1. The position of the one goes
from the first position of sc to the last position
# and is used for scanning. First this is done for graph 1, and than for
graph 2. Then it is checked if the following condition
# is fulfilled: (for all x) ((A[i,x] = 1) => (there exists one ore more
y) (M[x,y] * B[x,j] = 1)). The algorithm terminates if no more
# 1 can be changed in a 0
if (1 %in% collation(A[i,], sc) == FALSE){
refine4(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
lst[k] <- h
k <- k+1
refine4(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
}
refine4 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 4")
sc <- c(0,sc[-p_A])
h <- h+1
if (k != (rowSums(A)[i])+1){
refine3(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
refine5(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
}
refine5 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 5")
j <- 1
sc <- vector(mode = "numeric", length = p_B)
for (l in 1:p_B){
sc[l] <- 0
}
sc[1] <- 1
refine6(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
refine6 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 6")
if (1 %in% collation(M[i,], sc) == FALSE){
refine9(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
h <- 1
refine7(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
}
refine7 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 7")
x <- lst[h]
if (1 %in% collation(M[x,], B[,j]) == FALSE){
refine8(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
h <- h+1
if (h != (rowSums(A)[i])+1){
refine7(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
refine9(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
}
}
refine8 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 8")
not_sc <- vector(mode = "numeric", length = p_B)
for (n in 1:p_B){
if (sc[n] == 1){
not_sc[n] <- 0
}
else {
not_sc[n] <- 1
}
}
M[i,] <- collation(M[i,], not_sc)
elim <- elim+1
h <- h+1
refine9(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
refine9 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
#print("refine 9")
sc <- c(0,sc[-p_B])
j <- j+1
if (j != p_B+1){
refine6(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
refine10(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
}
refine10 <- function(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x){
print("refine 10")
if (1 %in% M[i,] == FALSE){
FAIL <- 1
print(M)
return(FAIL)
}
else {
i <- i+1
if (i != (p_A)+1){
refine2(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst, x)
}
else {
if (elim != 0){
refine1(M, A, B, p_A, p_B, FAIL, elim, i, j, k, sc, h, lst,
x)
}
else {
return(M)
}
}
}
}
I really don¹t now where the problem is. Hope that anybody can help me
solving it.
[[alternative HTML version deleted]]
------------------------------------------------------------------------
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
