R/perf.meas.R

In HEMDAG: Hierarchical Ensemble Methods for Directed Acyclic Graphs

##############################################
## Functions to evaluate HEMDAG performance ##
##############################################

#' @name auprc
#' @aliases auprc.single.class
#' @aliases auprc.single.over.classes
#' @title AUPRC measures
#' @description Compute the Area under the Precision Recall Curve (AUPRC) through \pkg{precrec} package.
#' @details The AUPRC (for a single class or for a set of classes) is computed either one-shot or averaged across stratified folds.
#' @details \code{auprc.single.class} computes the AUPRC just for a given class.
#' @details \code{auprc.single.over.classes} computes the AUPRC for a set of classes, returning also the averaged values across the classes.
#' @details For all those classes having zero annotations, the AUPRC is set to 0. These classes are discarded in the computing of the AUPRC
#' averaged across classes, both when the AUPRC is computed one-shot or averaged across stratified folds.
#' @details Names of rows and columns of \code{labels} and \code{predicted} matrix must be provided in the same order, otherwise a stop message is returned.
#' @param folds number of folds on which computing the AUPRC. If \code{folds=NULL} (\code{def.}), the AUPRC is computed one-shot,
#' otherwise the AUPRC is computed averaged across folds.
#' @param seed initialization seed for the random generator to create folds. Set \code{seed} only if \code{folds}\eqn{\neq}\code{NULL}.
#' If \code{seed=NULL} and \code{folds}\eqn{\neq}\code{NULL}, the AUPRC averaged across folds is computed without seed initialization.
#' @param labels vector of the true labels (0 negative, 1 positive examples).
#' @param target matrix with the target multilabel: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param scores a numeric vector of the values of the predicted labels (scores).
#' @param predicted a numeric matrix with predicted values (scores): rows correspond to examples and columns to classes.
#' @return \code{auprc.single.class} returns a numeric value corresponding to the AUPRC for the considered class;
#' \code{auprc.single.over.classes} returns a list with two elements:
#' \enumerate{
#'  \item average: the average AUPRC across classes;
#'  \item per.class: a named vector with AUPRC for each class. Names correspond to classes.
#' }
#' @export
#' @examples
#' data(labels);
#' data(scores);
#' data(graph);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' prc.single.class <- auprc.single.class(L[,3], S[,3], folds=5, seed=23);
#' prc.over.classes <- auprc.single.over.classes(L, S, folds=5, seed=23);
auprc.single.class <- function(labels, scores, folds=NULL, seed=NULL){
    if(is.matrix(labels) || is.matrix(scores))
        stop("labels or scores must be a vector");
    if(length(scores)!=length(labels))
        stop("length of true and predicted labels does not match");
    if(any(names(labels)!=names(scores)))
        stop("names of labels and scores are not in the same order");
    if(any((labels!=0) & (labels!=1)))
        stop("labels variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    if(!is.null(folds) && is.null(seed))
        warning("folds are generated without seed initialization");
    ## degenerate case when all labels are equals
    if((all(labels==0)) || (all(labels==1))){
        if(!is.null(folds)){
            prc.avg <- 0;
            prc.fold <- rep(0,folds);
            prc <- list(average=prc.avg, across.fold=prc.fold);
        }else{
            prc <- 0;
            names(prc) <- "one.shoot";
        }
        return(prc);
    }
    ## compute PRC averaged across folds
    if(!is.null(folds)){
        indices <- 1:length(labels);
        positives <- which(labels==1);
        foldIndex <- stratified.cv.data.single.class(indices, positives, kk=folds, seed=seed);
        testIndex <- mapply(c, foldIndex$fold.positives, foldIndex$fold.negatives, SIMPLIFY=FALSE); ## index of examples used for test set..
        fold.check <- unlist(lapply(testIndex,length));
        if(any(fold.check==0))
            stop("number of folds selected too high: some folds have no examples. Please reduce the number of folds");
        prc.fold <- numeric(folds);
        for(k in 1:folds){
            labels.test <- labels[testIndex[[k]]];
            scores.test <- scores[testIndex[[k]]];
            if(sum(labels.test) > 0){
                if(all(labels.test==1)){ ## degenerate case when all labels in the k fold are equals to 1
                    prc.fold[k] <- 0;
                }else{
                    res <- evalmod(scores=scores.test, labels=labels.test);
                    aucs <- auc(res);
                    prc <- subset(aucs, curvetypes=="PRC");
                    prc.fold[k] <- prc$aucs;
                }
            }else{
                prc.fold[k] <- 0;
            }
        }
        if(all(prc.fold==0)){
            prc.avg <- 0; ## degenerate case in which for each fold all the PRC computed by precrec are 0
        }else{
            classIndex <- which(prc.fold!=0);
            classPerfs <- prc.fold[classIndex];
            prc.avg <- mean(classPerfs);
        }
        prc <- list(average=prc.avg, across.fold=prc.fold);
        return(prc);
    }
    ## compute PRC one-shoot
    res <- evalmod(scores=scores, labels=labels);
    aucs <- auc(res);
    prc <- subset(aucs, curvetypes=="PRC")$aucs;
    names(prc) <- "one.shoot";
    return(prc);
}

#' @rdname auprc
#' @export
auprc.single.over.classes <- function(target, predicted, folds=NULL, seed=NULL){
    if(!is.matrix(target) || !is.matrix(predicted))
        stop("target or predicted must be a matrix");
    n.examples <- nrow(target);
    n.classes <- ncol(target);
    if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
        stop ("number of rows or columns do not match between target and predicted classes");
    if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
        stop("rows or columns names of target and predicted are not in the same order");
    if(any((target!=0) & (target!=1)))
        stop("target variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    ## AUPRC averaged across folds and classes
    if(!is.null(folds)){
        prc.class <- rep(0,ncol(predicted));
        names(prc.class) <- colnames(predicted);
        for(i in 1:ncol(predicted)){
            prc.class[i] <- auprc.single.class(target[,i], predicted[,i], folds=folds, seed=seed)$average;
        }
        if(all(prc.class==0)){
            prc.avg <- 0;
        }else{
            classIndex <- which(prc.class!=0);
            classPerfs <- prc.class[classIndex];
            prc.avg <- mean(classPerfs);
        }
        prc.res <- list(average=prc.avg, per.class=prc.class);
        return(prc.res);
    }
    ## AUPRC one-shot over classes
    ## if there are classes with zero annotations, we remove them..
    target.names <- colnames(target);
    sample.names <- rownames(target);
    class.ann <- apply(target,2,sum);
    class.noann <- which(class.ann==0);
    check <- length(class.noann)!=0;
    check.degen <- length(class.noann)!=ncol(target); ## degenerate case when all the classes have zero annotations
    if(check & check.degen){
        target <- target[,-class.noann];
        predicted <- predicted[,-class.noann];
    }
    ## degenerate case when target and predicted are vectors: just one class has an annotation. Might happen in cross validation..
    if(!is.matrix(target)){
        selected <- which(class.ann!=0);
        target <- matrix(target, nrow=n.examples, dimnames=list(sample.names, names(selected)));
    }
    if(!is.matrix(predicted)){
        selected <- which(class.ann!=0);
        predicted <- matrix(predicted, nrow=n.examples, dimnames=list(sample.names, names(selected)));
    }
    ## compute PRC considering only those class with non-zero annotations
    prc.class <- rep(0,ncol(predicted));
    names(prc.class) <- colnames(predicted);
    for(i in 1:ncol(predicted)){
        prc.class[i] <- auprc.single.class(target[,i], predicted[,i], folds=NULL, seed=NULL);
    }
    ## if there are classes with zero annotations, set the prc of those classes to zero and restore the start classes order
    if(check & check.degen){
        prc.class <- prc.class[target.names];
        prc.class[is.na(prc.class)] <- 0;
        names(prc.class) <- target.names;
    }
    if(all(prc.class==0)){
            prc.avg <- 0;
    }else{
        classIndex <- which(prc.class!=0);
        classPerfs <- prc.class[classIndex];
        prc.avg <- mean(classPerfs);
    }
    prc.res <- list(average=prc.avg, per.class=prc.class);
    return(prc.res);
}

#' @name auroc
#' @aliases auroc.single.class
#' @aliases auroc.single.over.classes
#' @title AUROC measures
#' @description Compute the Area under the ROC Curve (AUROC) through \pkg{precrec} package.
#' @details The AUROC (for a single class or for a set of classes) is computed either one-shot or averaged across stratified folds.
#' @details \code{auroc.single.class} computes the AUROC just for a given class.
#' @details \code{auroc.single.over.classes} computes the AUROC for a set of classes, including their average values across all the classes.
#' @details For all those classes having zero annotations, the AUROC is set to 0.5. These classes are included in the computing of the AUROC
#' averaged across classes, both when the AUROC is computed one-shot or averaged across stratified folds.
#' @details The AUROC is set to 0.5 to all those classes having zero annotations.
#' Names of rows and columns of \code{labels} and \code{predicted} must be provided in the same order, otherwise a stop message is returned.
#' @param folds number of folds on which computing the AUROC. If \code{folds=NULL} (\code{def.}), the AUROC is computed one-shot,
#' otherwise the AUROC is computed averaged across folds.
#' @param seed initialization seed for the random generator to create folds. Set \code{seed} only if \code{folds}\eqn{\neq}\code{NULL}.
#' If \code{seed=NULL} and \code{folds}\eqn{\neq}\code{NULL}, the AUROC averaged across folds is computed without seed initialization.
#' @param labels vector of the true labels (0 negative, 1 positive examples).
#' @param target annotation matrix: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param scores a numeric vector of the values of the predicted labels (scores).
#' @param predicted a numeric matrix with predicted values (scores): rows correspond to examples and columns to classes.
#' @return \code{auroc.single.class} returns a numeric value corresponding to the AUROC for the considered class;
#' \code{auprc.single.over.classes} returns a list with two elements:
#' \enumerate{
#'  \item average: the average AUROC across classes;
#'  \item per.class: a named vector with AUROC for each class. Names correspond to classes.
#' }
#' @export
#' @examples
#' data(labels);
#' data(scores);
#' data(graph);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' auc.single.class <- auroc.single.class(L[,3], S[,3], folds=5, seed=23);
#' auc.over.classes <- auroc.single.over.classes(L, S, folds=5, seed=23);
auroc.single.class <- function(labels, scores, folds=NULL, seed=NULL){
    if(is.matrix(labels) || is.matrix(scores))
        stop("labels or scores must be a vector");
    if(length(scores)!=length(labels))
        stop("length of true and predicted labels does not match");
    if(any(names(labels)!=names(scores)))
        stop("names of labels and scores are not in the same order");
    if(any((labels!=0) & (labels!=1)))
        stop("labels variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    if(!is.null(folds) && is.null(seed))
        warning("folds are generated without seed initialization");
    ## degenerate case when all labels are equals
    if((all(labels==0)) || (all(labels==1))){
        if(!is.null(folds)){
            auc.avg <- 0.5;
            auc.fold <- rep(0.5,folds);
            auc <- list(average=auc.avg, across.fold=auc.fold);
        }else{
            auc <- 0.5;
            names(auc) <- "one.shoot";
        }
        return(auc);
    }
    ## compute AUC averaged across folds
    if(!is.null(folds)){
        indices <- 1:length(labels);
        positives <- which(labels==1);
        foldIndex <- stratified.cv.data.single.class(indices, positives, kk=folds, seed=seed);
        testIndex <- mapply(c, foldIndex$fold.positives, foldIndex$fold.negatives, SIMPLIFY=FALSE); ## index of examples used for test set..
        fold.check <- unlist(lapply(testIndex,length));
        if(any(fold.check==0))
            stop("number of folds selected too high: some folds have no examples. Please reduce the number of folds");
        auc.fold <- numeric(folds);
        for(k in 1:folds){
            labels.test <- labels[testIndex[[k]]];
            scores.test <- scores[testIndex[[k]]];
            if(sum(labels.test) > 0){
                if(all(labels.test==1)){ ## degenerate case when all labels in the k fold are equals to 1
                    auc.fold[k] <- 0.5;
                }else{
                    res <- evalmod(scores=scores.test, labels=labels.test);
                    aucs <- auc(res);
                    auc <- subset(aucs, curvetypes=="ROC");
                    auc.fold[k] <- auc$aucs;
                }
            }else{
                auc.fold[k] <- 0.5;
            }
        }
        auc.avg <- mean(auc.fold);
        auc <- list(average=auc.avg, across.fold=auc.fold);
        return(auc);
    }
    ## compute AUC one-shoot
    res <- evalmod(scores=scores, labels=labels);
    aucs <- auc(res);
    auc <- subset(aucs, curvetypes=="ROC")$aucs;
    names(auc) <- "one.shoot";
    return(auc);
}

#' @rdname auroc
#' @export
auroc.single.over.classes <- function(target, predicted, folds=NULL, seed=NULL){
    if(!is.matrix(target) || !is.matrix(predicted))
        stop("target or predicted must be a matrix");
    n.examples <- nrow(target);
    n.classes <- ncol(target);
    if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
        stop ("number of rows or columns do not match between target and predicted classes");
    if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
        stop("rows or columns names of target and predicted are not in the same order");
    if(any((target!=0) & (target!=1)))
        stop("target variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    ## AUROC averaged across folds and classes
    if(!is.null(folds)){
        auc.class <- rep(0,ncol(predicted));
        names(auc.class) <- colnames(predicted);
        for(i in 1:ncol(predicted)){
            auc.class[i] <- auroc.single.class(target[,i], predicted[,i], folds=folds, seed=seed)$average;
        }
        auc.avg <- mean(auc.class);
        auc.res <- list(average=auc.avg, per.class=auc.class);
        return(auc.res);
    }
    ## if there are classes with zero annotations, we remove them..
    target.names <- colnames(target);
    sample.names <- rownames(target);
    class.ann <- apply(target,2,sum);
    class.noann <- which(class.ann==0);
    check <- length(class.noann)!=0;
    check.degen <- length(class.noann)!=ncol(target); ## degenerate case when all the classes have zero annotation
    if(check & check.degen){
        target <- target[,-class.noann];
        predicted <- predicted[,-class.noann];
    }
    ## degenerate case when target and predicted are vectors: just one class has an annotation. Might happen in cross validation..
    if(!is.matrix(target)){
        selected <- which(class.ann!=0);
        target <- matrix(target, nrow=n.examples, dimnames=list(sample.names, names(selected)));
    }
    if(!is.matrix(predicted)){
        selected <- which(class.ann!=0);
        predicted <- matrix(predicted, nrow=n.examples, dimnames=list(sample.names, names(selected)));
    }
    ## compute AUC considering only those class with non-zero annotations
    auc.class <- rep(0,ncol(predicted));
    names(auc.class) <- colnames(predicted);
    for(i in 1:ncol(predicted)){
        auc.class[i] <- auroc.single.class(target[,i],predicted[,i], folds=NULL, seed=NULL);
    }
    ## if there are classes with zero annotations, set the AUC of those classes to zero and restore the start classes order
    if(check & check.degen){
        auc.class <- auc.class[target.names];
        auc.class[is.na(auc.class)] <- 0.5;
        names(auc.class) <- target.names;
    }
    auc.avg <- mean(auc.class);
    auc.res <- list(average=auc.avg, per.class=auc.class);
    return(auc.res);
}

#' @name multilabel.F.measure
#' @aliases F.measure.multilabel
#' @title multilabel F-measure
#' @description Method for computing Precision, Recall, Specificity, Accuracy and F-measure for multiclass and multilabel classification.
#' @details Names of rows and columns of \code{target} and \code{predicted} matrix must be provided in the same order, otherwise a stop message is returned.
#' @param target matrix with the target multilabel: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param predicted a numeric matrix with discrete predicted values: rows correspond to examples and columns to classes.
#' \eqn{predicted[i,j]=1} if example \eqn{i} is predicted belonging to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param b.per.example a boolean value.
#' \itemize{
#'  \item \code{TRUE}: results are returned for each example;
#'  \item \code{FALSE}: only the average results are returned;
#' }
#' @return Two different outputs respect to the input parameter \code{b.per.example}:
#' \itemize{
#'  \item \code{b.per.example==FALSE}: a list with a single element average. A named vector with average precision (P), recall (R),
#'   specificity (S), F-measure (F), average F-measure (avF) and Accuracy (A) across examples. F is the F-measure computed as the
#'   harmonic mean between the average precision and recall; av.F is the F-measure computed as average across examples;
#'  \item \code{b.per.example==FALSE}: a list with two elements:
#'   \enumerate{
#'     \item average: a named vector with average precision (P), recall (R), specificity (S), F-measure (F), average F-measure (avF)
#'      and Accuracy (A) across examples;
#'     \item per.example: a named matrix with the Precision (P), Recall (R), Specificity (S), Accuracy (A), F-measure (F) and
#'      av.F-measure (av.F) for each example. Row names correspond to examples, column names correspond respectively to Precision (P), Recall (R),
#'      Specificity (S), Accuracy (A), F-measure (F) and av.F-measure (av.F);
#'   }
#' }
#' @examples
#' data(labels);
#' data(scores);
#' data(graph);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' S[S>0.7] <- 1;
#' S[S<0.7] <- 0;
#' fscore <- F.measure.multilabel(L,S);
#' @export
#' @docType methods
setGeneric("F.measure.multilabel",
    function(target, predicted, b.per.example=FALSE) standardGeneric("F.measure.multilabel"));
#' @rdname multilabel.F.measure
setMethod("F.measure.multilabel", signature(target="matrix", predicted="matrix"),
    function(target, predicted, b.per.example=FALSE){
        n.examples <- nrow(target);
        n.classes <- ncol(target);
        if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
            stop ("number of rows or columns do not match between target and predicted classes");
        if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
            stop("rows or columns names of target and predicted are not in the same order");
        if(any((target!=0) & (target!=1)) || any((predicted!=0) & (predicted!=1)))
            stop("target and predicted variables must take values 0 or 1");
        z <- target + predicted;
        TP <- apply(z, 1, function(x){
            return(sum(x==2));
            });
        TN <- apply(z, 1, function(x){
            return(sum(x==0));
        });
        z <- predicted - target;
        FP <- apply(z, 1, function(x){
            return(sum(x==1));
        });
        FN <- apply(z, 1, function(x){
            return(sum(x== -1));
        });
        rm(z);
        ## lines below were useful in debugging :)
        # n <- sum(TP)+sum(TN)+sum(FN)+sum(FP);
        # if(n != (n.examples*n.classes)){
        #     cat("n = ", n, "\n n.examples = ", n.examples, "\n n.classes = ", n.classes, "\n");
        #     cat (" sum(TP) = ", sum(TP), "\n sum(TN) = ", sum(TN), "\n sum(FN) = ", sum(FN), "\n sum(FP) = ", sum(FP), "\n");
        #     warning("something went wrong in F-measure");
        # }
        P <- TP+FP;
        P[which(P==0)] <- 1;  # to avoid division by 0 in precision
        sum.TP.FN <- TP+FN;
        sum.TN.FP <- TN+FP;
        sum.TP.FN[which(sum.TP.FN==0)] <- 1;  # to avoid division by 0 in recall
        sum.TN.FP[which(sum.TN.FP==0)] <- 1;  # to avoid division by 0 in specificity
        precision <- TP/P;
        recall <- TP/sum.TP.FN;
        specificity <- TN/sum.TN.FP;
        prec.rec <- precision+recall;
        prec.rec[which(prec.rec==0)] <- 1;  # to avoid division by 0 for f.measure
        f.measure <- (2*precision*recall)/prec.rec;
        accuracy <- (TP+TN)/n.classes;
        av.precision <- sum(precision)/n.examples;
        av.recall <- sum(recall)/n.examples;
        av.specificity <- sum(specificity)/n.examples;
        av.prec.rec <- av.precision+av.recall;
        if(av.prec.rec == 0){
            av.prec.rec <- 1;
        }
        overall.av.f.measure <- (2*av.precision*av.recall)/av.prec.rec;
        av.f.measure <- sum(f.measure)/n.examples;
        av.accuracy  <- sum(accuracy)/n.examples;
        average <- c(av.precision, av.recall, av.specificity, overall.av.f.measure, av.f.measure, av.accuracy);
        names(average) <- c("P", "R", "S", "F", "avF", "A");
       if(b.per.example){
            per.example <- cbind(precision, recall, specificity, f.measure, accuracy);
            colnames(per.example) <- c("P", "R", "S", "F","A");
            return (list(average=average, per.example=per.example))
        }else{
            return (list(average=average));
        }
    }
)

#' @title Best hierarchical F-score
#' @description Select the best hierarchical F-score by choosing an appropriate threshold in the scores.
#' @details All the examples having no positive annotations are discarded. The predicted scores matrix (\code{predicted}) is rounded
#' according to parameter \code{n.round} and all the values of \code{predicted} are divided by \code{max(predicted)}.
#' Then all the thresholds corresponding to all the different values included in \code{predicted} are attempted, and the threshold
#' leading to the maximum F-measure is selected.
#' @details Names of rows and columns of \code{target} and \code{predicted} matrix must be provided in the same order, otherwise a stop message is returned.
#' @param target matrix with the target multilabel: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param predicted a numeric matrix with continuous predicted values (scores): rows correspond to examples and columns to classes.
#' @param n.round number of rounding digits to be applied to predicted (\code{default=3}).
#' @param verbose a boolean value. If \code{TRUE} (def.) the number of iterations are printed on stdout.
#' @param b.per.example a boolean value.
#' \itemize{
#'  \item \code{TRUE}: results are returned for each example;
#'  \item \code{FALSE}: only the average results are returned;
#' }
#' @return Two different outputs respect to the input parameter \code{b.per.example}:
#' \itemize{
#'  \item \code{b.per.example==FALSE}: a list with a single element average. A named vector with 7 elements relative to the best result in terms
#'   of the F.measure: Precision (P), Recall (R), Specificity (S), F.measure (F), av.F.measure (av.F), Accuracy (A) and the best selected Threshold (T).
#'   F is the F-measure computed as the harmonic mean between the average precision and recall; av.F is the F-measure computed as the average across
#'   examples and T is the best selected threshold;
#'  \item \code{b.per.example==FALSE}: a list with two elements:
#'   \enumerate{
#'    \item average: a named vector with with 7 elements relative to the best result in terms of the F.measure: Precision (P), Recall (R),
#'     Specificity (S), F.measure (F), av.F.measure (av.F), Accuracy (A) and the best selected Threshold (T);
#'    \item per.example: a named matrix with the Precision (P), Recall (R), Specificity (S), Accuracy (A), F-measure (F), av.F-measure (av.F)
#'     and the best selected Threshold (T) for each example. Row names correspond to examples, column names correspond respectively
#'     to Precision (P), Recall (R), Specificity (S), Accuracy (A), F-measure (F), av.F-measure (av.F) and the best selected Threshold (T);
#'   }
#' }
#' @export
#' @examples
#' data(graph);
#' data(labels);
#' data(scores);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' fscore <- find.best.f(L, S, n.round=3, verbose=TRUE, b.per.example=TRUE);
find.best.f <- function(target, predicted, n.round=3, verbose=TRUE, b.per.example=FALSE){
    if(!is.matrix(target) || !is.matrix(predicted))
        stop("target or predicted must be a matrix");
    n.examples <- nrow(target);
    n.classes <- ncol(target);
    if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
        stop("number of rows or columns do not match between target and predicted classes");
    if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
        stop("rows or columns names of target and predicted are not in the same order");
    if(any((target!=0) & (target!=1)))
        stop("labels variable must take values 0 or 1");
    x <- apply(target,1,sum);
    selected <- which(x>0);
    ##  degenerate case when target is a full-zero matrix (all genes without annotations)
    if(length(selected)==0){
        selected <- which(x==0);
    }
    target <- target[selected,];
    ## degenerate case when target is a vector (just one annotated gene)
    if(!is.matrix(target)){
        target <- t(as.matrix(target));
        rownames(target) <- names(selected);
    }
    predicted <- predicted[selected,];
    ## degenerate case when predicted is a vector (just one annotated gene)
    if(!is.matrix(predicted)){
        predicted <- t(as.matrix(predicted));
        rownames(predicted) <- names(selected);
    }
    predicted <- predicted/max(predicted);
    predicted <- round(predicted,n.round);
    n.examples <- nrow(predicted);
    n.classes <- ncol(predicted);
    thresh <- unique(as.numeric(predicted));
    thresh <- sort(thresh);
    best.res <- best <- best.thresh <- 0;
    i <- 0;
    for(t in thresh){
        predicted.labels <- matrix(numeric(n.examples*n.classes), nrow=n.examples);
        predicted.labels[predicted>=t] <- 1;
        res <- F.measure.multilabel(target, predicted.labels, b.per.example);
        if(res$average["F"] > best){
            best <- res$average["F"];
            best.res <- res;
            best.thresh <- t;
        }
        i <- i+1;
        if(i%%100 == 0  && verbose){
            cat("iteration ", i,  "\n");
        }
    }
    ## degenerate case when target is a full-zero matrix: by.def F-score is zero
    if(!is.list(best.res)){
        best.res <- res;
    }
    if(b.per.example){
        best.res$average <- c(best.res$average, best.thresh);
        names(best.res$average)[7] <- "T";
        return(best.res);
    }else{
        best.res <- c(best.res$average, best.thresh);
        names(best.res)[7] <- "T";
        return(best.res);
    }
}

#' @name fmax
#' @title Compute Fmax
#' @description Compute the best hierarchical Fmax either one-shot or averaged across folds
#' @details Names of rows and columns of \code{target} and \code{predicted} matrix must be provided in the same order, otherwise a stop message is returned.
#' @param target matrix with the target multilabel: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param predicted a numeric matrix with predicted values (scores): rows correspond to examples and columns to classes.
#' @param n.round number of rounding digits to be applied to predicted (\code{default=3}).
#' @param verbose a boolean value. If \code{TRUE} (def.) the number of iterations are printed on stdout.
#' @param b.per.example a boolean value.
#' \itemize{
#'  \item \code{TRUE}: results are returned for each example;
#'  \item \code{FALSE}: only the average results are returned;
#' }
#' @param folds number of folds on which computing the Fmax If \code{folds=NULL} (\code{def.}), the Fmax is computed one-shot,
#' otherwise the Fmax is computed averaged across folds.
#' @param seed initialization seed for the random generator to create folds. Set \code{seed} only if \code{folds}\eqn{\neq}\code{NULL}.
#' If \code{seed=NULL} and \code{folds}\eqn{\neq}\code{NULL}, the Fmax averaged across folds is computed without seed initialization.
#' @return Two different outputs respect to the input parameter \code{b.per.example}:
#' \itemize{
#'  \item \code{b.per.example==FALSE}: a list with a single element average. A named vector with 7 elements relative to the best result in terms
#'   of the F.measure: Precision (P), Recall (R), Specificity (S), F.measure (F), av.F.measure (av.F), Accuracy (A) and the best selected Threshold (T).
#'   F is the F-measure computed as the harmonic mean between the average precision and recall; av.F is the F-measure computed as the average across
#'   examples and T is the best selected threshold;
#'  \item \code{b.per.example==FALSE}: a list with two elements:
#'   \enumerate{
#'    \item average: a named vector with with 7 elements relative to the best result in terms of the F.measure: Precision (P), Recall (R),
#'     Specificity (S), F.measure (F), av.F.measure (av.F), Accuracy (A) and the best selected Threshold (T);
#'    \item per.example: a named matrix with the Precision (P), Recall (R), Specificity (S), Accuracy (A), F-measure (F), av.F-measure (av.F)
#'     and the best selected Threshold (T) for each example. Row names correspond to examples, column names correspond respectively
#'     to Precision (P), Recall (R), Specificity (S), Accuracy (A), F-measure (F), av.F-measure (av.F) and the best selected Threshold (T);
#'   }
#' }
#' @export
#' @examples
#' data(graph);
#' data(labels);
#' data(scores);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' fmax <- compute.fmax(L, S, n.round=3, verbose=TRUE, b.per.example=TRUE, folds=5, seed=23);
compute.fmax <- function(target, predicted, n.round=3, verbose=TRUE, b.per.example=FALSE, folds=NULL, seed=NULL){
    if(!is.matrix(target) || !is.matrix(predicted))
        stop("target or predicted must be a matrix");
    n.examples <- nrow(target);
    n.classes <- ncol(target);
    if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
        stop("number of rows or columns do not match between target and predicted classes");
    if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
        stop("rows or columns names of target and predicted are not in the same order");
    if(any((target!=0) & (target!=1)))
        stop("target variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    if(!is.null(folds) && is.null(seed))
        warning("folds are generated without seed initialization");
    ## Fmax averaged across folds
    if(!is.null(folds)){
        testIndex <- unstratified.cv.data(predicted, kk=folds, seed=seed);
        avg.res.list <- list();
        res.per.example <- c();
        for(k in 1:folds){
            fold.res <- find.best.f(target[testIndex[[k]],], predicted[testIndex[[k]],], n.round=n.round, verbose=verbose, b.per.example=TRUE);
            avg.res.list[[k]] <- fold.res$average;
            res.per.example <- rbind(res.per.example, fold.res$per.example);
        }
        Fcv <- Reduce("+", avg.res.list)/folds;
        Fmeas <- apply(res.per.example,2,mean);
        names(Fmeas)[4] <- "avF";
        ## degenerate case when both precision and recall are zero..
        if((Fmeas[["P"]] && Fmeas[["R"]]) == 0){     ## sum(res.per.example[,"F"])==0
            Fmax <- 0;
        }else{
            Fmax <- 2 * (Fmeas[["P"]] * Fmeas[["R"]])/((Fmeas[["P"]] + Fmeas[["R"]]));
        }
        names(Fmax) <- "F";
        Fmax.avg <- append(Fmeas, Fmax, after=3);
        Fmax.avg <- append(Fmax.avg, Fcv["T"]);
        res <- list(average=Fmax.avg, per.example=res.per.example);
        if(b.per.example){
            return(res);
        }else{
            res <- res$average;
            return(res);
        }
    }
    ## Fmax one-shot
    res <- find.best.f(target=target, predicted=predicted, n.round=n.round, verbose=verbose, b.per.example=b.per.example);
    return(res);
}

#' @name pxr
#' @aliases precision.at.all.recall.levels.single.class
#' @aliases precision.at.given.recall.levels.over.classes
#' @title Precision-Recall curves
#' @description Compute the Precision-Recall (PxR) values through \pkg{precrec} package.
#' @details \code{precision.at.all.recall.levels.single.class} computes the precision at all recall levels just for a single class.
#' @details \code{precision.at.given.recall.levels.over.classes} computes the precision at fixed recall levels over classes.
#' @param labels vector of the true labels (0 negative, 1 positive examples).
#' @param target matrix with the target multilabel: rows correspond to examples and columns to classes.
#' \eqn{target[i,j]=1} if example \eqn{i} belongs to class \eqn{j}, \eqn{target[i,j]=0} otherwise.
#' @param scores a numeric vector of the values of the predicted labels (scores).
#' @param predicted a numeric matrix with predicted values (scores): rows correspond to examples and columns to classes.
#' @param folds number of folds on which computing the PXR. If \code{folds=NULL} (\code{def.}), the PXR is computed one-shot,
#' otherwise the PXR is computed averaged across folds.
#' @param seed initialization seed for the random generator to create folds. Set \code{seed} only if \code{folds}\eqn{\neq}\code{NULL}.
#' If \code{seed=NULL} and \code{folds}\eqn{\neq}\code{NULL}, the PXR averaged across folds is computed without seed initialization.
#' @param recall.levels a vector with the desired recall levels (\code{def:} \code{from:0.1}, \code{to:0.9}, \code{by:0.1}).
#' @return \code{precision.at.all.recall.levels.single.class} returns a two-columns matrix, representing a pair of precision and recall values.
#' The first column is the precision, the second the recall;
#' \code{precision.at.given.recall.levels.over.classes} returns a list with two elements:
#' \enumerate{
#'  \item average: a vector with the average precision at different recall levels across classes;
#'  \item fixed.recall: a matrix with the precision at different recall levels: rows are classes, columns precision at different recall levels;
#' }
#' @export
#' @examples
#' data(labels);
#' data(scores);
#' data(graph);
#' root <- root.node(g);
#' L <- L[,-which(colnames(L)==root)];
#' S <- S[,-which(colnames(S)==root)];
#' labels <- L[,1];
#' scores <- S[,1];
#' rec.levels <- seq(from=0.25, to=1, by=0.25);
#' pxr.single <- precision.at.all.recall.levels.single.class(labels, scores);
#' pxr <- precision.at.given.recall.levels.over.classes(L, S, folds=5, seed=23,
#'          recall.levels=rec.levels);
precision.at.all.recall.levels.single.class <- function(labels, scores){
    if(is.matrix(labels) || is.matrix(scores))
        stop("labels or scores must be a vector");
    if(length(scores)!=length(labels))
        stop("length of true and predicted labels does not match");
    if(any(names(labels)!=names(scores)))
        stop("names of labels and scores are not in the same order");
    if(any((labels!=0) & (labels!=1)))
        stop("labels variable must take values 0 or 1");
    if(all(labels==0) || all(labels==1)){
        recall <- seq(from=0.0, to=1, by=0.25);
        precision <- rep(0, length(recall));
        pxr <- cbind(precision=precision, recall=recall);
        return(pxr);
    }else{
        res <- evalmod(mode="basic", labels=labels, scores=scores);
        df <- data.frame(res, stringsAsFactors=TRUE);
        precision <- subset(df, df$type=="precision")$y;
        recall <- subset(df, df$type=="sensitivity")$y;
        pxr <- cbind(precision=precision, recall=recall);
        return(pxr);
    }
}

#' @rdname pxr
#' @export
precision.at.given.recall.levels.over.classes <- function(target, predicted, folds=NULL, seed=NULL, recall.levels=seq(from=0.1, to=1, by=0.1)){
    if(!is.matrix(target) || !is.matrix(predicted))
        stop("target or predicted must be a matrix");
    n.examples <- nrow(target);
    n.classes <- ncol(target);
    if((n.examples!=nrow(predicted)) || (n.classes!=ncol(predicted)))
        stop ("number of rows or columns do not match between target and predicted classes");
    if(any(colnames(target)!=colnames(predicted)) || any(rownames(target)!=rownames(predicted)))
        stop("rows or columns names of target and predicted are not in the same order");
    if(any((target!=0) & (target!=1)))
        stop("target variable must take values 0 or 1");
    if(is.null(folds) && !is.null(seed))
        seed <- NULL;
    if(!is.null(folds) && is.null(seed))
        warning("folds are generated without seed initialization");
    n.classes <- ncol(predicted);
    classes.names <- colnames(predicted);
    len.level <- length(recall.levels);
    pxr <- c();
    ## PXR cross-validated
    if(!is.null(folds)){
        for(i in 1:n.classes){
            labels <- target[,i];
            scores <- predicted[,i];
            df <- create.stratified.fold.df(labels=labels, scores=scores, folds=folds, seed=seed);
            nfold <- format_nfold(nfold_df=df,  score_cols="scores", lab_col="labels", fold_col="folds");
            prec2rec <- c(); ## storing the higher precisions at fixed recall level in the k_th fold
            pxr.fold <- c(); ## storing the higher precisions at fixed recall level average across the k_th folds
            for(k in 1:folds){
                for(j in 1:len.level){
                    ## degenerate case in which labels in the fold are all positive or all negative
                    if(all(nfold$labels[[k]]==1) || all(nfold$labels[[k]]==2)){
                        prec2rec[j] <- 0;
                    }else{
                        res <- evalmod(mode="basic", scores=nfold$scores[k], labels=nfold$labels[k], modnames="m1", dsids=k);
                        df <- data.frame(res, stringsAsFactors=TRUE);
                        precision <- subset(df, df$type=="precision")$y;
                        recall <- subset(df, df$type=="sensitivity")$y;
                        ## we take the higher precision value at the given recall level. NB: recall is monotone...
                        prec2rec[j] <- precision[which(recall - recall.levels[j]>=0)[1]];
                    }
                }
                pxr.fold <- c(pxr.fold, prec2rec);
            }
            if(all(pxr.fold==0)){
                prec2rec <- rep(0,len.level);
            }else{
                for(j in 1:len.level){
                    tmp <- pxr.fold[seq(j,len.level*folds,len.level)];
                    if(all(tmp==0)){
                        prec2rec[j] <- 0;
                    }else{
                        prec2rec[j] <- mean(tmp[which(tmp!=0)]);
                    }
                }
            }
            pxr <- rbind(pxr, prec2rec);
        }
        dimnames(pxr) <- list(classes.names, recall.levels);
        avgpxr <- apply(pxr, 2, mean);
        names(avgpxr) <- recall.levels;
        res <- list(average=avgpxr, fixed.recall=pxr);
        return(res);
    }
    # PXR one-shot
    for(i in 1:n.classes){
        labels <- target[,i];
        scores <- predicted[,i];
        prec2rec <- c();
        if(all(labels==0) || all(labels==1)){
            prec2rec <- rep(0,len.level);
        }else{
            for(j in 1:len.level){
                res <- evalmod(mode="basic", labels=labels, scores=scores);
                df <- data.frame(res, stringsAsFactors=TRUE);
                precision <- subset(df, df$type=="precision")$y;
                recall <- subset(df, df$type=="sensitivity")$y;
                ## we take the higher precision value at the given recall level. NB: recall is monotone...
                prec2rec[j] <- precision[which(recall - recall.levels[j]>=0)[1]];
            }
        }
        pxr <- rbind(pxr, prec2rec);
    }
    dimnames(pxr) <- list(classes.names, recall.levels);
    avgpxr <- apply(pxr, 2, mean);
    names(avgpxr) <- recall.levels;
    res <- list(average=avgpxr, fixed.recall=pxr);
    return(res);
}

#' @name stratified.cross.validation
#' @title Stratified cross validation
#' @description Generate data for the stratified cross-validation.
#' @details Folds are \emph{stratified}, i.e. contain the same amount of positive and negative examples.
#' @param labels labels matrix. Rows are genes and columns are classes. Let's denote \eqn{M} the labels matrix.
#' If \eqn{M[i,j]=1}, means that the gene \eqn{i} is annotated with the class \eqn{j}, otherwise \eqn{M[i,j]=0}.
#' @param examples indices or names of the examples. Can be either a vector of integers or a vector of names.
#' @param positives vector of integers or vector of names. The indices (or names) refer to the indices (or names) of 'positive' examples.
#' @param kk number of folds (\code{def. kk=5}).
#' @param seed seed of the random generator (\code{def. seed=NULL}). If is set to \code{NULL} no initialization is performed.
#' @return \code{stratified.cv.data.single.class} returns a list with 2 two component:
#' \itemize{
#'  \item fold.non.positives: a list with \eqn{k} components. Each component is a vector with the indices (or names) of the non-positive elements.
#'     Indexes (or names) refer to row numbers (or names) of a data matrix;
#'     \item fold.positives: a list with \eqn{k} components. Each component is a vector with the indices (or names) of the positive elements.
#'     Indexes (or names) refer to row numbers (or names) of a data matrix;
#' }
#' @examples
#' data(labels);
#' examples.index <- 1:nrow(L);
#' examples.name <- rownames(L);
#' positives <- which(L[,3]==1);
#' x <- stratified.cv.data.single.class(examples.index, positives, kk=5, seed=23);
#' y <- stratified.cv.data.single.class(examples.name, positives, kk=5, seed=23);
#' z <- stratified.cv.data.over.classes(L, examples.index, kk=5, seed=23);
#' k <- stratified.cv.data.over.classes(L, examples.name, kk=5, seed=23);
#' @export
stratified.cv.data.single.class <- function(examples, positives, kk=5, seed=NULL){
    set.seed(seed);
    if(is.numeric(examples) && length(names(positives))!=0)
        positives <- unname(positives);
    if(is.character(examples) && length(names(positives))!=0)
        positives <- names(positives);
    ## degenerate case when labels have only one positive
    if(length(positives)==1){
        positives <- positives;
    }else{
        positives <- sample(positives);
    }
    ## degenerate case when labels have only one negative
    negatives <- setdiff(examples,positives);
    if(length(negatives)==1){
        negatives <- negatives;
    }else{
        negatives <- sample(negatives);
    }
    n <- length(positives);
    m <- length(negatives);
    set.pos <- list();
    set.neg <- list();
    for (k in 1:kk) {
        ## folds of indices of positive examples
        last.pos <- (k * n) %/% kk;
        first.pos  <- ((k - 1) * n) %/% kk;
        size.pos <-  last.pos - first.pos;
        if(size.pos>1){subset.pos <- positives[1:size.pos];}
        if(size.pos==1){subset.pos <- positives[size.pos];}
        if(size.pos==0){subset.pos <- integer(0);}
        set.pos[[k]] <- subset.pos;
        positives <- setdiff(positives, subset.pos);
        ## folds of indices of negatives examples
        last.neg <- (k * m) %/% kk;
        first.neg  <- ((k - 1) * m) %/% kk;
        size.neg <-  last.neg - first.neg;
        if(size.neg>1){subset.neg <- negatives[1:size.neg];}
        if(size.neg==1){subset.neg <- negatives[size.neg];}
        if(size.neg==0){subset.neg <- integer(0);}
        set.neg[[k]] <- subset.neg;
        negatives <- setdiff(negatives, subset.neg);
    }
    return(list(fold.positives=set.pos, fold.negatives=set.neg));
}

#' @rdname stratified.cross.validation
#' @return \code{stratified.cv.data.over.classes} returns a list with \eqn{n} components, where \eqn{n} is the number of classes of the labels matrix.
#' Each component \eqn{n} is in turn a list with \eqn{k} elements, where \eqn{k} is the number of folds.
#' Each fold contains an equal amount of positives and negatives examples.
#' @export
stratified.cv.data.over.classes <- function(labels, examples, kk=5, seed=NULL){
    set.seed(seed);
    folds <- list();
    for(class in colnames(labels)){
        folds[[class]] <- list();
        positives <- which(labels[,class]==1);
        strfold <- stratified.cv.data.single.class(examples,positives, kk=kk, seed=seed);
        for(k in 1:kk){
            folds[[class]][[k]] <- list();
            names(folds[[class]])[k] <- paste0("fold",k);
            folds[[class]][[k]] <- append(strfold$fold.positives[[k]], strfold$fold.negatives[[k]]);
        }
    }
    return(folds);
}

#' @title DataFrame for stratified cross validation
#' @description Create a data frame for stratified cross-validation.
#' @details Folds are \emph{stratified}, i.e. contain the same amount of positive and negative examples.
#' @param labels vector of the true labels (0 negative, 1 positive).
#' @param scores a numeric vector of the values of the predicted labels.
#' @param seed initialization seed for the random generator to create folds (\code{def. seed=23}).
#' If \code{seed=NULL}, the stratified folds are generated without seed initialization.
#' @param folds number of folds of the cross validation (\code{def. folds=5}).
#' @return A data frame with three columns:
#' \itemize{
#'  \item \code{scores}: contains the predicted scores;
#'    \item \code{labels}: contains the labels as \code{pos} or \code{neg};
#'  \item \code{folds}: contains the index of the fold in which the example falls.
#'    The index can range from 1 to the number of folds.
#' }
#' @export
#' @examples
#' data(labels);
#' data(scores);
#' df <- create.stratified.fold.df(L[,3], S[,3], folds=5, seed=23);
create.stratified.fold.df <- function(labels, scores, folds=5, seed=23){
    if(is.matrix(labels) || is.matrix(scores))
        stop("labels or scores must be a vector");
    if(length(scores)!=length(labels))
        stop("length of true and predicted labels does not match");
    if(any((labels!=0) & (labels!=1)))
        stop("labels variable must take values 0 or 1");
    if(!(abs(folds - round(folds)) < .Machine$double.eps^0.5))
        stop("folds must be an integer number");
    if(is.null(seed))
        warning("folds are generated without seed initialization");
    indices <- 1:length(labels);
    positives <- which(labels==1);
    foldIndex <- stratified.cv.data.single.class(indices, positives, kk=folds, seed=seed);
    testIndex <- mapply(c, foldIndex$fold.positives, foldIndex$fold.negatives, SIMPLIFY=FALSE);
    fold.check <- unlist(lapply(testIndex,length));
    if(any(fold.check==0))
        stop("number of folds selected too high: some folds have no examples. Please reduce the number of folds");
    fold.len <- sapply(testIndex,length);
    tmp.list <- vector(mode="list", length=folds);
    for(k in 1:folds)
        tmp.list[[k]] <- rep(k, length(testIndex[[k]]));
    foldcol <- numeric(length(scores));
    labelschar <- ifelse(labels==1, "pos", "neg");
    df <- data.frame(scores, labelschar, foldcol, stringsAsFactors=TRUE);
    for(i in 1:length(tmp.list))
        df$foldcol[testIndex[[i]]] <- tmp.list[[i]];
    names(df) <- c("scores","labels","folds");
    return(df);
}