@@ -1,61 +1,73 @@
Loading
1 1
#' (faster) Linear SVM classifier for texts
2 2
#'
3 -
#' Fit a fast linear SVM classifier for texts, using the R interface to the
4 -
#' svmlin code by Vikas Sindhwani and S. Sathiya Keerthi for fast linear
5 -
#' transductive SVMs. This is passed through to \code{\link[RSSL]{svmlin}} as
6 -
#' implemented by the \pkg{RSSL} package.
3 +
#' Fit a fast linear SVM classifier for sparse text matrices, using svmlin C++
4 +
#' code written by Vikas Sindhwani and S. Sathiya Keerthi.  This method
5 +
#' implements the modified finite Newton L2-SVM method (L2-SVM-MFN) method
6 +
#' described in Sindhwani and Keerthi (2006). Currently,
7 +
#' `textmodel_svmlin()` only works for two-class problems.
7 8
#'
8 -
#' @description This function has been retained for testing purposes only;
9 -
#'   we recommend that you use \code{\link{textmodel_svm}} instead.  That
10 -
#'   function is more efficient, and implements prediction for more than
11 -
#'   two classes.
12 -
#' @param x the \link{dfm} on which the model will be fit.  Does not need to
13 -
#'   contain only the training documents.
9 +
#' @param x the [dfm] on which the model will be fit.  Does not need to contain
10 +
#'   only the training documents.
14 11
#' @param y vector of training labels associated with each document identified
15 -
#'   in \code{train}.  (These will be converted to factors if not already
16 -
#'   factors.)
17 -
#' @param intercept logical; if \code{TRUE}, add an intercept to the data
18 -
#' @param ... additional arguments passed to \code{\link[RSSL]{svmlin}}
19 -
#' @return
20 -
#' \code{textmodel_svmlin()} returns (for now) an object structured as a return
21 -
#' object from \code{\link[RSSL]{svmlin}}.
12 +
#'   in `train`.  (These will be converted to factors if not already factors.)
13 +
#' @param intercept logical; if `TRUE`, add an intercept to the data
14 +
#' @param lambda numeric; regularization parameter lambda (default 1)
15 +
#' @param cp numeric; Relative cost for "positive" examples (the second factor
16 +
#'   level)
17 +
#' @param cn numeric; Relative cost for "negative" examples (the first factor
18 +
#'   level)
19 +
#' @param scale	logical; if `TRUE`, normalize the feature counts
20 +
#' @param center logical; if `TRUE`, centre the feature counts
21 +
#' @return a fitted model object of class `textmodel_svmlin`
22 22
#' @references
23 -
#' Vikas Sindhwani and S. Sathiya Keerthi (2006).  Large Scale Semi-supervised
24 -
#' Linear SVMs. \emph{Proceedings of ACM SIGIR}.
23 +
#' Vikas Sindhwani and S. Sathiya Keerthi (2006).  [Large Scale Semi-supervised
24 +
#' Linear SVMs](https://vikas.sindhwani.org/sk_sigir06.pdf). *Proceedings of ACM
25 +
#' SIGIR*. August 6–11, 2006, Seattle.
25 26
#'
26 -
#' V. Sindhwani and S. Sathiya Keerthi (2006).  Newton Methods for Fast Solution of Semi-supervised
27 -
#' Linear SVMs. Book Chapter in \emph{Large Scale Kernel Machines}, MIT Press, 2006.
27 +
#' V. Sindhwani and S. Sathiya Keerthi (2006).  Newton Methods for Fast Solution
28 +
#' of Semi-supervised Linear SVMs. Book Chapter in *Large Scale Kernel
29 +
#' Machines*, MIT Press, 2006.
28 30
#'
29 -
#' @seealso \code{\link[RSSL]{svmlin}}, \code{text{textmodel_svm}}
31 +
#' @seealso [predict.textmodel_svmlin()]
30 32
#' @examples
31 33
#' # use Lenihan for govt class and Bruton for opposition
32 34
#' quanteda::docvars(data_corpus_irishbudget2010, "govtopp") <-
33 35
#'     c("Govt", "Opp", rep(NA, 12))
34 36
#' dfmat <- quanteda::dfm(data_corpus_irishbudget2010)
35 37
#'
36 -
#' tmod <- textmodel_svmlin(dfmat, y = quanteda::docvars(dfmat, "govtopp"),
37 -
#'                          pos_frac = 5/14)
38 +
#' tmod <- textmodel_svmlin(dfmat, y = quanteda::docvars(dfmat, "govtopp"))
38 39
#' predict(tmod)
39 -
#'
40 -
#' predict(textmodel_svmlin(dfmat, y = quanteda::docvars(dfmat, "govtopp"),
41 -
#'                          intercept = FALSE, pos_frac = 5/14))
42 40
#' @importFrom quanteda dfm_group as.dfm
43 41
#' @importFrom stats na.omit predict
44 42
#' @keywords textmodel
45 43
#' @export
46 -
textmodel_svmlin <- function(x, y, intercept = TRUE, ...) {
44 +
textmodel_svmlin <- function(x, y, intercept = TRUE, # x_u = NULL,
45 +
                             lambda = 1,
46 +
                             # algorithm = 1,
47 +
                             # lambda_u = 1, max_switch = 10000, # pos_frac = 0.5,
48 +
                             cp = 1, cn = 1,
49 +
                             scale = FALSE, center = FALSE) {
47 50
    UseMethod("textmodel_svmlin")
48 51
}
49 52
50 53
#' @export
51 -
textmodel_svmlin.default <- function(x, y, intercept = TRUE, ...) {
52 -
    stop(friendly_class_undefined_message(class(x), "textmodel_svmlin"))
54 +
textmodel_svmlin.default <-  function(x, y, intercept = TRUE, # x_u = NULL,
55 +
                                      lambda = 1,
56 +
                                      # algorithm = 1,
57 +
                                      # lambda_u = 1, max_switch = 10000, # pos_frac = 0.5,
58 +
                                      cp = 1, cn = 1,
59 +
                                      scale = FALSE, center = FALSE) {
60 +
  stop(friendly_class_undefined_message(class(x), "textmodel_svmlin"))
53 61
}
54 62
55 63
#' @export
56 -
#' @importFrom RSSL svmlin
57 64
#' @importFrom methods as
58 -
textmodel_svmlin.dfm <- function(x, y, intercept = TRUE, ...) {
65 +
textmodel_svmlin.dfm <-  function(x, y, intercept = TRUE, # x_u = NULL,
66 +
                                  lambda = 1,
67 +
                                  # algorithm = 1,
68 +
                                  # lambda_u = 1, max_switch = 10000, # pos_frac = 0.5,
69 +
                                  cp = 1, cn = 1,
70 +
                                  scale = FALSE, center = FALSE) {
59 71
    x <- as.dfm(x)
60 72
    if (!sum(x)) stop(message_error("dfm_empty"))
61 73
    call <- match.call()
@@ -67,19 +79,25 @@
Loading
67 79
    class <- y[!is.na(y)]
68 80
    temp <- dfm_group(temp, class, force = TRUE)
69 81
70 -
    svmlinfitted <- RSSL::svmlin(X = as(temp, "dgCMatrix"), X_u = NULL,
71 -
                                 y = factor(docnames(temp), levels = docnames(temp)),
72 -
                                 intercept = intercept,
73 -
                                 ...)
82 +
    svmlinfitted <- svmlin(X = as(temp, "dgCMatrix"), #X_u = x_u,
83 +
                           y = factor(docnames(temp), levels = docnames(temp)),
84 +
                           intercept = intercept,
85 +
                           algorithm = 1, lambda = lambda,
86 +
                           # lambda_u = lambda_u,
87 +
                           # max_switch = max_switch,
88 +
                           # pos_frac = pos_frac,
89 +
                           Cp = cp, Cn = cn, verbose = FALSE,
90 +
                           scale = scale, x_center = center)
91 +
74 92
    result <- list(
75 93
        x = x, y = y,
76 -
        weights = svmlinfitted@weights,
77 -
        algorithm = factor(svmlinfitted@algorithm, levels = 0:3,
78 -
                           labels = c("Regularized Least Squares Classification",
79 -
                                      "SVM",
80 -
                                      "Multi-switch Transductive SVM",
81 -
                                      "Deterministic Annealing Semi-supervised SVM")),
82 -
        classnames = svmlinfitted@classnames,
94 +
        weights = svmlinfitted$weights,
95 +
        # algorithm = factor(svmlinfitted$algorithm, levels = 0:3,
96 +
        #                    labels = c("Regularized Least Squares Classification",
97 +
        #                               "SVM",
98 +
        #                               "Multi-switch Transductive SVM",
99 +
        #                               "Deterministic Annealing Semi-supervised SVM")),
100 +
        # classnames = svmlinfitted$classnames,
83 101
        intercept = intercept,
84 102
        call = call
85 103
    )
@@ -94,18 +112,19 @@
Loading
94 112
95 113
#' Prediction from a fitted textmodel_svmlin object
96 114
#'
97 -
#' \code{predict.textmodel_svmlin()} implements class predictions from a fitted
115 +
#' `predict.textmodel_svmlin()` implements class predictions from a fitted
98 116
#' linear SVM model.
99 117
#' @param object a fitted linear SVM textmodel
100 118
#' @param newdata dfm on which prediction should be made
101 119
#' @param type the type of predicted values to be returned; see Value
102 -
#' @param force make newdata's feature set conformant to the model terms
120 +
#' @param force logical, if `TRUE`, make newdata's feature set conformant to the
121 +
#'   model terms
103 122
#' @param ... not used
104 -
#' @return \code{predict.textmodel_svmlin} returns either a vector of class
105 -
#'   predictions for each row of \code{newdata} (when \code{type = "class"}), or
106 -
#'   a document-by-class matrix of class probabilities (when \code{type =
107 -
#'   "probability"}).
108 -
#' @seealso \code{\link{textmodel_svmlin}}
123 +
#' @return `predict.textmodel_svmlin` returns either a vector of class
124 +
#'   predictions for each row of `newdata` (when `type = "class"`), or
125 +
#'   a document-by-class matrix of class probabilities (when `type =
126 +
#'   "probability"`).
127 +
#' @seealso [textmodel_svmlin()]
109 128
#' @importFrom stats predict
110 129
#' @method predict textmodel_svmlin
111 130
#' @keywords textmodel internal
@@ -132,8 +151,9 @@
Loading
132 151
    pred_y <- as.numeric(data %*% object$weights)
133 152
    names(pred_y) <- docnames(data)
134 153
154 +
    classnames <- levels(object$y)
135 155
    if (type == "class") {
136 -
        pred_y <- ifelse(pred_y < 0, object$classnames[1], object$classnames[2])
156 +
        pred_y <- ifelse(pred_y < 0, classnames[1], classnames[2])
137 157
    } else if (type == "probability") {
138 158
        stop("probability type not implemented yet")
139 159
    }
@@ -155,7 +175,7 @@
Loading
155 175
}
156 176
157 177
#' summary method for textmodel_svmlin objects
158 -
#' @param object output from \code{\link{textmodel_svmlin}}
178 +
#' @param object output from [textmodel_svmlin()]
159 179
#' @param n how many coefficients to print before truncating
160 180
#' @param ... additional arguments not used
161 181
#' @keywords textmodel internal
@@ -189,3 +209,107 @@
Loading
189 209
print.predict.textmodel_svmlin <- function(x, ...) {
190 210
    print(unclass(x))
191 211
}
212 +
213 +
214 +
215 +
# adapted from RSSL -----------------
216 +
217 +
svmlin <- function(X, y, X_u = NULL, algorithm = 1, lambda = 1, lambda_u = 1, max_switch = 10000,
218 +
                   pos_frac = 0.5, Cp = 1.0, Cn = 1.0, verbose = FALSE,
219 +
                   intercept = TRUE, scale = FALSE, x_center = FALSE) {
220 +
221 +
  use_Xu_for_scaling <- TRUE
222 +
  classnames <- levels(y)
223 +
224 +
    if (scale || x_center) {
225 +
        if (intercept) {
226 +
            cols <- 2:ncol(X) # do not scale the intercept column
227 +
        } else {
228 +
            cols <- 1:ncol(X)
229 +
        }
230 +
231 +
        if (!is.null(X_u) && use_Xu_for_scaling) {
232 +
            Xe <- rbind(X, X_u)
233 +
            scaling <- scaleMatrix(Xe[, cols, drop = FALSE], center = TRUE, scale = scale)
234 +
            X[, cols] <- predict(scaling, as.matrix(X[, cols, drop=FALSE]))
235 +
            X_u[, cols] <- predict(scaling, as.matrix(X_u[, cols, drop=FALSE]))
236 +
        } else {
237 +
            scaling <- scaleMatrix(X[, cols, drop = FALSE], center = TRUE, scale = scale)
238 +
            X[, cols] <- predict(scaling, as.matrix(X[, cols, drop=FALSE]))
239 +
            if (!is.null(X_u)) {
240 +
                X_u[, cols] <- predict(scaling,as.matrix(X_u[, cols, drop=FALSE]))
241 +
            }
242 +
        }
243 +
    } else {
244 +
        scaling = NULL
245 +
    }
246 +
247 +
    y <- as.numeric(y) * 2 - 3
248 +
249 +
    # Combine feature matrices, add intercept and transpose them to conform to the C++ datastructure
250 +
    if (is.null(X_u) || algorithm < 2) {
251 +
        if (intercept) {
252 +
            X <- cbind(X, 1)
253 +
            X_u <- cbind(X_u, 1)
254 +
        }
255 +
        Xall <- Matrix::t(X)
256 +
    } else {
257 +
        if (intercept) {
258 +
            X <- cbind(X, 1)
259 +
            X_u <- cbind(X_u, 1)
260 +
        }
261 +
        Xall <- Matrix::t(Matrix::rbind2(X, X_u))
262 +
        y <- c(y, rep(0,nrow(X_u)))
263 +
    }
264 +
265 +
    # Determine costs
266 +
    costs <- rep(1, ncol(Xall))
267 +
    if (algorithm < 1) {
268 +
        costs[y < 0] <- Cn
269 +
        costs[y > 0] <- Cp
270 +
    }
271 +
272 +
    res <- svmlin_rcpp(X = Xall, y = y, l = nrow(X), algorithm = algorithm,
273 +
                       lambda = lambda, lambda_u = lambda_u, max_switch = max_switch,
274 +
                       pos_frac = pos_frac, Cp = Cp, Cn = Cn, costs = costs,
275 +
                       verbose = verbose)
276 +
277 +
    list(classnames = classnames,
278 +
         weights = res$Weights,
279 +
         algorithm = algorithm,
280 +
         scaling = scaling,
281 +
         intercept = intercept)
282 +
}
283 +
284 +
setClass("scaleMatrix",
285 +
         representation(mean = "ANY", scale = "ANY"))
286 +
287 +
scaleMatrix <- function(x, center = TRUE, scale = TRUE) {
288 +
  if (center) {
289 +
    mean <- quanteda::colMeans(x)
290 +
    x <- sweep(x, 2, mean)
291 +
  } else {
292 +
    mean <- NULL
293 +
  }
294 +
295 +
  if (scale) {
296 +
    scale <- sqrt(colSums(sweep(x, 2, quanteda::colMeans(x)) ^ 2) / (nrow(x) - 1))
297 +
    x <- sweep(x, 2, scale, "/")
298 +
  } else {
299 +
    scale <- NULL
300 +
  }
301 +
302 +
  object <- new(Class="scaleMatrix", mean = mean, scale = scale)
303 +
  return(object)
304 +
}
305 +
306 +
setMethod("predict", signature=c("scaleMatrix"), function(object, newdata, ...) {
307 +
    if (!is.matrix(newdata)) stop("Incorrect newdata")
308 +
    if (!is.null(object@mean)) {
309 +
        newdata <- sweep(newdata, 2, object@mean)
310 +
    }
311 +
    if (!is.null(object@scale)) {
312 +
        newdata <- sweep(newdata, 2, object@scale, "/")
313 +
    }
314 +
    return(newdata)
315 +
})

@@ -0,0 +1,129 @@
Loading
1 +
/*    Copyright 2006 Vikas Sindhwani (vikass@cs.uchicago.edu)
2 +
 SVM-lin: Fast SVM Solvers for Supervised and Semi-supervised Learning
3 +
4 +
 This file is part of SVM-lin.
5 +
6 +
 SVM-lin is free software; you can redistribute it and/or modify
7 +
 it under the terms of the GNU General Public License as published by
8 +
 the Free Software Foundation; either version 2 of the License, or
9 +
 (at your option) any later version.
10 +
11 +
 SVM-lin is distributed in the hope that it will be useful,
12 +
 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 +
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 +
 GNU General Public License for more details.
15 +
16 +
 You should have received a copy of the GNU General Public License
17 +
 along with SVM-lin (see gpl.txt); if not, write to the Free Software
18 +
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
19 +
 */
20 +
21 +
#include <Rcpp.h>
22 +
#include <iostream>
23 +
#include <stdlib.h>
24 +
#include <stdio.h>
25 +
#include <ctype.h>
26 +
#include <cstring>
27 +
#include "ssl.h"
28 +
using namespace Rcpp;
29 +
30 +
struct options *Options = new options[1];
31 +
struct data *Data = new data[1];
32 +
struct vector_double *Weights = new vector_double[1];
33 +
struct vector_double *Outputs = new vector_double[1];
34 +
struct vector_double *Labels = new vector_double[1];
35 +
36 +
// [[Rcpp::export]]
37 +
List svmlin_rcpp(S4 X,
38 +
                 NumericVector y,
39 +
                 int l,
40 +
                 int algorithm,
41 +
                 double lambda,
42 +
                 double lambda_u,
43 +
                 int max_switch,
44 +
                 double pos_frac,
45 +
                 double Cp,
46 +
                 double Cn,
47 +
                 NumericVector costs,
48 +
                 bool verbose) {
49 +
  // Set options
50 +
  Options->algo = algorithm;
51 +
  Options->lambda=lambda;
52 +
  Options->lambda_u=lambda_u;
53 +
  Options->S=max_switch;
54 +
  Options->R=pos_frac;
55 +
  Options->epsilon=EPSILON;
56 +
  Options->cgitermax=CGITERMAX;
57 +
  Options->mfnitermax=MFNITERMAX;
58 +
  Options->Cp = Cp;
59 +
  Options->Cn = Cn;
60 +
  Options->verbose = verbose;
61 +
62 +
  NumericVector ycop( y.begin(), y.end() );
63 +
  NumericVector costcop( costs.begin(), costs.end() );
64 +
//   Rprintf("Step 1\n");
65 +
//   size_t size = ((DoubleVector)X.slot("x")).length()+((DoubleVector)Xu.slot("x")).length();
66 +
//   std::vector<double> vals(size);
67 +
//   std::vector<double>::iterator it = vals.begin();
68 +
//   vals.insert(it,((DoubleVector)X.slot("x")).begin(),((DoubleVector)Xu.slot("x")).end());
69 +
//   it = vals.begin()+((DoubleVector)X.slot("x")).length();
70 +
//   vals.insert(it,((DoubleVector)Xu.slot("x")).begin(),((DoubleVector)Xu.slot("x")).end());
71 +
//
72 +
//   Rprintf("Step 2\n");
73 +
//
74 +
//   size = ((IntegerVector)X.slot("i")).length()+((IntegerVector)Xu.slot("i")).length();
75 +
//   std::vector<int> colinds(size);
76 +
//   std::vector<int>::iterator it2 = colinds.begin();
77 +
//   colinds.insert(it2,((IntegerVector)X.slot("i")).begin(),((IntegerVector)X.slot("i")).end());
78 +
//   it2 = colinds.begin() + ((IntegerVector)X.slot("i")).length();
79 +
//   colinds.insert(it2,((IntegerVector)Xu.slot("i")).begin(),((IntegerVector)Xu.slot("i")).end());
80 +
//
81 +
//   size = ((IntegerVector)X.slot("p")).length()+((IntegerVector)Xu.slot("p")).length();
82 +
//   std::vector<int> rowpts(size);
83 +
//   it2 = rowpts.begin();
84 +
//   rowpts.insert(it2,((IntegerVector)X.slot("p")).begin(),((IntegerVector)X.slot("p")).end());
85 +
//   it2 = rowpts.begin() + ((IntegerVector)X.slot("p")).length();
86 +
//   rowpts.insert(it2,((IntegerVector)Xu.slot("p")).begin(),((IntegerVector)Xu.slot("p")).end());
87 +
88 +
  // R data to svmlin data structure
89 +
  Data->m=((IntegerVector)X.slot("Dim"))[1];
90 +
  Data->l=l;
91 +
  Data->u=Data->m-Data->l;
92 +
  Data->n=((IntegerVector)X.slot("Dim"))[0];
93 +
  Data->nz=((DoubleVector)X.slot("x")).size();
94 +
  Data->val=((DoubleVector)X.slot("x")).begin();
95 +
  Data->rowptr=((IntegerVector)X.slot("p")).begin();
96 +
  Data->colind=((IntegerVector)X.slot("i")).begin();
97 +
  Data->Y=ycop.begin();
98 +
  Data->C=costcop.begin();
99 +
100 +
  // TODO: load correct costs for unlabeled data.
101 +
  if (Options->verbose) {
102 +
    Rcout << "  Input Data Matrix Statistics:" << endl;
103 +
    Rcout << "      Examples: " << Data->m << endl;
104 +
    Rcout << "      Features: " << Data->n << " (including bias feature)" << endl;
105 +
    Rcout << "      Non-zeros:  " << Data->nz << " (including bias features)" << endl;
106 +
    Rcout << "      Average sparsity: " << Data->nz*1.0/Data->m << " non-zero features per example." << endl;
107 +
  }
108 +
    //   for (int i = 0; i<((DoubleVector)X.slot("x")).length();i++) {
109 +
//     Rprintf("val: %f \n",Data->val[i]);
110 +
//   }
111 +
//   for (int i = 0; i<((IntegerVector)X.slot("i")).length();i++) {
112 +
//     Rprintf("col: %d \n",Data->colind[i]);
113 +
//   }
114 +
//   for (int i = 0; i<((IntegerVector)X.slot("p")).length();i++) {
115 +
//     Rprintf("row: %d \n",Data->rowptr[i]);
116 +
//   }
117 +
118 +
119 +
  // Run
120 +
  ssl_train(Data,
121 +
            Options,
122 +
            Weights,
123 +
            Outputs);
124 +
//Clear(Data);
125 +
126 +
  return Rcpp::List::create(Rcpp::Named("Weights") = Rcpp::NumericVector(Weights->vec, Weights->vec+Weights->d),
127 +
                              Rcpp::Named("Outputs") = Rcpp::NumericVector(Outputs->vec, Outputs->vec+Outputs->d)
128 +
                              );
129 +
}

@@ -2,23 +2,23 @@
Loading
2 2
#'
3 3
#' Fit a fast linear SVM classifier for texts, using the
4 4
#' \pkg{LiblineaR} package.
5 -
#' @param x the \link{dfm} on which the model will be fit.  Does not need to
5 +
#' @param x the [dfm] on which the model will be fit.  Does not need to
6 6
#'   contain only the training documents.
7 7
#' @param y vector of training labels associated with each document identified
8 -
#'   in \code{train}.  (These will be converted to factors if not already
8 +
#'   in `train`.  (These will be converted to factors if not already
9 9
#'   factors.)
10 10
#' @param weight weights for different classes for imbalanced training sets,
11 -
#'   passed to \code{wi} in \code{\link[LiblineaR]{LiblineaR}}. \code{"uniform"}
12 -
#'   uses default; \code{"docfreq"} weights by the number of training examples,
13 -
#'   and \code{"termfreq"} by the relative sizes of the training classes in
11 +
#'   passed to `wi` in [LiblineaR::LiblineaR()]. `"uniform"`
12 +
#'   uses default; `"docfreq"` weights by the number of training examples,
13 +
#'   and `"termfreq"` by the relative sizes of the training classes in
14 14
#'   terms of their total lengths in tokens.
15 -
#' @param ... additional arguments passed to \code{\link[LiblineaR]{LiblineaR}}
15 +
#' @param ... additional arguments passed to [LiblineaR::LiblineaR()]
16 16
#' @references
17 17
#' R. E. Fan, K. W. Chang, C. J. Hsieh, X. R. Wang, and C. J. Lin. (2008)
18 18
#' LIBLINEAR: A Library for Large Linear Classification.
19 -
#' \emph{Journal of Machine Learning Research} 9: 1871-1874.
20 -
#' \url{http://www.csie.ntu.edu.tw/~cjlin/liblinear}.
21 -
#' @seealso \code{\link[LiblineaR]{LiblineaR}}
19 +
#' *Journal of Machine Learning Research* 9: 1871-1874.
20 +
#' <http://www.csie.ntu.edu.tw/~cjlin/liblinear>.
21 +
#' @seealso [LiblineaR::LiblineaR()] [predict.textmodel_svm()]
22 22
#' @examples
23 23
#' # use party leaders for govt and opposition classes
24 24
#' quanteda::docvars(data_corpus_irishbudget2010, "govtopp") <-
@@ -93,20 +93,22 @@
Loading
93 93
94 94
#' Prediction from a fitted textmodel_svm object
95 95
#'
96 -
#' \code{predict.textmodel_svm()} implements class predictions from a fitted
96 +
#' `predict.textmodel_svm()` implements class predictions from a fitted
97 97
#' linear SVM model.
98 98
#' @param object a fitted linear SVM textmodel
99 99
#' @param newdata dfm on which prediction should be made
100 100
#' @param type the type of predicted values to be returned; see Value
101 101
#' @param force make newdata's feature set conformant to the model terms
102 102
#' @param ... not used
103 -
#' @return \code{predict.textmodel_svm} returns either a vector of class
104 -
#'   predictions for each row of \code{newdata} (when \code{type = "class"}), or
105 -
#'   a document-by-class matrix of class probabilities (when \code{type =
106 -
#'   "probability"}).
107 -
#' @seealso \code{\link{textmodel_svm}}
103 +
#' @return `predict.textmodel_svm` returns either a vector of class
104 +
#'   predictions for each row of `newdata` (when `type = "class"`), or
105 +
#'   a document-by-class matrix of class probabilities (when `type =
106 +
#'   "probability"`).
107 +
#' @seealso [textmodel_svm()]
108 108
#' @keywords textmodel internal
109 109
#' @importFrom SparseM as.matrix.csr
110 +
#' @importFrom stats predict
111 +
#' @method predict textmodel_svm
110 112
#' @export
111 113
predict.textmodel_svm <- function(object, newdata = NULL,
112 114
                                  type = c("class", "probability"),
@@ -159,7 +161,7 @@
Loading
159 161
}
160 162
161 163
#' summary method for textmodel_svm objects
162 -
#' @param object output from \code{\link{textmodel_svm}}
164 +
#' @param object output from [textmodel_svm()]
163 165
#' @param n how many coefficients to print before truncating
164 166
#' @param ... additional arguments not used
165 167
#' @keywords textmodel internal
@@ -197,8 +199,8 @@
Loading
197 199
198 200
#' convert a dfm into a matrix.csr from SparseM package
199 201
#'
200 -
#' Utility to convert a dfm into a \link[SparseM]{matrix.csr} from the \pkg{SparseM} package.
201 -
#' @param x input \link{dfm}
202 +
#' Utility to convert a dfm into a [matrix.csr][SparseM::matrix.csr] from the \pkg{SparseM} package.
203 +
#' @param x input [dfm]
202 204
#' @importFrom SparseM as.matrix.csr
203 205
#' @importFrom methods new
204 206
#' @method as.matrix.csr dfm

@@ -0,0 +1,1263 @@
Loading
1 +
/*    Copyright 2006 Vikas Sindhwani (vikass@cs.uchicago.edu)
2 +
      SVM-lin: Fast SVM Solvers for Supervised and Semi-supervised Learning
3 +
4 +
      This file is part of SVM-lin.
5 +
6 +
      SVM-lin is free software; you can redistribute it and/or modify
7 +
      it under the terms of the GNU General Public License as published by
8 +
      the Free Software Foundation; either version 2 of the License, or
9 +
      (at your option) any later version.
10 +
11 +
      SVM-lin is distributed in the hope that it will be useful,
12 +
      but WITHOUT ANY WARRANTY; without even the implied warranty of
13 +
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 +
      GNU General Public License for more details.
15 +
16 +
      You should have received a copy of the GNU General Public License
17 +
      along with SVM-lin (see gpl.txt); if not, write to the Free Software
18 +
      Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
19 +
*/
20 +
21 +
#include <Rcpp.h>
22 +
#include <iostream>
23 +
#include <stdio.h>
24 +
#include <stdlib.h>
25 +
#include <algorithm>
26 +
#include <cmath>
27 +
#include <ctype.h>
28 +
#include "ssl.h"
29 +
30 +
#define VERBOSE 1
31 +
#define LOG2(x) 1.4426950408889634*log(x)
32 +
// for compatibility issues, not using log2
33 +
34 +
using namespace std;
35 +
36 +
void ssl_train(struct data *Data,
37 +
	       struct options *Options,
38 +
	       struct vector_double *Weights,
39 +
	       struct vector_double *Outputs)
40 +
{
41 +
  // initialize
42 +
  initialize(Weights,Data->n,0.0);
43 +
  initialize(Outputs,Data->m,0.0);
44 +
  vector_int    *Subset  = new vector_int[1];
45 +
  initialize(Subset,Data->m);
46 +
  // call the right algorithm
47 +
  int optimality = 0;
48 +
  switch(Options->algo)
49 +
    {
50 +
    case -1:
51 +
      if (Options->verbose) { Rcpp::Rcout << "Regularized Least Squares Regression (CGLS)\n" << endl; }
52 +
      optimality=CGLS(Data,Options,Subset,Weights,Outputs);
53 +
      break;
54 +
    case RLS:
55 +
      if (Options->verbose) { Rcpp::Rcout << "Regularized Least Squares Classification (CGLS)\n" << endl; }
56 +
      optimality=CGLS(Data,Options,Subset,Weights,Outputs);
57 +
      break;
58 +
    case SVM:
59 +
      if (Options->verbose) { Rcpp::Rcout << "Modified Finite Newton L2-SVM (L2-SVM-MFN)\n" << endl; }
60 +
      optimality=L2_SVM_MFN(Data,Options,Weights,Outputs,0);
61 +
      break;
62 +
    case TSVM:
63 +
      if (Options->verbose) { Rcpp::Rcout << "Transductive L2-SVM (TSVM)\n" << endl; }
64 +
      optimality=TSVM_MFN(Data,Options,Weights,Outputs);
65 +
      break;
66 +
    case DA_SVM:
67 +
      if (Options->verbose) { Rcpp::Rcout << "Deterministic Annealing Semi-supervised L2-SVM (DAS3VM)\n" << endl; }
68 +
      optimality=DA_S3VM(Data,Options,Weights,Outputs);
69 +
      break;
70 +
    default:
71 +
      ;
72 +
    }
73 +
  if (Options->verbose) { Rcpp::Rcout << "Optimality:" << optimality << endl; }
74 +
  return;
75 +
}
76 +
int CGLS(const struct data *Data,
77 +
	 const struct options *Options,
78 +
	 const struct vector_int *Subset,
79 +
	 struct vector_double *Weights,
80 +
	 struct vector_double *Outputs)
81 +
{
82 +
  if(VERBOSE_CGLS)
83 +
    if (Options->verbose) { Rcpp::Rcout << "CGLS starting..." << endl; }
84 +
  /* Disassemble the structures */
85 +
  timer tictoc;
86 +
  tictoc.restart();
87 +
  int active = Subset->d;
88 +
  int *J = Subset->vec;
89 +
  double *val = Data->val;
90 +
  int *row = Data->rowptr;
91 +
  int *col = Data->colind;
92 +
  double *Y = Data->Y;
93 +
  double *C = Data->C;
94 +
  int n  = Data->n;
95 +
  double lambda = Options->lambda;
96 +
  int cgitermax = Options->cgitermax;
97 +
  double epsilon = Options->epsilon;
98 +
  double *beta = Weights->vec;
99 +
  double *o  = Outputs->vec;
100 +
  // initialize z
101 +
  double *z = new double[active];
102 +
  double *q = new double[active];
103 +
  int ii=0;
104 +
  for(int i = active ; i-- ;){
105 +
    ii=J[i];
106 +
    z[i]  = C[ii]*(Y[ii] - o[ii]);
107 +
  }
108 +
  double *r = new double[n];
109 +
  for(int i = n ; i-- ;)
110 +
    r[i] = 0.0;
111 +
  for(int j=0; j < active; j++)
112 +
    {
113 +
      ii=J[j];
114 +
      for(int i=row[ii]; i < row[ii+1]; i++)
115 +
	r[col[i]]+=val[i]*z[j];
116 +
    }
117 +
  double *p = new double[n];
118 +
  double omega1 = 0.0;
119 +
  for(int i = n ; i-- ;)
120 +
    {
121 +
      r[i] -= lambda*beta[i];
122 +
      p[i] = r[i];
123 +
      omega1 += r[i]*r[i];
124 +
    }
125 +
  double omega_p = omega1;
126 +
  double omega_q = 0.0;
127 +
  double inv_omega2 = 1/omega1;
128 +
  double scale = 0.0;
129 +
  double omega_z=0.0;
130 +
  double gamma = 0.0;
131 +
  int cgiter = 0;
132 +
  int optimality = 0;
133 +
  double epsilon2 = epsilon*epsilon;
134 +
  // iterate
135 +
  while(cgiter < cgitermax)
136 +
    {
137 +
      cgiter++;
138 +
      omega_q=0.0;
139 +
      double t=0.0;
140 +
      int i,j;
141 +
      // #pragma omp parallel for private(i,j)
142 +
      for(i=0; i < active; i++)
143 +
	{
144 +
	  ii=J[i];
145 +
	  t=0.0;
146 +
	  for(j=row[ii]; j < row[ii+1]; j++)
147 +
	    t+=val[j]*p[col[j]];
148 +
	  q[i]=t;
149 +
	  omega_q += C[ii]*t*t;
150 +
	}
151 +
      gamma = omega1/(lambda*omega_p + omega_q);
152 +
      inv_omega2 = 1/omega1;
153 +
      for(int i = n ; i-- ;)
154 +
	{
155 +
	  r[i] = 0.0;
156 +
	  beta[i] += gamma*p[i];
157 +
	}
158 +
      omega_z=0.0;
159 +
      for(int i = active ; i-- ;)
160 +
	{
161 +
	  ii=J[i];
162 +
	  o[ii] += gamma*q[i];
163 +
	  z[i] -= gamma*C[ii]*q[i];
164 +
	  omega_z+=z[i]*z[i];
165 +
	}
166 +
      for(int j=0; j < active; j++)
167 +
	{
168 +
	  ii=J[j];
169 +
	  t=z[j];
170 +
	  for(int i=row[ii]; i < row[ii+1]; i++)
171 +
	    r[col[i]]+=val[i]*t;
172 +
	}
173 +
      omega1 = 0.0;
174 +
      for(int i = n ; i-- ;)
175 +
	{
176 +
	  r[i] -= lambda*beta[i];
177 +
	  omega1 += r[i]*r[i];
178 +
	}
179 +
      if(VERBOSE_CGLS)
180 +
        if (Options->verbose) { Rcpp::Rcout << "..." << cgiter << " ( " << omega1 << " )" ; }
181 +
      if(omega1 < epsilon2*omega_z)
182 +
	{
183 +
	  optimality=1;
184 +
	  break;
185 +
	}
186 +
      omega_p=0.0;
187 +
      scale=omega1*inv_omega2;
188 +
      for(int i = n ; i-- ;)
189 +
	{
190 +
	  p[i] = r[i] + p[i]*scale;
191 +
	  omega_p += p[i]*p[i];
192 +
	}
193 +
    }
194 +
  if(VERBOSE_CGLS)
195 +
    if (Options->verbose) { Rcpp::Rcout << "...Done." << endl; }
196 +
  tictoc.stop();
197 +
    if (Options->verbose) { Rcpp::Rcout << "CGLS converged in " << cgiter << " iteration(s) and " << tictoc.time() << " seconds." << endl; }
198 +
  delete[] z;
199 +
  delete[] q;
200 +
  delete[] r;
201 +
  delete[] p;
202 +
  return optimality;
203 +
}
204 +
int L2_SVM_MFN(const struct data *Data,
205 +
	       struct options *Options,
206 +
	       struct vector_double *Weights,
207 +
	       struct vector_double *Outputs,
208 +
	       int ini)
209 +
{
210 +
  /* Disassemble the structures */
211 +
  timer tictoc;
212 +
  tictoc.restart();
213 +
  double *val = Data->val;
214 +
  int *row = Data->rowptr;
215 +
  int *col = Data->colind;
216 +
  double *Y = Data->Y;
217 +
  double *C = Data->C;
218 +
  int n  = Data->n;
219 +
  int m  = Data->m;
220 +
  double lambda = Options->lambda;
221 +
  double epsilon;
222 +
  double *w = Weights->vec;
223 +
  double *o = Outputs->vec;
224 +
  double F_old = 0.0;
225 +
  double F = 0.0;
226 +
  double diff=0.0;
227 +
  vector_int *ActiveSubset = new vector_int[1];
228 +
  ActiveSubset->vec = new int[m];
229 +
  ActiveSubset->d = m;
230 +
  // initialize
231 +
  if(ini==0) {
232 +
    epsilon=BIG_EPSILON;
233 +
    Options->cgitermax=SMALL_CGITERMAX;
234 +
    Options->epsilon=BIG_EPSILON;
235 +
  }
236 +
  else {epsilon = Options->epsilon;}
237 +
  for(int i=0;i<n;i++) F+=w[i]*w[i];
238 +
  F=0.5*lambda*F;
239 +
  int active=0;
240 +
  int inactive=m-1; // l-1
241 +
  for(int i=0; i<m ; i++)
242 +
    {
243 +
      diff=1-Y[i]*o[i];
244 +
      if(diff>0)
245 +
	{
246 +
	  ActiveSubset->vec[active]=i;
247 +
	  active++;
248 +
	  F+=0.5*C[i]*diff*diff;
249 +
	}
250 +
      else
251 +
	{
252 +
	  ActiveSubset->vec[inactive]=i;
253 +
	  inactive--;
254 +
	}
255 +
    }
256 +
  ActiveSubset->d=active;
257 +
  int iter=0;
258 +
  int opt=0;
259 +
  int opt2=0;
260 +
  vector_double *Weights_bar = new vector_double[1];
261 +
  vector_double *Outputs_bar = new vector_double[1];
262 +
  double *w_bar = new double[n];
263 +
  double *o_bar = new double[m];
264 +
  Weights_bar->vec=w_bar;
265 +
  Outputs_bar->vec=o_bar;
266 +
  Weights_bar->d=n;
267 +
  Outputs_bar->d=m;
268 +
  double delta=0.0;
269 +
  double t=0.0;
270 +
  int ii = 0;
271 +
  while(iter<MFNITERMAX)
272 +
    {
273 +
      iter++;
274 +
    if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN Iteration# " << iter << " (" << active << " active examples, " << " objective_value = " << F << ")" << endl; }
275 +
      for(int i=n; i-- ;)
276 +
	w_bar[i]=w[i];
277 +
      for(int i=m; i-- ;)
278 +
	o_bar[i]=o[i];
279 +
280 +
      opt=CGLS(Data,Options,ActiveSubset,Weights_bar,Outputs_bar);
281 +
      for(int i=active; i < m; i++)
282 +
	{
283 +
	  ii=ActiveSubset->vec[i];
284 +
	  t=0.0;
285 +
	  for(int j=row[ii]; j < row[ii+1]; j++)
286 +
	    t+=val[j]*w_bar[col[j]];
287 +
	  o_bar[ii]=t;
288 +
	}
289 +
      if(ini==0) {Options->cgitermax=CGITERMAX; ini=1;};
290 +
      opt2=1;
291 +
      for(int i=0;i<m;i++)
292 +
	{
293 +
	  ii=ActiveSubset->vec[i];
294 +
	  if(i<active)
295 +
	    opt2=(opt2 && (Y[ii]*o_bar[ii]<=1+epsilon));
296 +
	  else
297 +
	    opt2=(opt2 && (Y[ii]*o_bar[ii]>=1-epsilon));
298 +
	  if(opt2==0) break;
299 +
	}
300 +
      if(opt && opt2) // l
301 +
	{
302 +
	  if(epsilon==BIG_EPSILON)
303 +
	    {
304 +
	      epsilon=EPSILON;
305 +
	      Options->epsilon=EPSILON;
306 +
	      if (Options->verbose) { Rcpp::Rcout << "  epsilon = " << BIG_EPSILON << " case converged (speedup heuristic 2). Continuing with epsilon=" <<  EPSILON << endl; }
307 +
	      continue;
308 +
	    }
309 +
	  else
310 +
	    {
311 +
	      for(int i=n; i-- ;)
312 +
		w[i]=w_bar[i];
313 +
	      for(int i=m; i-- ;)
314 +
		o[i]=o_bar[i];
315 +
	       delete[] ActiveSubset->vec;
316 +
	       delete[] ActiveSubset;
317 +
	       delete[] o_bar;
318 +
	       delete[] w_bar;
319 +
	       delete[] Weights_bar;
320 +
	       delete[] Outputs_bar;
321 +
	       tictoc.stop();
322 +
	       if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (optimality) in " << iter << " iteration(s) and "<< tictoc.time() << " seconds. \n" << endl; }
323 +
	       return 1;
324 +
	    }
325 +
	}
326 +
      if (Options->verbose) { Rcpp::Rcout << " " ; }
327 +
      delta=line_search(w,w_bar,lambda,o,o_bar,Y,C,n,m);
328 +
      if (Options->verbose) { Rcpp::Rcout << "LINE_SEARCH delta = " << delta << endl; }
329 +
      F_old=F;
330 +
      F=0.0;
331 +
      for(int i=n; i-- ;){
332 +
	w[i]+=delta*(w_bar[i]-w[i]);
333 +
	F+=w[i]*w[i];
334 +
      }
335 +
      F=0.5*lambda*F;
336 +
      active=0;
337 +
      inactive=m-1;
338 +
      for(int i=0; i<m ; i++)
339 +
	{
340 +
	  o[i]+=delta*(o_bar[i]-o[i]);
341 +
	  diff=1-Y[i]*o[i];
342 +
	  if(diff>0)
343 +
	    {
344 +
	      ActiveSubset->vec[active]=i;
345 +
	      active++;
346 +
	      F+=0.5*C[i]*diff*diff;
347 +
	    }
348 +
	  else
349 +
	    {
350 +
	      ActiveSubset->vec[inactive]=i;
351 +
	      inactive--;
352 +
	    }
353 +
	}
354 +
      ActiveSubset->d=active;
355 +
      if(fabs(F-F_old)<RELATIVE_STOP_EPS*fabs(F_old))
356 +
	{
357 +
        if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (rel. criterion) in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl; }
358 +
	  return 2;
359 +
	}
360 +
    }
361 +
  delete[] ActiveSubset->vec;
362 +
  delete[] ActiveSubset;
363 +
  delete[] o_bar;
364 +
  delete[] w_bar;
365 +
  delete[] Weights_bar;
366 +
  delete[] Outputs_bar;
367 +
  tictoc.stop();
368 +
  if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (max iter exceeded) in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl; }
369 +
  return 0;
370 +
}
371 +
372 +
double line_search(double *w,
373 +
                   double *w_bar,
374 +
                   double lambda,
375 +
                   double *o,
376 +
                   double *o_bar,
377 +
                   double *Y,
378 +
                   double *C,
379 +
                   int d, /* data dimensionality -- 'n' */
380 +
                   int l) /* number of examples */
381 +
{
382 +
   double omegaL = 0.0;
383 +
   double omegaR = 0.0;
384 +
   double diff=0.0;
385 +
   for(int i=d; i--; )
386 +
       {
387 +
         diff=w_bar[i]-w[i];
388 +
         omegaL+=w[i]*diff;
389 +
         omegaR+=w_bar[i]*diff;
390 +
   }
391 +
   omegaL=lambda*omegaL;
392 +
   omegaR=lambda*omegaR;
393 +
   double L=0.0;
394 +
   double R=0.0;
395 +
   int ii=0;
396 +
   for(int i=0;i<l;i++)
397 +
       {
398 +
         if(Y[i]*o[i]<1)
399 +
	   {
400 +
	     diff=C[i]*(o_bar[i]-o[i]);
401 +
	     L+=(o[i]-Y[i])*diff;
402 +
	     R+=(o_bar[i]-Y[i])*diff;
403 +
	   }
404 +
       }
405 +
   L+=omegaL;
406 +
   R+=omegaR;
407 +
   Delta* deltas=new Delta[l];
408 +
   int p=0;
409 +
   for(int i=0;i<l;i++)
410 +
     {
411 +
       diff=Y[i]*(o_bar[i]-o[i]);
412 +
413 +
       if(Y[i]*o[i]<1)
414 +
	 {
415 +
	   if(diff>0)
416 +
	     {
417 +
	       deltas[p].delta=(1-Y[i]*o[i])/diff;
418 +
	       deltas[p].index=i;
419 +
	       deltas[p].s=-1;
420 +
	       p++;
421 +
	     }
422 +
	 }
423 +
       else
424 +
	 {
425 +
	   if(diff<0)
426 +
	     {
427 +
	       deltas[p].delta=(1-Y[i]*o[i])/diff;
428 +
	       deltas[p].index=i;
429 +
	       deltas[p].s=1;
430 +
	       p++;
431 +
	     }
432 +
	 }
433 +
     }
434 +
   sort(deltas,deltas+p);
435 +
   double delta_prime=0.0;
436 +
   for(int i=0;i<p;i++)
437 +
     {
438 +
       delta_prime = L + deltas[i].delta*(R-L);
439 +
       if(delta_prime>=0)
440 +
	 break;
441 +
       ii=deltas[i].index;
442 +
       diff=(deltas[i].s)*C[ii]*(o_bar[ii]-o[ii]);
443 +
       L+=diff*(o[ii]-Y[ii]);
444 +
       R+=diff*(o_bar[ii]-Y[ii]);
445 +
     }
446 +
   delete [] deltas;
447 +
   return (-L/(R-L));
448 +
}
449 +
int TSVM_MFN(const struct data *Data,
450 +
	      struct options *Options,
451 +
	      struct vector_double *Weights,
452 +
	      struct vector_double *Outputs)
453 +
{
454 +
  /* Setup labeled-only examples and train L2_SVM_MFN */
455 +
  timer tictoc;
456 +
  tictoc.restart();
457 +
  struct data *Data_Labeled = new data[1];
458 +
  struct vector_double *Outputs_Labeled = new vector_double[1];
459 +
  initialize(Outputs_Labeled,Data->l,0.0);
460 +
  if (Options->verbose) { Rcpp::Rcout << "Initializing weights, unknown labels" << endl; }
461 +
  GetLabeledData(Data_Labeled,Data); /* gets labeled data and sets C=1/l */
462 +
  L2_SVM_MFN(Data_Labeled, Options, Weights,Outputs_Labeled,0);
463 +
  //Clear(Data_Labeled);
464 +
  /* Use this weight vector to classify R*u unlabeled examples as
465 +
       positive*/
466 +
  int p=0,q=0;
467 +
  double t=0.0;
468 +
  int *JU = new int[Data->u];
469 +
  double *ou = new double[Data->u];
470 +
  double lambda_0 = TSVM_LAMBDA_SMALL;
471 +
  for(int i=0;i<Data->m;i++)
472 +
    {
473 +
      if(Data->Y[i]==0.0)
474 +
	{
475 +
	  t=0.0;
476 +
	  for(int j=Data->rowptr[i]; j < Data->rowptr[i+1]; j++)
477 +
	    t+=Data->val[j]*Weights->vec[Data->colind[j]];
478 +
	  Outputs->vec[i]=t;
479 +
	  Data->C[i]=lambda_0*1.0/Data->u;
480 +
	  JU[q]=i;
481 +
	  ou[q]=t;
482 +
	  q++;
483 +
	}
484 +
      else
485 +
	{
486 +
	  Outputs->vec[i]=Outputs_Labeled->vec[p];
487 +
	  Data->C[i]=1.0/Data->l;
488 +
	  p++;
489 +
	}
490 +
    }
491 +
  nth_element(ou,ou+int((1-Options->R)*Data->u-1),ou+Data->u);
492 +
  double thresh=*(ou+int((1-Options->R)*Data->u)-1);
493 +
  delete [] ou;
494 +
  for(int i=0;i<Data->u;i++)
495 +
    {
496 +
      if(Outputs->vec[JU[i]]>thresh)
497 +
	Data->Y[JU[i]]=1.0;
498 +
      else
499 +
	Data->Y[JU[i]]=-1.0;
500 +
    }
501 +
  for(int i=0;i<Data->n;i++)
502 +
    Weights->vec[i]=0.0;
503 +
  for(int i=0;i<Data->m;i++)
504 +
    Outputs->vec[i]=0.0;
505 +
  L2_SVM_MFN(Data,Options,Weights,Outputs,0);
506 +
  int num_switches=0;
507 +
  int s=0;
508 +
  int last_round=0;
509 +
  while(lambda_0 <= Options->lambda_u)
510 +
    {
511 +
      int iter2=0;
512 +
      while(1){
513 +
	s=switch_labels(Data->Y,Outputs->vec,JU,Data->u,Options->S);
514 +
	if(s==0) break;
515 +
	iter2++;
516 +
	if (Options->verbose) { Rcpp::Rcout << "****** lambda_0 = " << lambda_0 << " iteration = " << iter2 << " ************************************"  << endl; }
517 +
	if (Options->verbose) { Rcpp::Rcout << "Optimizing unknown labels. switched " << s << " labels." << endl; }
518 +
	num_switches+=s;
519 +
	if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
520 +
	L2_SVM_MFN(Data,Options,Weights,Outputs,1);
521 +
      }
522 +
      if(last_round==1) break;
523 +
      lambda_0=TSVM_ANNEALING_RATE*lambda_0;
524 +
      if(lambda_0 >= Options->lambda_u) {lambda_0 = Options->lambda_u; last_round=1;}
525 +
      for(int i=0;i<Data->u;i++)
526 +
	Data->C[JU[i]]=lambda_0*1.0/Data->u;
527 +
      if (Options->verbose) { Rcpp::Rcout << endl << "****** lambda0 increased to " << lambda_0*100/Options->lambda_u << "%" << " of lambda_u = " << Options->lambda_u << " ************************" << endl; }
528 +
      if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
529 +
      L2_SVM_MFN(Data,Options,Weights,Outputs,1);
530 +
    }
531 +
  if (Options->verbose) { Rcpp::Rcout <<  "Total Number of Switches = " << num_switches << endl; }
532 +
  /* reset labels */
533 +
  for(int i=0;i<Data->u;i++) Data->Y[JU[i]] = 0.0;
534 +
  double F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
535 +
  if (Options->verbose) { Rcpp::Rcout << "Objective Value = " << F << endl; }
536 +
  delete [] JU;
537 +
  tictoc.stop();
538 +
  if (Options->verbose) { Rcpp::Rcout << "Transductive L2_SVM_MFN took " << tictoc.time() << " seconds." << endl; }
539 +
  return num_switches;
540 +
}
541 +
int switch_labels(double* Y, double* o, int* JU, int u, int S)
542 +
{
543 +
  int npos=0;
544 +
  int nneg=0;
545 +
  for(int i=0;i<u;i++)
546 +
    {
547 +
      if((Y[JU[i]]>0) && (o[JU[i]]<1.0)) npos++;
548 +
      if((Y[JU[i]]<0) && (-o[JU[i]]<1.0)) nneg++;
549 +
    }
550 +
  Delta* positive=new Delta[npos];
551 +
  Delta* negative=new Delta[nneg];
552 +
  int p=0;
553 +
  int n=0;
554 +
  int ii=0;
555 +
  for(int i=0;i<u;i++)
556 +
    {
557 +
      ii=JU[i];
558 +
      if((Y[ii]>0.0) && (o[ii]<1.0)) {
559 +
	positive[p].delta=o[ii];
560 +
	positive[p].index=ii;
561 +
	positive[p].s=0;
562 +
	p++;};
563 +
      if((Y[ii]<0.0) && (-o[ii]<1.0))
564 +
	{
565 +
	  negative[n].delta=-o[ii];
566 +
	  negative[n].index=ii;
567 +
	  negative[n].s=0;
568 +
	  n++;};
569 +
    }
570 +
  sort(positive,positive+npos);
571 +
  sort(negative,negative+nneg);
572 +
  int s=-1;
573 +
  //  cout << "Switching Labels: ";
574 +
  while(1)
575 +
    {
576 +
      s++;
577 +
      if((s>=S) || (positive[s].delta>=-negative[s].delta) || (s>=npos) || (s>=nneg))
578 +
	break;
579 +
      Y[positive[s].index]=-1.0;
580 +
      Y[negative[s].index]= 1.0;
581 +
      //  cout << positive[s].index << " " << negative[s].index << "...";
582 +
    }
583 +
  // cout << "...switched " << s  << " labels." << endl;
584 +
  delete [] positive;
585 +
  delete [] negative;
586 +
  return s;
587 +
}
588 +
589 +
int DA_S3VM(struct data *Data,
590 +
	   struct options *Options,
591 +
	   struct vector_double *Weights,
592 +
	   struct vector_double *Outputs)
593 +
{
594 +
  timer tictoc;
595 +
  tictoc.restart();
596 +
  double T = DA_INIT_TEMP*Options->lambda_u;
597 +
  int iter1 = 0, iter2 =0;
598 +
  double *p = new double[Data->u];
599 +
  double *q = new double[Data->u];
600 +
  double *g = new double[Data->u];
601 +
  double F,F_min;
602 +
  double *w_min = new double[Data->n];
603 +
  double *o_min = new double[Data->m];
604 +
  double *w = Weights->vec;
605 +
  double *o = Outputs->vec;
606 +
  double kl_divergence = 1.0;
607 +
  /*initialize */
608 +
  if (Options->verbose) { Rcpp::Rcout << "Initializing weights, p" << endl; }
609 +
  for(int i=0;i<Data->u; i++)
610 +
    p[i] = Options->R;
611 +
  /* record which examples are unlabeled */
612 +
  int *JU = new int[Data->u];
613 +
  int j=0;
614 +
  for(int i=0;i<Data->m;i++)
615 +
    {
616 +
      if(Data->Y[i]==0.0)
617 +
	{JU[j]=i;j++;}
618 +
    }
619 +
  double H = entropy(p,Data->u);
620 +
  optimize_w(Data,p,Options,Weights,Outputs,0);
621 +
  F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
622 +
  F_min = F;
623 +
  for(int i=0;i<Weights->d;i++)
624 +
    w_min[i]=w[i];
625 +
  for(int i=0;i<Outputs->d;i++)
626 +
    o_min[i]=o[i];
627 +
  while((iter1 < DA_OUTER_ITERMAX) && (H > Options->epsilon))
628 +
    {
629 +
      iter1++;
630 +
      iter2=0;
631 +
      kl_divergence=1.0;
632 +
      while((iter2 < DA_INNER_ITERMAX) && (kl_divergence > Options->epsilon))
633 +
	{
634 +
	  iter2++;
635 +
	  for(int i=0;i<Data->u;i++)
636 +
	    {
637 +
	      q[i]=p[i];
638 +
	      g[i] = Options->lambda_u*((o[JU[i]] > 1 ? 0 : (1 - o[JU[i]])*(1 - o[JU[i]])) - (o[JU[i]]< -1 ? 0 : (1 + o[JU[i]])*(1 + o[JU[i]])));
639 +
	    }
640 +
	  if (Options->verbose) { Rcpp::Rcout << "Optimizing p. " ; }
641 +
	  optimize_p(g,Data->u,T,Options->R,p);
642 +
	  kl_divergence=KL(p,q,Data->u);
643 +
	  if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
644 +
	  optimize_w(Data,p,Options,Weights,Outputs,1);
645 +
	  F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
646 +
	  if(F < F_min)
647 +
	    {
648 +
	      F_min = F;
649 +
	      for(int i=0;i<Weights->d;i++)
650 +
		w_min[i]=w[i];
651 +
	      for(int i=0;i<Outputs->d;i++)
652 +
		o_min[i]=o[i];
653 +
	    }
654 +
	  if (Options->verbose) { Rcpp::Rcout << "***** outer_iter = " << iter1 << "  T = " << T << "  inner_iter = " << iter2 << "  kl = "<< kl_divergence <<"  cost = "<<F<<" *****"<<endl; }
655 +
	}
656 +
      H = entropy(p,Data->u);
657 +
      if (Options->verbose) { Rcpp::Rcout << "***** Finished outer_iter = " << iter1 << "T = " << T << " Entropy = " << H <<endl; }
658 +
      T = T/DA_ANNEALING_RATE;
659 +
    }
660 +
  for(int i=0;i<Weights->d;i++)
661 +
    w[i]=w_min[i];
662 +
  for(int i=0;i<Outputs->d;i++)
663 +
    o[i]=o_min[i];
664 +
  /* may want to reset the original Y */
665 +
  delete [] p;
666 +
  delete [] q;
667 +
  delete [] g;
668 +
  delete [] JU;
669 +
  delete [] w_min;
670 +
  delete [] o_min;
671 +
  tictoc.stop();
672 +
  if (Options->verbose) { Rcpp::Rcout << endl << "(min) Objective Value = " << F_min << endl; }
673 +
  if (Options->verbose) { Rcpp::Rcout << "DA-SVM took " << tictoc.time() << " seconds." << endl; }
674 +
  return 1;
675 +
}
676 +
int optimize_w(const struct data *Data,
677 +
	       const double *p,
678 +
	       struct options *Options,
679 +
	       struct vector_double *Weights,
680 +
	       struct vector_double *Outputs,
681 +
	       int ini)
682 +
{
683 +
  timer tictoc;
684 +
  tictoc.restart();
685 +
  double *val = Data->val;
686 +
  int *row = Data->rowptr;
687 +
  int *col = Data->colind;
688 +
  int n  = Data->n;
689 +
  int m  = Data->m;
690 +
  int u  = Data->u;
691 +
  double lambda = Options->lambda;
692 +
  double epsilon;
693 +
  double *w = Weights->vec;
694 +
  double *o = new double[m+u];
695 +
  double *Y = new double[m+u];
696 +
  double *C = new double[m+u];
697 +
  int *labeled_indices = new int[m];
698 +
  double F_old = 0.0;
699 +
  double F = 0.0;
700 +
  double diff=0.0;
701 +
  double lambda_u_by_u = Options->lambda_u/u;
702 +
  vector_int *ActiveSubset = new vector_int[1];
703 +
  ActiveSubset->vec = new int[m];
704 +
  ActiveSubset->d = m;
705 +
  // initialize
706 +
  if(ini==0)
707 +
    {
708 +
      epsilon=BIG_EPSILON;
709 +
      Options->cgitermax=SMALL_CGITERMAX;
710 +
      Options->epsilon=BIG_EPSILON;}
711 +
  else {epsilon = Options->epsilon;}
712 +
713 +
  for(int i=0;i<n;i++) F+=w[i]*w[i];
714 +
  F=lambda*F;
715 +
  int active=0;
716 +
  int inactive=m-1; // l-1
717 +
  int j = 0;
718 +
  double temp1;
719 +
  double temp2;
720 +
  for(int i=0; i<m ; i++)
721 +
    {
722 +
      o[i]=Outputs->vec[i];
723 +
      if(Data->Y[i]==0.0)
724 +
	{
725 +
	  labeled_indices[i]=0;
726 +
	  o[m+j]=o[i];
727 +
	  Y[i]=1;
728 +
	  Y[m+j]=-1;
729 +
	  C[i]=lambda_u_by_u*p[j];
730 +
	  C[m+j]=lambda_u_by_u*(1-p[j]);
731 +
	  ActiveSubset->vec[active]=i;
732 +
	  active++;
733 +
	  diff = 1 - fabs(o[i]);
734 +
	  if(diff>0)
735 +
	    {
736 +
	      Data->Y[i] = 2*p[j]-1;
737 +
	      Data->C[i] = lambda_u_by_u;
738 +
	      temp1 = (1 - o[i]);
739 +
	      temp2 = (1 + o[i]);
740 +
	      F+=lambda_u_by_u*(p[j]*temp1*temp1 + (1-p[j])*temp2*temp2);
741 +
	    }
742 +
	  else
743 +
	    {
744 +
	      if(o[i]>0)
745 +
		{
746 +
		  Data->Y[i] = -1.0;
747 +
		  Data->C[i] = C[m+j];
748 +
		}
749 +
	      else
750 +
		{
751 +
		  Data->Y[i] = 1.0;
752 +
		  Data->C[i] = C[i];
753 +
		}
754 +
	      temp1 = (1-Data->Y[i]*o[i]);
755 +
	      F+= Data->C[i]*temp1*temp1;
756 +
	    }
757 +
	  j++;
758 +
	}
759 +
      else
760 +
	{
761 +
	  labeled_indices[i]=1;
762 +
	  Y[i]=Data->Y[i];
763 +
	  C[i]=1.0/Data->l;
764 +
	  Data->C[i]=1.0/Data->l;
765 +
	  diff=1-Data->Y[i]*o[i];
766 +
	  if(diff>0)
767 +
	    {
768 +
	      ActiveSubset->vec[active]=i;
769 +
	      active++;
770 +
	      F+=Data->C[i]*diff*diff;
771 +
	    }
772 +
	  else
773 +
	    {
774 +
	      ActiveSubset->vec[inactive]=i;
775 +
	      inactive--;
776 +
	    }
777 +
	}
778 +
    }
779 +
  F=0.5*F;
780 +
  ActiveSubset->d=active;
781 +
  int iter=0;
782 +
  int opt=0;
783 +
  int opt2=0;
784 +
  vector_double *Weights_bar = new vector_double[1];
785 +
  vector_double *Outputs_bar = new vector_double[1];
786 +
  double *w_bar = new double[n];
787 +
  double *o_bar = new double[m+u];
788 +
  Weights_bar->vec=w_bar;
789 +
  Outputs_bar->vec=o_bar;
790 +
  Weights_bar->d=n;
791 +
  Outputs_bar->d=m; /* read only the top m ; bottom u will be copies */
792 +
  double delta=0.0;
793 +
  double t=0.0;
794 +
  int ii = 0;
795 +
  while(iter<MFNITERMAX)
796 +
    {
797 +
      iter++;
798 +
      if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN Iteration# " << iter << " (" << active << " active examples, " << " objective_value = " << F << ")" << endl; }
799 +
      for(int i=n; i-- ;)
800 +
	w_bar[i]=w[i];
801 +
802 +
      for(int i=m+u; i-- ;)
803 +
	o_bar[i]=o[i];
804 +
      if (Options->verbose) { Rcpp::Rcout << "_" ; }
805 +
      opt=CGLS(Data,Options,ActiveSubset,Weights_bar,Outputs_bar);
806 +
      for(int i=active; i < m; i++)
807 +
	{
808 +
	  ii=ActiveSubset->vec[i];
809 +
	  t=0.0;
810 +
	  for(int j=row[ii]; j < row[ii+1]; j++)
811 +
	    t+=val[j]*w_bar[col[j]];
812 +
	  o_bar[ii]=t;
813 +
	}
814 +
      // make o_bar consistent in the bottom half
815 +
      int j=0;
816 +
      for(int i=0; i<m;i++)
817 +
	{
818 +
	  if(labeled_indices[i]==0)
819 +
	    {o_bar[m+j]=o_bar[i]; j++;};
820 +
	}
821 +
      if(ini==0) {Options->cgitermax=CGITERMAX; ini=1;};
822 +
      opt2=1;
823 +
      for(int i=0; i < m ;i++)
824 +
	{
825 +
	  ii=ActiveSubset->vec[i];
826 +
	  if(i<active)
827 +
	    {
828 +
	      if(labeled_indices[ii]==1)
829 +
		opt2=(opt2 && (Data->Y[ii]*o_bar[ii]<=1+epsilon));
830 +
	      else
831 +
		{
832 +
		  if(fabs(o[ii])<1)
833 +
		    opt2=(opt2 && (fabs(o_bar[ii])<=1+epsilon));
834 +
		  else
835 +
		    opt2=(opt2 && (fabs(o_bar[ii])>=1-epsilon));
836 +
		}
837 +
	    }
838 +
	  else
839 +
	    opt2=(opt2 && (Data->Y[ii]*o_bar[ii]>=1-epsilon));
840 +
	  if(opt2==0) break;
841 +
	}
842 +
      if(opt && opt2) // l
843 +
	{
844 +
	  if(epsilon==BIG_EPSILON)
845 +
	    {
846 +
	      epsilon=EPSILON;
847 +
	      Options->epsilon=EPSILON;
848 +
	      if (Options->verbose) { Rcpp::Rcout << "  epsilon = " << BIG_EPSILON << " case converged (speedup heuristic 2). Continuing with epsilon=" <<  EPSILON << endl; }
849 +
	      continue;
850 +
	    }
851 +
	  else
852 +
	    {
853 +
	      for(int i=n; i-- ;)
854 +
		w[i]=w_bar[i];
855 +
	      for(int i=m; i-- ;)
856 +
		Outputs->vec[i]=o_bar[i];
857 +
	      for(int i=m; i-- ;)
858 +
		{
859 +
		  if(labeled_indices[i]==0)
860 +
		    Data->Y[i]=0.0;
861 +
		}
862 +
	       delete[] ActiveSubset->vec;
863 +
	       delete[] ActiveSubset;
864 +
	       delete[] o_bar;
865 +
	       delete[] w_bar;
866 +
	       delete[] o;
867 +
	       delete[] Weights_bar;
868 +
	       delete[] Outputs_bar;
869 +
	       delete[] Y;
870 +
	       delete[] C;
871 +
	       delete[] labeled_indices;
872 +
	       tictoc.stop();
873 +
	       if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged in " << iter << " iteration(s) and "<< tictoc.time() << " seconds. \n" << endl; }
874 +
	       return 1;
875 +
	    }
876 +
	}
877 +
878 +
      delta=line_search(w,w_bar,lambda,o,o_bar,Y,C,n,m+u);
879 +
      if (Options->verbose) { Rcpp::Rcout << " LINE_SEARCH delta = " << delta << endl; }
880 +
      F_old=F;
881 +
      F=0.0;
882 +
      for(int i=0;i<n;i++) {w[i]+=delta*(w_bar[i]-w[i]);  F+=w[i]*w[i];}
883 +
      F=lambda*F;
884 +
      j=0;
885 +
      active=0;
886 +
      inactive=m-1;
887 +
      for(int i=0; i<m ; i++)
888 +
	{
889 +
	   o[i]+=delta*(o_bar[i]-o[i]);
890 +
	   if(labeled_indices[i]==0)
891 +
	    {
892 +
	      o[m+j]=o[i];
893 +
	      ActiveSubset->vec[active]=i;
894 +
	      active++;
895 +
	      diff = 1 - fabs(o[i]);
896 +
	      if(diff>0)
897 +
		{
898 +
		  Data->Y[i] = 2*p[j]-1;
899 +
		  Data->C[i] = lambda_u_by_u;
900 +
		  temp1 = (1 - o[i]);
901 +
		  temp2 = (1 + o[i]);
902 +
		  F+=lambda_u_by_u*(p[j]*temp1*temp1 + (1-p[j])*temp2*temp2);
903 +
		}
904 +
	      else
905 +
		{
906 +
		  if(o[i]>0)
907 +
		    {
908 +
		      Data->Y[i] = -1;
909 +
		      Data->C[i] = C[m+j];
910 +
		    }
911 +
		  else
912 +
		    {
913 +
		      Data->Y[i] = +1;
914 +
		      Data->C[i] = C[i];
915 +
		    }
916 +
		  temp1=(1-Data->Y[i]*o[i]);
917 +
		  F+= Data->C[i]*temp1*temp1;
918 +
		}
919 +
	      j++;
920 +
	    }
921 +
	  else
922 +
	    {
923 +
	      diff=1-Data->Y[i]*o[i];
924 +
	      if(diff>0)
925 +
		{
926 +
		  ActiveSubset->vec[active]=i;
927 +
		  active++;
928 +
		  F+=Data->C[i]*diff*diff;
929 +
		}
930 +
	      else
931 +
		{
932 +
		  ActiveSubset->vec[inactive]=i;
933 +
		  inactive--;
934 +
		}
935 +
	    }
936 +
	}
937 +
      F=0.5*F;
938 +
      ActiveSubset->d=active;
939 +
      if(fabs(F-F_old)<EPSILON)
940 +
	break;
941 +
    }
942 +
  for(int i=m; i-- ;)
943 +
    {
944 +
      Outputs->vec[i]=o[i];
945 +
      if(labeled_indices[i]==0)
946 +
	Data->Y[i]=0.0;
947 +
    }
948 +
  delete[] ActiveSubset->vec;
949 +
  delete[] labeled_indices;
950 +
  delete[] ActiveSubset;
951 +
  delete[] o_bar;
952 +
  delete[] w_bar;
953 +
  delete[] o;
954 +
  delete[] Weights_bar;
955 +
  delete[] Outputs_bar;
956 +
  delete[] Y;
957 +
  delete[] C;
958 +
  tictoc.stop();
959 +
  if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl;}
960 +
  return 0;
961 +
}
962 +
void optimize_p(const double* g, int u, double T, double r, double* p)
963 +
{
964 +
  int iter=0;
965 +
  double epsilon=1e-10;
966 +
  int maxiter=500;
967 +
  double nu_minus=g[0];
968 +
  double nu_plus=g[0];
969 +
  for(int i=0;i<u;i++)
970 +
    {
971 +
      if(g[i]<nu_minus) nu_minus=g[i];
972 +
      if(g[i]>nu_plus) nu_plus=g[i];
973 +
    };
974 +
975 +
  double b=T*log((1-r)/r);
976 +
  nu_minus-=b;
977 +
  nu_plus-=b;
978 +
  double nu=(nu_plus+nu_minus)/2;
979 +
  double Bnu=0.0;
980 +
  double BnuPrime=0.0;
981 +
  double s=0.0;
982 +
  double tmp=0.0;
983 +
  for(int i=0;i<u;i++)
984 +
    {
985 +
      s=exp((g[i]-nu)/T);
986 +
      if(!(std::isinf(s)))
987 +
	{
988 +
	  tmp=1.0/(1.0+s);
989 +
	  Bnu+=tmp;
990 +
	  BnuPrime+=s*tmp*tmp;
991 +
	}
992 +
    }
993 +
  Bnu=Bnu/u;
994 +
  Bnu-=r;
995 +
  BnuPrime=BnuPrime/(T*u);
996 +
  double nuHat=0.0;
997 +
  while((fabs(Bnu)>epsilon) && (iter < maxiter))
998 +
    {
999 +
      iter++;
1000 +
      if(fabs(BnuPrime)>0.0)
1001 +
        nuHat=nu-Bnu/BnuPrime;
1002 +
      if((fabs(BnuPrime) > 0.0) | (nuHat>nu_plus)  | (nuHat < nu_minus))
1003 +
        nu=(nu_minus+nu_plus)/2.0;
1004 +
      else
1005 +
        nu=nuHat;
1006 +
      Bnu=0.0;
1007 +
      BnuPrime=0.0;
1008 +
      for(int i=0;i<u;i++)
1009 +
	{
1010 +
	  s=exp((g[i]-nu)/T);
1011 +
	  if(!(std::isinf(s)))
1012 +
	    {
1013 +
	      tmp=1.0/(1.0+s);
1014 +
	      Bnu+=tmp;
1015 +
	      BnuPrime+=s*tmp*tmp;
1016 +
	    }
1017 +
	}
1018 +
      Bnu=Bnu/u;
1019 +
      Bnu-=r;
1020 +
      BnuPrime=BnuPrime/(T*u);
1021 +
      if(Bnu<0)
1022 +
        nu_minus=nu;
1023 +
      else
1024 +
        nu_plus=nu;
1025 +
      if(fabs(nu_minus-nu_plus)<epsilon)
1026 +
	break;
1027 +
    }
1028 +
  if(fabs(Bnu)>epsilon)
1029 +
    Rcpp::Rcout << "Warning (Root): root not found to required precision" << endl;
1030 +
1031 +
  for(int i=0;i<u;i++)
1032 +
    {
1033 +
      s=exp((g[i]-nu)/T);
1034 +
      if(std::isinf(s)) p[i]=0.0;
1035 +
      else p[i]=1.0/(1.0+s);
1036 +
    }
1037 +
  //Rcpp::Rcout << " root (nu) = " << nu << " B(nu) = " << Bnu << endl;
1038 +
}
1039 +
double transductive_cost(double normWeights,double *Y, double *Outputs, int m, double lambda,double lambda_u)
1040 +
{
1041 +
  double F1=0.0,F2=0.0, o=0.0, y=0.0;
1042 +
  int u=0,l=0;
1043 +
  for(int i=0;i<m;i++)
1044 +
    {
1045 +
      o=Outputs[i];
1046 +
      y=Y[i];
1047 +
      if(y==0.0)
1048 +
	{F1 += fabs(o) > 1 ? 0 : (1 - fabs(o))*(1 - fabs(o)); u++;}
1049 +
      else
1050 +
	{F2 += y*o > 1 ? 0 : (1-y*o)*(1-y*o); l++;}
1051 +
    }
1052 +
  double F;
1053 +
  F = 0.5*(lambda*normWeights + lambda_u*F1/u + F2/l);
1054 +
  return F;
1055 +
}
1056 +
1057 +
double entropy(const double *p, int u)
1058 +
{
1059 +
  double h=0.0;
1060 +
  double q=0.0;
1061 +
  for(int i=0;i<u;i++)
1062 +
    {
1063 +
      q=p[i];
1064 +
      if(q>0 && q<1)
1065 +
	h+= -(q*LOG2(q) + (1-q)*LOG2(1-q));
1066 +
    }
1067 +
  return h/u;
1068 +
}
1069 +
double KL(const double *p, const double *q, int u)
1070 +
{
1071 +
  double h=0.0;
1072 +
  double p1=0.0;
1073 +
  double q1=0.0;
1074 +
  double g=0.0;
1075 +
  for(int i=0;i<u;i++)
1076 +
    {
1077 +
      p1=p[i];
1078 +
      q1=q[i];
1079 +
      if(p1>1-1e-8) p1-=1e-8;
1080 +
      if(p1<1-1e-8) p1+=1e-8;
1081 +
      if(q1>1-1e-8) q1-=1e-8;
1082 +
      if(q1<1-1e-8) q1+=1e-8;
1083 +
      g= (p1*LOG2(p1/q1) + (1-p1)*LOG2((1-p1)/(1-q1)));
1084 +
      if(fabs(g)<1e-12 || std::isnan(g)) g=0.0;
1085 +
      h+=g;
1086 +
    }
1087 +
    return h/u;
1088 +
}
1089 +
/* Prediction */
1090 +
/* the test file can potentially be huge..dont need to load it.*/
1091 +
void ssl_predict(char *inputs_file_name, const struct vector_double *Weights, struct vector_double *Outputs)
1092 +
{
1093 +
  FILE *fpin;
1094 +
  double *w = Weights->vec;
1095 +
  int n = Weights->d;
1096 +
  fpin = fopen(inputs_file_name, "r");
1097 +
  int m=0;
1098 +
  if (fpin == NULL)
1099 +
    {
1100 +
    Rcpp::stop("Cannot open input file\n");
1101 +
1102 +
    }
1103 +
 /* read and predict */
1104 +
  while (1)
1105 +
    {
1106 +
      int c = fgetc(fpin);
1107 +
      switch(c)
1108 +
	{
1109 +
	case '\n':
1110 +
	  ++m;
1111 +
	  break;
1112 +
	case EOF:
1113 +
	  goto out;
1114 +
	default:
1115 +
	  ;
1116 +
	}
1117 +
    }
1118 +
 out:
1119 +
  initialize(Outputs,m,0.0);
1120 +
  rewind(fpin);
1121 +
  int colind;
1122 +
  double val;
1123 +
  double t=0.0;
1124 +
  for(int i=0;i<m;i++)
1125 +
    {
1126 +
      t=0.0;
1127 +
      while(1)
1128 +
	{
1129 +
	  int c;
1130 +
	  do {
1131 +
	    c = getc(fpin);
1132 +
	    if(c=='\n') goto out2;
1133 +
	  } while(isspace(c));
1134 +
	  ungetc(c,fpin);
1135 +
	  int ret;
1136 +
	  ret = fscanf(fpin,"%d:%lf",&colind,&val);
1137 +
	  if (ret==-1) { Rcpp::Rcout << "EOF" << endl; }
1138 +
	  colind--;
1139 +
	  if(colind<n)
1140 +
	    t+=val*w[colind];
1141 +
	}
1142 +
    out2:
1143 +
      Outputs->vec[i]=t + w[n-1];
1144 +
    }
1145 +
}
1146 +
/* Evaluation -- confusion matrix, accuracy, precision, recall, prbep */
1147 +
/* roc area */
1148 +
void ssl_evaluate(struct vector_double *Outputs, struct vector_double *TrueLabels)
1149 +
{
1150 +
  double accuracy=0.0;
1151 +
  for(int i=0;i<Outputs->d;i++)
1152 +
    accuracy+=(Outputs->vec[i]*TrueLabels->vec[i])>0;
1153 +
  Rcpp::Rcout << "Accuracy = " << accuracy*100.0/Outputs->d << " %" << endl;
1154 +
}
1155 +
1156 +
/********************** UTILITIES ********************/
1157 +
double norm_square(const vector_double *A)
1158 +
{
1159 +
  double x=0.0, t=0.0;
1160 +
  for(int i=0;i<A->d;i++)
1161 +
    {
1162 +
      t=A->vec[i];
1163 +
      x+=t*t;
1164 +
    }
1165 +
  return x;
1166 +
}
1167 +
void initialize(struct vector_double *A, int k, double a)
1168 +
 {
1169 +
  double *vec = new double[k];
1170 +
  for(int i=0;i<k;i++)
1171 +
    vec[i]=a;
1172 +
  A->vec = vec;
1173 +
  A->d   = k;
1174 +
  return;
1175 +
}
1176 +
void initialize(struct vector_int *A, int k)
1177 +
{
1178 +
  int *vec = new int[k];
1179 +
  for(int i=0;i<k;i++)
1180 +
    vec[i]=i;
1181 +
  A->vec = vec;
1182 +
  A->d   = k;
1183 +
  return;
1184 +
}
1185 +
void Write(const char *file_name,
1186 +
	   const struct vector_double *somevector)
1187 +
{
1188 +
  FILE* fp = fopen(file_name,"w");
1189 +
  for(int i=0;i<somevector->d;i++)
1190 +
    fprintf(fp,"%g\n",somevector->vec[i]);
1191 +
  return;
1192 +
}
1193 +
void GetLabeledData(struct data *D, const struct data *Data)
1194 +
{
1195 +
  int *J = new int[Data->l];
1196 +
  D->C   = new double[Data->l];
1197 +
  D->Y   = new double[Data->l];
1198 +
  int nz=0;
1199 +
  int k=0;
1200 +
  int rowptrs_=Data->l;
1201 +
  for(int i=0;i<Data->m;i++)
1202 +
    {
1203 +
      if(Data->Y[i]!=0.0)
1204 +
	{
1205 +
	  J[k]=i;
1206 +
	  D->Y[k]=Data->Y[i];
1207 +
	  D->C[k]=1.0/Data->l;
1208 +
	  nz+=(Data->rowptr[i+1] - Data->rowptr[i]);
1209 +
	  k++;
1210 +
	}
1211 +
    }
1212 +
  D->val    = new double[nz];
1213 +
  D->colind = new int[nz];
1214 +
  D->rowptr = new int[rowptrs_+1];
1215 +
  nz=0;
1216 +
  for(int i=0;i<Data->l;i++)
1217 +
    {
1218 +
      D->rowptr[i]=nz;
1219 +
      for(int j=Data->rowptr[J[i]];j<Data->rowptr[J[i]+1];j++)
1220 +
	{
1221 +
	  D->val[nz] = Data->val[j];
1222 +
	  D->colind[nz] = Data->colind[j];
1223 +
	  nz++;
1224 +
	}
1225 +
    }
1226 +
  D->rowptr[rowptrs_]=nz;
1227 +
  D->nz=nz;
1228 +
  D->l=Data->l;
1229 +
  D->m=Data->l;
1230 +
  D->n=Data->n;
1231 +
  D->u=0;
1232 +
  delete [] J;
1233 +
}
1234 +
void SetData(struct data *a, int m,int n, int l,int u, int nz, double *VAL, int *R, int *C, double *Y, double *COSTS)
1235 +
{
1236 +
  a->m=m;
1237 +
  a->u=u;
1238 +
  a->l=m-u;
1239 +
  a->n=n;
1240 +
  a->nz=nz;
1241 +
  a->val=VAL;
1242 +
  a->rowptr=R;
1243 +
  a->colind=C;
1244 +
  a->Y=Y;
1245 +
  a->C=COSTS;
1246 +
  return;
1247 +
}
1248 +
void Clear(struct data *a)
1249 +
{
1250 +
  delete [] a->val;
1251 +
  delete [] a->rowptr;
1252 +
  delete [] a->colind;
1253 +
  delete [] a->Y;
1254 +
  delete [] a->C;
1255 +
  delete [] a;
1256 +
  return;
1257 +
}
1258 +
void Clear(struct vector_double *c)
1259 +
{ delete[] c->vec; delete [] c; return;}
1260 +
void Clear(struct vector_int *c)
1261 +
{ delete[] c->vec; delete [] c; return;}
1262 +
void Clear(struct options *opt)
1263 +
{ delete[] opt; delete [] opt; return;}

@@ -145,6 +145,8 @@
Loading
145 145
#' predict(tmod, se.fit = TRUE, interval = "confidence")
146 146
#' predict(tmod, se.fit = TRUE, interval = "confidence", rescaling = "lbg")
147 147
#' @keywords textmodel internal
148 +
#' @importFrom stats predict
149 +
#' @method predict textmodel_wordscores
148 150
#' @export
149 151
#' @importFrom stats qnorm median sd
150 152
predict.textmodel_wordscores <- function(object,

@@ -120,6 +120,7 @@
Loading
120 120
#' @seealso [textmodel_lr()]
121 121
#' @keywords textmodel internal
122 122
#' @importFrom stats predict
123 +
#' @method predict textmodel_lr
123 124
#' @export
124 125
predict.textmodel_lr <- function(object, newdata = NULL,
125 126
                                 type = c("class", "probability"),

@@ -110,6 +110,8 @@
Loading
110 110
#' @param ... unused
111 111
#' @return `predict()` returns a predicted [textmodel_lsa] object, projecting the patterns onto
112 112
#' new data.
113 +
#' @importFrom stats predict
114 +
#' @method predict textmodel_lsa
113 115
#' @keywords textmodel internal
114 116
#' @export
115 117
predict.textmodel_lsa <- function(object, newdata = NULL, ...) {

@@ -201,6 +201,8 @@
Loading
201 201
#' (tmod <- textmodel_nb(data_dfm_lbgexample, y = c("A", "A", "B", "C", "C", NA)))
202 202
#' predict(tmod)
203 203
#' predict(tmod, type = "logposterior")
204 +
#' @importFrom stats predict
205 +
#' @method predict textmodel_nb
204 206
#' @keywords textmodel internal
205 207
#' @export
206 208
predict.textmodel_nb <- function(object, newdata = NULL,

@@ -230,6 +230,8 @@
Loading
230 230
#' out-of-sample data.
231 231
#' @param object a fitted wordfish model
232 232
#' @inheritParams predict.textmodel_wordscores
233 +
#' @importFrom stats predict
234 +
#' @method predict textmodel_wordfish
233 235
#' @keywords textmodel internal
234 236
#' @export
235 237
predict.textmodel_wordfish <- function(object,

@@ -452,6 +452,8 @@
Loading
452 452
#' @return `predict()` returns a list of predicted affinity textmodel
453 453
#'   quantities.
454 454
#' @importFrom methods new
455 +
#' @importFrom stats predict
456 +
#' @method predict textmodel_affinity
455 457
#' @keywords textmodel internal
456 458
#' @seealso [influence.predict.textmodel_affinity()] for methods of
457 459
#'   computing the influence of particular features from a predicted

@@ -0,0 +1,200 @@
Loading
1 +
/*    Copyright 2006 Vikas Sindhwani (vikass@cs.uchicago.edu)
2 +
      SVM-lin: Fast SVM Solvers for Supervised and Semi-supervised Learning
3 +
4 +
      This file is part of SVM-lin.
5 +
6 +
      SVM-lin is free software; you can redistribute it and/or modify
7 +
      it under the terms of the GNU General Public License as published by
8 +
      the Free Software Foundation; either version 2 of the License, or
9 +
      (at your option) any later version.
10 +
11 +
      SVM-lin is distributed in the hope that it will be useful,
12 +
      but WITHOUT ANY WARRANTY; without even the implied warranty of
13 +
      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 +
      GNU General Public License for more details.
15 +
16 +
      You should have received a copy of the GNU General Public License
17 +
      along with SVM-lin (see gpl.txt); if not, write to the Free Software
18 +
      Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
19 +
*/
20 +
#ifndef _svmlin_H
21 +
#define _svmlin_H
22 +
#include <vector>
23 +
#include <ctime>
24 +
25 +
using namespace std;
26 +
27 +
/* OPTIMIZATION CONSTANTS */
28 +
#define CGITERMAX 10000 /* maximum number of CGLS iterations */
29 +
#define SMALL_CGITERMAX 10 /* for heuristic 1 in reference [2] */
30 +
#define EPSILON   1e-6 /* most tolerances are set to this value */
31 +
#define BIG_EPSILON 0.01 /* for heuristic 2 in reference [2] */
32 +
#define RELATIVE_STOP_EPS 1e-9 /* for L2-SVM-MFN relative stopping criterion */
33 +
#define MFNITERMAX 50 /* maximum number of MFN iterations */
34 +
#define TSVM_ANNEALING_RATE 1.5 /* rate at which lambda_u is increased in TSVM */
35 +
#define TSVM_LAMBDA_SMALL 1e-5 /* lambda_u starts from this value */
36 +
#define DA_ANNEALING_RATE 1.5 /* annealing rate for DA */
37 +
#define DA_INIT_TEMP 10 /* initial temperature relative to lambda_u */
38 +
#define DA_INNER_ITERMAX 100 /* maximum fixed temperature iterations for DA */
39 +
#define DA_OUTER_ITERMAX 30 /* maximum number of outer loops for DA */
40 +
41 +
#define VERBOSE_CGLS 0
42 +
43 +
/* Data: Input examples are stored in sparse (Compressed Row Storage) format */
44 +
struct data
45 +
{
46 +
  int m; /* number of examples */
47 +
  int l; /* number of labeled examples */
48 +
  int u; /* number of unlabeled examples l+u = m */
49 +
  int n; /* number of features */
50 +
  int nz; /* number of non-zeros */
51 +
  double *val; /* data values (nz elements) [CRS format] */
52 +
  int *rowptr; /* n+1 vector [CRS format] */
53 +
  int *colind; /* nz elements [CRS format] */
54 +
  double *Y;   /* labels */
55 +
  double *C;   /* cost associated with each example */
56 +
};
57 +
58 +
struct vector_double /* defines a vector of doubles */
59 +
{
60 +
  int d; /* number of elements */
61 +
  double *vec; /* ptr to vector elements*/
62 +
};
63 +
64 +
65 +
66 +
struct vector_int /* defines a vector of ints for index subsets */
67 +
{
68 +
  int d; /* number of elements */
69 +
  int *vec; /* ptr to vector elements */
70 +
};
71 +
72 +
enum { RLS, SVM, TSVM, DA_SVM }; /* currently implemented algorithms */
73 +
74 +
struct options
75 +
{
76 +
  /* user options */
77 +
  int algo; /* 1 to 4 for RLS,SVM,TSVM,DASVM */
78 +
  double lambda; /* regularization parameter */
79 +
  double lambda_u; /* regularization parameter over unlabeled examples */
80 +
  int S; /* maximum number of TSVM switches per fixed-weight label optimization */
81 +
  double R; /* expected fraction of unlabeled examples in positive class */
82 +
  double Cp; /* cost for positive examples */
83 +
  double Cn; /* cost for negative examples */
84 +
  /*  internal optimization options */
85 +
  double epsilon; /* all tolerances */
86 +
  int cgitermax;  /* max iterations for CGLS */
87 +
  int mfnitermax; /* max iterations for L2_SVM_MFN */
88 +
  bool verbose; /* Should the output be verbose? */
89 +
};
90 +
91 +
class timer { /* to output run time */
92 +
protected:
93 +
  double start, finish;
94 +
public:
95 +
  vector<double> times;
96 +
  void record() {
97 +
    times.push_back(time());
98 +
  }
99 +
  void reset_vectors() {
100 +
    times.erase(times.begin(), times.end());
101 +
  }
102 +
  void restart() { start = clock(); }
103 +
  void stop() { finish = clock(); }
104 +
  double time() const { return ((double)(finish - start))/CLOCKS_PER_SEC; }
105 +
};
106 +
class Delta { /* used in line search */
107 +
 public:
108 +
   Delta() {delta=0.0; index=0;s=0;};
109 +
   double delta;
110 +
   int index;
111 +
   int s;
112 +
};
113 +
inline bool operator<(const Delta& a , const Delta& b) { return (a.delta < b.delta);}
114 +
115 +
void initialize(struct vector_double *A, int k, double a);
116 +
/* initializes a vector_double to be of length k, all elements set to a */
117 +
void initialize(struct vector_int *A, int k);
118 +
/* initializes a vector_int to be of length k, elements set to 1,2..k. */
119 +
void SetData(struct data *Data, int m,int n, int l,int u, int nz,
120 +
	     double *VAL, int *R, int *C, double *Y, double *COSTS); /* sets data fields */
121 +
void GetLabeledData(struct data *Data_Labeled, const struct data *Data);
122 +
/* extracts labeled data from Data and copies it into Data_Labeled */
123 +
void Write(const char *file_name, const struct vector_double *somevector);
124 +
/* writes a vector into filename, one element per line */
125 +
void Clear(struct data *a); /* deletes a */
126 +
void Clear(struct vector_double *a); /* deletes a */
127 +
void Clear(struct vector_int *a); /* deletes a */
128 +
double norm_square(const vector_double *A); /* returns squared length of A */
129 +
130 +
/* ssl_train: takes data, options, uninitialized weight and output
131 +
   vector_doubles, routes it to the algorithm */
132 +
/* the learnt weight vector and the outputs it gives on the data matrix are saved */
133 +
void ssl_train(struct data *Data,
134 +
	       struct options *Options,
135 +
	       struct vector_double *W, /* weight vector */
136 +
	       struct vector_double *O); /* output vector */
137 +
138 +
/* Main svmlin Subroutines */
139 +
/*ssl_predict: reads test inputs from input_file_name, a weight vector, and an
140 +
 uninitialized outputs vector. Performs */
141 +
void ssl_predict(char *inputs_file_name, const struct vector_double *Weights,
142 +
		 struct vector_double *Outputs);
143 +
/* ssl_evaluate: if test labels are given in the vector True, and predictions in vector Output,
144 +
   this code prints out various performance statistics. Currently only accuracy. */
145 +
void ssl_evaluate(struct vector_double *Outputs,struct vector_double *True);
146 +
147 +
/* svmlin algorithms and their subroutines */
148 +
149 +
/* Conjugate Gradient for Sparse Linear Least Squares Problems */
150 +
/* Solves: min_w 0.5*Options->lamda*w'*w + 0.5*sum_{i in Subset} Data->C[i] (Y[i]- w' x_i)^2 */
151 +
/* over a subset of examples x_i specified by vector_int Subset */
152 +
int CGLS(const struct data *Data,
153 +
	 const struct options *Options,
154 +
	 const struct vector_int *Subset,
155 +
	 struct vector_double *Weights,
156 +
	 struct vector_double *Outputs);
157 +
158 +
/* Linear Modified Finite Newton L2-SVM*/
159 +
/* Solves: min_w 0.5*Options->lamda*w'*w + 0.5*sum_i Data->C[i] max(0,1 - Y[i] w' x_i)^2 */
160 +
int L2_SVM_MFN(const struct data *Data,
161 +
	       struct options *Options,
162 +
	       struct vector_double *Weights,
163 +
	       struct vector_double *Outputs,
164 +
	       int ini); /* use ini=0 if no good starting guess for Weights, else 1 */
165 +
double line_search(double *w,
166 +
                   double *w_bar,
167 +
                   double lambda,
168 +
                   double *o,
169 +
                   double *o_bar,
170 +
                   double *Y,
171 +
                   double *C,
172 +
                   int d,
173 +
                   int l);
174 +
175 +
/* Transductive L2-SVM */
176 +
/* Solves : min_(w, Y[i],i in UNlabeled) 0.5*Options->lamda*w'*w + 0.5*(1/Data->l)*sum_{i in labeled} max(0,1 - Y[i] w' x_i)^2 + 0.5*(Options->lambda_u/Data->u)*sum_{i in UNlabeled} max(0,1 - Y[i] w' x_i)^2
177 +
 subject to: (1/Data->u)*sum_{i in UNlabeled} max(0,Y[i]) = Options->R */
178 +
int   TSVM_MFN(const struct data *Data,
179 +
	      struct options *Options,
180 +
	      struct vector_double *Weights,
181 +
	      struct vector_double *Outputs);
182 +
int switch_labels(double* Y, double* o, int* JU, int u, int S);
183 +
184 +
/* Deterministic Annealing*/
185 +
int DA_S3VM(struct data *Data,
186 +
	   struct options *Options,
187 +
	   struct vector_double *Weights,
188 +
	   struct vector_double *Outputs);
189 +
void optimize_p(const double* g, int u, double T, double r, double*p);
190 +
int optimize_w(const struct data *Data,
191 +
	       const  double *p,
192 +
	       struct options *Options,
193 +
	       struct vector_double *Weights,
194 +
	       struct vector_double *Outputs,
195 +
	       int ini);
196 +
double transductive_cost(double normWeights,double *Y, double *Outputs, int m, double lambda,double lambda_u);
197 +
double entropy(const  double *p, int u);
198 +
double KL(const  double *p, const  double *q, int u); /* KL-divergence */
199 +
200 +
#endif
Files Coverage
R 84.11%
src 47.32%
Project Totals (20 files) 65.74%

No yaml found.

Create your codecov.yml to customize your Codecov experience

Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading