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 1
void ssl_train(struct data *Data,
37
	       struct options *Options,
38
	       struct vector_double *Weights,
39
	       struct vector_double *Outputs)
40
{
41
  // initialize
42 1
  initialize(Weights,Data->n,0.0);
43 1
  initialize(Outputs,Data->m,0.0);
44 1
  vector_int    *Subset  = new vector_int[1];
45 1
  initialize(Subset,Data->m);
46
  // call the right algorithm
47 1
  int optimality = 0;
48 1
  switch(Options->algo)
49
    {
50
    case -1:
51 0
      if (Options->verbose) { Rcpp::Rcout << "Regularized Least Squares Regression (CGLS)\n" << endl; }
52 0
      optimality=CGLS(Data,Options,Subset,Weights,Outputs);
53 0
      break;
54
    case RLS:
55 0
      if (Options->verbose) { Rcpp::Rcout << "Regularized Least Squares Classification (CGLS)\n" << endl; }
56 0
      optimality=CGLS(Data,Options,Subset,Weights,Outputs);
57 0
      break;
58
    case SVM:
59 1
      if (Options->verbose) { Rcpp::Rcout << "Modified Finite Newton L2-SVM (L2-SVM-MFN)\n" << endl; }
60 1
      optimality=L2_SVM_MFN(Data,Options,Weights,Outputs,0);
61 1
      break;
62
    case TSVM:
63 0
      if (Options->verbose) { Rcpp::Rcout << "Transductive L2-SVM (TSVM)\n" << endl; }
64 0
      optimality=TSVM_MFN(Data,Options,Weights,Outputs);
65 0
      break;
66
    case DA_SVM:
67 0
      if (Options->verbose) { Rcpp::Rcout << "Deterministic Annealing Semi-supervised L2-SVM (DAS3VM)\n" << endl; }
68 0
      optimality=DA_S3VM(Data,Options,Weights,Outputs);
69 0
      break;
70
    default:
71
      ;
72
    }
73 1
  if (Options->verbose) { Rcpp::Rcout << "Optimality:" << optimality << endl; }
74 1
  return;
75
}
76 1
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 1
  timer tictoc;
86 1
  tictoc.restart();
87 1
  int active = Subset->d;
88 1
  int *J = Subset->vec;
89 1
  double *val = Data->val;
90 1
  int *row = Data->rowptr;
91 1
  int *col = Data->colind;
92 1
  double *Y = Data->Y;
93 1
  double *C = Data->C;
94 1
  int n  = Data->n;
95 1
  double lambda = Options->lambda;
96 1
  int cgitermax = Options->cgitermax;
97 1
  double epsilon = Options->epsilon;
98 1
  double *beta = Weights->vec;
99 1
  double *o  = Outputs->vec;
100
  // initialize z
101 1
  double *z = new double[active];
102 1
  double *q = new double[active];
103 1
  int ii=0;
104 1
  for(int i = active ; i-- ;){
105 1
    ii=J[i];
106 1
    z[i]  = C[ii]*(Y[ii] - o[ii]);
107
  }
108 1
  double *r = new double[n];
109 1
  for(int i = n ; i-- ;)
110 1
    r[i] = 0.0;
111 1
  for(int j=0; j < active; j++)
112
    {
113 1
      ii=J[j];
114 1
      for(int i=row[ii]; i < row[ii+1]; i++)
115 1
	r[col[i]]+=val[i]*z[j];
116
    }
117 1
  double *p = new double[n];
118 1
  double omega1 = 0.0;
119 1
  for(int i = n ; i-- ;)
120
    {
121 1
      r[i] -= lambda*beta[i];
122 1
      p[i] = r[i];
123 1
      omega1 += r[i]*r[i];
124
    }
125 1
  double omega_p = omega1;
126 1
  double omega_q = 0.0;
127 1
  double inv_omega2 = 1/omega1;
128 1
  double scale = 0.0;
129 1
  double omega_z=0.0;
130 1
  double gamma = 0.0;
131 1
  int cgiter = 0;
132 1
  int optimality = 0;
133 1
  double epsilon2 = epsilon*epsilon;
134
  // iterate
135 1
  while(cgiter < cgitermax)
136
    {
137 1
      cgiter++;
138 1
      omega_q=0.0;
139 1
      double t=0.0;
140
      int i,j;
141
      // #pragma omp parallel for private(i,j)
142 1
      for(i=0; i < active; i++)
143
	{
144 1
	  ii=J[i];
145 1
	  t=0.0;
146 1
	  for(j=row[ii]; j < row[ii+1]; j++)
147 1
	    t+=val[j]*p[col[j]];
148 1
	  q[i]=t;
149 1
	  omega_q += C[ii]*t*t;
150
	}
151 1
      gamma = omega1/(lambda*omega_p + omega_q);
152 1
      inv_omega2 = 1/omega1;
153 1
      for(int i = n ; i-- ;)
154
	{
155 1
	  r[i] = 0.0;
156 1
	  beta[i] += gamma*p[i];
157
	}
158 1
      omega_z=0.0;
159 1
      for(int i = active ; i-- ;)
160
	{
161 1
	  ii=J[i];
162 1
	  o[ii] += gamma*q[i];
163 1
	  z[i] -= gamma*C[ii]*q[i];
164 1
	  omega_z+=z[i]*z[i];
165
	}
166 1
      for(int j=0; j < active; j++)
167
	{
168 1
	  ii=J[j];
169 1
	  t=z[j];
170 1
	  for(int i=row[ii]; i < row[ii+1]; i++)
171 1
	    r[col[i]]+=val[i]*t;
172
	}
173 1
      omega1 = 0.0;
174 1
      for(int i = n ; i-- ;)
175
	{
176 1
	  r[i] -= lambda*beta[i];
177 1
	  omega1 += r[i]*r[i];
178
	}
179
      if(VERBOSE_CGLS)
180
        if (Options->verbose) { Rcpp::Rcout << "..." << cgiter << " ( " << omega1 << " )" ; }
181 1
      if(omega1 < epsilon2*omega_z)
182
	{
183 1
	  optimality=1;
184 1
	  break;
185
	}
186 1
      omega_p=0.0;
187 1
      scale=omega1*inv_omega2;
188 1
      for(int i = n ; i-- ;)
189
	{
190 1
	  p[i] = r[i] + p[i]*scale;
191 1
	  omega_p += p[i]*p[i];
192
	}
193
    }
194
  if(VERBOSE_CGLS)
195
    if (Options->verbose) { Rcpp::Rcout << "...Done." << endl; }
196 1
  tictoc.stop();
197 1
    if (Options->verbose) { Rcpp::Rcout << "CGLS converged in " << cgiter << " iteration(s) and " << tictoc.time() << " seconds." << endl; }
198 1
  delete[] z;
199 1
  delete[] q;
200 1
  delete[] r;
201 1
  delete[] p;
202 1
  return optimality;
203
}
204 1
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 1
  timer tictoc;
212 1
  tictoc.restart();
213 1
  double *val = Data->val;
214 1
  int *row = Data->rowptr;
215 1
  int *col = Data->colind;
216 1
  double *Y = Data->Y;
217 1
  double *C = Data->C;
218 1
  int n  = Data->n;
219 1
  int m  = Data->m;
220 1
  double lambda = Options->lambda;
221
  double epsilon;
222 1
  double *w = Weights->vec;
223 1
  double *o = Outputs->vec;
224 1
  double F_old = 0.0;
225 1
  double F = 0.0;
226 1
  double diff=0.0;
227 1
  vector_int *ActiveSubset = new vector_int[1];
228 1
  ActiveSubset->vec = new int[m];
229 1
  ActiveSubset->d = m;
230
  // initialize
231 1
  if(ini==0) {
232 1
    epsilon=BIG_EPSILON;
233 1
    Options->cgitermax=SMALL_CGITERMAX;
234 1
    Options->epsilon=BIG_EPSILON;
235
  }
236 0
  else {epsilon = Options->epsilon;}
237 1
  for(int i=0;i<n;i++) F+=w[i]*w[i];
238 1
  F=0.5*lambda*F;
239 1
  int active=0;
240 1
  int inactive=m-1; // l-1
241 1
  for(int i=0; i<m ; i++)
242
    {
243 1
      diff=1-Y[i]*o[i];
244 1
      if(diff>0)
245
	{
246 1
	  ActiveSubset->vec[active]=i;
247 1
	  active++;
248 1
	  F+=0.5*C[i]*diff*diff;
249
	}
250
      else
251
	{
252 0
	  ActiveSubset->vec[inactive]=i;
253 0
	  inactive--;
254
	}
255
    }
256 1
  ActiveSubset->d=active;
257 1
  int iter=0;
258 1
  int opt=0;
259 1
  int opt2=0;
260 1
  vector_double *Weights_bar = new vector_double[1];
261 1
  vector_double *Outputs_bar = new vector_double[1];
262 1
  double *w_bar = new double[n];
263 1
  double *o_bar = new double[m];
264 1
  Weights_bar->vec=w_bar;
265 1
  Outputs_bar->vec=o_bar;
266 1
  Weights_bar->d=n;
267 1
  Outputs_bar->d=m;
268 1
  double delta=0.0;
269 1
  double t=0.0;
270 1
  int ii = 0;
271 1
  while(iter<MFNITERMAX)
272
    {
273 1
      iter++;
274 1
    if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN Iteration# " << iter << " (" << active << " active examples, " << " objective_value = " << F << ")" << endl; }
275 1
      for(int i=n; i-- ;)
276 1
	w_bar[i]=w[i];
277 1
      for(int i=m; i-- ;)
278 1
	o_bar[i]=o[i];
279

280 1
      opt=CGLS(Data,Options,ActiveSubset,Weights_bar,Outputs_bar);
281 1
      for(int i=active; i < m; i++)
282
	{
283 0
	  ii=ActiveSubset->vec[i];
284 0
	  t=0.0;
285 0
	  for(int j=row[ii]; j < row[ii+1]; j++)
286 0
	    t+=val[j]*w_bar[col[j]];
287 0
	  o_bar[ii]=t;
288
	}
289 1
      if(ini==0) {Options->cgitermax=CGITERMAX; ini=1;};
290 1
      opt2=1;
291 1
      for(int i=0;i<m;i++)
292
	{
293 1
	  ii=ActiveSubset->vec[i];
294 1
	  if(i<active)
295 1
	    opt2=(opt2 && (Y[ii]*o_bar[ii]<=1+epsilon));
296
	  else
297 0
	    opt2=(opt2 && (Y[ii]*o_bar[ii]>=1-epsilon));
298 1
	  if(opt2==0) break;
299
	}
300 1
      if(opt && opt2) // l
301
	{
302 1
	  if(epsilon==BIG_EPSILON)
303
	    {
304 1
	      epsilon=EPSILON;
305 1
	      Options->epsilon=EPSILON;
306 1
	      if (Options->verbose) { Rcpp::Rcout << "  epsilon = " << BIG_EPSILON << " case converged (speedup heuristic 2). Continuing with epsilon=" <<  EPSILON << endl; }
307 1
	      continue;
308
	    }
309
	  else
310
	    {
311 1
	      for(int i=n; i-- ;)
312 1
		w[i]=w_bar[i];
313 1
	      for(int i=m; i-- ;)
314 1
		o[i]=o_bar[i];
315 1
	       delete[] ActiveSubset->vec;
316 1
	       delete[] ActiveSubset;
317 1
	       delete[] o_bar;
318 1
	       delete[] w_bar;
319 1
	       delete[] Weights_bar;
320 1
	       delete[] Outputs_bar;
321 1
	       tictoc.stop();
322 1
	       if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (optimality) in " << iter << " iteration(s) and "<< tictoc.time() << " seconds. \n" << endl; }
323 1
	       return 1;
324
	    }
325
	}
326 0
      if (Options->verbose) { Rcpp::Rcout << " " ; }
327 0
      delta=line_search(w,w_bar,lambda,o,o_bar,Y,C,n,m);
328 0
      if (Options->verbose) { Rcpp::Rcout << "LINE_SEARCH delta = " << delta << endl; }
329 0
      F_old=F;
330 0
      F=0.0;
331 0
      for(int i=n; i-- ;){
332 0
	w[i]+=delta*(w_bar[i]-w[i]);
333 0
	F+=w[i]*w[i];
334
      }
335 0
      F=0.5*lambda*F;
336 0
      active=0;
337 0
      inactive=m-1;
338 0
      for(int i=0; i<m ; i++)
339
	{
340 0
	  o[i]+=delta*(o_bar[i]-o[i]);
341 0
	  diff=1-Y[i]*o[i];
342 0
	  if(diff>0)
343
	    {
344 0
	      ActiveSubset->vec[active]=i;
345 0
	      active++;
346 0
	      F+=0.5*C[i]*diff*diff;
347
	    }
348
	  else
349
	    {
350 0
	      ActiveSubset->vec[inactive]=i;
351 0
	      inactive--;
352
	    }
353
	}
354 0
      ActiveSubset->d=active;
355 0
      if(fabs(F-F_old)<RELATIVE_STOP_EPS*fabs(F_old))
356
	{
357 0
        if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (rel. criterion) in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl; }
358 0
	  return 2;
359
	}
360
    }
361 0
  delete[] ActiveSubset->vec;
362 0
  delete[] ActiveSubset;
363 0
  delete[] o_bar;
364 0
  delete[] w_bar;
365 0
  delete[] Weights_bar;
366 0
  delete[] Outputs_bar;
367 0
  tictoc.stop();
368 0
  if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged (max iter exceeded) in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl; }
369 0
  return 0;
370
}
371

372 0
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 0
   double omegaL = 0.0;
383 0
   double omegaR = 0.0;
384 0
   double diff=0.0;
385 0
   for(int i=d; i--; )
386
       {
387 0
         diff=w_bar[i]-w[i];
388 0
         omegaL+=w[i]*diff;
389 0
         omegaR+=w_bar[i]*diff;
390
   }
391 0
   omegaL=lambda*omegaL;
392 0
   omegaR=lambda*omegaR;
393 0
   double L=0.0;
394 0
   double R=0.0;
395 0
   int ii=0;
396 0
   for(int i=0;i<l;i++)
397
       {
398 0
         if(Y[i]*o[i]<1)
399
	   {
400 0
	     diff=C[i]*(o_bar[i]-o[i]);
401 0
	     L+=(o[i]-Y[i])*diff;
402 0
	     R+=(o_bar[i]-Y[i])*diff;
403
	   }
404
       }
405 0
   L+=omegaL;
406 0
   R+=omegaR;
407 0
   Delta* deltas=new Delta[l];
408 0
   int p=0;
409 0
   for(int i=0;i<l;i++)
410
     {
411 0
       diff=Y[i]*(o_bar[i]-o[i]);
412

413 0
       if(Y[i]*o[i]<1)
414
	 {
415 0
	   if(diff>0)
416
	     {
417 0
	       deltas[p].delta=(1-Y[i]*o[i])/diff;
418 0
	       deltas[p].index=i;
419 0
	       deltas[p].s=-1;
420 0
	       p++;
421
	     }
422
	 }
423
       else
424
	 {
425 0
	   if(diff<0)
426
	     {
427 0
	       deltas[p].delta=(1-Y[i]*o[i])/diff;
428 0
	       deltas[p].index=i;
429 0
	       deltas[p].s=1;
430 0
	       p++;
431
	     }
432
	 }
433
     }
434 0
   sort(deltas,deltas+p);
435 0
   double delta_prime=0.0;
436 0
   for(int i=0;i<p;i++)
437
     {
438 0
       delta_prime = L + deltas[i].delta*(R-L);
439 0
       if(delta_prime>=0)
440 0
	 break;
441 0
       ii=deltas[i].index;
442 0
       diff=(deltas[i].s)*C[ii]*(o_bar[ii]-o[ii]);
443 0
       L+=diff*(o[ii]-Y[ii]);
444 0
       R+=diff*(o_bar[ii]-Y[ii]);
445
     }
446 0
   delete [] deltas;
447 0
   return (-L/(R-L));
448
}
449 0
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 0
  timer tictoc;
456 0
  tictoc.restart();
457 0
  struct data *Data_Labeled = new data[1];
458 0
  struct vector_double *Outputs_Labeled = new vector_double[1];
459 0
  initialize(Outputs_Labeled,Data->l,0.0);
460 0
  if (Options->verbose) { Rcpp::Rcout << "Initializing weights, unknown labels" << endl; }
461 0
  GetLabeledData(Data_Labeled,Data); /* gets labeled data and sets C=1/l */
462 0
  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 0
  int p=0,q=0;
467 0
  double t=0.0;
468 0
  int *JU = new int[Data->u];
469 0
  double *ou = new double[Data->u];
470 0
  double lambda_0 = TSVM_LAMBDA_SMALL;
471 0
  for(int i=0;i<Data->m;i++)
472
    {
473 0
      if(Data->Y[i]==0.0)
474
	{
475 0
	  t=0.0;
476 0
	  for(int j=Data->rowptr[i]; j < Data->rowptr[i+1]; j++)
477 0
	    t+=Data->val[j]*Weights->vec[Data->colind[j]];
478 0
	  Outputs->vec[i]=t;
479 0
	  Data->C[i]=lambda_0*1.0/Data->u;
480 0
	  JU[q]=i;
481 0
	  ou[q]=t;
482 0
	  q++;
483
	}
484
      else
485
	{
486 0
	  Outputs->vec[i]=Outputs_Labeled->vec[p];
487 0
	  Data->C[i]=1.0/Data->l;
488 0
	  p++;
489
	}
490
    }
491 0
  nth_element(ou,ou+int((1-Options->R)*Data->u-1),ou+Data->u);
492 0
  double thresh=*(ou+int((1-Options->R)*Data->u)-1);
493 0
  delete [] ou;
494 0
  for(int i=0;i<Data->u;i++)
495
    {
496 0
      if(Outputs->vec[JU[i]]>thresh)
497 0
	Data->Y[JU[i]]=1.0;
498
      else
499 0
	Data->Y[JU[i]]=-1.0;
500
    }
501 0
  for(int i=0;i<Data->n;i++)
502 0
    Weights->vec[i]=0.0;
503 0
  for(int i=0;i<Data->m;i++)
504 0
    Outputs->vec[i]=0.0;
505 0
  L2_SVM_MFN(Data,Options,Weights,Outputs,0);
506 0
  int num_switches=0;
507 0
  int s=0;
508 0
  int last_round=0;
509 0
  while(lambda_0 <= Options->lambda_u)
510
    {
511 0
      int iter2=0;
512 0
      while(1){
513 0
	s=switch_labels(Data->Y,Outputs->vec,JU,Data->u,Options->S);
514 0
	if(s==0) break;
515 0
	iter2++;
516 0
	if (Options->verbose) { Rcpp::Rcout << "****** lambda_0 = " << lambda_0 << " iteration = " << iter2 << " ************************************"  << endl; }
517 0
	if (Options->verbose) { Rcpp::Rcout << "Optimizing unknown labels. switched " << s << " labels." << endl; }
518 0
	num_switches+=s;
519 0
	if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
520 0
	L2_SVM_MFN(Data,Options,Weights,Outputs,1);
521
      }
522 0
      if(last_round==1) break;
523 0
      lambda_0=TSVM_ANNEALING_RATE*lambda_0;
524 0
      if(lambda_0 >= Options->lambda_u) {lambda_0 = Options->lambda_u; last_round=1;}
525 0
      for(int i=0;i<Data->u;i++)
526 0
	Data->C[JU[i]]=lambda_0*1.0/Data->u;
527 0
      if (Options->verbose) { Rcpp::Rcout << endl << "****** lambda0 increased to " << lambda_0*100/Options->lambda_u << "%" << " of lambda_u = " << Options->lambda_u << " ************************" << endl; }
528 0
      if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
529 0
      L2_SVM_MFN(Data,Options,Weights,Outputs,1);
530
    }
531 0
  if (Options->verbose) { Rcpp::Rcout <<  "Total Number of Switches = " << num_switches << endl; }
532
  /* reset labels */
533 0
  for(int i=0;i<Data->u;i++) Data->Y[JU[i]] = 0.0;
534 0
  double F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
535 0
  if (Options->verbose) { Rcpp::Rcout << "Objective Value = " << F << endl; }
536 0
  delete [] JU;
537 0
  tictoc.stop();
538 0
  if (Options->verbose) { Rcpp::Rcout << "Transductive L2_SVM_MFN took " << tictoc.time() << " seconds." << endl; }
539 0
  return num_switches;
540
}
541 0
int switch_labels(double* Y, double* o, int* JU, int u, int S)
542
{
543 0
  int npos=0;
544 0
  int nneg=0;
545 0
  for(int i=0;i<u;i++)
546
    {
547 0
      if((Y[JU[i]]>0) && (o[JU[i]]<1.0)) npos++;
548 0
      if((Y[JU[i]]<0) && (-o[JU[i]]<1.0)) nneg++;
549
    }
550 0
  Delta* positive=new Delta[npos];
551 0
  Delta* negative=new Delta[nneg];
552 0
  int p=0;
553 0
  int n=0;
554 0
  int ii=0;
555 0
  for(int i=0;i<u;i++)
556
    {
557 0
      ii=JU[i];
558 0
      if((Y[ii]>0.0) && (o[ii]<1.0)) {
559 0
	positive[p].delta=o[ii];
560 0
	positive[p].index=ii;
561 0
	positive[p].s=0;
562 0
	p++;};
563 0
      if((Y[ii]<0.0) && (-o[ii]<1.0))
564
	{
565 0
	  negative[n].delta=-o[ii];
566 0
	  negative[n].index=ii;
567 0
	  negative[n].s=0;
568 0
	  n++;};
569
    }
570 0
  sort(positive,positive+npos);
571 0
  sort(negative,negative+nneg);
572 0
  int s=-1;
573
  //  cout << "Switching Labels: ";
574 0
  while(1)
575
    {
576 0
      s++;
577 0
      if((s>=S) || (positive[s].delta>=-negative[s].delta) || (s>=npos) || (s>=nneg))
578 0
	break;
579 0
      Y[positive[s].index]=-1.0;
580 0
      Y[negative[s].index]= 1.0;
581
      //  cout << positive[s].index << " " << negative[s].index << "...";
582
    }
583
  // cout << "...switched " << s  << " labels." << endl;
584 0
  delete [] positive;
585 0
  delete [] negative;
586 0
  return s;
587
}
588

589 0
int DA_S3VM(struct data *Data,
590
	   struct options *Options,
591
	   struct vector_double *Weights,
592
	   struct vector_double *Outputs)
593
{
594 0
  timer tictoc;
595 0
  tictoc.restart();
596 0
  double T = DA_INIT_TEMP*Options->lambda_u;
597 0
  int iter1 = 0, iter2 =0;
598 0
  double *p = new double[Data->u];
599 0
  double *q = new double[Data->u];
600 0
  double *g = new double[Data->u];
601
  double F,F_min;
602 0
  double *w_min = new double[Data->n];
603 0
  double *o_min = new double[Data->m];
604 0
  double *w = Weights->vec;
605 0
  double *o = Outputs->vec;
606 0
  double kl_divergence = 1.0;
607
  /*initialize */
608 0
  if (Options->verbose) { Rcpp::Rcout << "Initializing weights, p" << endl; }
609 0
  for(int i=0;i<Data->u; i++)
610 0
    p[i] = Options->R;
611
  /* record which examples are unlabeled */
612 0
  int *JU = new int[Data->u];
613 0
  int j=0;
614 0
  for(int i=0;i<Data->m;i++)
615
    {
616 0
      if(Data->Y[i]==0.0)
617 0
	{JU[j]=i;j++;}
618
    }
619 0
  double H = entropy(p,Data->u);
620 0
  optimize_w(Data,p,Options,Weights,Outputs,0);
621 0
  F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
622 0
  F_min = F;
623 0
  for(int i=0;i<Weights->d;i++)
624 0
    w_min[i]=w[i];
625 0
  for(int i=0;i<Outputs->d;i++)
626 0
    o_min[i]=o[i];
627 0
  while((iter1 < DA_OUTER_ITERMAX) && (H > Options->epsilon))
628
    {
629 0
      iter1++;
630 0
      iter2=0;
631 0
      kl_divergence=1.0;
632 0
      while((iter2 < DA_INNER_ITERMAX) && (kl_divergence > Options->epsilon))
633
	{
634 0
	  iter2++;
635 0
	  for(int i=0;i<Data->u;i++)
636
	    {
637 0
	      q[i]=p[i];
638 0
	      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 0
	  if (Options->verbose) { Rcpp::Rcout << "Optimizing p. " ; }
641 0
	  optimize_p(g,Data->u,T,Options->R,p);
642 0
	  kl_divergence=KL(p,q,Data->u);
643 0
	  if (Options->verbose) { Rcpp::Rcout << "Optimizing weights" << endl; }
644 0
	  optimize_w(Data,p,Options,Weights,Outputs,1);
645 0
	  F = transductive_cost(norm_square(Weights),Data->Y,Outputs->vec,Outputs->d,Options->lambda,Options->lambda_u);
646 0
	  if(F < F_min)
647
	    {
648 0
	      F_min = F;
649 0
	      for(int i=0;i<Weights->d;i++)
650 0
		w_min[i]=w[i];
651 0
	      for(int i=0;i<Outputs->d;i++)
652 0
		o_min[i]=o[i];
653
	    }
654 0
	  if (Options->verbose) { Rcpp::Rcout << "***** outer_iter = " << iter1 << "  T = " << T << "  inner_iter = " << iter2 << "  kl = "<< kl_divergence <<"  cost = "<<F<<" *****"<<endl; }
655
	}
656 0
      H = entropy(p,Data->u);
657 0
      if (Options->verbose) { Rcpp::Rcout << "***** Finished outer_iter = " << iter1 << "T = " << T << " Entropy = " << H <<endl; }
658 0
      T = T/DA_ANNEALING_RATE;
659
    }
660 0
  for(int i=0;i<Weights->d;i++)
661 0
    w[i]=w_min[i];
662 0
  for(int i=0;i<Outputs->d;i++)
663 0
    o[i]=o_min[i];
664
  /* may want to reset the original Y */
665 0
  delete [] p;
666 0
  delete [] q;
667 0
  delete [] g;
668 0
  delete [] JU;
669 0
  delete [] w_min;
670 0
  delete [] o_min;
671 0
  tictoc.stop();
672 0
  if (Options->verbose) { Rcpp::Rcout << endl << "(min) Objective Value = " << F_min << endl; }
673 0
  if (Options->verbose) { Rcpp::Rcout << "DA-SVM took " << tictoc.time() << " seconds." << endl; }
674
  return 1;
675
}
676 0
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 0
  timer tictoc;
684 0
  tictoc.restart();
685 0
  double *val = Data->val;
686 0
  int *row = Data->rowptr;
687 0
  int *col = Data->colind;
688 0
  int n  = Data->n;
689 0
  int m  = Data->m;
690 0
  int u  = Data->u;
691 0
  double lambda = Options->lambda;
692
  double epsilon;
693 0
  double *w = Weights->vec;
694 0
  double *o = new double[m+u];
695 0
  double *Y = new double[m+u];
696 0
  double *C = new double[m+u];
697 0
  int *labeled_indices = new int[m];
698 0
  double F_old = 0.0;
699 0
  double F = 0.0;
700 0
  double diff=0.0;
701 0
  double lambda_u_by_u = Options->lambda_u/u;
702 0
  vector_int *ActiveSubset = new vector_int[1];
703 0
  ActiveSubset->vec = new int[m];
704 0
  ActiveSubset->d = m;
705
  // initialize
706 0
  if(ini==0)
707
    {
708 0
      epsilon=BIG_EPSILON;
709 0
      Options->cgitermax=SMALL_CGITERMAX;
710 0
      Options->epsilon=BIG_EPSILON;}
711 0
  else {epsilon = Options->epsilon;}
712

713 0
  for(int i=0;i<n;i++) F+=w[i]*w[i];
714 0
  F=lambda*F;
715 0
  int active=0;
716 0
  int inactive=m-1; // l-1
717 0
  int j = 0;
718
  double temp1;
719
  double temp2;
720 0
  for(int i=0; i<m ; i++)
721
    {
722 0
      o[i]=Outputs->vec[i];
723 0
      if(Data->Y[i]==0.0)
724
	{
725 0
	  labeled_indices[i]=0;
726 0
	  o[m+j]=o[i];
727 0
	  Y[i]=1;
728 0
	  Y[m+j]=-1;
729 0
	  C[i]=lambda_u_by_u*p[j];
730 0
	  C[m+j]=lambda_u_by_u*(1-p[j]);
731 0
	  ActiveSubset->vec[active]=i;
732 0
	  active++;
733 0
	  diff = 1 - fabs(o[i]);
734 0
	  if(diff>0)
735
	    {
736 0
	      Data->Y[i] = 2*p[j]-1;
737 0
	      Data->C[i] = lambda_u_by_u;
738 0
	      temp1 = (1 - o[i]);
739 0
	      temp2 = (1 + o[i]);
740 0
	      F+=lambda_u_by_u*(p[j]*temp1*temp1 + (1-p[j])*temp2*temp2);
741
	    }
742
	  else
743
	    {
744 0
	      if(o[i]>0)
745
		{
746 0
		  Data->Y[i] = -1.0;
747 0
		  Data->C[i] = C[m+j];
748
		}
749
	      else
750
		{
751 0
		  Data->Y[i] = 1.0;
752 0
		  Data->C[i] = C[i];
753
		}
754 0
	      temp1 = (1-Data->Y[i]*o[i]);
755 0
	      F+= Data->C[i]*temp1*temp1;
756
	    }
757 0
	  j++;
758
	}
759
      else
760
	{
761 0
	  labeled_indices[i]=1;
762 0
	  Y[i]=Data->Y[i];
763 0
	  C[i]=1.0/Data->l;
764 0
	  Data->C[i]=1.0/Data->l;
765 0
	  diff=1-Data->Y[i]*o[i];
766 0
	  if(diff>0)
767
	    {
768 0
	      ActiveSubset->vec[active]=i;
769 0
	      active++;
770 0
	      F+=Data->C[i]*diff*diff;
771
	    }
772
	  else
773
	    {
774 0
	      ActiveSubset->vec[inactive]=i;
775 0
	      inactive--;
776
	    }
777
	}
778
    }
779 0
  F=0.5*F;
780 0
  ActiveSubset->d=active;
781 0
  int iter=0;
782 0
  int opt=0;
783 0
  int opt2=0;
784 0
  vector_double *Weights_bar = new vector_double[1];
785 0
  vector_double *Outputs_bar = new vector_double[1];
786 0
  double *w_bar = new double[n];
787 0
  double *o_bar = new double[m+u];
788 0
  Weights_bar->vec=w_bar;
789 0
  Outputs_bar->vec=o_bar;
790 0
  Weights_bar->d=n;
791 0
  Outputs_bar->d=m; /* read only the top m ; bottom u will be copies */
792 0
  double delta=0.0;
793 0
  double t=0.0;
794 0
  int ii = 0;
795 0
  while(iter<MFNITERMAX)
796
    {
797 0
      iter++;
798 0
      if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN Iteration# " << iter << " (" << active << " active examples, " << " objective_value = " << F << ")" << endl; }
799 0
      for(int i=n; i-- ;)
800 0
	w_bar[i]=w[i];
801

802 0
      for(int i=m+u; i-- ;)
803 0
	o_bar[i]=o[i];
804 0
      if (Options->verbose) { Rcpp::Rcout << "_" ; }
805 0
      opt=CGLS(Data,Options,ActiveSubset,Weights_bar,Outputs_bar);
806 0
      for(int i=active; i < m; i++)
807
	{
808 0
	  ii=ActiveSubset->vec[i];
809 0
	  t=0.0;
810 0
	  for(int j=row[ii]; j < row[ii+1]; j++)
811 0
	    t+=val[j]*w_bar[col[j]];
812 0
	  o_bar[ii]=t;
813
	}
814
      // make o_bar consistent in the bottom half
815 0
      int j=0;
816 0
      for(int i=0; i<m;i++)
817
	{
818 0
	  if(labeled_indices[i]==0)
819 0
	    {o_bar[m+j]=o_bar[i]; j++;};
820
	}
821 0
      if(ini==0) {Options->cgitermax=CGITERMAX; ini=1;};
822 0
      opt2=1;
823 0
      for(int i=0; i < m ;i++)
824
	{
825 0
	  ii=ActiveSubset->vec[i];
826 0
	  if(i<active)
827
	    {
828 0
	      if(labeled_indices[ii]==1)
829 0
		opt2=(opt2 && (Data->Y[ii]*o_bar[ii]<=1+epsilon));
830
	      else
831
		{
832 0
		  if(fabs(o[ii])<1)
833 0
		    opt2=(opt2 && (fabs(o_bar[ii])<=1+epsilon));
834
		  else
835 0
		    opt2=(opt2 && (fabs(o_bar[ii])>=1-epsilon));
836
		}
837
	    }
838
	  else
839 0
	    opt2=(opt2 && (Data->Y[ii]*o_bar[ii]>=1-epsilon));
840 0
	  if(opt2==0) break;
841
	}
842 0
      if(opt && opt2) // l
843
	{
844 0
	  if(epsilon==BIG_EPSILON)
845
	    {
846 0
	      epsilon=EPSILON;
847 0
	      Options->epsilon=EPSILON;
848 0
	      if (Options->verbose) { Rcpp::Rcout << "  epsilon = " << BIG_EPSILON << " case converged (speedup heuristic 2). Continuing with epsilon=" <<  EPSILON << endl; }
849 0
	      continue;
850
	    }
851
	  else
852
	    {
853 0
	      for(int i=n; i-- ;)
854 0
		w[i]=w_bar[i];
855 0
	      for(int i=m; i-- ;)
856 0
		Outputs->vec[i]=o_bar[i];
857 0
	      for(int i=m; i-- ;)
858
		{
859 0
		  if(labeled_indices[i]==0)
860 0
		    Data->Y[i]=0.0;
861
		}
862 0
	       delete[] ActiveSubset->vec;
863 0
	       delete[] ActiveSubset;
864 0
	       delete[] o_bar;
865 0
	       delete[] w_bar;
866 0
	       delete[] o;
867 0
	       delete[] Weights_bar;
868 0
	       delete[] Outputs_bar;
869 0
	       delete[] Y;
870 0
	       delete[] C;
871 0
	       delete[] labeled_indices;
872 0
	       tictoc.stop();
873 0
	       if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged in " << iter << " iteration(s) and "<< tictoc.time() << " seconds. \n" << endl; }
874 0
	       return 1;
875
	    }
876
	}
877

878 0
      delta=line_search(w,w_bar,lambda,o,o_bar,Y,C,n,m+u);
879 0
      if (Options->verbose) { Rcpp::Rcout << " LINE_SEARCH delta = " << delta << endl; }
880 0
      F_old=F;
881 0
      F=0.0;
882 0
      for(int i=0;i<n;i++) {w[i]+=delta*(w_bar[i]-w[i]);  F+=w[i]*w[i];}
883 0
      F=lambda*F;
884 0
      j=0;
885 0
      active=0;
886 0
      inactive=m-1;
887 0
      for(int i=0; i<m ; i++)
888
	{
889 0
	   o[i]+=delta*(o_bar[i]-o[i]);
890 0
	   if(labeled_indices[i]==0)
891
	    {
892 0
	      o[m+j]=o[i];
893 0
	      ActiveSubset->vec[active]=i;
894 0
	      active++;
895 0
	      diff = 1 - fabs(o[i]);
896 0
	      if(diff>0)
897
		{
898 0
		  Data->Y[i] = 2*p[j]-1;
899 0
		  Data->C[i] = lambda_u_by_u;
900 0
		  temp1 = (1 - o[i]);
901 0
		  temp2 = (1 + o[i]);
902 0
		  F+=lambda_u_by_u*(p[j]*temp1*temp1 + (1-p[j])*temp2*temp2);
903
		}
904
	      else
905
		{
906 0
		  if(o[i]>0)
907
		    {
908 0
		      Data->Y[i] = -1;
909 0
		      Data->C[i] = C[m+j];
910
		    }
911
		  else
912
		    {
913 0
		      Data->Y[i] = +1;
914 0
		      Data->C[i] = C[i];
915
		    }
916 0
		  temp1=(1-Data->Y[i]*o[i]);
917 0
		  F+= Data->C[i]*temp1*temp1;
918
		}
919 0
	      j++;
920
	    }
921
	  else
922
	    {
923 0
	      diff=1-Data->Y[i]*o[i];
924 0
	      if(diff>0)
925
		{
926 0
		  ActiveSubset->vec[active]=i;
927 0
		  active++;
928 0
		  F+=Data->C[i]*diff*diff;
929
		}
930
	      else
931
		{
932 0
		  ActiveSubset->vec[inactive]=i;
933 0
		  inactive--;
934
		}
935
	    }
936
	}
937 0
      F=0.5*F;
938 0
      ActiveSubset->d=active;
939 0
      if(fabs(F-F_old)<EPSILON)
940 0
	break;
941
    }
942 0
  for(int i=m; i-- ;)
943
    {
944 0
      Outputs->vec[i]=o[i];
945 0
      if(labeled_indices[i]==0)
946 0
	Data->Y[i]=0.0;
947
    }
948 0
  delete[] ActiveSubset->vec;
949 0
  delete[] labeled_indices;
950 0
  delete[] ActiveSubset;
951 0
  delete[] o_bar;
952 0
  delete[] w_bar;
953 0
  delete[] o;
954 0
  delete[] Weights_bar;
955 0
  delete[] Outputs_bar;
956 0
  delete[] Y;
957 0
  delete[] C;
958 0
  tictoc.stop();
959 0
  if (Options->verbose) { Rcpp::Rcout << "L2_SVM_MFN converged in " << iter << " iterations and "<< tictoc.time() << " seconds. \n" << endl;}
960 0
  return 0;
961
}
962 0
void optimize_p(const double* g, int u, double T, double r, double* p)
963
{
964 0
  int iter=0;
965 0
  double epsilon=1e-10;
966 0
  int maxiter=500;
967 0
  double nu_minus=g[0];
968 0
  double nu_plus=g[0];
969 0
  for(int i=0;i<u;i++)
970
    {
971 0
      if(g[i]<nu_minus) nu_minus=g[i];
972 0
      if(g[i]>nu_plus) nu_plus=g[i];
973
    };
974

975 0
  double b=T*log((1-r)/r);
976 0
  nu_minus-=b;
977 0
  nu_plus-=b;
978 0
  double nu=(nu_plus+nu_minus)/2;
979 0
  double Bnu=0.0;
980 0
  double BnuPrime=0.0;
981 0
  double s=0.0;
982 0
  double tmp=0.0;
983 0
  for(int i=0;i<u;i++)
984
    {
985 0
      s=exp((g[i]-nu)/T);
986 0
      if(!(std::isinf(s)))
987
	{
988 0
	  tmp=1.0/(1.0+s);
989 0
	  Bnu+=tmp;
990 0
	  BnuPrime+=s*tmp*tmp;
991
	}
992
    }
993 0
  Bnu=Bnu/u;
994 0
  Bnu-=r;
995 0
  BnuPrime=BnuPrime/(T*u);
996 0
  double nuHat=0.0;
997 0
  while((fabs(Bnu)>epsilon) && (iter < maxiter))
998
    {
999 0
      iter++;
1000 0
      if(fabs(BnuPrime)>0.0)
1001 0
        nuHat=nu-Bnu/BnuPrime;
1002 0
      if((fabs(BnuPrime) > 0.0) | (nuHat>nu_plus)  | (nuHat < nu_minus))
1003 0
        nu=(nu_minus+nu_plus)/2.0;
1004
      else
1005 0
        nu=nuHat;
1006 0
      Bnu=0.0;
1007 0
      BnuPrime=0.0;
1008 0
      for(int i=0;i<u;i++)
1009
	{
1010 0
	  s=exp((g[i]-nu)/T);
1011 0
	  if(!(std::isinf(s)))
1012
	    {
1013 0
	      tmp=1.0/(1.0+s);
1014 0
	      Bnu+=tmp;
1015 0
	      BnuPrime+=s*tmp*tmp;
1016
	    }
1017
	}
1018 0
      Bnu=Bnu/u;
1019 0
      Bnu-=r;
1020 0
      BnuPrime=BnuPrime/(T*u);
1021 0
      if(Bnu<0)
1022 0
        nu_minus=nu;
1023
      else
1024 0
        nu_plus=nu;
1025 0
      if(fabs(nu_minus-nu_plus)<epsilon)
1026 0
	break;
1027
    }
1028 0
  if(fabs(Bnu)>epsilon)
1029 0
    Rcpp::Rcout << "Warning (Root): root not found to required precision" << endl;
1030

1031 0
  for(int i=0;i<u;i++)
1032
    {
1033 0
      s=exp((g[i]-nu)/T);
1034 0
      if(std::isinf(s)) p[i]=0.0;
1035 0
      else p[i]=1.0/(1.0+s);
1036
    }
1037
  //Rcpp::Rcout << " root (nu) = " << nu << " B(nu) = " << Bnu << endl;
1038
}
1039 0
double transductive_cost(double normWeights,double *Y, double *Outputs, int m, double lambda,double lambda_u)
1040
{
1041 0
  double F1=0.0,F2=0.0, o=0.0, y=0.0;
1042 0
  int u=0,l=0;
1043 0
  for(int i=0;i<m;i++)
1044
    {
1045 0
      o=Outputs[i];
1046 0
      y=Y[i];
1047 0
      if(y==0.0)
1048 0
	{F1 += fabs(o) > 1 ? 0 : (1 - fabs(o))*(1 - fabs(o)); u++;}
1049
      else
1050 0
	{F2 += y*o > 1 ? 0 : (1-y*o)*(1-y*o); l++;}
1051
    }
1052
  double F;
1053 0
  F = 0.5*(lambda*normWeights + lambda_u*F1/u + F2/l);
1054 0
  return F;
1055
}
1056

1057 0
double entropy(const double *p, int u)
1058
{
1059 0
  double h=0.0;
1060 0
  double q=0.0;
1061 0
  for(int i=0;i<u;i++)
1062
    {
1063 0
      q=p[i];
1064 0
      if(q>0 && q<1)
1065 0
	h+= -(q*LOG2(q) + (1-q)*LOG2(1-q));
1066
    }
1067 0
  return h/u;
1068
}
1069 0
double KL(const double *p, const double *q, int u)
1070
{
1071 0
  double h=0.0;
1072 0
  double p1=0.0;
1073 0
  double q1=0.0;
1074 0
  double g=0.0;
1075 0
  for(int i=0;i<u;i++)
1076
    {
1077 0
      p1=p[i];
1078 0
      q1=q[i];
1079 0
      if(p1>1-1e-8) p1-=1e-8;
1080 0
      if(p1<1-1e-8) p1+=1e-8;
1081 0
      if(q1>1-1e-8) q1-=1e-8;
1082 0
      if(q1<1-1e-8) q1+=1e-8;
1083 0
      g= (p1*LOG2(p1/q1) + (1-p1)*LOG2((1-p1)/(1-q1)));
1084 0
      if(fabs(g)<1e-12 || std::isnan(g)) g=0.0;
1085 0
      h+=g;
1086
    }
1087 0
    return h/u;
1088
}
1089
/* Prediction */
1090
/* the test file can potentially be huge..dont need to load it.*/
1091 0
void ssl_predict(char *inputs_file_name, const struct vector_double *Weights, struct vector_double *Outputs)
1092
{
1093
  FILE *fpin;
1094 0
  double *w = Weights->vec;
1095 0
  int n = Weights->d;
1096 0
  fpin = fopen(inputs_file_name, "r");
1097 0
  int m=0;
1098 0
  if (fpin == NULL)
1099
    {
1100 0
    Rcpp::stop("Cannot open input file\n");
1101

1102
    }
1103
 /* read and predict */
1104 0
  while (1)
1105
    {
1106 0
      int c = fgetc(fpin);
1107 0
      switch(c)
1108
	{
1109
	case '\n':
1110 0
	  ++m;
1111 0
	  break;
1112
	case EOF:
1113 0
	  goto out;
1114
	default:
1115
	  ;
1116
	}
1117
    }
1118
 out:
1119 0
  initialize(Outputs,m,0.0);
1120 0
  rewind(fpin);
1121
  int colind;
1122
  double val;
1123 0
  double t=0.0;
1124 0
  for(int i=0;i<m;i++)
1125
    {
1126 0
      t=0.0;
1127 0
      while(1)
1128
	{
1129
	  int c;
1130 0
	  do {
1131 0
	    c = getc(fpin);
1132 0
	    if(c=='\n') goto out2;
1133 0
	  } while(isspace(c));
1134 0
	  ungetc(c,fpin);
1135
	  int ret;
1136 0
	  ret = fscanf(fpin,"%d:%lf",&colind,&val);
1137 0
	  if (ret==-1) { Rcpp::Rcout << "EOF" << endl; }
1138 0
	  colind--;
1139 0
	  if(colind<n)
1140 0
	    t+=val*w[colind];
1141
	}
1142
    out2:
1143 0
      Outputs->vec[i]=t + w[n-1];
1144
    }
1145
}
1146
/* Evaluation -- confusion matrix, accuracy, precision, recall, prbep */
1147
/* roc area */
1148 0
void ssl_evaluate(struct vector_double *Outputs, struct vector_double *TrueLabels)
1149
{
1150 0
  double accuracy=0.0;
1151 0
  for(int i=0;i<Outputs->d;i++)
1152 0
    accuracy+=(Outputs->vec[i]*TrueLabels->vec[i])>0;
1153 0
  Rcpp::Rcout << "Accuracy = " << accuracy*100.0/Outputs->d << " %" << endl;
1154
}
1155

1156
/********************** UTILITIES ********************/
1157 0
double norm_square(const vector_double *A)
1158
{
1159 0
  double x=0.0, t=0.0;
1160 0
  for(int i=0;i<A->d;i++)
1161
    {
1162 0
      t=A->vec[i];
1163 0
      x+=t*t;
1164
    }
1165 0
  return x;
1166
}
1167 1
void initialize(struct vector_double *A, int k, double a)
1168
 {
1169 1
  double *vec = new double[k];
1170 1
  for(int i=0;i<k;i++)
1171 1
    vec[i]=a;
1172 1
  A->vec = vec;
1173 1
  A->d   = k;
1174 1
  return;
1175
}
1176 1
void initialize(struct vector_int *A, int k)
1177
{
1178 1
  int *vec = new int[k];
1179 1
  for(int i=0;i<k;i++)
1180 1
    vec[i]=i;
1181 1
  A->vec = vec;
1182 1
  A->d   = k;
1183 1
  return;
1184
}
1185 0
void Write(const char *file_name,
1186
	   const struct vector_double *somevector)
1187
{
1188 0
  FILE* fp = fopen(file_name,"w");
1189 0
  for(int i=0;i<somevector->d;i++)
1190 0
    fprintf(fp,"%g\n",somevector->vec[i]);
1191 0
  return;
1192
}
1193 0
void GetLabeledData(struct data *D, const struct data *Data)
1194
{
1195 0
  int *J = new int[Data->l];
1196 0
  D->C   = new double[Data->l];
1197 0
  D->Y   = new double[Data->l];
1198 0
  int nz=0;
1199 0
  int k=0;
1200 0
  int rowptrs_=Data->l;
1201 0
  for(int i=0;i<Data->m;i++)
1202
    {
1203 0
      if(Data->Y[i]!=0.0)
1204
	{
1205 0
	  J[k]=i;
1206 0
	  D->Y[k]=Data->Y[i];
1207 0
	  D->C[k]=1.0/Data->l;
1208 0
	  nz+=(Data->rowptr[i+1] - Data->rowptr[i]);
1209 0
	  k++;
1210
	}
1211
    }
1212 0
  D->val    = new double[nz];
1213 0
  D->colind = new int[nz];
1214 0
  D->rowptr = new int[rowptrs_+1];
1215 0
  nz=0;
1216 0
  for(int i=0;i<Data->l;i++)
1217
    {
1218 0
      D->rowptr[i]=nz;
1219 0
      for(int j=Data->rowptr[J[i]];j<Data->rowptr[J[i]+1];j++)
1220
	{
1221 0
	  D->val[nz] = Data->val[j];
1222 0
	  D->colind[nz] = Data->colind[j];
1223 0
	  nz++;
1224
	}
1225
    }
1226 0
  D->rowptr[rowptrs_]=nz;
1227 0
  D->nz=nz;
1228 0
  D->l=Data->l;
1229 0
  D->m=Data->l;
1230 0
  D->n=Data->n;
1231 0
  D->u=0;
1232 0
  delete [] J;
1233
}
1234 0
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 0
  a->m=m;
1237 0
  a->u=u;
1238 0
  a->l=m-u;
1239 0
  a->n=n;
1240 0
  a->nz=nz;
1241 0
  a->val=VAL;
1242 0
  a->rowptr=R;
1243 0
  a->colind=C;
1244 0
  a->Y=Y;
1245 0
  a->C=COSTS;
1246 0
  return;
1247
}
1248 0
void Clear(struct data *a)
1249
{
1250 0
  delete [] a->val;
1251 0
  delete [] a->rowptr;
1252 0
  delete [] a->colind;
1253 0
  delete [] a->Y;
1254 0
  delete [] a->C;
1255 0
  delete [] a;
1256 0
  return;
1257
}
1258 0
void Clear(struct vector_double *c)
1259 0
{ delete[] c->vec; delete [] c; return;}
1260 0
void Clear(struct vector_int *c)
1261 0
{ delete[] c->vec; delete [] c; return;}
1262 0
void Clear(struct options *opt)
1263 0
{ delete[] opt; delete [] opt; return;}

Read our documentation on viewing source code .

Loading