## ------------------------------------------------------------------------------------------------|
#     Copyright (C) 2012 - 2018  Reza Mohammadi                                                    |
#                                                                                                  |
#     This file is part of BDgraph package.                                                        |
#                                                                                                  |
#     BDgraph is free software: you can redistribute it and/or modify it under                     |
#     the terms of the GNU General Public License as published by the Free                         |
#     Software Foundation; see <https://cran.r-project.org/web/licenses/GPL-3>.                    |
#                                                                                                  |
#     Maintainer: Reza Mohammadi <a.mohammadi@uva.nl>                                              |
## ------------------------------------------------------------------------------------------------|
#     This function reports below measures to assess the performance of estimated graphs:          |
## ------------------------------------------------------------------------------------------------|
# True positive:  number of correctly estimated links.                                             |
# True negative:  number of true non-existing links which is correctly estimated.                  |
# False positive: number of links which are not in the true graph, but are incorrectly estimated.  |
# False negative: number of links which they are in the true graph, but are not estimated.         |
# F1-score:       weighted average of the positive predictive and true positive rate.              |
# Specificity:    Specificity value reaches.                                                       |
# Sensitivity:    Sensitivity value reaches.                                                       |
# MCC:            Matthews Correlation Coefficients (MCC).                                         |
## ------------------------------------------------------------------------------------------------|

compare = function( sim.obj, bdgraph.obj, bdgraph.obj2 = NULL, bdgraph.obj3 = NULL, 
                    colnames = NULL, vis = FALSE ) 
{
	if( is.matrix( sim.obj ) )
	{
	    if( ( sum( sim.obj == 0 ) + sum( sim.obj == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'sim.obj' must be 0 or 1" )
	    G    = sim.obj
	}
    
	if( is.matrix( bdgraph.obj ) )
	{
	    if( ( sum( bdgraph.obj == 0 ) + sum( bdgraph.obj == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj' must be 0 or 1" )
	    est  = bdgraph.obj
	}
	
	if( !is.matrix( sim.obj )      && ( class( sim.obj )      == "sim"     ) )  G    <- sim.obj $ G 

	if( !is.matrix( bdgraph.obj )  && ( class( bdgraph.obj )  == "bdgraph" ) )  est  <- BDgraph::select( bdgraph.obj ) 
	if( !is.matrix( bdgraph.obj )  && ( class( bdgraph.obj )  == "ssgraph" ) )  est  <- BDgraph::select( bdgraph.obj ) 
	if( class( bdgraph.obj )  == "select"  )  est  <- bdgraph.obj $ 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 )
	
	if( sum( dim( G ) == dim( est ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj' have non-conforming size" )
	
	G[ lower.tri( G, diag = TRUE ) ]     = 0
	est[ lower.tri( est, diag = TRUE ) ] = 0
	   	   
	result = matrix( 1, 8, 2 )
	result[ , 2 ] = compute_measures( G = G, est_G = est )
   
    if( !is.null( bdgraph.obj2 ) )
    {
		if( is.matrix( bdgraph.obj2 ) )
		{
		    if( ( sum( bdgraph.obj2 == 0 ) + sum( bdgraph.obj2 == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj2' must be 0 or 1" )
		    est2 = bdgraph.obj2
		}
		
		if( !is.matrix( bdgraph.obj2 ) && ( class( bdgraph.obj2 ) == "bdgraph" ) )  est2 <- BDgraph::select( bdgraph.obj2 ) 
		if( !is.matrix( bdgraph.obj2 ) && ( class( bdgraph.obj2 ) == "ssgraph" ) )  est2 <- BDgraph::select( bdgraph.obj2 ) 
		if( class( bdgraph.obj2 ) == "select"  )  est2 <- bdgraph.obj2 $ refit

		est2 = as.matrix( est2 )
		
		if( sum( dim( G ) == dim( est2 ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj2' have non-conforming size" )
		
		est2[ lower.tri( est2, diag = TRUE ) ] = 0

		result = cbind( result, compute_measures( G = G, est_G = est2 ) )
	} 
	
    if( !is.null( bdgraph.obj3 ) )
    { 
		if( is.matrix( bdgraph.obj3 ) ) 
		{
		    if( ( sum( bdgraph.obj3 == 0 ) + sum( bdgraph.obj3 == 1 ) ) != ( p ^ 2 ) ) stop( "Element of 'bdgraph.obj3' must be 0 or 1" )
		    est3 = bdgraph.obj3
		}
		
		if( !is.matrix( bdgraph.obj3 ) && ( class( bdgraph.obj3 ) == "bdgraph" ) )  est3 <- BDgraph::select( bdgraph.obj3 ) 
		if( !is.matrix( bdgraph.obj3 ) && ( class( bdgraph.obj3 ) == "ssgraph" ) )  est3 <- BDgraph::select( bdgraph.obj3 ) 
		if( class( bdgraph.obj3 ) == "select"  )  est3 <- bdgraph.obj3 $ refit
		  
		est3 = as.matrix( est3 ) 
		
		if( sum( dim( G ) == dim( est3 ) ) != 2 ) stop( "'sim.obj' and 'bdgraph.obj3' have non-conforming size" )
		
		est3[ lower.tri( est3, diag = TRUE ) ] = 0

		result = cbind( result, compute_measures( G = G, est_G = est3 ) )
	} 
	
	result[ c( 3, 4 ), 1 ]    = 0
	result[ 1, 1 ]            = sum( G )
	result[ 2, 1 ]            = p * ( p - 1 ) / 2 - result[ 1, 1 ]  
	result[ is.na( result ) ] = 0

	if( is.null( colnames ) ) 
	{
		colnames = c( "True", "estimate" )
		if( !is.null( bdgraph.obj2 ) ) colnames = c( colnames, "estimate2" )
		if( !is.null( bdgraph.obj3 ) ) colnames = c( colnames, "estimate3" )
	}
		
	colnames( result ) <- colnames

	rownames( result ) <- c( "true positive", "true negative", "false positive", "false negative", 
                             "F1-score", "specificity", "sensitivity", "MCC" )
				
   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

		row_plot = ifelse( is.null( bdgraph.obj2 ), 1, 2 )
		op       = par( mfrow = c( row_plot, 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( bdgraph.obj2 ) )
		{
			est2_igraph <- graph.adjacency( as.matrix(est2), mode = "undirected", diag = FALSE )
			plot.igraph( est2_igraph, layout = layout.circle, main = colnames[3], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )			
		}
		
		if( !is.null( bdgraph.obj3 ) )
		{ 
			est3_igraph <- graph.adjacency( as.matrix(est3), mode = "undirected", diag = FALSE ) 
			plot.igraph( est3_igraph, layout = layout.circle, main = colnames[4], vertex.color = "white", vertex.size = sizev, vertex.label.color = 'black' )			
		}
		
		par( op )
   }
   
	return( round( result, 3 ) )
}
    
## ------------------------------------------------------------------------------------------------|
#    To compare measures the performance of estimated graphs based on true graph
## ------------------------------------------------------------------------------------------------|
compute_measures = function( G, est_G ) 
{
	upper_G     = G[     upper.tri( G     ) ]
	upper_est_G = est_G[ upper.tri( est_G ) ]
		
	tp = sum( ( upper_G != 0 ) * ( upper_est_G != 0 ) ) 
	tn = sum( ( upper_G == 0 ) * ( upper_est_G == 0 ) )
	fp = sum( ( upper_G == 0 ) * ( upper_est_G != 0 ) ) 
	fn = sum( ( upper_G != 0 ) * ( upper_est_G == 0 ) )
		
	# harmonic mean of precision and recall, called F-measure or balanced F-score
	F1score = ( 2 * tp ) / ( 2 * tp + fp + fn )

	specificity  = tn / ( tn + fp )
	sensitivity  = tp / ( tp + fn )
	# Matthews Correlation Coefficients (MCC)
	mcc          = ( ( tp * tn ) - ( fp * fn ) ) / ( sqrt( ( tp + fp ) * ( tp + fn ) ) * sqrt( ( tn + fp ) * ( tn + fn ) ) )
	
	return( c( tp, tn, fp, fn, F1score, specificity, sensitivity, mcc ) )
}