# March 2016	  
# April 2016: all high level function modified: some input parameters of find.best.f function added
# April 2016: all high level function modified: normaliation method added
# July 2016: added PRC computed by precrec package (only in holdout)
# July 2016: added Do.tpr.'tpr-variant'.holdout

#*****************************************************************************************#
# libraries and source files to be loaded: 
 library(PerfMeas)  			   ## compute AUROC and precision at fixed recall rate
 library(precrec)				   ## compute AUPRC
 library(preprocessCore)		   ## Qnorm
 source("flat.score.norm.R")  	   ## Maxnorm
 source("graph.utils.R")  		   ## graph utility functions 
 source("F-hier.R")				   ## compute Kiritchenko-like multi-label F-scores
 source("Do.flat.normalization.R") ## high level function to compute Maxnorm e Qnorm
#*****************************************************************************************#

# High level function to compute hierarchical correction according to TPR-threshold algorithm. 
# It perform a k fold cross-validation to find the best threshold maximizing on F.max measure.
# INPUT: 
# threshold: range of threshold values to be tested in order to find the best threshold. The denser the range is, the higher the probability to find the best
# 			 theshold, but the execution time will be increasing (def: from:0.1, to:0.9, by:0.1 step).
# kk:  number of folds of the cross validation (def: 5).
# seed:  intialization seed for the random generator to create folds (def:0). If NULL (default) no initialization is performed.
# norm: boolean value: 1.TRUE means that the flat scores matrix has been already normalized in according to a normalization method (def);
#       			   2.FALSE means that the flat scores matrix has NOT been normalized yet.
# norm.type: this variable can assume three values: MaxNorm, Qnorm, NONE. We have two case respect to norm:
#			 1.if norm==FALSE, two kind of normalizations are possible:	1. MaxNorm: each score is divided w.r.t. the max of each class 
#																		2. Qnorm: quantile normalization is applied. Library preprocessCore is used.
#			 2.if norm==TRUE, set norm.type=="NONE" (def);
# flat.file: name of flat scores matrix already normalized or to be normalized in according to norm.type (without rda extension); 
# ann.file: name of the target labels (without rda extension). It must be  an .rda file containing the label matrix of the examples (def: ann.file)
# dag.file: name of the graph that represents the hierarchy of the classes. 
# flat.dir: relative path to folder where flat normalized scores matrix is stored
# ann.dir: relative path to folder where annotation matrix is stored
# dag.dir: relative path to folder where graph is stored
# flat.norm.dir: 1.if norm=FALSE, relative path where flat normalized scores matrix is strored;
#				 2.if norm=TRUE, the flat scores matrix is already normalized, than it is set to NULL (def)
# n.round: 	number of rounding digits to be applied to the hierarchical scores matrix (def. 3). 
#			It's used for choosing the best threshold on the basis of the best F.measure (see f.criterion parameter).
# f.criterion: character. Type of F-measure to be used to select the best F.measure. There are 2 possibilities: 
#              1. "F" (default) corresponds to the harmonic mean between the average precision and recall; 
#              2. "avF" corresponds to the per-example F-score averaged across all the examples.
# hierScore.dir: relative path to folder where the matrix with the scores of the classes corrected according to TPR-threshold algorithm is stored 
# macro.dir: relative path to folder where the macro-centric measures (i.e. AUC and PxR across classes) are stored
# Fmeas.dir: relative path to folder where example-centric measures (i.e. Precision, Recall, Specificity, F-measure, Accuracy across example) are stored
# OUTPUT:
# 4 rda files stored in the corresponding output directory:
# - Matrix with examples on rows and classes on colums representing the hierarchical scores of the classes computed with TPR-threshold algorithm.
#	Stored in hierScore.dir folder
# - Example-centric measures computed through find.best.f from F-hier.R file. Stored in macro.dir folder
# - AUC (average and per classes) computed through AUC.single.over.classes from package PerfMeas. Stored in macro.dir
# - Precision at fixed recall levels (average and per classes) computed through precision.at.multiple.recall.level.over.classes from package PerfMeas.
# 	Stored in macro.dir folder 
Do.tpr.threshold.cv <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), kk=5, seed=NULL, norm=TRUE, norm.type= "NONE",
									flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
									flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir,	flat.norm.dir=NULL, 
									n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
								){
	## Loading Data ############
	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	ann.no.root <- hpo.ann[,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## Hierarchical Correction #########################
	## TPR-Threshold correction using k-fold cross-validation to find best threshold 

	## splitting data in 5 fold. Each fold is a sub matrix. Rows are a random sub-sampling of a set of genes; columns are all HPO terms
	folds <- do.unstratified.cv.data(S.flat, kk=kk, seed=seed);

	## storing the best F.max and the best threshold for each of k training set
	training.top.Fmax <- vector(mode="list", length=kk);
	names(training.top.Fmax) <- paste0(rep("fold",kk), 1:kk);

	## storing hierarchical scores for each of k test set
	test.hier.score <- vector(mode="list", length=kk);
	names(test.hier.score) <- paste0(rep("fold",kk), 1:kk);

	## storing all the best protein-centric measures for each of k test set 
	best.avg.meas.test <- vector(mode="list", length=kk);
	names(best.avg.meas.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average macro AUC for each of k test set 
	AUC.average.test <- vector(mode="list", length=kk);
	names(AUC.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class macro AUC for each of k test set 
	AUC.class.test <- vector(mode="list", length=kk);
	names(AUC.class.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average PxR for each of k test set 
	PXR.average.test <- vector(mode="list", length=kk);
	names(PXR.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class PxR for each of k test set 
	PXR.class.test <- vector(mode="list", length=kk);
	names(PXR.class.test) <- paste0(rep("fold",kk), 1:kk);

	## variable for hosting the k sub-matrix assembled 
	S.tpr <- c();

	## Let's start k-fold crossing validation for choosing best threshold maximizing on F.max...
	for(k in 1:kk){
		test <- S.flat[folds[[k]],];		# test set: 1 out of k folds (e.g.: if k=5, testing set is 1/5 of data)
		training <- S.flat[!rownames(S.flat) %in% rownames(test),];		# training set: (k-1) out of k folds (e.g.: if k=5, training set is 4/5 of data)
		gc();

		top.Fmax <- 0;
		best.Fmaxt <- 0;
		for(t in threshold){
			training.hier.score <- tpr.threshold(training, hpo, root=root, t=t);	## training set hierarchical correction...
			pred.training <- training.hier.score;
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxt <- t;
				training.top.Fmax[[k]] <- c(best.Fmax=top.Fmax, best.thres=best.Fmaxt);
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, sep="\t", "\n");
			}
		}
		test.hier.score[[k]] <- tpr.threshold(test, hpo, root=root, t=training.top.Fmax[[k]][2]);	## testing set hierarchical correction...
		pred.test <- test.hier.score[[k]];
		target.test <- hpo.ann[rownames(pred.test),colnames(pred.test)];
		test.avg.meas <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## testing set evaluating...
		best.avg.meas.test[[k]] <- test.avg.meas;

		## AUC (average and per.class) computed with PerfMeas package
		AUC.average.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$average;	## storing average AUC of each k testing set
		AUC.class.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$per.class;	## storing AUC per.class of each k testing set

		## PxR at fixed recall levels (average and per.class) computed with PerfMeas package
		PXR.average.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$avgPXR;	## storing average PxR of each k testing set
		PXR.class.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$PXR;	## storing PxR per.class of each k testing set

		## assembling of all the hierarchical scores of each k sub-matrix to build the full matrix 
		S.tpr <- rbind(S.tpr,test.hier.score[[k]]);
	}
	## put the rows (i.e. genes) of assembled k sub-matrix in the same order of the full beginning matrix
	S.tpr <- S.tpr[rownames(S.flat),];

	## Averaging Performances Measures across k testing set ###################
	## averaging protein-centric measures across k testing sets 
	F.meas <- Reduce("+", best.avg.meas.test)/kk;

	## averaging macro-AUC (average and per.class) across k testing sets 
	AUC.average.tpr.over.test <- Reduce("+", AUC.average.test)/kk;
	AUC.class.tpr.over.test <- Reduce("+", AUC.class.test)/kk;
	AUC.tpr <- list(average=AUC.average.tpr.over.test, per.class=AUC.class.tpr.over.test);

	## averaging PxR (average and per.class) across k testing sets 
	PXR.average.tpr.over.test <- Reduce("+", PXR.average.test)/kk;
	PXR.class.tpr.over.test <- Reduce("+", PXR.class.test)/kk;
	PXR.tpr <- list(average=PXR.average.tpr.over.test, per.class=PXR.class.tpr.over.test);

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
	}
}

# TPR-threshold Hold-Out version
# High level function to correct the scores with a hierarchy according to TPR-threshold algorithm performing a classical holdout procedure
# All input parameters are the same of the above corresponding function that implement a kfcv, except the following:
# ind.test.set: vector of integer. Indices refer to the examples of the adjancency matrix to be used in the test set. 
# ind.dir: relative path to folder where ind.test.set is stored
# An output file is added: PRC (average and per.class) computed by precrec package. It assumed to store in macro.dir
Do.tpr.threshold.holdout <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), kk=5, seed=NULL, norm=TRUE, norm.type= "NONE",
										flat.file=flat.file, ann.file=ann.file, dag.file=dag.file, ind.test.set=ind.test.set,
										ind.dir=ind.dir,flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir, flat.norm.dir=NULL, 
										n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", 
										Fmeas.dir="Fmeas.dir/"
									){
	## Loading Data ############
	## loading examples indices of the test set
	ind.set <- paste0(ind.dir, ind.test.set, ".rda");
	ind.test <- get(load(ind.set));

	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## scores flat matrix shrinked to test and training test respectively
	S.test <- S.flat[ind.test,];
	S.training <- S.flat[-ind.test,];
	S.flat <- S.test;

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table and shrinking annotation table to tet set
	ann.no.root <- hpo.ann[ind.test,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold correction using k-fold cross-validation to find best threshold on training set
	folds <- do.unstratified.cv.data(S.training, kk=kk, seed=seed);
			
	## Let's start k-fold crossing validation for choosing best threshold maximizing on the basis of F.max value...
	for(k in 1:kk){
		training <- S.training[folds[[k]],];	
		
		top.Fmax <- 0;
		best.Fmaxt <- 0;
		for(t in threshold){
			pred.training <- tpr.threshold(training, hpo, root=root, t=t);	## training set hierarchical correction...
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxt <- t;
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, sep="\t", "\n");
			}
		}
	}
	S.tpr <- pred.test <- tpr.threshold(S.test, hpo, root=root, t=best.Fmaxt);	## testing set hierarchical correction..
	target.test <- ann.no.root;
	
	## AUC (average and per.class) computed with PerfMeas package on the test set
	AUC.tpr <- AUC.single.over.classes(target.test, pred.test, hpo, root)[c("average","per.class")];	

	## PxR at fixed recall levels (average and per.class) computed with PerfMeas package on the test set
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(target.test, pred.test);

	## F.measure: Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprT.rda"), compress=TRUE);
	}
}


# High level function to compute hierarchical correction according to TPR-weighted-threshold-free algorithm. 
# It perform a k fold cross-validation to find the best weight maximizing on F.max measure.
# INPUT: 
# weight: range of weight values to be tested in order to find the best weight. The denser the range is, the higher the probability to find the best
# 		  weigth, but the execution time will be increasing (def: from:0.1, to:1, by:0.1 step).
# kk:  number of folds of the cross validation (def: 5).
# seed:  intialization seed for the random generator to create folds (def:0). If NULL (default) no initialization is performed.
# norm: boolean value: 1.TRUE means that the flat scores matrix has been already normalized in according to a normalization method (def);
#       			   2.FALSE means that the flat scores matrix has NOT been normalized yet.
# norm.type: this variable can assume three values: MaxNorm, Qnorm, NONE. We have two case respect to norm:
#			 1.if norm==FALSE, two kind of normalizations are possible:	1. MaxNorm: each score is divided w.r.t. the max of each class 
#																		2. Qnorm: quantile normalization is applied. Library preprocessCore is used.
#			 2.if norm==TRUE, set norm.type=="NONE" (def);
# flat.file: name of flat scores matrix already normalized or to be normalized in according to norm.type (without rda extension); 
# ann.file: name of the target labels (without rda extension). It must be  an .rda file containing the label matrix of the examples (def: ann.file)
# dag.file: name of the graph that represents the hierarchy of the classes. 
# flat.dir: relative path to folder where flat normalized scores matrix is stored
# ann.dir: relative path to folder where annotation matrix is stored
# dag.dir: relative path to folder where graph is stored
# flat.norm.dir: 1.if norm=FALSE, relative path where flat normalized scores matrix is strored;
#				 2.if norm=TRUE, the flat scores matrix is already normalized, than it is set to NULL (def)
# n.round: 	number of rounding digits to be applied to the hierarchical scores matrix (def. 3). 
#			It's used for choosing the best threshold on the basis of the best F.measure (see f.criterion parameter).
# f.criterion: character. Type of F-measure to be used to select the best F.measure. There are 2 possibilities: 
#              1. "F" (default) corresponds to the harmonic mean between the average precision and recall; 
#              2. "avF" corresponds to the per-example F-score averaged across all the examples.
# hierScore.dir: relative path to folder where the matrix with the scores of the classes corrected in according to TPR-weighted-threshold-free algorithm is stored 
# macro.dir: relative path to folder where the macro-centric measures (i.e. AUC and PxR across classes) are stored
# Fmeas.dir: relative path to folder where example-centric measures (i.e. Precision, Recall, Specificity, F-measure, Accuracy across example) are stored
# OUTPUT:
# 4 rda files stored in the rispective output directory:
# - Matrix with examples on rows and classes on colums representing the hierarchical scores of the classes computed with TPR-weighted-threshold-free algorithm.
#	Stored in hierScore.dir folder
# - Example-centric measures computed through find.best.f from F-hier.R file. Stored in macro.dir folder
# - AUC (average and per classes) computed through AUC.single.over.classes from package PerfMeas. Stored in macro.dir
# - Precision at fixed recall levels (average and per classes) computed through precision.at.multiple.recall.level.over.classes from package PerfMeas.
# 	Stored in macro.dir folder 
Do.tpr.weighted.threshold.free.cv <- function(	weight=seq(from=0.1, to=1, by=0.1), kk=5, seed=NULL, norm=TRUE, norm.type= "NONE",
												flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
												flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir,	flat.norm.dir=NULL, 
												n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
											){
	## Loading Data ############
	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	root <- root.node(hpo);
	ann.no.root <- hpo.ann[,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## Hierarchical Correction #########################
	## TPR-Weighted-Threshold-Free correction using k-fold cross-validation to find best weight 

	## splitting data in 5 fold. Each fold is a sub matrix. Rows are a random sub-sampling of a set of genes; columns are all HPO terms
	folds <- do.unstratified.cv.data(S.flat, kk=kk, seed=seed);

	## storing the best F.max and the best threshold for each of k training set
	training.top.Fmax <- vector(mode="list", length=kk);
	names(training.top.Fmax) <- paste0(rep("fold",kk), 1:kk);

	## storing hierarchical scores for each of k test set
	test.hier.score <- vector(mode="list", length=kk);
	names(test.hier.score) <- paste0(rep("fold",kk), 1:kk);

	## storing all the best protein-centric measures for each of k test set 
	best.avg.meas.test <- vector(mode="list", length=kk);
	names(best.avg.meas.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average macro AUC for each of k test set 
	AUC.average.test <- vector(mode="list", length=kk);
	names(AUC.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class macro AUC for each of k test set 
	AUC.class.test <- vector(mode="list", length=kk);
	names(AUC.class.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average PxR for each of k test set 
	PXR.average.test <- vector(mode="list", length=kk);
	names(PXR.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class PxR for each of k test set 
	PXR.class.test <- vector(mode="list", length=kk);
	names(PXR.class.test) <- paste0(rep("fold",kk), 1:kk);

	## variable for hosting the k sub-matrix assembled 
	S.tpr <- c();

	## Let's start k-fold crossing validation for choosing best threshold maximizing on F.max...
	for(k in 1:kk){
		test <- S.flat[folds[[k]],];		# test set: 1 out of k folds (e.g.: if k=5, testing set is 1/5 of data)
		training <- S.flat[!rownames(S.flat) %in% rownames(test),];		# training set: (k-1) out of k folds (e.g.: if k=5, training set is 4/5 of data)
		gc();

		top.Fmax <- 0;
		best.Fmaxw <- 0;
		for(w in weight){
			training.hier.score <- tpr.weighted.threshold.free(training, hpo, root=root, w=w);	## training set hierarchical correction...
			pred.training <- training.hier.score;
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxw <- w;
				training.top.Fmax[[k]] <- c(best.Fmax=top.Fmax, best.weigh=best.Fmaxw);
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best weight:",best.Fmaxw, sep="\t", "\n");
			}
		}
		test.hier.score[[k]] <- tpr.weighted.threshold.free(test, hpo, root=root, w=training.top.Fmax[[k]][2]);	## testing set hierarchical correction...
		pred.test <- test.hier.score[[k]];
		target.test <- hpo.ann[rownames(pred.test),colnames(pred.test)];
		test.avg.meas <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## testing set evaluating...
		best.avg.meas.test[[k]] <- test.avg.meas;

		## Hierarchical AUC (average and per.class) computed with PerfMeas package
		AUC.average.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$average;	## storing average AUC of each k testing set
		AUC.class.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$per.class;	## storing AUC per.class of each k testing set

		## Hierarchical PxR at fixed recall levels (average and per.class) computed with PerfMeas package
		PXR.average.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$avgPXR;	## storing average PxR of each k testing set
		PXR.class.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$PXR;	## storing PxR per.class of each k testing set

		## assembling of all the hierarchical scores of each k sub-matrix to build the full matrix 
		S.tpr <- rbind(S.tpr,test.hier.score[[k]]);
	}
	## put the rows (i.e. genes) of assembled k sub-matrix in the same order of the full beginning matrix
	S.tpr <- S.tpr[rownames(S.flat),];

	## Averaging Performances Measures across k testing set ###################
	## averaging protein-centric measures across k testing sets 
	F.meas <- Reduce("+", best.avg.meas.test)/kk;

	## averaging macro-AUC (average and per.class) across k testing sets 
	AUC.average.tpr.over.test <- Reduce("+", AUC.average.test)/kk;
	AUC.class.tpr.over.test <- Reduce("+", AUC.class.test)/kk;
	AUC.tpr <- list(average=AUC.average.tpr.over.test, per.class=AUC.class.tpr.over.test);

	## averaging PxR (average and per.class) across k testing sets 
	PXR.average.tpr.over.test <- Reduce("+", PXR.average.test)/kk;
	PXR.class.tpr.over.test <- Reduce("+", PXR.class.test)/kk;
	PXR.tpr <- list(average=PXR.average.tpr.over.test, per.class=PXR.class.tpr.over.test);

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.",flat.file, ".hierScores.tprW.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.",norm.type,".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.",norm.type,".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", norm.type,".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
	}
}

# High level function to correct the scores with a hierarchy according to TPR-weighted-threshold-free algorithm performing a classical holdout procedure
# All input parameters are the same of the above corresponding function that implement a kfcv, except the following:
# ind.test.set: vector of integer. Indices refer to the examples of the adjancency matrix to be used in the test set. 
# ind.dir: relative path to folder where ind.test.set is stored
# An output file is added: PRC (average and per.class) computed by precrec package. It assumed to store in macro.dir
Do.tpr.weighted.threshold.free.holdout <- function(	weight=seq(from=0.1, to=1, by=0.1), kk=5, seed=NULL, norm=TRUE, norm.type= "NONE",
													flat.file=flat.file, ann.file=ann.file, dag.file=dag.file, ind.test.set=ind.test.set,
													ind.dir=ind.dir, flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir,	 
													flat.norm.dir=NULL, n.round=3, f.criterion ="F", 
													hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
												){
	## Loading Data ############
	## loading examples indices of the test set
	ind.set <- paste0(ind.dir, ind.test.set, ".rda");
	ind.test <- get(load(ind.set));

	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 																		 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## scores flat matrix shrinked to test and training test respectively
	S.test <- S.flat[ind.test,];
	S.training <- S.flat[-ind.test,];
	S.flat <- S.test;

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table and shrinking annotation table to tet set
	ann.no.root <- hpo.ann[ind.test,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold correction using k-fold cross-validation to find best threshold on training set
	folds <- do.unstratified.cv.data(S.training, kk=kk, seed=seed);
			
	## Let's start k-fold crossing validation for choosing best threshold maximizing on the basis of F.max value...
	for(k in 1:kk){
		training <- S.training[folds[[k]],];	
		
		top.Fmax <- 0;
		best.Fmaxw <- 0;
		for(w in weight){
			pred.training <- tpr.weighted.threshold.free(training, hpo, root=root, w=w);	## training set hierarchical correction...
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxw <- w;
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxw, sep="\t", "\n");
			}
		}
	}
	S.tpr <- pred.test <- tpr.weighted.threshold.free(S.test, hpo, root=root, w=best.Fmaxw);	## testing set hierarchical correction..
	target.test <- ann.no.root;
	
	## AUC (average and per.class) computed with PerfMeas package on the test set
	AUC.tpr <- AUC.single.over.classes(target.test, pred.test, hpo, root)[c("average","per.class")];	

	## PxR at fixed recall levels (average and per.class) computed with PerfMeas package on the test set
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(target.test, pred.test);

	## F.measure: Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprW.rda"), compress=TRUE);
	}
}


# High level function to compute hierarchical correction according to TPR-weighted-threshold algorithm. 
# It performs a k fold cross-validation to find the best threshold and weight maximizing the F.max measure.
# INPUT: 
# threshold: range of threshold values to be tested in order to find the best threshold. The denser the range is, the higher the probability to find the best
# 			 theshold, but the execution time will be increasing (def: from:0.1, to:0.9, by:0.1 step).
# weight: range of weight values to be tested in order to find the best weight. The denser the range is, the higher the probability to find the best
# 		  weigth, but the execution time will be increasing (def: from:0.1, to:1, by:0.1 step).
# kk:  number of folds of the cross validation (def: 5).
# seed:  intialization seed for the random generator to create folds (def:0). If NULL (default) no initialization is performed.
# norm: boolean value: 1.TRUE means that the flat scores matrix has been already normalized in according to a normalization method (def);
#       			   2.FALSE means that the flat scores matrix has NOT been normalized yet.
# norm.type: this variable can assume three values: MaxNorm, Qnorm, NONE. We have two case respect to norm:
#			 1.if norm==FALSE, two kind of normalizations are possible:	1. MaxNorm: each score is divided w.r.t. the max of each class 
#																		2. Qnorm: quantile normalization is applied. Library preprocessCore is used.
#			 2.if norm==TRUE, set norm.type=="NONE" (def);
# flat.file: name of flat scores matrix already normalized or to be normalized in according to norm.type (without rda extension); 
# ann.file: name of the target labels (without rda extension). It must be  an .rda file containing the label matrix of the examples (def: ann.file)
# dag.file: name of the graph that represents the hierarchy of the classes. 
# flat.dir: relative path to folder where flat normalized scores matrix is stored
# ann.dir: relative path to folder where annotation matrix is stored
# dag.dir: relative path to folder where graph is stored
# flat.norm.dir: 1.if norm=FALSE, relative path where flat normalized scores matrix is strored;
#				 2.if norm=TRUE, the flat scores matrix is already normalized, than it is set to NULL (def)
# n.round: 	number of rounding digits to be applied to the hierarchical scores matrix (def. 3). 
#			It's used for choosing the best threshold on the basis of the best F.measure (see f.criterion parameter).
# f.criterion: character. Type of F-measure to be used to select the best F.measure. There are 2 possibilities: 
#              1. "F" (default) corresponds to the harmonic mean between the average precision and recall; 
#              2. "avF" corresponds to the per-example F-score averaged across all the examples.
# hierScore.dir: relative path to folder where the matrix with the scores of the classes corrected in according to TPR-weighted-threshold algorithm is stored 
# macro.dir: relative path to folder where the macro-centric measures (i.e. AUC and PxR across classes) are stored
# Fmeas.dir: relative path to folder where example-centric measures (i.e. Precision, Recall, Specificity, F-measure, Accuracy across example) are stored
# OUTPUT:
# 4 rda files stored in the rispective output directory:
# - Matrix with examples on rows and classes on colums representing the hierarchical scores of the classes computed with TPR-weighted-threshold algorithm.
#	Stored in hierScore.dir folder
# - Example-centric measures computed through find.best.f from F-hier.R file. Stored in macro.dir folder
# - AUC (average and per classes) computed through AUC.single.over.classes from package PerfMeas. Stored in macro.dir
# - Precision at fixed recall levels (average and per classes) computed through precision.at.multiple.recall.level.over.classes from package PerfMeas.
# 	Stored in macro.dir folder 
Do.tpr.weighted.threshold.cv <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), weight=seq(from=0.1, to=1, by=0.1), 
											kk=5, seed=NULL,  norm=TRUE, norm.type= "NONE",
											flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
											flat.dir=flat.dir, ann.dir=ann.dir, dag.dir=dag.dir, flat.norm.dir=NULL,
											n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
										){
	## Loading Data ############
	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	root <- root.node(hpo);
	ann.no.root <- hpo.ann[,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## Hierarchical Correction #########################
	## TPR-Weighted-Threshold correction using k-fold cross-validation to find best threshold and weight 

	## splitting data in 5 fold. Each fold is a sub matrix. Rows are a random sub-sampling of a set of genes; columns are all HPO terms
	folds <- do.unstratified.cv.data(S.flat, kk=kk, seed=seed);

	## storing the best F.max and the best threshold and weight for each of k training set
	training.top.Fmax <- vector(mode="list", length=kk);
	names(training.top.Fmax) <- paste0(rep("fold",kk), 1:kk);

	## storing hierarchical scores for each of k test set
	test.hier.score <- vector(mode="list", length=kk);
	names(test.hier.score) <- paste0(rep("fold",kk), 1:kk);

	## storing all the best protein-centric measures for each of k test set 
	best.avg.meas.test <- vector(mode="list", length=kk);
	names(best.avg.meas.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average macro AUC for each of k test set 
	AUC.average.test <- vector(mode="list", length=kk);
	names(AUC.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class macro AUC for each of k test set 
	AUC.class.test <- vector(mode="list", length=kk);
	names(AUC.class.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average PxR for each of k test set 
	PXR.average.test <- vector(mode="list", length=kk);
	names(PXR.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class PxR for each of k test set 
	PXR.class.test <- vector(mode="list", length=kk);
	names(PXR.class.test) <- paste0(rep("fold",kk), 1:kk);

	## variable for hosting the k sub-matrix assembled 
	S.tpr <- c();

	## Let's start k-fold crossing validation for choosing best threshold and weight maximizing on F.max...
	for(k in 1:kk){
		test <- S.flat[folds[[k]],];	# test set: 1 out of k folds (e.g.: if k=5, testing set is 1/5 of data)
		training <- S.flat[!rownames(S.flat) %in% rownames(test),];		# training set: (k-1) out of k folds (e.g.: if k=5, training set is 4/5 of data)
		gc();

		top.Fmax <- 0;
		best.Fmaxt <- 0;
		best.Fmaxw <- 0;
		for(t in threshold){
			for(w in weight){
				training.hier.score <- tpr.weighted.threshold(training, hpo, root=root, w=w, t=t);	## training set hierarchical correction...
				pred.training <- training.hier.score;
				target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
				training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## training evaluating... 
				training.Fmax <- training.avg.meas[4];	## F.max maximization...
				if(training.Fmax > top.Fmax){
					top.Fmax <- training.Fmax;
					best.Fmaxt <- t;
					best.Fmaxw <- w;
					training.top.Fmax[[k]] <- c(best.Fmax=top.Fmax, best.thres=best.Fmaxt, best.weight=best.Fmaxw);
					cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, "best weight:", best.Fmaxw,sep="\t", "\n");
				}
			}
		}
		test.hier.score[[k]] <- tpr.weighted.threshold(test, hpo, root=root, t=training.top.Fmax[[k]][2], w=training.top.Fmax[[k]][3]);	## testing set hierarchical correction...
		pred.test <- test.hier.score[[k]];
		target.test <- hpo.ann[rownames(pred.test),colnames(pred.test)];
		test.avg.meas <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## testing set evaluating...
		best.avg.meas.test[[k]] <- test.avg.meas;

		## Hierarchical AUC (average and per.class) computed with PerfMeas package
		AUC.average.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$average;	## storing average AUC of each k testing set
		AUC.class.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$per.class;	## storing AUC per.class of each k testing set

		## Hierarchical PxR at fixed recall levels (average and per.class) computed with PerfMeas package
		PXR.average.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$avgPXR;	## storing average PxR of each k testing set
		PXR.class.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$PXR;	## storing PxR per.class of each k testing set

		## assembling of all the hierarchical scores of each k sub-matrix to build the full matrix 
		S.tpr <- rbind(S.tpr,test.hier.score[[k]]);
	}
	## put the rows (i.e. genes) of assembled k sub-matrix in the same order of the full beginning matrix
	S.tpr <- S.tpr[rownames(S.flat),];

	## Averaging Performances Measures across k testing set ###################
	## averaging protein-centric measures across k testing sets 
	F.meas <- Reduce("+", best.avg.meas.test)/kk;

	## averaging macro-AUC (average and per.class) across k testing sets 
	AUC.average.tpr.over.test <- Reduce("+", AUC.average.test)/kk;
	AUC.class.tpr.over.test <- Reduce("+", AUC.class.test)/kk;
	AUC.tpr <- list(average=AUC.average.tpr.over.test, per.class=AUC.class.tpr.over.test);

	## averaging PxR (average and per.class) across k testing sets 
	PXR.average.tpr.over.test <- Reduce("+", PXR.average.test)/kk;
	PXR.class.tpr.over.test <- Reduce("+", PXR.class.test)/kk;
	PXR.tpr <- list(average=PXR.average.tpr.over.test, per.class=PXR.class.tpr.over.test);

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type,".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type,".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", norm.type,".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
	}
}

# Hold-Out version
# High level function to correct the scores with a hierarchy according to TPR-weighted-threshold algorithm performing a classical holdout procedure
# All input parameters are the same of the corresponding function above, except the following:
# ind.test.set: vector of integer. Indices refer to the examples of the adjancency matrix to be used in the test set. 
# ind.dir: relative path to folder where ind.test.set is stored
# An output file is added: PRC (average and per.class) computed by precrec package. Stored in macro.dir
Do.tpr.weighted.threshold.holdout <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), weight=seq(from=0.1, to=1, by=0.1),
												kk=5, seed=NULL, norm=TRUE, norm.type= "NONE", flat.file=flat.file, ann.file=ann.file, 
												dag.file=dag.file, ind.test.set=ind.test.set, ind.dir=ind.dir, flat.dir=flat.dir, 
												ann.dir=ann.dir,dag.dir=dag.dir, flat.norm.dir=NULL, n.round=3, f.criterion ="F", 
												hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"												
											){
	## Loading Data ############
	## loading examples indices of the test set
	ind.set <- paste0(ind.dir, ind.test.set, ".rda");
	ind.test <- get(load(ind.set));

	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 																		 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## scores flat matrix shrinked to test and training test respectively
	S.test <- S.flat[ind.test,];
	S.training <- S.flat[-ind.test,];
	S.flat <- S.test;

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table and shrinking annotation table to tet set
	ann.no.root <- hpo.ann[ind.test,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold correction using k-fold cross-validation to find best threshold on training set
	folds <- do.unstratified.cv.data(S.training, kk=kk, seed=seed);
			
	## Let's start k-fold crossing validation for choosing best threshold maximizing on the basis of F.max value...
	for(k in 1:kk){
		training <- S.training[folds[[k]],];	
		
		top.Fmax <- 0;
		best.Fmaxt <- 0;
		best.Fmaxw <- 0;
		for(t in threshold){
			for(w in weight){
				pred.training <- tpr.weighted.threshold(training, hpo, root=root, w=w, t=t); ## training set hierarchical correction...
				target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
				training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	## training evaluating... 
				training.Fmax <- training.avg.meas[4];	## F.max maximization...
				if(training.Fmax > top.Fmax){
					top.Fmax <- training.Fmax;
					best.Fmaxt <- t;
					best.Fmaxw <- w;
					cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, "best weight:",best.Fmaxw, sep="\t","\n");
				}
			}
		}
	}
	S.tpr <- pred.test <- tpr.weighted.threshold(S.test, hpo, root=root, t=best.Fmaxt, w=best.Fmaxw);	## testing set hierarchical correction..
	target.test <- ann.no.root;
	
	## AUC (average and per.class) computed with PerfMeas package on the test set
	AUC.tpr <- AUC.single.over.classes(target.test, pred.test, hpo, root)[c("average","per.class")];	

	## PxR at fixed recall levels (average and per.class) computed with PerfMeas package on the test set
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(target.test, pred.test);

	## F.measure: Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprWT.rda"), compress=TRUE);
	}
}


# High level function to compute hierarchical correction according to TPR-threshold-free algorithm. 
# INPUT: 
# norm: boolean value: 1.TRUE means that the flat scores matrix has been already normalized in according to a normalization method (def);
#       			   2.FALSE means that the flat scores matrix has NOT been normalized yet.
# norm.type: this variable can assume three values: MaxNorm, Qnorm, NONE. We have two case respect to norm:
#			 1.if norm==FALSE, two kind of normalizations are possible:	1. MaxNorm: each score is divided w.r.t. the max of each class 
#																		2. Qnorm: quantile normalization is applied. Library preprocessCore is used.
#			 2.if norm==TRUE, set norm.type=="NONE" (def);
# flat.file: name of flat scores matrix already normalized or to be normalized in according to norm.type (without rda extension);
# ann.file: name of the target labels (without rda extension). It must be  an .rda file containing the label matrix of the examples (def: ann.file)
# dag.file: name of the graph that represents the hierarchy of the classes (def dag.file) 
# flat.dir: relative path to folder where flat scores matrix (already normalized or to normalize) is stored (def flat.dir)
# ann.dir: relative path to folder where annotation matrix is stored
# dag.dir: relative path to folder where graph is stored
# flat.norm.dir: 1.if norm=FALSE, relative path where flat normalized scores matrix is strored;
#				 2.if norm=TRUE, the flat scores matrix is already normalized, than it is set to NULL (def)
# n.round: 	number of rounding digits to be applied to the hierarchical scores matrix (def. 3). 
#			It's used for choosing the best threshold on the basis of the best F.measure (see f.criterion parameter).
# f.criterion: character. Type of F-measure to be used to select the best F.measure. There are 2 possibilities: 
#              1. "F" (default) corresponds to the harmonic mean between the average precision and recall; 
#              2. "avF" corresponds to the per-example F-score averaged across all the examples.
# hierScore.dir: relative path to folder where the matrix with the scores of the classes corrected in according to TPR-threshold-free algorithm is stored 
# macro.dir: relative path to folder where the macro-centric measures (i.e. AUC and PxR across classes) are stored
# Fmeas.dir: relative path to folder where example-centric measures (i.e. Precision, Recall, Specificity, F-measure, Accuracy across example) are stored
# OUTPUT:
# 5 rda files stored in the rispective output directory:
# - Matrix with examples on rows and classes on colums representing the hierarchical scores of the classes computed with TPR-threshold-free algorithm.
#	Stored in hierScore.dir folder
# - Example-centric measures computed through find.best.f from F-hier.R file. Stored in macro.dir folder
# - AUC (average and per classes) computed through AUC.single.over.classes from package PerfMeas. Stored in macro.dir
# - Precision at fixed recall levels (average and per classes) computed through precision.at.multiple.recall.level.over.classes from package PerfMeas.
# 	Stored in macro.dir folder 
# - PRC (average and per.class) computed by precrec package. Stored in macro.dir
Do.tpr.threshold.free <- function(	norm=TRUE, norm.type= "NONE", flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
									flat.dir=flat.dir, ann.dir=ann.dir, dag.dir=dag.dir, flat.norm.dir=NULL,
									n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
								){
	## Loading Data ############
	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	root <- root.node(hpo);
	ann.no.root <- target <- hpo.ann[,-which(colnames(hpo.ann)==root)];

	## Computing Flat Perfromances
	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold-free correction 
	S.tpr <- pred <- tpr.threshold.free(S.flat, hpo, root=root);

	## Computing Hierarchical Performances
	## Hierarchical AUC (average and per.class) computed by PerfMeas package
	AUC.tpr <- AUC.single.over.classes(ann.no.root, S.tpr, hpo, root=root)[c("average","per.class")];

	## Hierarchical PxR at fixed recall levels (average and per.class) computed by PerfMeas package
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.tpr);

	## Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target, pred, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type,".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type,".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
	}
}

# Hold-Out version
# High level function to correct the scores with a hierarchy according to TPR-threshold-free algorithm performing a classical holdout procedure
# All input parameters are the same of the corresponding function above, except the following:
# ind.test.set: vector of integer. Indices refer to the examples of the adjancency matrix to be used in the test set. 
# ind.dir: relative path to folder where ind.test.set is stored
# An output file is added: PRC (average and per.class) computed by precrec package. Stored in macro.dir
Do.tpr.threshold.free.holdout <- function(	norm=TRUE, norm.type= "NONE", flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
											flat.dir=flat.dir, ann.dir=ann.dir, dag.dir=dag.dir, flat.norm.dir=NULL,
											ind.test.set=ind.test.set, ind.dir=ind.dir, n.round=3, f.criterion ="F",											 
											hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
										){
	## Loading Data ############
	# loading examples indices of the test set
	ind.set <- paste0(ind.dir, ind.test.set, ".rda");
	ind.test <- get(load(ind.set));

	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 																		 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## shrinking the size of S.flat to the examples of test set
	S.flat <- S.flat[ind.test,];

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	root <- root.node(hpo);
	ann.no.root <- target <- hpo.ann[ind.test,-which(colnames(hpo.ann)==root)];

	## Computing Flat Perfromances
	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold-free correction 
	S.tpr <- pred <- tpr.threshold.free(S.flat, hpo, root=root);

	## Computing Hierarchical Performances
	## Hierarchical AUC (average and per.class) computed by PerfMeas package
	AUC.tpr <- AUC.single.over.classes(ann.no.root, S.tpr, hpo, root=root)[c("average","per.class")];

	## Hierarchical PxR at fixed recall levels (average and per.class) computed by PerfMeas package
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.tpr);

	## Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target, pred, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprTF.rda"), compress=TRUE);
	}
}


##**************************************************##
## TPR version with DESCENDANTS instead of children ##
##**************************************************##

# High level function to compute hierarchical correction according to TPR-threshold-descendants algorithm. 
# It perform a k fold cross-validation to find the best threshold maximizing on F.max measure.
# INPUT: 
# threshold: range of threshold values to be tested in order to find the best threshold. The denser the range is, the higher the probability to find the best
# 			 theshold, but the execution time will be increasing (def: from:0.1, to:0.9, by:0.1 step).
# kk:  number of folds of the cross validation (def: 5).
# seed:  intialization seed for the random generator to create folds (def:0). If NULL (default) no initialization is performed.
# flat.norm.file: name of the flat scores matrix already normalized in according to a normalization method (without rda extension). 
# ann.file: name of the target labels (without rda extension). It must be  an .rda file containing the label matrix of the examples (def: ann.file)
# dag.file: name of the graph that represents the hierarchy of the classes. 
# n.round: 	number of rounding digits to be applied to the hierarchical scores matrix (def. 3). 
#			It's used for choosing the best threshold on the basis of the best F.measure (see f.criterion parameter).
# f.criterion: character. Type of F-measure to be used to select the best F.measure. There are 2 possibilities: 
#              1. "F" (default) corresponds to the harmonic mean between the average precision and recall; 
#              2. "avF" corresponds to the per-example F-score averaged across all the examples.
# flat.dir: relative path to folder where flat normalized scores matrix is stored
# ann.dir: relative path to folder where annotation matrix is stored
# dag.dir: relative path to folder where graph is stored
# hierScore.dir: relative path to folder where the matrix with the scores of the classes corrected in according to TPR-threshold-descendants algorithm is stored 
# macro.dir: relative path to folder where the macro-centric measures (i.e. AUC and PxR across classes) are stored
# Fmeas.dir: relative path to folder where example-centric measures (i.e. Precision, Recall, Specificity, F-measure, Accuracy across example) are stored
# OUTPUT:
# 4 rda files stored in the rispective output directory:
# - Matrix with examples on rows and classes on colums representing the hierarchical scores of the classes computed with TPR-threshold-descendants algorithm.
#	Stored in hierScore.dir folder
# - Example-centric measures computed through find.best.f from F-hier.R file. Stored in macro.dir folder
# - AUC (average and per classes) computed through AUC.single.over.classes from package PerfMeas. Stored in macro.dir
# - Precision at fixed recall levels (average and per classes) computed through precision.at.multiple.recall.level.over.classes from package PerfMeas.
# 	Stored in macro.dir folder 
Do.tpr.threshold.descendants.cv <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), kk=5, seed=NULL,norm=TRUE, norm.type= "NONE",
												flat.file=flat.file, ann.file=ann.file, dag.file=dag.file,
												flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir,	flat.norm.dir=NULL,											
												n.round=3, f.criterion ="F",	 
												hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", Fmeas.dir="Fmeas.dir/"
											){
	## Loading Data ############
	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 	
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table 
	root <- root.node(hpo);
	ann.no.root <- hpo.ann[,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## Hierarchical Correction #########################
	## TPR-threshold-descendants correction using k-fold cross-validation to find best threshold 

	## splitting data in 5 fold. Each fold is a sub matrix. Rows are a random sub-sampling of a set of genes; columns are all HPO terms
	folds <- do.unstratified.cv.data(S.flat, kk=kk, seed=seed);

	## storing the best F.max and the best threshold for each of k training set
	training.top.Fmax <- vector(mode="list", length=kk);
	names(training.top.Fmax) <- paste0(rep("fold",kk), 1:kk);

	## storing hierarchical scores for each of k test set
	test.hier.score <- vector(mode="list", length=kk);
	names(test.hier.score) <- paste0(rep("fold",kk), 1:kk);

	## storing all the best protein-centric measures for each of k test set 
	best.avg.meas.test <- vector(mode="list", length=kk);
	names(best.avg.meas.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average macro AUC for each of k test set 
	AUC.average.test <- vector(mode="list", length=kk);
	names(AUC.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class macro AUC for each of k test set 
	AUC.class.test <- vector(mode="list", length=kk);
	names(AUC.class.test) <- paste0(rep("fold",kk), 1:kk);

	## storing average PxR for each of k test set 
	PXR.average.test <- vector(mode="list", length=kk);
	names(PXR.average.test) <- paste0(rep("fold",kk), 1:kk);

	## storing per class PxR for each of k test set 
	PXR.class.test <- vector(mode="list", length=kk);
	names(PXR.class.test) <- paste0(rep("fold",kk), 1:kk);

	## variable for hosting the k sub-matrix assembled 
	S.tpr <- c();

	## Let's start k-fold crossing validation for choosing best threshold maximizing on F.max...
	for(k in 1:kk){
		test <- S.flat[folds[[k]],];		# test set: 1 out of k folds (e.g.: if k=5, testing set is 1/5 of data)
		training <- S.flat[!rownames(S.flat) %in% rownames(test),];		# training set: (k-1) out of k folds (e.g.: if k=5, training set is 4/5 of data)
		gc();

		top.Fmax <- 0;
		best.Fmaxt <- 0;
		for(t in threshold){
			training.hier.score <- tpr.threshold.descendants(training, hpo, root=root, t=t);	## training set hierarchical correction...
			pred.training <- training.hier.score;
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxt <- t;
				training.top.Fmax[[k]] <- c(best.Fmax=top.Fmax, best.thres=best.Fmaxt);
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, sep="\t", "\n");
			}
		}
		test.hier.score[[k]] <- tpr.threshold.descendants(test, hpo, root=root, t=training.top.Fmax[[k]][2]);	## testing set hierarchical correction...
		pred.test <- test.hier.score[[k]];
		target.test <- hpo.ann[rownames(pred.test),colnames(pred.test)];
		test.avg.meas <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion =f.criterion, verbose=FALSE);;	## testing set evaluating...
		best.avg.meas.test[[k]] <- test.avg.meas;

		## Hierarchical AUC (average and per.class) computed with PerfMeas package
		AUC.average.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$average;	## storing average AUC of each k testing set
		AUC.class.test[[k]] <- AUC.single.over.classes(target.test,pred.test,hpo,root)$per.class;	## storing AUC per.class of each k testing set

		## Hierarchical PxR at fixed recall levels (average and per.class) computed with PerfMeas package
		PXR.average.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$avgPXR;	## storing average PxR of each k testing set
		PXR.class.test[[k]] <- precision.at.multiple.recall.level.over.classes(target.test,pred.test)$PXR;	## storing PxR per.class of each k testing set

		## assembling of all the hierarchical scores of each k sub-matrix to build the full matrix 
		S.tpr <- rbind(S.tpr,test.hier.score[[k]]);
	}
	## put the rows (i.e. genes) of assembled k sub-matrix in the same order of the full beginning matrix
	S.tpr <- S.tpr[rownames(S.flat),];

	## Averaging Performances Measures across k testing set ###################
	## averaging protein-centric measures across k testing sets 
	F.meas <- Reduce("+", best.avg.meas.test)/kk;

	## averaging macro-AUC (average and per.class) across k testing sets 
	AUC.average.tpr.over.test <- Reduce("+", AUC.average.test)/kk;
	AUC.class.tpr.over.test <- Reduce("+", AUC.class.test)/kk;
	AUC.tpr <- list(average=AUC.average.tpr.over.test, per.class=AUC.class.tpr.over.test);

	## averaging PxR (average and per.class) across k testing sets 
	PXR.average.tpr.over.test <- Reduce("+", PXR.average.test)/kk;
	PXR.class.tpr.over.test <- Reduce("+", PXR.class.test)/kk;
	PXR.tpr <- list(average=PXR.average.tpr.over.test, per.class=PXR.class.tpr.over.test);

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.",  flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type,".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type,".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(F.meas, file=paste0(Fmeas.dir, "PCM.", norm.type,".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
	}
}

# High level function to correct the scores with a hierarchy according to TPR-threshold-descendants algorithm performing a classical holdout procedure
# All input parameters are the same of the above corresponding function that implement a kfcv, except the following:
# ind.test.set: vector of integer. Indices refer to the examples of the adjancency matrix to be used in the test set. 
# ind.dir: relative path to folder where ind.test.set is stored
# An output file is added: PRC (average and per.class) computed by precrec package. It assumed to store in macro.dir
Do.tpr.threshold.descendants.holdout <- function(	threshold=seq(from=0.1, to=0.9, by=0.1), kk=5, seed=NULL, norm=TRUE, norm.type= "NONE",
													flat.file=flat.file, ann.file=ann.file, dag.file=dag.file, ind.test.set=ind.test.set,
													ind.dir=ind.dir,flat.dir=flat.dir, ann.dir=ann.dir,dag.dir=dag.dir, flat.norm.dir=NULL, 
													n.round=3, f.criterion ="F", hierScore.dir="hierScore.dir/", macro.dir="macro.dir/", 
													Fmeas.dir="Fmeas.dir/"
												){
	## Loading Data ############
	## loading examples indices of the test set
	ind.set <- paste0(ind.dir, ind.test.set, ".rda");
	ind.test <- get(load(ind.set));

	## loading hpo dag
	dag.path <- paste0(dag.dir, dag.file,".rda");
	hpo <- get(load(dag.path));
	
	##root node
	root <- root.node(hpo);

	## loading flat scores matrix relative to a specific subontology
	flat.path <- paste0(flat.dir, flat.file,".rda");
	if(norm==TRUE){
		S.flat <- get(load(flat.path));
		gc();	##in order to save ram memory..

		## removing root node from flat norm matrix if it exists
		if(root %in% colnames(S.flat)){
			S.flat <- S.flat[,-which(colnames(S.flat)==root)];
		}
	}else{
		Do.FLAT.scores.normalization(	norm.type= norm.type, flat.file=flat.file, ann.file=ann.file, 
						 				dag.file=dag.file, flat.dir=flat.dir, ann.dir=ann.dir, 
						 				dag.dir=dag.dir, flat.norm.dir=flat.norm.dir
									)
		flat.path <- paste0(flat.norm.dir, norm.type,".",flat.file,".rda");
		S.flat <- get(load(flat.path));
	}

	## scores flat matrix shrinked to test and training test respectively
	S.test <- S.flat[ind.test,];
	S.training <- S.flat[-ind.test,];
	S.flat <- S.test;

	## loading annotation matrix
	ann.path <- paste0(ann.dir, ann.file,".rda");
	hpo.ann <- get(load(ann.path));
	gc();

	## removing root node from annotation table and shrinking annotation table to tet set
	ann.no.root <- hpo.ann[ind.test,-which(colnames(hpo.ann)==root)];

	## FLAT AUC computed by PerfMeas package
	AUC.flat <- AUC.single.over.classes(ann.no.root, S.flat, hpo, root=root)[c("average","per.class")];

	# FLAT PxRs computed by PerfMeas pacakge
	PXR.flat <- precision.at.multiple.recall.level.over.classes(ann.no.root, S.flat);

	## F.measure: Computing Flat Examples-Measures 
	FMM.flat <- find.best.f(ann.no.root, S.flat, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);

	## FLAT PRC computed by precrec package (more precise and accurate than PerfMeas)
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.flat <- join_scores(S.flat);
	res.flat <- evalmod(scores=scores.flat, labels=labels, dsids=1:ncol(S.flat), modnames=colnames(S.flat));
	meas.flat <- auc(res.flat);
	prc.flat <- subset(meas.flat, curvetypes == "PRC");
	prc.mean.flat <- mean(prc.flat$aucs);
	PRC.class <- prc.flat$aucs;
	names(PRC.class) <- prc.flat$modnames;
	PRC.flat <- list(average=prc.mean.flat, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas 

	## Hierarchical Correction #########################
	## TPR-Threshold correction using k-fold cross-validation to find best threshold on training set
	folds <- do.unstratified.cv.data(S.training, kk=kk, seed=seed);
			
	## Let's start k-fold crossing validation for choosing best threshold maximizing on the basis of F.max value...
	for(k in 1:kk){
		training <- S.training[folds[[k]],];	
		
		top.Fmax <- 0;
		best.Fmaxt <- 0;
		for(t in threshold){
			pred.training <- tpr.threshold.descendants(training, hpo, root=root, t=t);	## training set hierarchical correction...
			target.training <- hpo.ann[rownames(pred.training),colnames(pred.training)];
			training.avg.meas <- find.best.f(target.training, pred.training, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	## training evaluating... 
			training.Fmax <- training.avg.meas[4];	## F.max maximization...
			if(training.Fmax > top.Fmax){
				top.Fmax <- training.Fmax;
				best.Fmaxt <- t;
				cat("training fold:",k, "better Fmax.avg found:",top.Fmax, "best threshold:",best.Fmaxt, sep="\t", "\n");
			}
		}
	}
	S.tpr <- pred.test <- tpr.threshold.descendants(S.test, hpo, root=root, t=best.Fmaxt);	## testing set hierarchical correction..
	target.test <- ann.no.root;
	
	## AUC (average and per.class) computed with PerfMeas package on the test set
	AUC.tpr <- AUC.single.over.classes(target.test, pred.test, hpo, root)[c("average","per.class")];	

	## PxR at fixed recall levels (average and per.class) computed with PerfMeas package on the test set
	PXR.tpr <- precision.at.multiple.recall.level.over.classes(target.test, pred.test);

	## F.measure: Computing Hierarchical Examples-Measures 
	FMM.tpr <- find.best.f(target.test, pred.test, n.round=n.round, f.criterion=f.criterion, verbose=FALSE);	

	## Hierarchical PRC (average and per.class) computed by precrec package
	labels <- join_labels(ann.no.root);	##N.B.: if there are some terms with NO annotations, the function return a *stop* mex
	scores.tpr <- join_scores(S.tpr);
	res.tpr <- evalmod(scores=scores.tpr, labels=labels, dsids=1:ncol(S.tpr), modnames=colnames(S.tpr));
	meas.tpr <- auc(res.tpr);
	prc.tpr <- subset(meas.tpr, curvetypes == "PRC");
	prc.mean.tpr <- mean(prc.tpr$aucs);
	PRC.class <- prc.tpr$aucs;
	names(PRC.class) <- prc.tpr$modnames;
	PRC.tpr <- list(average=prc.mean.tpr, per.class=PRC.class); #saving PRC result in the same format of package PerfMeas

	## Storing Results #########
	if(norm==TRUE){
		save(S.tpr, file=paste0(hierScore.dir, flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
	}else{
		save(S.tpr, file=paste0(hierScore.dir, norm.type,".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(AUC.flat, AUC.tpr, file=paste0(macro.dir, "AUC.", norm.type, ".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PXR.flat, PXR.tpr, file=paste0(macro.dir, "PXR.", norm.type, ".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(FMM.flat, FMM.tpr, file=paste0(Fmeas.dir, "PCM.", norm.type, ".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
		save(PRC.flat, PRC.tpr, file=paste0(macro.dir, "PRC.", norm.type, ".", flat.file, ".hierScores.tprTdesc.rda"), compress=TRUE);
	}
}