quanteda / quanteda.textmodels

Compare e591752 ... +4 ... c468e51

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.

Showing 3 of 5 files from the diff.
Newly tracked file
src/ssl.h created.
Newly tracked file
src/ssl.cpp created.
Newly tracked file
src/svmlin.cpp created.
Other files ignored by Codecov

@@ -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

@@ -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;}

@@ -0,0 +1,128 @@
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 +
#include <Rcpp.h>
21 +
#include <iostream>
22 +
#include <stdlib.h>
23 +
#include <stdio.h>
24 +
#include <ctype.h>
25 +
#include <cstring>
26 +
#include "ssl.h"
27 +
using namespace Rcpp;
28 +
29 +
struct options *Options = new options[1];
30 +
struct data *Data = new data[1];
31 +
struct vector_double *Weights = new vector_double[1];
32 +
struct vector_double *Outputs = new vector_double[1];
33 +
struct vector_double *Labels = new vector_double[1];
34 +
35 +
// [[Rcpp::export]]
36 +
List svmlin_rcpp(S4 X,
37 +
                 NumericVector y,
38 +
                 int l,
39 +
                 int algorithm,
40 +
                 double lambda,
41 +
                 double lambda_u,
42 +
                 int max_switch,
43 +
                 double pos_frac,
44 +
                 double Cp,
45 +
                 double Cn,
46 +
                 NumericVector costs,
47 +
                 bool verbose) {
48 +
  // Set options
49 +
  Options->algo = algorithm;
50 +
  Options->lambda=lambda;
51 +
  Options->lambda_u=lambda_u;
52 +
  Options->S=max_switch;
53 +
  Options->R=pos_frac;
54 +
  Options->epsilon=EPSILON;
55 +
  Options->cgitermax=CGITERMAX;
56 +
  Options->mfnitermax=MFNITERMAX;
57 +
  Options->Cp = Cp;
58 +
  Options->Cn = Cn;
59 +
  Options->verbose = verbose;
60 +
61 +
  NumericVector ycop( y.begin(), y.end() );
62 +
  NumericVector costcop( costs.begin(), costs.end() );
63 +
//   Rprintf("Step 1\n");
64 +
//   size_t size = ((DoubleVector)X.slot("x")).length()+((DoubleVector)Xu.slot("x")).length();
65 +
//   std::vector<double> vals(size);
66 +
//   std::vector<double>::iterator it = vals.begin();
67 +
//   vals.insert(it,((DoubleVector)X.slot("x")).begin(),((DoubleVector)Xu.slot("x")).end());
68 +
//   it = vals.begin()+((DoubleVector)X.slot("x")).length();
69 +
//   vals.insert(it,((DoubleVector)Xu.slot("x")).begin(),((DoubleVector)Xu.slot("x")).end());
70 +
//
71 +
//   Rprintf("Step 2\n");
72 +
//
73 +
//   size = ((IntegerVector)X.slot("i")).length()+((IntegerVector)Xu.slot("i")).length();
74 +
//   std::vector<int> colinds(size);
75 +
//   std::vector<int>::iterator it2 = colinds.begin();
76 +
//   colinds.insert(it2,((IntegerVector)X.slot("i")).begin(),((IntegerVector)X.slot("i")).end());
77 +
//   it2 = colinds.begin() + ((IntegerVector)X.slot("i")).length();
78 +
//   colinds.insert(it2,((IntegerVector)Xu.slot("i")).begin(),((IntegerVector)Xu.slot("i")).end());
79 +
//
80 +
//   size = ((IntegerVector)X.slot("p")).length()+((IntegerVector)Xu.slot("p")).length();
81 +
//   std::vector<int> rowpts(size);
82 +
//   it2 = rowpts.begin();
83 +
//   rowpts.insert(it2,((IntegerVector)X.slot("p")).begin(),((IntegerVector)X.slot("p")).end());
84 +
//   it2 = rowpts.begin() + ((IntegerVector)X.slot("p")).length();
85 +
//   rowpts.insert(it2,((IntegerVector)Xu.slot("p")).begin(),((IntegerVector)Xu.slot("p")).end());
86 +
87 +
  // R data to svmlin data structure
88 +
  Data->m=((IntegerVector)X.slot("Dim"))[1];
89 +
  Data->l=l;
90 +
  Data->u=Data->m-Data->l;
91 +
  Data->n=((IntegerVector)X.slot("Dim"))[0];
92 +
  Data->nz=((DoubleVector)X.slot("x")).size();
93 +
  Data->val=((DoubleVector)X.slot("x")).begin();
94 +
  Data->rowptr=((IntegerVector)X.slot("p")).begin();
95 +
  Data->colind=((IntegerVector)X.slot("i")).begin();
96 +
  Data->Y=ycop.begin();
97 +
  Data->C=costcop.begin();
98 +
99 +
  // TODO: load correct costs for unlabeled data.
100 +
  if (Options->verbose) {
101 +
    Rcout << "  Input Data Matrix Statistics:" << endl;
102 +
    Rcout << "      Examples: " << Data->m << endl;
103 +
    Rcout << "      Features: " << Data->n << " (including bias feature)" << endl;
104 +
    Rcout << "      Non-zeros:  " << Data->nz << " (including bias features)" << endl;
105 +
    Rcout << "      Average sparsity: " << Data->nz*1.0/Data->m << " non-zero features per example." << endl;
106 +
  }
107 +
    //   for (int i = 0; i<((DoubleVector)X.slot("x")).length();i++) {
108 +
//     Rprintf("val: %f \n",Data->val[i]);
109 +
//   }
110 +
//   for (int i = 0; i<((IntegerVector)X.slot("i")).length();i++) {
111 +
//     Rprintf("col: %d \n",Data->colind[i]);
112 +
//   }
113 +
//   for (int i = 0; i<((IntegerVector)X.slot("p")).length();i++) {
114 +
//     Rprintf("row: %d \n",Data->rowptr[i]);
115 +
//   }
116 +
117 +
118 +
  // Run
119 +
  ssl_train(Data,
120 +
            Options,
121 +
            Weights,
122 +
            Outputs);
123 +
//Clear(Data);
124 +
125 +
  return Rcpp::List::create(Rcpp::Named("Weights") = std::vector<double>(Weights->vec, Weights->vec+Weights->d),
126 +
                              Rcpp::Named("Outputs") = std::vector<double>(Outputs->vec, Outputs->vec+Outputs->d)
127 +
                              );
128 +
}

Learn more Showing 3 files with coverage changes found.

New file src/ssl.h
New
Loading file...
New file src/ssl.cpp
New
Loading file...
New file src/svmlin.cpp
New
Loading file...
Files Coverage
R 85.31%
src -64.58% 29.51%
Project Totals (19 files) 56.03%
Loading