"MTcoherence.fun" = function(Dataset, plot = T)
{
	options(warn = -1)
	# Dataset: theta, var.theta, index with duplication in multitreatment and 
	#two treatments as.characters # if multiple treatment studies, they must contain ALL comparison 
	#WARNING: do not use numbers as names for the treatments if you have more than 9!!!
	#-----Get the data
theta<-Dataset[, 1]
var.theta<-Dataset[, 2]
id<-Dataset[, 3]
treatA<-Dataset[, 4]
treatB<-Dataset[, 5]
no.rows <- length(theta)
	#-----Convert treatments to characters
	if(data.class(treatA) != "character" | data.class(treatB) != 
		"character") {
		treatA <- as.character(treatA)
		treatB <- as.character(treatB)
	}
	#-----Sort the treatments so that horizontaly treatA>treatB
	inlex.order <- as.numeric(apply(Dataset[, 4:5], 1, order)[1,  ] == 1)
	treat.A <- treatA
	treat.A[inlex.order == 0] <- treatB[inlex.order == 0]
	treat.B <- treatB
	treat.B[inlex.order == 0] <- treatA[inlex.order == 0]
	treatA <- treat.A
	treatB <- treat.B
	theta1 <- theta
	theta1[inlex.order == 0] <-  - theta1[inlex.order == 0]
	theta <- theta1
	#---- Buid up comparisons
	comp <- c()
	for(i in 1:no.rows) {
		comp <- append(comp, paste(sort(c(treatA[i], treatB[i])),
			sep = "", collapse = ""))
	}
	#---- Restructure Dataset
	Dataset <- cbind.data.frame(theta1, var.theta, id, comp)
	treatments <- unique(c(treatA, treatB))
	av.treat <- unique(comp)
	#---- Report a little 
	cat("\n", 
		"\n  *-----  Evaluating the coherence of the network ------*",
		"\n")
	cat("\n  Nr of treatments: ", length(treatments))
	#---- Construct the network
	# Builts all possible triangles
	pos.l <- .possible.loops(treatments)[1]
	cat("\n  Nr of all possible first order loops (triangles): ", length(
		pos.l[[1]]))
	#---- Available loops
	availability <- apply(pos.l[[1]], 2, match, av.treat)
	availability <- !is.na(apply(availability, 2, sum))
	pos.l <- pos.l[[1]][, availability]
	pos.l <- lapply(apply(pos.l, 2, list), unlist)
	no.loops <- length(pos.l)
	loops <- length(pos.l)
	cat("\n  Nr of available first order loops: ", length(pos.l), "\n",
		"\n")
	#---- Evaluation of the loops
	EVALUATE.LOOP <- c()
	LOOPS <- c()

	for(i in 1:no.loops) {
		cat("\n", i, ": Evaluation of the loop", names(pos.l[i]))
		loop <- c()
		loop.dimension <- length(pos.l[[i]])
		table.comp <- table(match(comp, pos.l[[i]]))
		names(table.comp) <- pos.l[[i]]
		cat("\n", "Direct comparisons in the loop:", "\n")
		print(table.comp)
		duplicat <- table(id[!is.na(match(comp, pos.l[[i]]))]) > 2
		# searching for multitreatment studies reporting more than 2 studies in the loop
		if(sum(duplicat) > 0) {
			duplication <- names(duplicat[duplicat])
			DatasetOnuse <- Dataset
			leaveout <- c()
			for(k in 1:length(duplication)) {
				most.rep.comp <- names(table.comp[table.comp ==
					max(table.comp)])
				#which is the most popular contrast? to take it out?
				cat("\n The study with id=", duplication[
					k], sep = "")
				cat(" has more than two treatments in this loop"
					)
				if(length(most.rep.comp) > 1) {
					most.rep.comp <- as.character(
						most.rep.comp[1])
				}
				leaveout <- append(leaveout, c(1:no.rows)[
					id == duplication[k] & comp == 
					most.rep.comp])
				cat(" (out is row nr ", leaveout[k], 
					" for comparison ", as.character(
					comp[leaveout[k]]), ")", sep = "")
			}
			DatasetOnuse <- DatasetOnuse[ - leaveout,  ]
		}
		else {
			DatasetOnuse <- Dataset
		}
		for(j in 1:loop.dimension) {
			cat("\n   Meta-analysis for the", pos.l[[i]][j], "comparison"
				)
			dta <- DatasetOnuse[DatasetOnuse$comp == pos.l[[i]][
				j], 1:2]
			if(dim(dta)[1] == 1) {
				metaanalysis <- c(Value = dta[, 1], var = sqrt(
					dta[, 2]))
			}
			else {
				if(dim(dta)[1] >= 1) {
				o<-metagen(dta[,1],sqrt(dta[,2]))
				metaanalysis <- c(o$TE.random, o$seTE.random)}
				}
				
			cat("\n   mean(se)=", round(metaanalysis[1], 3))
			cat("(", round(metaanalysis[2], 3), ")", sep = "")
			loop <- rbind(loop, metaanalysis)
		}
		dimnames(loop)[1] <- list(pos.l[[i]])
		loop[, 2] <- loop[, 2]^2
		leftLoop <- c(sum(loop[1:(loop.dimension - 1), 1]), sum(loop[
			1:(loop.dimension - 1), 2]))
		cat("\n   Indirect summary estimate for the", names(table.comp)[
			3], "comparison")
		cat("\n   Mean(se)=", round(leftLoop[1], 3))
		cat("(", round(sqrt(leftLoop[2]), 3), ")", sep = "")
		EvaluateLoop <- c(leftLoop[1] - loop[loop.dimension, 1], 
			leftLoop[2] + loop[loop.dimension, 2])
		cat("\n", "\n   Incoherence within the loop")
		cat(":  Mean(se)=", round(EvaluateLoop[1], 3))
		cat("(", round(sqrt(EvaluateLoop[2]), 3), ")", sep = "", "\n",
			"\n")
		LOOPS <- append(LOOPS, list(loop))
		EVALUATE.LOOP <- rbind(EVALUATE.LOOP, EvaluateLoop)
	}
	names(LOOPS) <- names(pos.l)
	dimnames(EVALUATE.LOOP)[1] <- list(names(pos.l))
	#---- Plot
	if(plot) {
		plot(metagen(abs(EVALUATE.LOOP[,1]),sqrt(EVALUATE.LOOP[,2])), studlab=names(LOOPS),xlab="Incoherence")
	}
	out <- cbind.data.frame(names(LOOPS),abs(EVALUATE.LOOP[,1]),sqrt(EVALUATE.LOOP[,2]))
names(out)<-c("Loop","ES.Inc","seES.Inc")
	invisible()
	options(warn = 1)
	out
}





".possible.loops" = function(treatments)
{
	#Given the names of the treatments, it returns all possible loops to evaluate for concordence
	#and the 'price' to pay for any loop in terms of needed associations. 
	# Note that the bigest loop is not returned
	# Mind you! notation follows a circle
	t <- length(treatments)
	z <- list()
	u <- c()
	for(i in 1:t) {
		u <- cbind(u, (rbind(treatments[i], combn(2, x = treatments[
			 - i]))))
	}
	u <- apply(u, 2, sort)
	price.in.vectors <- apply(u, 2, .needed.vectors)
	dimnames(price.in.vectors) <- list(NULL, apply(u, 2, paste, collapse
		 = ""))
	price.in.vectors <- price.in.vectors[, unique(apply(u, 2, paste, 
		collapse = ""))]
	z <- append(z, list(price.in.vectors))
	z
}



"combn" = function(m, x, fun = NULL, simplify = TRUE, ...)
{
	if(length(m) > 1) {
		warning(paste("Argument m has", length(m), 
			"elements: only the first used"))
		m <- m[1]
	}
	if(m < 0)
		stop("m < 0")
	if(m == 0)
		return(if(simplify) vector(mode(x), 0) else list())
	if(is.numeric(x) && length(x) == 1 && x > 0 && trunc(x) == x)
		x <- seq(x)
	n <- length(x)
	if(n < m)
		stop("n < m")
	e <- 0
	h <- m
	a <- 1:m
	nofun <- is.null(fun)
	count <- nCm(n, m, 0.1)
	out <- vector("list", count)
	out[[1]] <- if(nofun) x[a] else fun(x[a], ...)
	if(simplify) {
		dim.use <- NULL
		if(nofun) {
			if(count > 1)
				dim.use <- c(m, count)
		}
		else {
			out1 <- out[[1]]
			d <- dim(out1)
			if(count > 1) {
				if(length(d) > 1)
					dim.use <- c(d, count)
				else if(length(out1) > 1)
					dim.use <- c(length(out1), count)
			}
			else if(length(d) > 1)
				dim.use <- d
		}
	}
	i <- 2
	nmmp1 <- n - m + 1
	mp1 <- m + 1
	while(a[1] != nmmp1) {
		if(e < n - h) {
			h <- 1
			e <- a[m]
			j <- 1
		}
		else {
			h <- h + 1
			e <- a[mp1 - h]
			j <- 1:h
		}
		a[m - h + j] <- e + j
		out[[i]] <- if(nofun) x[a] else fun(x[a], ...)
		i <- i + 1
	}
	if(simplify) {
		if(is.null(dim.use))
			out <- unlist(out)
		else out <- array(unlist(out), dim.use)
	}
	out
}


"nCm" = function(n, m, tol = 1e-008)
{
	len <- max(length(n), length(m))
	out <- numeric(len)
	n <- rep(n, length = len)
	m <- rep(m, length = len)
	mint <- (trunc(m) == m)
	out[!mint] <- NA
	out[m == 0] <- 1
	whichm <- (mint & m > 0)
	whichn <- (n < 0)
	which <- (whichm & whichn)
	if(any(which)) {
		nnow <- n[which]
		mnow <- m[which]
		out[which] <- ((-1)^mnow) * Recall(mnow - nnow - 1, mnow)
	}
	whichn <- (n > 0)
	nint <- (trunc(n) == n)
	which <- (whichm & whichn & !nint & n < m)
	if(any(which)) {
		nnow <- n[which]
		mnow <- m[which]
		foo <- function(j, nn, mm)
		{
			n <- nn[j]
			m <- mm[j]
			iseq <- seq(n - m + 1, n)
			negs <- sum(iseq < 0)
			((-1)^negs) * exp(sum(log(abs(iseq))) - lgamma(m +
				1))
		}
		out[which] <- unlist(lapply(seq(along = nnow), foo, nn = nnow,
			mm = mnow))
	}
	which <- (whichm & whichn & n >= m)
	nnow <- n[which]
	mnow <- m[which]
	out[which] <- exp(lgamma(nnow + 1) - lgamma(mnow + 1) - lgamma(nnow -
		mnow + 1))
	nna <- !is.na(out)
	outnow <- out[nna]
	rout <- round(outnow)
	smalldif <- abs(rout - outnow) < tol
	outnow[smalldif] <- rout[smalldif]
	out[nna] <- outnow
	out
}



".needed.vectors" = function(loop)
{
	# function for the possible loops. It calculates the vectors needed to evaluate a loop
	#loops are just a loop<-c("a","b","c")
	leng <- length(loop)
	need.vec <- apply(apply(rbind(loop, c(loop[-1], loop[1])), 2, sort),
		2, paste, collapse = "")
	need.vec
}