@@ -1,8 +1,8 @@

 Package: kebabs

 Type: Package

 Title: Kernel-Based Analysis Of Biological Sequences

-Version: 1.1.5

-Date: 2014-01-12

+Version: 1.1.6

+Date: 2015-01-21

 Author: Johannes Palme

 Maintainer: Ulrich Bodenhofer <bodenhofer@bioinf.jku.at>

 Description: The package provides functionality for kernel-based analysis of

@@ -54,4 +54,4 @@ Imports: methods, stats, Rcpp (>= 0.11.2), Matrix, XVector (>= 0.5.8),

 LinkingTo: IRanges, XVector, Biostrings, Rcpp, S4Vectors

 Suggests: SparseM, apcluster, Biobase, BiocGenerics

 biocViews: SupportVectorMachine, Classification, Clustering, Regression

-Packaged: 2015-01-12 13:37:03 UTC; jo

+Packaged: 2015-01-21 11:10:57 UTC; jo

@@ -136,6 +136,53 @@ ersTransposedAsdgCMatrix <- function(from)

     return(dgC)

 }

+distWeightKernelToString <- function(distWeight)

+{

+    dwLength <- length(distWeight)

+    distWeightChar <- ""

+    

+    if (is.numeric(distWeight))

+    {

+        if (dwLength == 1)

+            distWeightChar <- distWeight

+        else if (dwLength < 7)

+        {

+            distWeightChar <-

+                paste("c(",paste(distWeight,collapse=","),")", sep="")

+        }

+        else

+        {

+            distWeightChar <-

+                paste("c(", paste(distWeight[1:3],collapse=","), " ... ",

+                      paste(distWeight[dwLength - 2:dwLength],collapse=","),

+                      ")", sep="")

+        }

+    }

+    else if (typeof(distWeight) == "closure")

+    {

+        func <- deparse(distWeight)[2]

+        index <- grep("(", strsplit(func, split="")[[1]], fixed=TRUE,

+                      value=FALSE)

+        

+        if (length(index) > 0)

+            distWeightChar <- substr(func, 1, index[1] - 1)

+        {

+            if (distWeightChar %in% c("linWeight", "gaussWeight",

+                                      "expWeight"))

+            {

+                sigma <- get("sigma", envir = environment(distWeight))

+                distWeightChar <- paste(distWeightChar, "(sigma=", sigma, ")",

+                                        sep="")

+            }

+            else if (distWeightChar == "swdWeight")

+                distWeightChar <- paste(distWeightChar, "()", sep="")

+            else

+                distWeightChar <- paste(distWeightChar, "( . . . )", sep="")

+        }

+    }

+    distWeightChar

+}

+

 #' @rdname sequenceKernel

 #' @aliases

 #' seqKernelAsChar

@@ -182,8 +229,8 @@ seqKernelAsChar <- function(from)

     if (length(kernelParameters(from)$distWeight) > 0)

     {

-        kChar <- paste(kChar, ", distWeight=",

-                       kernelParameters(from)$distWeight, sep="")

+        dwString <- distWeightKernelToString(kernelParameters(from)$distWeight)

+        kChar <- paste(kChar, ", distWeight=", dwString, sep="")

     }

     if (!kernelParameters(from)$normalized)

@@ -194,9 +194,21 @@ gappyPairKernel <- function(k=1, m=1, r=1, annSpec=FALSE, distWeight=numeric(0),

     if (!isTRUEorFALSE(annSpec))

         stop("'annSpec' must be TRUE or FALSE\n")

-    if (length(distWeight) > 0 &&

-        !(is.numeric(distWeight) || is.function(distWeight)))

-        stop("'distWeight' must be a numeric vector or a function\n")

+    if (length(distWeight) > 0)

+    {

+        if (!(is.numeric(distWeight) || is.function(distWeight)))

+            stop("'distWeight' must be a numeric vector or a function\n")

+    

+        if (is.function(distWeight))

+        {

+            func <- deparse(distWeight)[2]

+            index <- grep("(", strsplit(func, split="")[[1]], fixed=TRUE,

+                          value=FALSE)

+        

+            if (length(index) < 1)

+                stop("Missing parentheses in 'distWeight'\n")

+        }

+    }

     if (presence && length(distWeight) > 0)

     {

@@ -164,16 +164,14 @@ kbsvm.seqs <- function(x, y, kernel=NULL, pkg="auto", svm="C-svc",

                         if (ncol(exRepSV) == 0)

                             stop("explicit representaton is not available\n")

+                    }

+                    else

+                        exRepSV <- NULL

-                        model@featureWeights <-

-                            getFeatureWeights(model=model,

+                    model@featureWeights <-

+                        getFeatureWeights(model=model,

                                           exrep=exRepSV,

                                           weightLimit=model@svmInfo@weightLimit)

-                    }

-                    else

-                        stop("feature weights cannot be computed because\n",

-                             "        the selected kernel does not support\n",

-                             "        an explicit representation\n")

                     model@b <- getSVMSlotValue("b", model)

                 }

@@ -511,20 +509,31 @@ kbsvm.seqs <- function(x, y, kernel=NULL, pkg="auto", svm="C-svc",

 #' Training of biological sequences with a sequence kernel\cr

 #'

 #' Training is performed via the method \code{kbsvm} for classification and

-#' regression tasks. The user passes sequence data in \code{XStringSet} format,

-#' the label vector, a sequence kernel object and the requested SVM along with

-#' SVM parameters and receives the training results in the form of a KeBABS

-#' model object of class \code{\linkS4class{KBModel}}. The accessor

+#' regression tasks. The user passes sequence data, the response vector, a

+#' sequence kernel object and the requested SVM along with SVM parameters

+#' to \code{kbsvm} and receives the training results in the form of a

+#' KeBABS model object of class \code{\linkS4class{KBModel}}. The accessor

 #' \code{svmModel} allows to retrieve the SVM specific model from the KeBABS

-#' model object. However for regular operation a detailed look into the SVM

-#' specific model is usually not necessary. (When repeat regions are coded as

-#' lowercase characters and should be excluded from the analysis the sequence

-#' data can be passed as \code{\linkS4class{BioVector}} which also supports

-#' lowercase characters instead of \code{\linkS4class{XStringSet}} format.

-#' Please note that the classes derived from \code{\linkS4class{XStringSet}}

-#' are much more powerful than the \code{\linkS4class{BioVector}} derived

-#' classes and should be used in all cases where lowercase characters are not

-#' needed). Apart from SVM training \code{kbsvm} can be also used for cross

+#' model object. However, for regular operation a detailed look into the SVM

+#' specific model is usually not necessary.

+#'

+#' The standard data format for sequences in KeBABS are the

+#' \code{XStringSet}-derived classes \code{\linkS4class{DNAStringSet}},

+#' \code{\linkS4class{RNAStringSet}} and \code{\linkS4class{AAStringSet}}.

+#' (When repeat regions are coded as lowercase characters and should be

+#' excluded from the analysis the sequence data can be passed as

+#' \code{\linkS4class{BioVector}} which also supports lowercase characters

+#' instead of \code{\linkS4class{XStringSet}} format. Please note that the

+#' classes derived from \code{\linkS4class{XStringSet}} are much more

+#' powerful than the \code{\linkS4class{BioVector}} derived classes and

+#' should be used in all cases where lowercase characters are not needed).

+#'

+#' Instead of sequences also a precomputed explicit representation or

+#' a precomputed kernel matrix can be used for training. Examples for

+#' training with kernel matrix and explicit representation can be found on

+#' the help page for the prediction method \code{\link{predict}}.

+#'

+#' Apart from SVM training \code{kbsvm} can be also used for cross

 #' validation (see \link{crossValidation} and parameters \code{cross} and

 #' \code{noCross}), grid search for SVM- and kernel-parameter values (see

 #' \link{gridSearch}) and model selection (see \link{modelSelection} and

@@ -653,9 +662,9 @@ kbsvm.seqs <- function(x, y, kernel=NULL, pkg="auto", svm="C-svc",

 #' with the parameter \code{weightLimit} which defines the cutoff for

 #' small feature weights not stored in the model.

 #'

-#' Hint: For multiclass prediction is currently not performed via feature

-#' weights but native in the SVM. This means that for multiclass a pruned

-#' model cannot be used for prediction.\cr\cr

+#' Hint: For training with a precomputed kernel matrix feature weights are

+#' not available. For multiclass prediction is currently not performed via

+#' feature weights but native in the SVM.\cr\cr

 #'

 #' Cross validation, grid search and model selection\cr

 #'

@@ -754,6 +754,18 @@ predict.KBModel <- function(object, x, predictionType="response", sel=NULL,

 #' do not support the generation of an explicit representation, e.g. the

 #' position dependent kernel variants.\cr\cr

 #'

+#' Prediction with precomputed kernel matrix

+#'

+#' When training was performed with a precomputed kernel matrix also in

+#' prediction a precomputed kernel matrix must be passed to the \code{predict}

+#' method. In contrast to the quadratic and symmetric kernel matrix used

+#' in training the kernel matrix for prediction is rectangular and contains

+#' the similarities of test samples (rows) against support vectors (columns).

+#' support vector indices can be read from the model with the accessor SVindex.

+#' Please not that these indices refer to the sample subset used in training.

+#' An example for training and prediction via precomputed kernel matrix is

+#' shown below.

+#'

 #' Generation of prediction profiles

 #'

 #' The parameter \code{predProfiles} controls whether prediction profiles

@@ -791,15 +803,14 @@ predict.KBModel <- function(object, x, predictionType="response", sel=NULL,

 #'

 #' ## run training with explicit representation

 #' model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=gappy,

-#'                pkg="LiblineaR", svm="C-svc", cost=10, explicit="yes",

-#'                featureWeights="yes")

+#'                pkg="LiblineaR", svm="C-svc", cost=10)

 #'

 #' ## show feature weights in KeBABS model

 #' featureWeights(model)[1:8]

 #'

 #' ## predict the test sequences

 #' pred <- predict(model, enhancerFB[test])

-#' evaluatePrediction(pred, yFB[test], allLabels=c(-1,1))

+#' evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))

 #' pred[1:10]

 #'

 #' ## output decision values instead

@@ -807,11 +818,30 @@ predict.KBModel <- function(object, x, predictionType="response", sel=NULL,

 #' pred[1:10]

 #'

 #' \dontrun{

-#' ## computing of probability model via Platt scaling during training

+#' ## example for training and prediction via precomputed kernel matrix

+#'

+#' ## compute quadratic kernel matrix of training samples

+#' kmtrain <- getKernelMatrix(gappy, x=enhancerFB, selx=train)

+#'

+#' ## train model with kernel matrix

+#' model <- kbsvm(x=kmtrain, y=yFB[train], kernel=gappy,

+#'                pkg="kernlab", svm="C-svc", cost=10)

+#'

+#' ## compute rectangular kernel matrix of test samples versus

+#' ## support vectors

+#' kmtest <- getKernelMatrix(gappy, x=enhancerFB, y=enhancerFB,

+#'                           selx=test, sely=train[SVindex(model)])

+#'

+#' ## predict with kernel matrix

+#' pred <- predict(model, kmtest)

+#' evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))

+#'

+#' ## example for probability model generation during training

+#'

+#' ## compute probability model via Platt scaling during training

 #' ## and predict class membership probabilities

 #' model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=gappy,

-#'                pkg="e1071", svm="C-svc", cost=10, explicit="yes",

-#'                probModel=TRUE)

+#'                pkg="e1071", svm="C-svc", cost=10, probModel=TRUE)

 #'

 #' ## show parameters of the fitted probability model which are the parameters

 #' ## probA and probB for the fitted sigmoid function in case of classification

@@ -358,26 +358,8 @@ setMethod("show", signature(object="SpectrumKernel"),

             cat(paste(", annSpec=TRUE"))

         if (length(object@distWeight) > 0)

         {

-            if (is.numeric(object@distWeight))

-            {

-                if (length(object@distWeight) == 1)

-                {

-                    cat(paste(", distWeight=",

-                              object@distWeight, sep=""))

-                }

-                else

-                {

-                    cat(paste(", distWeight=",

-                              paste("c(",paste(object@distWeight,

-                                               collapse=","),")",

-                                    sep=""), sep=""))

-                }

-            }

-            else

-            {

-                cat(", \ndistWeight=")

-                cat(format(object@distWeight))

-            }

+             dwString <- distWeightKernelToString(object@distWeight)

+             cat(", distWeight=", dwString, sep="")

         }

         if (object@normalized == FALSE)

             cat(", normalized=FALSE")

@@ -170,9 +170,21 @@ spectrumKernel <- function(k=3, r=1, annSpec=FALSE, distWeight=numeric(0),

     if (!isTRUEorFALSE(annSpec))

         stop("'annSpec' must be TRUE or FALSE\n")

-    if (length(distWeight) > 0 &&

-        !(is.numeric(distWeight) || is.function(distWeight)))

-        stop("'distWeight' must be a numeric vector or a function\n")

+    if (length(distWeight) > 0)

+    {

+        if (!(is.numeric(distWeight) || is.function(distWeight)))

+            stop("'distWeight' must be a numeric vector or a function\n")

+        

+        if (is.function(distWeight))

+        {

+            func <- deparse(distWeight)[2]

+            index <- grep("(", strsplit(func, split="")[[1]], fixed=TRUE,

+                          value=FALSE)

+            

+            if (length(index) < 1)

+                stop("Missing parentheses in 'distWeight'\n")

+        }

+    }

     if (presence && length(distWeight) > 0)

     {

@@ -1,6 +1,19 @@

 Change history of package kebabs:

 ====================================

+Version 1.1.6

+- error correction for training with position specific kernel and

+  computation of feature weights

+- error correction in coercion of kernel to character for distance

+  weighting

+- error correction in spectrum, gappy pair and motif kernel

+  for kernel matrix - last feature was missing in kernel value

+  in rare situations

+- correction of Windows build problem in linearKernel

+- build warnings on Windows removed

+- minor changes in help pages

+- minor changes in vignette

+

 Version 1.1.5

 - new method heatmap to display heatmap of prediction profiles

 - extension of function linearKernel to optionally return

@@ -291,30 +291,43 @@ Training of biological sequences with a sequence kernel\cr

 Training is performed via the method \code{kbsvm} for

 classification and regression tasks. The user passes

-sequence data in \code{XStringSet} format, the label

-vector, a sequence kernel object and the requested SVM

-along with SVM parameters and receives the training results

-in the form of a KeBABS model object of class

+sequence data, the response vector, a sequence kernel

+object and the requested SVM along with SVM parameters to

+\code{kbsvm} and receives the training results in the form

+of a KeBABS model object of class

 \code{\linkS4class{KBModel}}. The accessor \code{svmModel}

 allows to retrieve the SVM specific model from the KeBABS

-model object. However for regular operation a detailed look

-into the SVM specific model is usually not necessary. (When

-repeat regions are coded as lowercase characters and should

-be excluded from the analysis the sequence data can be

-passed as \code{\linkS4class{BioVector}} which also

-supports lowercase characters instead of

+model object. However, for regular operation a detailed

+look into the SVM specific model is usually not necessary.

+

+The standard data format for sequences in KeBABS are the

+\code{XStringSet}-derived classes

+\code{\linkS4class{DNAStringSet}},

+\code{\linkS4class{RNAStringSet}} and

+\code{\linkS4class{AAStringSet}}. (When repeat regions are

+coded as lowercase characters and should be excluded from

+the analysis the sequence data can be passed as

+\code{\linkS4class{BioVector}} which also supports

+lowercase characters instead of

 \code{\linkS4class{XStringSet}} format. Please note that

 the classes derived from \code{\linkS4class{XStringSet}}

 are much more powerful than the

 \code{\linkS4class{BioVector}} derived classes and should

 be used in all cases where lowercase characters are not

-needed). Apart from SVM training \code{kbsvm} can be also

-used for cross validation (see \link{crossValidation} and

-parameters \code{cross} and \code{noCross}), grid search

-for SVM- and kernel-parameter values (see

-\link{gridSearch}) and model selection (see

-\link{modelSelection} and parameters \code{nestedCross} and

-\code{noNestedCross}).\cr\cr

+needed).

+

+Instead of sequences also a precomputed explicit

+representation or a precomputed kernel matrix can be used

+for training. Examples for training with kernel matrix and

+explicit representation can be found on the help page for

+the prediction method \code{\link{predict}}.

+

+Apart from SVM training \code{kbsvm} can be also used for

+cross validation (see \link{crossValidation} and parameters

+\code{cross} and \code{noCross}), grid search for SVM- and

+kernel-parameter values (see \link{gridSearch}) and model

+selection (see \link{modelSelection} and parameters

+\code{nestedCross} and \code{noNestedCross}).\cr\cr

 Package and SVM selection\cr

@@ -428,10 +441,10 @@ Pruning of feature weights can be achieved with the

 parameter \code{weightLimit} which defines the cutoff for

 small feature weights not stored in the model.

-Hint: For multiclass prediction is currently not performed

-via feature weights but native in the SVM. This means that

-for multiclass a pruned model cannot be used for

-prediction.\cr\cr

+Hint: For training with a precomputed kernel matrix feature

+weights are not available. For multiclass prediction is

+currently not performed via feature weights but native in

+the SVM.\cr\cr

 Cross validation, grid search and model selection\cr

@@ -107,6 +107,20 @@ variants which do not support the generation of an explicit

 representation, e.g. the position dependent kernel

 variants.\cr\cr

+Prediction with precomputed kernel matrix

+

+When training was performed with a precomputed kernel

+matrix also in prediction a precomputed kernel matrix must

+be passed to the \code{predict} method. In contrast to the

+quadratic and symmetric kernel matrix used in training the

+kernel matrix for prediction is rectangular and contains

+the similarities of test samples (rows) against support

+vectors (columns). support vector indices can be read from

+the model with the accessor SVindex. Please not that these

+indices refer to the sample subset used in training. An

+example for training and prediction via precomputed kernel

+matrix is shown below.

+

 Generation of prediction profiles

 The parameter \code{predProfiles} controls whether

@@ -132,15 +146,14 @@ gappy

 ## run training with explicit representation

 model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=gappy,

-               pkg="LiblineaR", svm="C-svc", cost=10, explicit="yes",

-               featureWeights="yes")

+               pkg="LiblineaR", svm="C-svc", cost=10)

 ## show feature weights in KeBABS model

 featureWeights(model)[1:8]

 ## predict the test sequences

 pred <- predict(model, enhancerFB[test])

-evaluatePrediction(pred, yFB[test], allLabels=c(-1,1))

+evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))

 pred[1:10]

 ## output decision values instead

@@ -148,11 +161,30 @@ pred <- predict(model, enhancerFB[test], predictionType="decision")

 pred[1:10]

 \dontrun{

-## computing of probability model via Platt scaling during training

+## example for training and prediction via precomputed kernel matrix

+

+## compute quadratic kernel matrix of training samples

+kmtrain <- getKernelMatrix(gappy, x=enhancerFB, selx=train)

+

+## train model with kernel matrix

+model <- kbsvm(x=kmtrain, y=yFB[train], kernel=gappy,

+               pkg="kernlab", svm="C-svc", cost=10)

+

+## compute rectangular kernel matrix of test samples versus

+## support vectors

+kmtest <- getKernelMatrix(gappy, x=enhancerFB, y=enhancerFB,

+                          selx=test, sely=train[SVindex(model)])

+

+## predict with kernel matrix

+pred <- predict(model, kmtest)

+evaluatePrediction(pred, yFB[test], allLabels=unique(yFB))

+

+## example for probability model generation during training

+

+## compute probability model via Platt scaling during training

 ## and predict class membership probabilities

 model <- kbsvm(x=enhancerFB[train], y=yFB[train], kernel=gappy,

-               pkg="e1071", svm="C-svc", cost=10, explicit="yes",

-               probModel=TRUE)

+               pkg="e1071", svm="C-svc", cost=10, probModel=TRUE)

 ## show parameters of the fitted probability model which are the parameters

 ## probA and probB for the fitted sigmoid function in case of classification

@@ -900,13 +900,14 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                       int sizeX, int sizeY, IntegerVector selX, IntegerVector selY,

                       ByteStringVector annCharset, ByteStringVector annX, ByteStringVector annY,

                       bool unmapped, int k, int m, bool normalized, bool symmetric, bool presence,

-                      bool reverseComplement, int maxSeqLength, struct alphaInfo *alphaInf)

+                      bool reverseComplement, int maxSeqLength, uint64_t dimFeatureSpace,

+                      struct alphaInfo *alphaInf)

 {

     T featureIndex, annotIndex;

     int i, j, l, iX, iY, elemIndex, noSamples = sizeX;

     int freeNode, nodeLimit, currBlock, currIndex, newBlock, maxBlockIndex;

     int maxNoOfNodes, seqnchar, currStack, stack[8*k];

-    uint64_t maxNodesPerSequence, numAnnPow2K = 0;

+    uint64_t maxNodesPerSequence, maxNumFeatures, numAnnPow2K = 0;

     double kv;

     bool printWarning = TRUE;

     const char *seqptr, *annptr;

@@ -928,14 +929,17 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

         noSamples += sizeY;

     // add one for the sentinel

-    maxSeqLength += 1;

-

+    if (dimFeatureSpace < maxSeqLength)

+        maxNumFeatures = dimFeatureSpace + 1;

+    else

+        maxNumFeatures = (maxSeqLength - 2*k - 0.5*m + 1) * (m + 1) + 1;

+    

     // allocate arrays for sparse feature vectors with 32 or 64 bit index

     // store only unnormalized k-mer count to avoid double space usage

     // add one for the sentinel

     featVectorArrays featVectors;

-    featVectors.x   = (int32_t *) R_alloc(noSamples * ((maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1), sizeof(int32_t));

-    featVectors.idx = (T *) R_alloc(noSamples * ((maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1), sizeof(T));

+    featVectors.x   = (int32_t *) R_alloc(noSamples * maxNumFeatures, sizeof(int32_t));

+    featVectors.idx = (T *) R_alloc(noSamples * maxNumFeatures, sizeof(T));

     double *normValues   = (double *) R_alloc(noSamples, sizeof(double));

@@ -1016,12 +1020,13 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             currBlock = 0;

             currIndex = 0;

             currStack = -1;

-            elemIndex = i * ((maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1);

+            elemIndex = i * maxNumFeatures;

             featureIndex = 0;

             annotIndex = 0;

             // for situation with no features set sentinel

             featVectors.idx[elemIndex] = maxUnSignedIndex;

+            featVectors.x[elemIndex] = MAXINT32;

             while (currStack >= 0 || currIndex <= maxBlockIndex)

             {

@@ -1117,10 +1122,11 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             }

             featVectors.idx[elemIndex] = maxUnSignedIndex;

+            featVectors.x[elemIndex] = MAXINT32;

         }

         computeKernelMatrix(maxUnSignedIndex, featVectors.idx, featVectors.x, km, normValues,

-                            (maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1, sizeX, sizeY, normalized, symmetric);

+                            maxNumFeatures, sizeX, sizeY, normalized, symmetric);

     }

     else // not symmetric

     {

@@ -1144,10 +1150,6 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             }

             else

             {

-                // this statement causes problems in the optimization level -O2

-//                iY = selY[i - sizeX];

-//                Rprintf("i: %d  sizeX: %d  i-sizeX: %d\n", i, sizeX, i - sizeX);

-

                 indexY++;

                 iY = selY[indexY];

@@ -1162,8 +1164,8 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             //create new tree for sequence

             kv = createTreeGappy((const char *) seqptr, seqnchar, annptr, k, m, &annotationIndexMap,

-                                 presence, reverseComplement, pTree, maxNoOfNodes, &freeNode, &nullBlock, &printWarning,

-                                 alphaInf);

+                                 presence, reverseComplement, pTree, maxNoOfNodes, &freeNode,

+                                 &nullBlock, &printWarning, alphaInf);

             if (kv == NA_REAL)

             {

@@ -1185,12 +1187,13 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             currBlock = 0;

             currIndex = 0;

             currStack = -1;

-            elemIndex = i * ((maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1);

+            elemIndex = i * maxNumFeatures;

             featureIndex = 0;

             annotIndex = 0;

             // for situation with no features set sentinel

             featVectors.idx[elemIndex] = maxUnSignedIndex;

+            featVectors.x[elemIndex] = MAXINT32;

             while (currStack >= 0 || currIndex <= maxBlockIndex)

             {

@@ -1286,10 +1289,11 @@ void getKMStdAnnGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             }

             featVectors.idx[elemIndex] = maxUnSignedIndex;

+            featVectors.x[elemIndex] = MAXINT32;

         }

         computeKernelMatrix(maxUnSignedIndex, featVectors.idx, featVectors.x, km, normValues,

-                            (maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1, sizeX, sizeY, normalized, symmetric);

+                            maxNumFeatures, sizeX, sizeY, normalized, symmetric);

     }

     return;

@@ -1303,9 +1307,9 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                        NumericVector distWeight, int maxSeqLength, struct alphaInfo *alphaInf)

 {

     T prevIndex;

-    uint32_t km1, inew, invalidCharPos, elemIndex, tempIndex1, tempIndex2, fIndex;

-    uint64_t featureIndex, fDimArray, numAlphaPowerK, numAlphaPowerK1;

-    int i, j, l, offset, noSamples = sizeX;

+    uint32_t km1, inew, elemIndex, tempIndex1, tempIndex2, fIndex;

+    uint64_t featureIndex, fDimArray, invalidCharPos, numAlphaPowerK, numAlphaPowerK1, *featVectorsStart;

+    int i, j, l, offset, maxNumPatternsPerPos, maxFeaturesPerSample, noSamples = sizeX;

     int patLength, seqnchar, index, iold, iX, iY, upper;

     const char *seqptr;

     double kv, *normValues;

@@ -1321,19 +1325,25 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

         noSamples += sizeY;

     km1 = k + m + 1;

-    T *oldIndex = (T *) R_alloc(k, sizeof(uint64_t));

-    T *newIndex = (T *) R_alloc(km1, sizeof(uint64_t));

+    uint64_t *oldIndex = (uint64_t *) R_alloc(k, sizeof(uint64_t));

+    uint64_t *newIndex = (uint64_t *) R_alloc(km1, sizeof(uint64_t));

     numAlphaPowerK1 = ipow64(alphaInf->numAlphabetChars, k - 1);

     numAlphaPowerK = numAlphaPowerK1 * alphaInf->numAlphabetChars;

     // allocate arrays for feature vectors with 8, 16, 32 or 64 bit index

-    fDimArray = (maxSeqLength - 2*k + 0.5*m) * (m + 1) + 1;

+    fDimArray = (maxSeqLength - 2*k - 0.5*m + 1) * (m + 1);

     featVectorArrays featVectors;

     featVectors.x   = (int32_t *) R_alloc(noSamples * fDimArray, sizeof(int32_t));

     featVectors.idx = (T *) R_alloc(noSamples * fDimArray, sizeof(T));

+    featVectorsStart = (uint64_t *) R_alloc(noSamples + 1, sizeof(uint64_t));

     normValues = (double *) R_alloc(noSamples, sizeof(double));

+    maxNumPatternsPerPos = m + 1;

+    featVectorsStart[0] = 0;

+    maxFeaturesPerSample = 0;

+    elemIndex = 0;

+    

     // walk along sequence and create sparse feature vector

     for (i = 0; i < noSamples; i++)

     {

@@ -1362,18 +1372,17 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

         patLength = 0;

         featureIndex = 0;

-        elemIndex = i * fDimArray;

         iold = 0;

         inew = 0;

         kv = 0;

-        invalidCharPos = maxUnSignedIndex;

+        invalidCharPos = MAXUINT64;

         for (l = 0; l < seqnchar - 2 * k + 1; l++)

         {

-            if (invalidCharPos != maxUnSignedIndex && l <= (int)invalidCharPos)

+            if (invalidCharPos != MAXUINT64 && l <= (int)invalidCharPos)

                 continue;

             else

-                invalidCharPos = maxUnSignedIndex;

+                invalidCharPos = MAXUINT64;

             // look ahead and fill newIndex buffer to process

             // each of the m + 1 patterns at a position immediately

@@ -1386,9 +1395,9 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             {

                 for (j = 0; j < upper; j++)

                 {

-                    if (invalidCharPos != maxUnSignedIndex)

+                    if (invalidCharPos != MAXUINT64)

                     {

-                        newIndex[inew++] = maxUnSignedIndex;

+                        newIndex[inew++] = MAXUINT64;

                         if (inew == km1)

                             inew = 0;

@@ -1438,7 +1447,7 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                     else

                     {

                         //store max index to buffer for not valid pattern

-                        newIndex[inew++] = maxUnSignedIndex;

+                        newIndex[inew++] = MAXUINT64;

                         if (inew == km1)

                             inew = 0;

@@ -1452,7 +1461,7 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                 }

             }

-            if (l + 2 * k + m - 1 < seqnchar && invalidCharPos == maxUnSignedIndex)

+            if (l + 2 * k + m - 1 < seqnchar && invalidCharPos == MAXUINT64)

             {

                 //calc next look ahead feature index and store to buffer

                 index = alphaInf->seqIndexMap[(int) seqptr[l + 2 * k + m - 1]];

@@ -1475,7 +1484,7 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                 else

                 {

                     //store max index to buffer for not valid pattern

-                    newIndex[inew++] = maxUnSignedIndex;

+                    newIndex[inew++] = MAXUINT64;

                     if (inew == km1)

                         inew = 0;

@@ -1499,13 +1508,13 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

             else

                 upper = seqnchar - 2 * k - l;

-            if (invalidCharPos != maxUnSignedIndex && upper > (int) invalidCharPos - l - 2 * k)

+            if (invalidCharPos != MAXUINT64 && upper > (int) invalidCharPos - l - 2 * k)

                 upper = invalidCharPos - l - 2 * k;

             for (j = 0; j <= upper; j++)

             {

-                if (newIndex[inew] == maxUnSignedIndex ||

-                    newIndex[(inew + k + j) % km1] == maxUnSignedIndex)

+                if (newIndex[inew] == MAXUINT64 ||

+                    newIndex[(inew + k + j) % km1] == MAXUINT64)

                     continue;

                 if (reverseComplement)

@@ -1537,11 +1546,14 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

                 kv++;

             }

-            if (invalidCharPos != maxUnSignedIndex)

+            if (invalidCharPos != MAXUINT64)

                 inew = 0;

         }

-        featVectors.idx[elemIndex] = maxUnSignedIndex;

+        featVectorsStart[i + 1] = elemIndex;

+

+        if (maxFeaturesPerSample < (featVectorsStart[i + 1] - featVectorsStart[i]))

+            maxFeaturesPerSample = featVectorsStart[i + 1] - featVectorsStart[i];

         // for dist weighting the kernel value is determined in computeKernelMatrixPos

         if (distWeight.length() == 0)

@@ -1553,8 +1565,9 @@ void getKMPosDistGappy(T maxUnSignedIndex, NumericMatrix km, ByteStringVector x,

         }

     }

-    computeKernelMatrixPos(maxUnSignedIndex, featVectors.idx, featVectors.x, km, normValues, fDimArray,

-                           m + 1, sizeX, sizeY, normalized, symmetric, FALSE, distWeight);

+    computeKernelMatrixPos(maxUnSignedIndex, featVectors.idx, featVectors.x, featVectorsStart, km,

+                           normValues, maxFeaturesPerSample, maxNumPatternsPerPos, sizeX, sizeY,

+                           normalized, symmetric, FALSE, distWeight);

     return;

 }

@@ -3464,7 +3477,7 @@ RcppExport SEXP gappyPairKernelMatrixC(SEXP xR, SEXP yR, SEXP selXR, SEXP selYR,

             {

                 getKMStdAnnGappy(maxUIndex8, km, x, y, sizeX, sizeY, selX, selY, annCharset, annX, annY,

                                  unmapped, k, m, normalized, symmetric, presence, reverseComplement,

-                                 maxSeqLength, &alphaInf);

+                                 maxSeqLength, dimFeatureSpace, &alphaInf);

             }

             break;

@@ -3482,7 +3495,7 @@ RcppExport SEXP gappyPairKernelMatrixC(SEXP xR, SEXP yR, SEXP selXR, SEXP selYR,

             {

                 getKMStdAnnGappy(maxUIndex16, km, x, y, sizeX, sizeY, selX, selY, annCharset, annX, annY,

                                  unmapped, k, m, normalized, symmetric, presence, reverseComplement,

-                                 maxSeqLength, &alphaInf);

+                                 maxSeqLength, dimFeatureSpace, &alphaInf);

             }

             break;

@@ -3501,7 +3514,7 @@ RcppExport SEXP gappyPairKernelMatrixC(SEXP xR, SEXP yR, SEXP selXR, SEXP selYR,

             {

                 getKMStdAnnGappy(maxUIndex32, km, x, y, sizeX, sizeY, selX, selY, annCharset, annX, annY,

                                  unmapped, k, m, normalized, symmetric, presence, reverseComplement,

-                                 maxSeqLength, &alphaInf);

+                                 maxSeqLength, dimFeatureSpace, &alphaInf);

             }

             break;

@@ -3519,7 +3532,7 @@ RcppExport SEXP gappyPairKernelMatrixC(SEXP xR, SEXP yR, SEXP selXR, SEXP selYR,

             {

                 getKMStdAnnGappy(maxUIndex64, km, x, y, sizeX, sizeY, selX, selY, annCharset, annX, annY,

                                  unmapped, k, m, normalized, symmetric, presence, reverseComplement,

-                                 maxSeqLength, &alphaInf);

+                                 maxSeqLength, dimFeatureSpace, &alphaInf);

             }

             break;

@@ -3536,7 +3549,7 @@ uint64_t * featureNamesToIndexGappyPair(SEXP featureNames, int numFeatures, Byte

                                         struct alphaInfo *alphaInf)

 {

     int i, j, l, numDots;

-    uint64_t featureIndex, *featIndex, annIndex, numAnnPow2K, tempIndex, fIndex;

+    uint64_t featureIndex, *featIndex, annIndex, tempIndex, fIndex, numAnnPow2K = 0;

     const char *pattern;

     featIndex = (uint64_t *) R_alloc(numFeatures, sizeof(uint64_t));

@@ -3600,7 +3613,7 @@ RcppExport void featuresToHashmapGappyPair(NumericMatrix featureWeights, int svm

                                            IntegerVector annotationIndexMap)

 {

     int i, j, numFeatures, result, numDots;

-    uint64_t featureIndex, annotIndex, numAnnPow2K;

+    uint64_t featureIndex, annotIndex, numAnnPow2K = 0;

     const char *pattern;

     SEXP dimNames, colNames;

     khiter_t iter;

@@ -18,6 +18,7 @@

 #define CHAR_UNMAPPED        FALSE

 #define MAXUINT8             0xFF

 #define MAXUINT16            0xFFFF

+#define MAXINT32             0x7FFFFFFF

 #define MAXUINT32            0xFFFFFFFF

 #define MAXUINT64            0xFFFFFFFFFFFFFFFF

 #define MAXINDEX32           MAXUINT32 - 1

@@ -237,11 +237,11 @@ void computeKernelMatrix(T maxUnSignedIndex, T *featVectorIndex, int32_t *featVe

                 while (ix < endx && iy < endy)

                 {

-                    if (featVectorIndex[ix] == maxUnSignedIndex ||

-                        featVectorIndex[iy] == maxUnSignedIndex)

-                    {

+                    if ((featVectorIndex[ix] == maxUnSignedIndex &&

+                         featVectorValue[ix] == MAXINT32)||

+                        (featVectorIndex[iy] == maxUnSignedIndex &&

+                         featVectorValue[iy] == MAXINT32))

                         break;

-                    }

                     if (featVectorIndex[ix] < featVectorIndex[iy])

                         ix++;

@@ -285,8 +285,10 @@ void computeKernelMatrix(T maxUnSignedIndex, T *featVectorIndex, int32_t *featVe

                 while (ix < endx && iy < endy)

                 {

-                    if (featVectorIndex[ix] == maxUnSignedIndex ||

-                        featVectorIndex[iy] == maxUnSignedIndex)

+                    if ((featVectorIndex[ix] == maxUnSignedIndex &&

+                         featVectorValue[ix] == MAXINT32)||

+                        (featVectorIndex[iy] == maxUnSignedIndex &&

+                         featVectorValue[iy] == MAXINT32))

                         break;

                     if (featVectorIndex[ix] < featVectorIndex[iy])

@@ -317,9 +319,9 @@ void computeKernelMatrix(T maxUnSignedIndex, T *featVectorIndex, int32_t *featVe

 template<typename T>

 void computeKernelMatrixPos(T maxUnSignedIndex, T *featVectorIndex, int32_t *featVectorValue,

-                            Rcpp::NumericMatrix km, double *normValues, uint32_t fDim, int maxNumPatterns,

-                            int sizeX, int sizeY, bool normalized, bool symmetric, bool computePosition,

-                            Rcpp::NumericVector distWeight)

+                            uint64_t *featVectorStart, Rcpp::NumericMatrix km, double *normValues,

+                            uint32_t fDim, int maxNumPatterns, int sizeX, int sizeY, bool normalized,

+                            bool symmetric, bool computePosition, Rcpp::NumericVector distWeight)

 {

     uint32_t endx, endy, ix, iy, prevIndex;

     int i, j, j1, j2, posX, posY, distWeightLength, startIndex, numSamples, yOffset, numPatterns1, numPatterns2;

@@ -360,10 +362,10 @@ void computeKernelMatrixPos(T maxUnSignedIndex, T *featVectorIndex, int32_t *fea

             for (j = startIndex; j < sizeY; j++)

             {

-                ix = i * fDim;

-                iy = (yOffset + j) * fDim;

-                endx = ix + fDim;

-                endy = iy + fDim;

+                ix = featVectorStart[i];

+                iy = featVectorStart[yOffset + j];

+                endx = featVectorStart[i + 1];

+                endy = featVectorStart[yOffset + j + 1];

                 prevIndex = 0;

                 kv = 0;

@@ -378,10 +380,6 @@ void computeKernelMatrixPos(T maxUnSignedIndex, T *featVectorIndex, int32_t *fea

                     while (ix < endx && iy < endy)

                     {

-                        if (featVectorIndex[ix] == maxUnSignedIndex ||

-                            featVectorIndex[iy] == maxUnSignedIndex)

-                            break;

-

                         if (featVectorIndex[ix] == featVectorIndex[iy])

                             kv++;

@@ -393,10 +391,6 @@ void computeKernelMatrixPos(T maxUnSignedIndex, T *featVectorIndex, int32_t *fea

                 {

                     while (ix < endx && iy < endy)

                     {