## This program 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, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.

# DISTANCE BETWEEN TWO NETWORKS, FOLLOWING APPROACH DEFINED BY Andrade et al (2008) 
# Andrade RFS, Miranda JGV, Pinho STR, Lobão TP (2008b) Measuring distances between complex networks. Phys Lett A 372: 5265–5269
# http://www.sciencedirect.com/science/article/pii/S0375960108009158
# AUTHOR: CHARLES NOVAES DE SANTANA, MAY 24th 2012
# 

library(igraph)

neighborhood_matrix<-function(net){
        nverts<-length(V(net));
        shortest_paths<-as.matrix(shortest.paths(net,mode="in",weights=E(net)$weight));
        for(i in 1:nverts){
                        inf_paths<-which(is.infinite(shortest_paths[i,]))
                        shortest_paths[i,inf_paths]<-0;
        }       
        return(shortest_paths);
}    

distance_between_networks<-function(net1,net2){
	m1<-neighborhood_matrix(net1);
	d1<-diameter(net1,directed=TRUE,weights=E(net1)$weight);
	m2<-neighborhood_matrix(net2);
	d2<-diameter(net2,directed=TRUE,weights=E(net2)$weight);
	nverts1<-length(V(net1));
	nverts2<-length(V(net2));
	distance<-0;
	if(nverts1 == nverts2){
		acummulator<-0;
		for(i in 1:nverts1){
			for(j in 1:nverts2){
				if((d1!=0)&&(d2!=0)){
					acummulator<-acummulator+( ( (m1[i,j]/d1)-(m2[i,j]/d2) )^2 )
				}
				if(d1==0){
					acummulator<-acummulator+( ( 0-(m2[i,j]/d2) )^2 )
				}	
				if(d2==0){
					acummulator<-acummulator+( ( (m1[i,j]/d1)-0 )^2 )
				}
				if((d1==0)&&(d2==0)){
					acummulator<-0;
				}
			}
		}
		distance<-ifelse(acummulator,acummulator/(nverts1*(nverts1-1)),0);
		distance<-sqrt(distance);	
	}
	if(nverts1 != nverts2)
	{
		distance<-NaN;
	}
	return(distance);
}