# To compare the result
ROC = function ( G, est ) 
{
	tp   = sum( ( G != 0 ) * ( est != 0 ) )
	fp   = sum( ( G == 0 ) * ( est != 0 ) )
	fn   = sum( ( G != 0 ) * ( est == 0 ) )
	tn_M = ( G == 0 ) * ( est == 0 )
	tn   = sum( tn_M[ upper.tri( tn_M == 1 ) ] )
	
	# Precision is the probability that a randomly selected link is relevant 
	Precision <- tp / ( tp + fp ) 
	
	# Recall is the probability that a randomly selected relevant link 	
	Recall <- tp / ( tp + fn ) # also called TPR
	
	FPR <- fp / ( fp + tn ) # False positive rate
	
	Accuracy <- ( tp + tn ) / ( tp + tn + fp + fn )
	
	# harmonic mean of precision and recall, called F-measure or balanced F-score
	F1score <- ( 2 * tp ) / ( 2 * tp + fp + fn )
	
	return( c( tp, tn, fp, fn, Recall, FPR, Accuracy, F1score, Precision ) )
}
   
# To compare the result according to the true graph
compare = function ( G, est, est2 = NULL, est3 = NULL, colnames = NULL, vis = FALSE ) 
{
	if( class(G)    == "sim" )      G    <- G $ G 
	if( class(G)    == "bdgraph" ){ es   <- select(G); G <- est $ G ; est <- es }
	if( class(est)  == "bdgraph" )  est  <- select( est ) 
	if( class(est2) == "bdgraph" )  est2 <- select( est2 ) 
	if( class(est3) == "bdgraph" )  est3 <- select( est3 ) 
	if( class(est)  == "select"  )  est  <- est $ refit
	if( class(est2) == "select"  )  est2 <- est2 $ refit
	if( class(est3) == "select"  )  est3 <- est3 $ refit

	G   = as.matrix(G)        # G is the adjacency matrix of true graph 
	est = as.matrix(est)      # est is the adjacency matrix of estimated graph 
	p   = nrow( G )
	G[ lower.tri( G, diag = TRUE ) ]     = 0
	est[ lower.tri( est, diag = TRUE ) ] = 0
   	   
	if( is.null( est2 ) & is.null( est3 ) )
	{
		result = matrix( 1, 9, 2 )
		result[ , 2 ] = ROC( G = G, est = est )
		if( is.null( colnames ) ) colnames <- c( "True graph", "estimate" )
	}

	if( !is.null( est2 ) & is.null( est3 ) )
	{
		est2 = as.matrix( est2 )       
		est2[ lower.tri( est2, diag = TRUE ) ] = 0

		result = matrix( 1, 9, 3 )
		result[ , 2 ] = ROC( G = G, est = est )
		result[ , 3 ] = ROC( G = G, est = est2 )
		if( is.null( colnames ) ) colnames <- c( "True graph", "estimate", "estimate2" )
	} 
	
	if( !is.null( est3 ) & !is.null( est3 ) )
	{
		est2 = as.matrix( est2 )       
		est2[ lower.tri( est2, diag = TRUE ) ] = 0
		est3 = as.matrix( est3 )       
		est3[ lower.tri( est3, diag = TRUE ) ] = 0

		result = matrix( 1, 9, 4 )
		result[ , 2 ] = ROC( G = G, est = est )
		result[ , 3 ] = ROC( G = G, est = est2 )
		result[ , 4 ] = ROC( G = G, est = est3 )
		if( is.null( colnames ) ) colnames <- c( "True graph", "estimate", "estimate2", "estimate3" )
	} 
	
   if( vis )
   {
		G_igraph   <- graph.adjacency( G,   mode = "undirected", diag = FALSE )
		est_igraph <- graph.adjacency( est, mode = "undirected", diag = FALSE )
		if ( p < 20 ) sizev = 15 else sizev = 2

		if( is.null( est2 ) )
		{
			op <- par( mfrow = c( 1, 2 ), pty = "s", omi = c( 0.3, 0.3, 0.3, 0.3 ), mai = c( 0.3, 0.3, 0.3, 0.3 ) )
			plot.igraph( G_igraph,   layout = layout.circle, main = colnames[1], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
			plot.igraph( est_igraph, layout = layout.circle, main = colnames[2], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
		}
		 
		if( !is.null( est2 ) & is.null( est3 ) )
		{
			op   <- par( mfrow = c( 2, 2 ), pty = "s", omi = c( 0.3, 0.3, 0.3, 0.3 ), mai = c( 0.3, 0.3, 0.3, 0.3 ) )
			est2_igraph <- graph.adjacency( as.matrix(est2), mode = "undirected", diag = FALSE )
			plot.igraph( G_igraph,    layout = layout.circle, main = colnames[1], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
			plot.igraph( est_igraph,  layout = layout.circle, main = colnames[2], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
			plot.igraph( est2_igraph, layout = layout.circle, main = colnames[3], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )			
		}
		
		if( !is.null( est2 ) & !is.null( est3 ) )
		{
			op   <- par( mfrow = c( 2, 2 ), pty = "s", omi = c( 0.3, 0.3, 0.3, 0.3 ), mai = c( 0.3, 0.3, 0.3, 0.3 ) )
			est2_igraph <- graph.adjacency( as.matrix(est2), mode = "undirected", diag = FALSE ) 
			est3_igraph <- graph.adjacency( as.matrix(est3), mode = "undirected", diag = FALSE ) 
			plot.igraph( G_igraph,    layout = layout.circle, main = colnames[1], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
			plot.igraph( est_igraph,  layout = layout.circle, main = colnames[2], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )
			plot.igraph( est2_igraph, layout = layout.circle, main = colnames[3], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )			
			plot.igraph( est3_igraph, layout = layout.circle, main = colnames[4], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )			
		}
		
		par(op)
   }
   
	result[ c( 3, 4, 6 ),1 ] = 0
	result[ 1, 1 ] = sum( G )
	result[ 2, 1 ] = p * ( p - 1 ) / 2 - result[ 1, 1 ]  

	result[ is.na( result ) ] = 0

	colnames( result ) <- colnames
	rownames( result ) <- c("true positive", "true negative", "false positive", "false negative", 
                "true positive rate", "false positive rate", "accuracy", "balanced F-score", "positive predictive")
				
	return( round( result, 3 ) )
}