##----------------------------------------------------------##
## Practical session: reconstruction of undirected networks ##
##                                                          ##
## Jean-Philippe Vert, 10/11/2010                           ##
##----------------------------------------------------------##

# This file is available at
# http://cbio.ensmp.fr/~jvert/svn/tutorials/practical/netinference/netinference.R
# The corresponding note is available at
# http://cbio.ensmp.fr/~jvert/svn/tutorials/practical/netinference/netinference_notes.pdf


## Data preparation ##
##------------------##

# Download the data from
# http://cbio.ensmp.fr/~jvert/svn/tutorials/practical/netinference/ppimetabo.tar.gz
# then tar xvzf ppimetabo.tar.gz

# Here we just focus on metabolic pathway and integrated data

# Load network
admat <- as.matrix(read.table('metabolic/gold_admat.txt'))
n <- dim(admat)[1]

# Load data (kernel)
K <- as.matrix(read.table('metabolic/Kmat_intg.txt'))

# We extract the list of interactions
xpos <- which(admat>0,arr.ind=TRUE)
xpos <- xpos[xpos[,1]>xpos[,2],]
npos <- dim(xpos)[1]

# Generate the same number of negative pairs
ineg <- sample(seq(admat)[-which(admat!=0)] , npos)
xneg <- cbind((ineg-1)%%n+1 , (ineg-1)%/%n+1)

# Create the full dataset
x <- rbind(xpos,xneg)
y <- c(rep(1,npos),rep(-1,npos))



## Engineer TPPK pairwise kernels ##
##--------------------------------##

# We construct a functions that computes the TPPK pairwise kernel

tppk <- function(preK=matrix())
# preK is a kernel matrix between individual. 
{
	rval <- function(x, y = NULL) {
		## we assume that x and y are just indices to be evaluated
		if (is.null(y)) {
			preK[x[1],x[1]]*preK[x[2],x[2]] + preK[x[1],x[2]]*preK[x[1],x[2]]
		} else {
			preK[x[1],y[1]]*preK[x[2],y[2]] + preK[x[1],y[2]]*preK[x[2],y[1]]
		}
	}
	return(new("kernel", .Data=rval, kpar=list(preK = preK)))
}


## Run SVM with TPPK kernel ##
##--------------------------##

# Make a TPPK kernel from the baseline kernel
pairK <- tppk(K)

# We make experiments on a small sample
ntrain <- 1000
ntest <- 1000
itrain <- sample(length(y),ntrain)
itest <- sample(seq(length(y))[-itrain],ntest)

m <- ksvm(x[itrain,],y[itrain],type="C-svc",kernel=pairK,scale=c(),C=1)
ypred <- predict(m,x[itest,],type="decision")

pred <- prediction(ypred,y[itest])
roc <- performance(pred, measure = "tpr", x.measure = "fpr")
pre <- performance(pred, measure = "prec", x.measure = "rec") 
plot(roc)
plot(pre)