开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 多国语言处理 > 查看源码
SVM.cpp
资源名称：TextClassify.rar　[点击查看]
上传用户：sanxfzhen
上传日期：2014-12-28
资源大小：2324k
文件大小：133k
源码类别：

多国语言处理

开发平台：

Visual C++

SVM.cpp：源码内容
  1. // SVM.cpp: implementation of the CSVM class.
  2. //
  3. //////////////////////////////////////////////////////////////////////
  4. #include "stdafx.h"
  5. #include "SVM.h"
  6. #include "Message.h"
  8. # include <stdio.h>
  9. # include <ctype.h>
  10. # include <math.h>
  11. # include <string.h>
  12. # include <stdlib.h>
  13. # include <time.h> 
  14. # include <float.h>
  16. #ifdef _DEBUG
  17. #undef THIS_FILE
  18. static char THIS_FILE[]=__FILE__;
  19. #define new DEBUG_NEW
  20. #endif
  22. extern CCompute_Prompt com_pro;
  23. extern CCompute_Param com_param;
  24. extern CCompute_Result com_result;
  26. # define VERSION       "V3.50"
  27. # define VERSION_DATE  "01.11.00"
  29. # define PRIMAL_OPTIMAL      1
  30. # define DUAL_OPTIMAL        2
  31. # define MAXITER_EXCEEDED    3
  32. # define NAN_SOLUTION        4
  33. # define ONLY_ONE_VARIABLE   5
  35. # define LARGEROUND          0
  36. # define SMALLROUND          1
  38. # define DEF_PRECISION          1E-5
  39. # define DEF_MAX_ITERATIONS     200
  40. # define DEF_LINDEP_SENSITIVITY 1E-8
  41. # define EPSILON_HIDEO          1E-20
  42. # define EPSILON_EQ             1E-5
  45. CSVM theSVM;
  46. //////////////////////////////////////////////////////////////////////
  47. // Construction/Destruction
  48. //////////////////////////////////////////////////////////////////////
  50. CSVM::CSVM()
  51. {
  52. precision_violations=0;
  53. opt_precision=DEF_PRECISION;
  54. maxiter=DEF_MAX_ITERATIONS;
  55. lindep_sensitivity=DEF_LINDEP_SENSITIVITY;
  56. smallroundcount=0;
  57. }
  59. CSVM::~CSVM()
  60. {
  62. }
  64. /***************************svm_learn_main****************************/
  65. void CSVM::set_learn_parameters(LEARN_PARM* learn_parm,KERNEL_PARM* kernel_parm)
  66. {
  67.   learn_parm->biased_hyperplane=com_param.biased_Hyperplane;
  68.   learn_parm->remove_inconsistent=com_param.remove_inconsitant;
  69.   learn_parm->skip_final_opt_check=com_param.final_test;
  70.   learn_parm->svm_maxqpsize=com_param.maximum_size;
  71.   learn_parm->svm_newvarsinqp=com_param.new_variable;
  72.   learn_parm->svm_iter_to_shrink=com_param.iteration_time;
  73.   learn_parm->svm_c=com_param.C;
  74.   learn_parm->transduction_posratio=com_param.fraction;
  75.   learn_parm->svm_costratio=com_param.cost_factor;
  76.   learn_parm->svm_costratio_unlab=1.0;
  77.   learn_parm->svm_unlabbound=1E-5;
  78.   learn_parm->epsilon_crit=0.001;
  79.   learn_parm->epsilon_a=1E-15;
  80.   learn_parm->compute_loo=com_param.loo;
  81.   learn_parm->rho=com_param.rho;
  82.   learn_parm->xa_depth=com_param.search_depth;
  84.   kernel_parm->kernel_type=com_param.kernel_type;
  85.   kernel_parm->poly_degree=com_param.poly_degree;
  86.   kernel_parm->rbf_gamma=com_param.rbf_gamma;
  87.   kernel_parm->coef_lin=com_param.poly_s;
  88.   kernel_parm->coef_const=com_param.poly_c;
  89.   //strcpy(kernel_parm->custom,com_param.);
  90. }
  93. int CSVM::svm_learn_main (int pos_label)
  94. {  
  95.   DOC *docs;  /* training examples */
  96.   long *label,max_docs,max_words_doc;
  97.   long totwords,totdoc,ll,i;
  98.   KERNEL_CACHE kernel_cache;
  99.   LEARN_PARM learn_parm;
  100.   KERNEL_PARM kernel_parm;
  101.   MODEL model;
  102.   char docfile[200];
  103.   char modelfile[200];
  104.   
  105.   if (com_pro.show_action)
  106.   printm("begin to compute");
  107.   if (com_pro.show_action)
  108. printm("Scanning examples...");
  110.   set_learn_parameters(&learn_parm,&kernel_parm);
  112.   strcpy(docfile,com_param.trainfile);
  113.   strcpy(modelfile,com_param.modelfile);
  114. //  kernel_cache_size=com_param.cache_size;
  117.   nol_ll(docfile,&max_docs,&max_words_doc,&ll); /* scan size of input file */
  118.   max_words_doc+=2;
  119.   ll+=2;
  120.   max_docs+=2;
  123.   docs = (DOC *)my_malloc(sizeof(DOC)*max_docs);
  124.   label = (long *)my_malloc(sizeof(long)*max_docs);
  126.   
  127.   read_documents(docfile,docs,label,max_words_doc,ll,&totwords,&totdoc,pos_label);
  129.   if(kernel_parm.kernel_type == LINEAR) /* don't need the cache */
  130.   { 
  131.     svm_learn(docs,label,totdoc,totwords,&learn_parm,&kernel_parm,NULL,&model);
  132.   }
  133.   else 
  134.   {
  135.     kernel_cache_init(&kernel_cache,totdoc,com_param.cache_size);
  136.     svm_learn(docs,label,totdoc,totwords,&learn_parm,&kernel_parm,&kernel_cache,&model);
  137.     kernel_cache_cleanup(&kernel_cache);
  138.   }
  140.   write_model(modelfile,&model);
  142.   free(model.supvec);
  143.   free(model.alpha);
  144.   free(model.index);
  145.   for(i=0;i<totdoc;i++) 
  146.     free(docs[i].words);
  147.   free(docs);
  148.   free(label);
  149.   if (com_pro.show_action)
  150. printm("Cease to compute");
  151.   return(0);
  152. }
  153. /********************************svm_learn_main****************************/
  155. /********************************svm_classify****************************/
  156. double CSVM::svm_classify(DOC &doc)
  157. {
  158. long llsv,max_sv,max_words_sv,pred_format,j;
  159. double sim=0;
  160. MODEL model; 
  161. char modelfile[MAX_PATH];
  163. strcpy(modelfile,com_param.modelfile);
  164. pred_format=0;
  165. // scan size of model file
  166. nol_ll(modelfile,&max_sv,&max_words_sv,&llsv);
  167. max_words_sv+=2;
  168. llsv+=2;
  169. model.supvec = (DOC **)my_malloc(sizeof(DOC *)*max_sv);
  170. model.alpha = (double *)my_malloc(sizeof(double)*max_sv);
  171. read_model(modelfile,&model,max_words_sv,llsv);
  173. doc.twonorm_sq=sprod_ss(doc.words,doc.words);
  174. if(model.kernel_parm.kernel_type == 0) 
  175. {   // linear kernel, compute weight vector
  176. add_weight_vector_to_linear_model(&model);
  177. for(j=0;(doc.words[j]).wnum != 0;j++) 
  178. {  
  179. // Check if feature numbers are not larger than in    
  180. // model. Remove feature if necessary.
  181. if((doc.words[j]).wnum>model.totwords) 
  182. (doc.words[j]).wnum=0;              
  183. }                                        
  184. sim=classify_example_linear(&model,&doc);
  185. }
  186. // non-linear kernel
  187. else sim=classify_example(&model,&doc);
  189. free(model.supvec);
  190. free(model.alpha);
  191. //linear kernel
  192. if(model.kernel_parm.kernel_type == 0) free(model.lin_weights);
  193. return sim;
  194. }
  196. int CSVM::svm_classify(int post_label, double* weight)
  197. {
  198. DOC doc;   /* test example */
  199. long max_docs,max_words_doc,lld,llsv;
  200. long max_sv,max_words_sv,totdoc=0,doc_label;
  201. long wnum,pred_format;
  202. long j;
  203. char *line; 
  204. FILE *docfl;
  205. MODEL model; 
  206. char docfile[MAX_PATH];
  207. char modelfile[MAX_PATH];
  209. strcpy(docfile,com_param.classifyfile);
  210. strcpy(modelfile,com_param.modelfile);
  211. pred_format=0;
  213. nol_ll(docfile,&max_docs,&max_words_doc,&lld); /* scan size of input file */
  214. max_words_doc+=2;
  215. lld+=2;
  216. nol_ll(modelfile,&max_sv,&max_words_sv,&llsv); /* scan size of model file */
  217. max_words_sv+=2;
  218. llsv+=2;
  219. model.supvec = (DOC **)my_malloc(sizeof(DOC *)*max_sv);
  220. model.alpha = (double *)my_malloc(sizeof(double)*max_sv);
  221. line = (char *)my_malloc(sizeof(char)*lld);
  222. doc.words = (SVM_WORD *)my_malloc(sizeof(SVM_WORD)*(max_words_doc+10));
  223. read_model(modelfile,&model,max_words_sv,llsv);
  224. /* linear kernel */
  225. /* compute weight vector */
  226. if(model.kernel_parm.kernel_type == 0)
  227. { 
  228. add_weight_vector_to_linear_model(&model);
  229. }
  230. if (com_pro.show_action)
  231. printm("Classifying test examples..");
  233. if ((docfl = fopen (docfile, "r")) == NULL)
  234. {
  235. printe(docfile);
  236. return -1;
  237. }
  238. //chen 10.9.2001
  240. while((!feof(docfl)) && fgets(line,(int)lld,docfl)) 
  241. {
  242. if(line[0] == '#') continue;  /* line contains comments */
  243. parse_document(line,&doc,&doc_label,&wnum,max_words_doc);
  244. if(model.kernel_parm.kernel_type == 0) 
  245. {   /* linear kernel */
  246. for(j=0;(doc.words[j]).wnum != 0;j++) 
  247. {  /* Check if feature numbers are not larger than in    
  248.  model. Remove feature if necessary.  */
  249. if((doc.words[j]).wnum>model.totwords) 
  250. (doc.words[j]).wnum=0;              
  251. }                                        
  252. weight[totdoc]=classify_example_linear(&model,&doc);
  253. }
  254. else /* non-linear kernel */                        
  255. weight[totdoc]=classify_example(&model,&doc);
  256. totdoc++;
  257. }  
  258. free(line);
  259. free(doc.words);
  260. free(model.supvec);
  261. free(model.alpha);
  262. if(model.kernel_parm.kernel_type == 0) 
  263. { /* linear kernel */
  264. free(model.lin_weights);
  265. }
  266. return(0);
  267. }
  268. /********************************svm_classify****************************/
  271. /********************************svm_common****************************/
  272. double CSVM::classify_example(MODEL *model,DOC *ex)
  273. {
  274. register long i;
  275. register double dist;
  277. dist=0;
  278. for(i=1;i<model->sv_num;i++) {  
  279. dist+=kernel(&model->kernel_parm,model->supvec[i],ex)*model->alpha[i];
  280. }
  281. return(dist-model->b);
  282. }
  284. /*    classifies example for linear kernel 
  285. important: the model must have the linear weight vector computed 
  286. important: the feature numbers in the example to classify must 
  287. not be larger than the weight vector!               */
  288. double CSVM::classify_example_linear(MODEL *model,DOC *ex)
  289. {
  290. return((double)(sprod_ns(model->lin_weights,ex->words)-model->b));
  291. }
  293. /* calculate the kernel function */
  294. CFLOAT CSVM::kernel(KERNEL_PARM *kernel_parm,DOC *a,DOC*b)
  295. {
  296. com_result.kernel_cache_statistic++;
  297. switch(kernel_parm->kernel_type)
  298. {
  299.     case 0: /* linear */ 
  300. return((CFLOAT)sprod_ss(a->words,b->words)); 
  301.     case 1: /* polynomial */
  302. return((CFLOAT)pow(kernel_parm->coef_lin*sprod_ss(a->words,b->words)+kernel_parm->coef_const,(double)kernel_parm->poly_degree)); 
  303.     case 2: /* radial basis function */
  304. return((CFLOAT)exp(-kernel_parm->rbf_gamma*(a->twonorm_sq-2*sprod_ss(a->words,b->words)+b->twonorm_sq)));
  305.     case 3: /* sigmoid neural net */
  306. return((CFLOAT)tanh(kernel_parm->coef_lin*sprod_ss(a->words,b->words)+kernel_parm->coef_const)); 
  307.     case 4: /* custom-kernel supplied in file kernel.h*/
  308. return((CFLOAT)custom_kernel(kernel_parm,a,b)); 
  309. //chen .test sum of 
  310. //return((CFLOAT)pow(kernel_parm->coef_lin*sprod_ss(a->words,b->words)+kernel_parm->coef_const,(double)kernel_parm->poly_degree)+exp(-kernel_parm->rbf_gamma*(a->twonorm_sq-2*sprod_ss(a->words,b->words)+b->twonorm_sq))); 
  311.     default: sprintf(temstr,"Error: Unknown kernel function");
  312. printm(temstr);
  313. return (-1);
  314. }
  315. }
  317. /* compute the inner product of two sparse vectors */
  318. double CSVM::sprod_ss(SVM_WORD *a,SVM_WORD*b)
  319. {
  320.     register FVAL sum=0;
  321.     register SVM_WORD *ai,*bj;
  322.     ai=a;
  323.     bj=b;
  324.     while (ai->wnum && bj->wnum) {
  325. if(ai->wnum > bj->wnum) {
  326. bj++;
  327. }
  328. else if (ai->wnum < bj->wnum) {
  329. ai++;
  330. }
  331. else {
  332. sum+=ai->weight * bj->weight;
  333. ai++;
  334. bj++;
  335. }
  336.     }
  337.     return((double)sum);
  338. }
  342. /* compute the inner product of two sparse vectors,b right shit 1 bit */
  343. double CSVM::sprod_ss1(SVM_WORD *a,SVM_WORD*b,int offset)
  344. {
  345.     register FVAL sum=0;
  346.     register SVM_WORD *ai,*bj;
  347.     ai=a;
  348.     bj=b;
  349.     while (ai->wnum && bj->wnum) {
  350. if(ai->wnum > bj->wnum+offset) {
  351. bj++;
  352. }
  353. else if (ai->wnum < bj->wnum+offset) {
  354. ai++;
  355. }
  356. else 
  357. {
  358. int np=(ai->wnum-1)%16+1+offset;
  359. if (np>0 && np<17) 
  360. sum+=ai->weight * bj->weight;
  361. ai++;
  362. bj++;
  363. }
  364.     }
  365.     return((double)sum);
  366. }
  368. double CSVM::sprod_ss2(SVM_WORD *a,SVM_WORD*b,int offset)
  369. {
  370.     register FVAL sum=0;
  371.     register SVM_WORD *ai,*bj;
  372.     ai=a;
  373.     bj=b;
  374.     while (ai->wnum && bj->wnum) {
  375. if(ai->wnum > bj->wnum+offset) {
  376. bj++;
  377. }
  378. else if (ai->wnum < bj->wnum+offset) {
  379. ai++;
  380. }
  381. else 
  382. {
  383. int np=ai->wnum+offset;
  384. if (np>0 && np<257) 
  385. sum+=ai->weight * bj->weight;
  386. ai++;
  387. bj++;
  388. }
  389.     }
  390.     return((double)sum);
  391. }
  393. /* compute length of weight vector */
  394. double CSVM::model_length_s(MODEL *model,KERNEL_PARM *kernel_parm)
  395. {
  396. register long i,j;
  397. register double sum=0,alphai;
  398. register DOC *supveci;
  400. for(i=1;i<model->sv_num;i++) {  
  401. alphai=model->alpha[i];
  402. supveci=model->supvec[i];
  403. for(j=1;j<model->sv_num;j++) {
  404. sum+=alphai*model->alpha[j]
  405. *kernel(kernel_parm,supveci,model->supvec[j]);
  406. }
  407. }
  408. return(sqrt(sum));
  409. }
  411. void CSVM::clear_vector_n(double *vec,long n)
  412. {
  413. register long i;
  414. for(i=0;i<=n;i++) vec[i]=0;
  415. }
  417. void CSVM::add_vector_ns(double *vec_n,SVM_WORD *vec_s,double faktor)
  418. {
  419. register SVM_WORD *ai;
  420. ai=vec_s;
  421. while (ai->wnum) {
  422. vec_n[ai->wnum]+=(faktor*ai->weight);
  423. ai++;
  424. }
  425. }
  427. double CSVM::sprod_ns(double *vec_n,SVM_WORD *vec_s)
  428. {
  429. register double sum=0;
  430. register SVM_WORD *ai;
  431. ai=vec_s;
  432. while (ai->wnum) {
  433. sum+=(vec_n[ai->wnum]*ai->weight);
  434. ai++;
  435. }
  436. return(sum);
  437. }
  439. /* compute weight vector in linear case and add to model*/
  440. void CSVM::add_weight_vector_to_linear_model(MODEL *model)
  441. {
  442. long i;
  444. model->lin_weights=(double *)my_malloc(sizeof(double)*(model->totwords+1));
  445. clear_vector_n(model->lin_weights,model->totwords);
  446. for(i=1;i<model->sv_num;i++) {
  447. add_vector_ns(model->lin_weights,(model->supvec[i])->words,
  448. model->alpha[i]);
  449. }
  450. }
  452. int CSVM::read_model(char *modelfile,MODEL *model,long max_words,long ll)
  453. {
  454. FILE *modelfl;
  455. long j,i;
  456. char *line;
  457. SVM_WORD *words;
  458. register long wpos;
  459. long wnum,pos;
  460. double weight;
  461. char version_buffer[100];
  463. if (com_pro.show_action)
  464. {
  466. sprintf(temstr,"Reading model..."); 
  467. printm(temstr); 
  468. }
  470. words = (SVM_WORD *)my_malloc(sizeof(SVM_WORD)*(max_words+10));
  471. line = (char *)my_malloc(sizeof(char)*ll);
  473. if ((modelfl = fopen (modelfile, "r")) == NULL)
  474. {
  475. printe (modelfile);  
  476. return -1;
  477. }
  479. fscanf(modelfl,"SVM-light Version %sn",version_buffer);
  481. if(strcmp(version_buffer,VERSION)) 
  482. {
  483. printe ("Version of model-file does not match version of svm_classify!"); 
  484. return -1;
  485. }
  488. fscanf(modelfl,"%ld # kernel typen",&(model->kernel_parm.kernel_type));
  489. fscanf(modelfl,"%ld # kernel parameter -d n",&(model->kernel_parm.poly_degree));
  490. fscanf(modelfl,"%lf # kernel parameter -g n",&(model->kernel_parm.rbf_gamma));
  491. fscanf(modelfl,"%lf # kernel parameter -s n",&(model->kernel_parm.coef_lin));
  492. fscanf(modelfl,"%lf # kernel parameter -r n",&(model->kernel_parm.coef_const));
  493. fscanf(modelfl,"%s # kernel parameter -u n",&(model->kernel_parm.custom));
  494. fscanf(modelfl,"%ld # highest feature index n",&(model->totwords));
  495. fscanf(modelfl,"%ld # number of training documents n",&(model->totdoc));
  496. fscanf(modelfl,"%ld # number of support vectors plus 1 n",&(model->sv_num));
  497. fscanf(modelfl,"%lf # threshold b n",&(model->b));
  499. for(i=1;i<model->sv_num;i++) 
  500. {
  501. fgets(line,(int)ll,modelfl);
  502. pos=0;
  503. wpos=0;
  504. sscanf(line,"%lf",&model->alpha[i]);
  505. while(line[++pos]>' ');
  506. while((sscanf(line+pos,"%ld:%lf",&wnum,&weight) != EOF) && (wpos<max_words)) 
  507. {
  508. while(line[++pos]>' ');
  509. words[wpos].wnum=wnum;
  510. words[wpos].weight=weight; 
  511. wpos++;
  512. } 
  513. model->supvec[i] = (DOC *)my_malloc(sizeof(DOC));
  514. (model->supvec[i])->words = (SVM_WORD *)my_malloc(sizeof(SVM_WORD)*(wpos+1));
  515. for(j=0;j<wpos;j++)
  516. {
  517. (model->supvec[i])->words[j]=words[j]; 
  518. }
  519. ((model->supvec[i])->words[wpos]).wnum=0;
  520. (model->supvec[i])->twonorm_sq = sprod_ss((model->supvec[i])->words,(model->supvec[i])->words);
  521. (model->supvec[i])->docnum = -1;
  522. }
  523. fclose(modelfl);
  524. free(line);
  525. free(words);
  526. if (com_pro.show_readfile)
  527. {
  528. sprintf(temstr, "OK. (%d support vectors read)",(int)(model->sv_num-1));
  529. printm(temstr);
  530. }
  531. }
  532. /*read the data from text documents*/
  533. int CSVM::read_documents(char *docfile,
  534. DOC  *docs,
  535. long *label,
  536. long max_words_doc,
  537. long ll,
  538. long *totwords,
  539. long *totdoc,
  540. int post_label)
  541. {
  542. char *line;
  543. DOC doc;
  544. long dnum=0,wpos,i,dpos=0,dneg=0,dunlab=0;
  545. long doc_label;
  546. FILE *docfl;
  548. line = (char *)my_malloc(sizeof(char)*ll);
  549. if ((docfl = fopen (docfile, "r")) == NULL)
  550. { 
  551. printe (docfile);  
  552. return -1;
  553. }
  555. doc.words = (SVM_WORD *)my_malloc(sizeof(SVM_WORD)*(max_words_doc+10));
  556. if (com_pro.show_readfile)
  557. {
  558. sprintf(temstr,"Reading examples into memory..."); 
  559. printm(temstr);
  560. }
  562. dnum=0;
  563. (*totwords)=0;
  564. while((!feof(docfl)) && fgets(line,(int)ll,docfl)) {
  565. if(line[0] == '#') continue;  /* line contains comments */
  566. if(!parse_document(line,&doc,&doc_label,&wpos,max_words_doc)) 
  567. {
  568. sprintf(temstr,"Parsing error in line %ld!",dnum);
  569. printm(temstr);
  570. }
  572. if(doc_label==0)
  573. {
  574. label[dnum]=0;
  575. dunlab++;
  576. }
  577. else if(doc_label==post_label)
  578. {
  579. label[dnum]=1;
  580. dpos++;
  581. }
  582. else
  583. {
  584. label[dnum]=-1;
  585. dneg++;
  586. }
  588. if((wpos>1) && ((doc.words[wpos-2]).wnum>(*totwords))) 
  589. (*totwords)=(doc.words[wpos-2]).wnum;
  590. docs[dnum].words = (SVM_WORD *)my_malloc(sizeof(SVM_WORD)*wpos);
  591. docs[dnum].docnum=dnum;
  592. for(i=0;i<wpos;i++) 
  593. docs[dnum].words[i]=doc.words[i];
  594. docs[dnum].twonorm_sq=doc.twonorm_sq;
  595. dnum++;  
  596. if((dnum % 100) == 0&&com_pro.show_readfile) 
  597. {
  598. sprintf(temstr,"read %ld..",dnum); 
  599. printm(temstr);
  600. }
  601. } 
  603. fclose(docfl);
  604. free(line);
  605. free(doc.words);
  606. if (com_pro.show_action)
  607. {
  608. sprintf(temstr, "OK. (%ld examples read)", dnum);
  609. printm(temstr);
  610. sprintf(temstr,"%ld positive, %ld negative, and %ld unlabeled examples.",dpos,dneg,dunlab); 
  611. printm(temstr);
  612. }
  613. (*totdoc)=dnum;
  614. }
  615. /*Parse one  line of data file */
  616. int CSVM::parse_document(char *line,DOC *doc,long *label,long*numwords,long max_words_doc)
  617. {
  618. register long wpos,pos;
  619. long wnum;
  620. double weight;
  622. pos=0;
  623. while(line[pos]) 
  624. {      /* cut off comments */
  625. if(line[pos] == '#') 
  626. {
  627. line[pos]=0;
  628. }
  629. else 
  630. {
  631. pos++;
  632. }
  633. }
  634. wpos=0;
  635. if((sscanf(line,"%ld",label)) == EOF) return(0);
  636. pos=0;
  637. while(line[pos]==' ') pos++;
  638. while(line[pos]>' ') pos++;
  639. while((sscanf(line+pos,"%ld:%lf",&wnum,&weight) != EOF) &&  (wpos<max_words_doc))
  640. {
  641. while(line[pos++]==' ');
  642. while(line[++pos]>' ');
  643. if(wnum<=0) 
  644. { 
  645. printe ("Feature numbers must be larger or equal to 1!!!"); 
  646. sprintf(temstr,"LINE: %s",line);
  647. printm(temstr);
  648. return (0);
  649.  
  650. }
  651. if((wpos>0) && ((doc->words[wpos-1]).wnum >= wnum))
  652. { 
  653. printe ("Features must be in increasing order!!!"); 
  654. sprintf(temstr,"LINE: %s",line);
  655. printm(temstr);
  656. return (0);
  657.  
  658. }
  659. (doc->words[wpos]).wnum=wnum;
  660. (doc->words[wpos]).weight=weight; 
  661. wpos++;
  662. }
  663. (doc->words[wpos]).wnum=0;
  664. (*numwords)=wpos+1;
  665. doc->docnum=-1;
  666. doc->twonorm_sq=sprod_ss(doc->words,doc->words);
  667. return(1);
  668. }
  669. /* grep through file and count number of lines, 
  670. maximum number of spaces per line, and 
  671. longest line. */
  672. int CSVM::nol_ll(char *file,long *nol,long *wol,long *ll) 
  673. {
  674. FILE *fl;
  675. int ic;
  676. char c;
  677. long current_length,current_wol;
  679. if ((fl = fopen (file, "r")) == NULL)
  680. { 
  681. printe (file);   
  682. return -1;
  683. }
  684. current_length=0;
  685. current_wol=0;
  686. (*ll)=0;
  687. (*nol)=1;
  688. (*wol)=0;
  689. while((ic=getc(fl)) != EOF) 
  690. {
  691. c=(char)ic;
  692. current_length++;
  693. if(c == ' ') 
  694. {
  695. current_wol++;
  696. }
  697. if(c == 'n') 
  698. {
  699. (*nol)++;
  700. if(current_length>(*ll)) 
  701. {
  702. (*ll)=current_length;
  703. }
  704. if(current_wol>(*wol)) 
  705. {
  706. (*wol)=current_wol;
  707. }
  708. current_length=0;
  709. current_wol=0;
  710. }
  711. }
  712. fclose(fl);
  713. }
  715. long CSVM::minl(long a,long b)
  716. {
  717. if(a<b)
  718. return(a);
  719. else
  720. return(b);
  721. }
  723. long CSVM::maxl(long a,long b)
  724. {
  725. if(a>b)
  726. return(a);
  727. else
  728. return(b);
  729. }
  731. long CSVM::get_runtime() 
  732. {
  733. clock_t start;
  734. start = clock();
  735. return((long)((double)start*100.0/(double)CLOCKS_PER_SEC));
  736. }
  739. int CSVM::isnan(double a)
  740. {
  741. return(_isnan(a));
  742. }
  744. void * CSVM::my_malloc(long size) 
  745. {
  746. void *ptr;
  747. ptr=(void *)malloc(size);
  748. if(!ptr)
  749. { 
  750. printe ("Out of memory!"); 
  751. return (NULL);
  752. }
  753. return(ptr);
  754. }
  755. //print error on screen
  757. void CSVM::printe(char* pInfo)
  758. {
  759. CMessage::PrintError(pInfo);
  760. }
  762. //print message on screen
  763. void CSVM::printm(char* pInfo)
  764. {
  765. CMessage::PrintInfo(pInfo);
  766. }
  767. //custom kernel function
  768. /////////////////////////////////////////////////////////////////
  769. //chen 2001.09.14
  770. /////////////////////////////////////////////////////////////////
  772. double CSVM::ktt(int ta,int tb,double pa[],double pb[])
  773. {
  774. int ya,yb;
  775. ya=((ta-1)%13)+1;
  776. yb=((tb-1)%13)+1;
  778. if (ya<13&&yb<13) 
  779. return pa[ta]*pa[ta+1]*pb[tb]*pb[tb+1];
  780. else return 0.0;
  782. }
  784. double CSVM::kt(int t,double ta[],double tb[])
  785. {
  787. int x,y;
  788. double sum=0.0;
  789. x=((t-1)/16)+1;
  790. y=((t-1)%16)+1;
  791. if (x>1) 
  792. sum+=ta[t]*tb[t]*ta[t-16]*tb[t-16];
  793. if (x>2) 
  794. sum+=ta[t]*tb[t]*ta[t-32]*tb[t-32];
  795. if (x>3) 
  796. sum+=ta[t]*tb[t]*ta[t-48]*tb[t-48];
  797. if (x<14) 
  798. sum+=ta[t]*tb[t]*ta[t-48]*tb[t-48];
  799. if (x<15) 
  800. sum+=ta[t]*tb[t]*ta[t-32]*tb[t-32];
  801. if (x<16)
  802. sum+=ta[t]*tb[t]*ta[t+16]*tb[t+16];
  803. if (y>3) 
  804. sum+=ta[t]*tb[t]*ta[t-3]*tb[t-3];
  805. if (y>2) 
  806. sum+=ta[t]*tb[t]*ta[t-2]*tb[t-2];
  807. if (y>1) 
  808. sum+=ta[t]*tb[t]*ta[t-1]*tb[t-1];
  809. if (y<14) 
  810. sum+=ta[t]*tb[t]*ta[t+3]*tb[t+3];
  811. if (y<15) 
  812. sum+=ta[t]*tb[t]*ta[t+2]*tb[t+2];
  813. if (y<16) 
  814. sum+=ta[t]*tb[t]*ta[t+1]*tb[t+1];
  815. return sum;
  816. }
  818. double CSVM::fi(double* tt)
  819. {
  820. int x,y;
  821. double sum=0.0;
  822. for (int t=1;t<=52;t++)
  823. {
  824. x=((t-1)/13)+1;
  825. y=((t-1)%13)+1;
  826. if (y<13)
  827. sum+=tt[t]*tt[t+1];
  828. }
  829. return sum;
  830. }
  831. double CSVM::fs(double ta[])
  832. {
  833. double sum=0.0;
  834. int x,y;
  835. for (int i=1;i<256;i++)
  836. {
  837. x=((i-1)/16)+1;
  838. y=((i-1)%16)+1;
  840. if (x<16)
  841. sum+=ta[i]*ta[i+16];
  842. if (y<16) 
  843. sum+=ta[i]*ta[i+1];
  844. }
  845. return sum;
  846. }
  847. double CSVM::sumofword(DOC* a)
  848. {
  849. double sum=0.0;
  850. SVM_WORD* pwa=a->words;
  851. while (pwa->wnum)
  852. {
  853. sum+=pwa->weight;
  854. pwa++;
  855. }
  856. return sum;
  860. }
  861. double CSVM::custom_kernel(KERNEL_PARM *kernel_parm,DOC *a,DOC*b)
  862. {
  863. double sum=0;
  864.     SVM_WORD *ai,*bj;
  865.     ai=a->words;
  866.     bj=b->words;
  867.     while (ai->wnum || bj->wnum) 
  868. {
  869. if(ai->wnum == bj->wnum) 
  870. {
  871. sum+=(fabs(ai->weight-bj->weight))*(fabs(ai->weight-bj->weight));
  872. ai++;bj++;
  873. }
  874. else if ((ai!=0) &&(ai->wnum<bj->wnum || bj->wnum==0))
  875. {
  876. sum+=fabs(ai->weight)*fabs(ai->weight);
  877. ai++;
  878. }
  879. else if ((bj!=0) &&(bj->wnum<ai->wnum|| ai->wnum==0))
  880. {
  881. sum+=fabs(bj->weight)*fabs(bj->weight);
  882. bj++;
  883. }
  884.     }
  885. //   case 1: /* polynomial *///
  886. //return((CFLOAT)pow(kernel_parm->coef_lin*sprod_ss(a->words,b->words)+kernel_parm->coef_const,(double)kernel_parm->poly_degree)); 
  887.     //case 2: /* radial basis function */
  888. // return((CFLOAT)exp(-kernel_parm->rbf_gamma*(a->twonorm_sq-2*sprod_ss(a->words,b->words)+b->twonorm_sq)));
  889.     //case 3: /* sigmoid neural net */
  890. // return((CFLOAT)tanh(kernel_parm->coef_lin*sprod_ss(a->words,b->words)+kernel_parm->coef_const)); 
  891.     //case 4: /* custom-kernel supplied in file kernel.h*/
  892. // return((CFLOAT)custom_kernel(kernel_parm,a,b)); 
  894. /*
  895.     SVM_WORD *ai,*bj;
  896.     ai=a->words;
  897.     bj=b->words;
  898. double suma=0.0,sumb=0.0;
  900.     while (ai->wnum ) 
  901. {
  902. suma+=ai->weight;
  903. ai++;
  904. }
  905.     while (bj->wnum ) 
  906. {
  907. sumb+=bj->weight;
  908. bj++;
  909. }
  910. */
  911. // double K_rbf=exp(-0.001*(a->twonorm_sq-2*sprod_ss(a->words,b->words)+b->twonorm_sq));
  912. double K_Laplace=exp(-0.0001*sum);
  913. double K_poly=pow(sprod_ss(a->words,b->words)+20,3);
  914. // double K_linear=sprod_ss(a->words,b->words);
  915. //double sum;
  916. //double sum=suma*sumb;
  917. //sum=K_rbf*K_poly;
  918. //double sum=fabs(pro*pro+pro-tan(pro));
  919. return K_Laplace*K_poly;
  921. }
  922. /********************************svm_common****************************/
  925. /********************************svm_hideo****************************/
  926. double *CSVM::optimize_qp(QP *qp,
  927. double *epsilon_crit,
  928. long nx, /* Maximum number of variables in QP */
  929. double *threshold, 
  930. LEARN_PARM *learn_parm)
  931. /* start the optimizer and return the optimal values */
  932. /* The HIDEO optimizer does not necessarily fully solve the problem. */
  933. /* Since it requires a strictly positive definite hessian, the solution */
  934. /* is restricted to a linear independent subset in case the matrix is */
  935. /* only semi-definite. */
  936. {
  937.   long i;
  938.   int result;
  939.   double eq;
  941.   if(!primal) { /* allocate memory at first call */
  942.     primal=(double *)my_malloc(sizeof(double)*nx);
  943.     dual=(double *)my_malloc(sizeof(double)*((nx+1)*2));
  944.     nonoptimal=(long *)my_malloc(sizeof(long)*(nx));
  945.     buffer=(double *)my_malloc(sizeof(double)*((nx+1)*2*(nx+1)*2+
  946.        nx*nx+2*(nx+1)*2+2*nx+1+2*nx+
  947.        nx+nx+nx*nx));
  948.     (*threshold)=0;
  949.   }
  951.     eq=qp->opt_ce0[0];
  952.     for(i=0;i<qp->opt_n;i++) {
  953.       eq+=qp->opt_xinit[i]*qp->opt_ce[i];
  954.     }
  955.     
  956.   
  958.   result=optimize_hildreth_despo(qp->opt_n,qp->opt_m,
  959.  opt_precision,(*epsilon_crit),
  960.  learn_parm->epsilon_a,maxiter,
  961.  /* (long)PRIMAL_OPTIMAL, */
  962.  (long)0, (long)0,
  963.  lindep_sensitivity,
  964.  qp->opt_g,qp->opt_g0,qp->opt_ce,qp->opt_ce0,
  965.  qp->opt_low,qp->opt_up,primal,qp->opt_xinit,
  966.  dual,nonoptimal,buffer);
  968.   if(learn_parm->totwords < learn_parm->svm_maxqpsize) { 
  969.     /* larger working sets will be linear dependent anyway */
  970.     learn_parm->svm_maxqpsize=maxl(learn_parm->totwords,(long)2);
  971.   }
  973.   if(result == NAN_SOLUTION) {
  974.     lindep_sensitivity*=2;  /* throw out linear dependent examples more */
  975.                             /* generously */
  976.     if(learn_parm->svm_maxqpsize>2) {
  977.       learn_parm->svm_maxqpsize--;  /* decrease size of qp-subproblems */
  978.     }
  979.     precision_violations++;
  980.   }
  982.   /* take one round of only two variable to get unstuck */
  983.   if(result != PRIMAL_OPTIMAL) {
  985.     smallroundcount++;
  987.     result=optimize_hildreth_despo(qp->opt_n,qp->opt_m,
  988.    opt_precision,(*epsilon_crit),
  989.    learn_parm->epsilon_a,(long)maxiter,
  990.    (long)PRIMAL_OPTIMAL,(long)SMALLROUND,
  991.    lindep_sensitivity,
  992.    qp->opt_g,qp->opt_g0,qp->opt_ce,qp->opt_ce0,
  993.    qp->opt_low,qp->opt_up,primal,qp->opt_xinit,
  994.    dual,nonoptimal,buffer);
  995.     }
  997.     if(result != PRIMAL_OPTIMAL) {
  998.       if(result != ONLY_ONE_VARIABLE) 
  999. precision_violations++;
  1000.       if(result == MAXITER_EXCEEDED) 
  1001. maxiter+=100;
  1002.       if(result == NAN_SOLUTION) {
  1003. lindep_sensitivity*=2;  /* throw out linear dependent examples more */
  1004.                         /* generously */
  1005. /* results not valid, so return inital values */
  1006. for(i=0;i<qp->opt_n;i++) {
  1007.   primal[i]=qp->opt_xinit[i];
  1008. }
  1009.       }
  1010.     }
  1013.   if(precision_violations > 50) {
  1014.     precision_violations=0;
  1015.     (*epsilon_crit)*=10.0; 
  1016.     }
  1019.   if((qp->opt_m>0) && (result != NAN_SOLUTION))
  1020.     (*threshold)=dual[1]-dual[0];
  1021.   else
  1022.     (*threshold)=0;
  1024.   return(primal);
  1025. }
  1029. int CSVM::optimize_hildreth_despo(
  1030.      long   n,            /* number of variables */
  1031.      long   m,            /* number of linear equality constraints [0,1] */
  1032.      double precision,    /* solve at least to this dual precision */
  1033.      double epsilon_crit, /* stop, if KT-Conditions approx fulfilled */
  1034.      double epsilon_a,    /* precision of alphas at bounds */
  1035.      long   maxiter,      /* stop after this many iterations */
  1036.      long   goal,         /* keep going until goal fulfilled */
  1037.      long   smallround,   /* use only two variables of steepest descent */
  1038.      double lindep_sensitivity, /* epsilon for detecting linear dependent ex */
  1039.      double *g,           /* hessian of objective */
  1040.      double *g0,          /* linear part of objective */
  1041.      double *ce,double *ce0,     /* linear equality constraints */
  1042.      double *low,double *up,     /* box constraints */
  1043.      double *primal,      /* primal variables */
  1044.      double *init,        /* initial values of primal */
  1045.      double *dual,        /* dual variables */
  1046.      long   *lin_dependent,
  1047.      double *buffer)
  1048. {
  1049.   long i,j,k,from,to,n_indep,changed;
  1050.   double sum,bmin=0,bmax=0;
  1051.   double *d,*d0,*ig,*dual_old,*temp,*start;       
  1052.   double *g0_new,*g_new,*ce_new,*ce0_new,*low_new,*up_new;
  1053.   double add,t;
  1054.   int result;
  1055.   //double obj_before,obj_after; 
  1056.   long b1,b2;
  1058.   g0_new=&(buffer[0]);    /* claim regions of buffer */
  1059.   d=&(buffer[n]);
  1060.   d0=&(buffer[n+(n+m)*2*(n+m)*2]);
  1061.   ce_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2]);
  1062.   ce0_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n]);
  1063.   ig=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m]);
  1064.   dual_old=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n]);
  1065.   low_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2]);
  1066.   up_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n]);
  1067.   start=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n]);
  1068.   g_new=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n+n]);
  1069.   temp=&(buffer[n+(n+m)*2*(n+m)*2+(n+m)*2+n+m+n*n+(n+m)*2+n+n+n+n*n]);
  1071.   b1=-1;
  1072.   b2=-1;
  1073.   for(i=0;i<n;i++) {   /* get variables with steepest feasible descent */
  1074.     sum=g0[i];         
  1075.     for(j=0;j<n;j++) 
  1076.       sum+=init[j]*g[i*n+j];
  1077.     sum=sum*ce[i];
  1078.     if(((b1==-1) || (sum<bmin)) 
  1079.        && (!((init[i]<=(low[i]+epsilon_a)) && (ce[i]<0.0)))
  1080.        && (!((init[i]>=( up[i]-epsilon_a)) && (ce[i]>0.0)))
  1081.        ) {
  1082.       bmin=sum;
  1083.       b1=i;
  1084.     }
  1085.     if(((b2==-1) || (sum>bmax)) 
  1086.        && (!((init[i]<=(low[i]+epsilon_a)) && (ce[i]>0.0)))
  1087.        && (!((init[i]>=( up[i]-epsilon_a)) && (ce[i]<0.0)))
  1088.        ) {
  1089.       bmax=sum;
  1090.       b2=i;
  1091.     }
  1092.   }
  1093.   
  1094.   for(i=0;i<n;i++) {
  1095.     start[i]=init[i];
  1096.   }
  1098.   /* in case both examples are equal */
  1099.   add=0;
  1100.   changed=0;
  1101.   if((-g[b1*n+b2] == g[b1*n+b1]) && (-g[b1*n+b2] == g[b2*n+b2])) {
  1102.     changed=1;
  1103.     if(ce[b1] == ce[b2]) { /* distribute evenly */
  1104.       start[b1]=(init[b1]+init[b2])/2.0;
  1105.       start[b2]=(init[b1]+init[b2])/2.0;
  1106.       if(start[b2] > up[b2]) {
  1107. t=start[b2]-up[b2];
  1108. start[b2]=up[b2];
  1109. start[b1]+=t;
  1110.       }
  1111.       if(start[b1] > up[b1]) {
  1112. t=start[b1]-up[b1];
  1113. start[b1]=up[b1];
  1114. start[b2]+=t;
  1115.       }
  1116.     }
  1117.     else { /* set to upper bound */
  1118.       t=up[b1]-init[b1];
  1119.       if((up[b2]-init[b2]) < t) {
  1120. t=up[b2]-init[b2];
  1121.       }
  1122.       start[b1]=init[b1]+t;
  1123.       start[b2]=init[b2]+t;
  1124.     }
  1125.   }
  1126.   /* if we have a biased hyperplane, then adding a constant to the */
  1127.   /* hessian does not change the solution. So that is done for examples */
  1128.   /* with zero diagonal entry, since HIDEO cannot handle them. */
  1129.   else if((fabs(g[b1*n+b1]) < lindep_sensitivity) 
  1130.   || (fabs(g[b2*n+b2]) < lindep_sensitivity)) {
  1131.     add+=0.093274;
  1132.   }    
  1133.   /* in case both examples are linear dependent */
  1134.   else if(fabs(g[b1*n+b1]/g[b1*n+b2] - g[b1*n+b2]/g[b2*n+b2])
  1135.   < lindep_sensitivity) { 
  1136.     add+=0.078274;
  1137.   }
  1140. /* sprintf(temstr,"b1=%ld,b2=%ldn",b1,b2); */
  1142.   lcopy_matrix(g,n,d);
  1143.   if((m==1) && (add>0.0)) {
  1144.     for(j=0;j<n;j++) {
  1145.       for(k=0;k<n;k++) {
  1146. d[j*n+k]+=add*ce[j]*ce[k];
  1147.       }
  1148.     }
  1149.   }
  1150.   else {
  1151.     add=0.0;
  1152.   }
  1154.   if(n>2) {                    /* switch, so that variables are better mixed */
  1155.     lswitchrk_matrix(d,n,b1,(long)0); 
  1156.     if(b2 == 0) 
  1157.       lswitchrk_matrix(d,n,b1,(long)1); 
  1158.     else
  1159.       lswitchrk_matrix(d,n,b2,(long)1); 
  1160.   }
  1161.   if(smallround == SMALLROUND) {
  1162.     for(i=2;i<n;i++) {
  1163.       lin_dependent[i]=1;
  1164.     }
  1165.     lin_dependent[0]=0;
  1166.     lin_dependent[1]=0;
  1167.   }
  1168.   else {
  1169.     for(i=0;i<n;i++) {
  1170.       lin_dependent[i]=0;
  1171.     }
  1172.   }
  1173.   linvert_matrix(d,n,ig,lindep_sensitivity,lin_dependent);
  1174.   if(n>2) {                    /* now switch back */
  1175.     if(b2 == 0) {
  1176.       lswitchrk_matrix(ig,n,b1,(long)1); 
  1177.       i=lin_dependent[1];  
  1178.       lin_dependent[1]=lin_dependent[b1];
  1179.       lin_dependent[b1]=i;
  1180.     }
  1181.     else {
  1182.       lswitchrk_matrix(ig,n,b2,(long)1); 
  1183.       i=lin_dependent[1];  
  1184.       lin_dependent[1]=lin_dependent[b2];
  1185.       lin_dependent[b2]=i;
  1186.     }
  1187.     lswitchrk_matrix(ig,n,b1,(long)0); 
  1188.     i=lin_dependent[0];  
  1189.     lin_dependent[0]=lin_dependent[b1];
  1190.     lin_dependent[b1]=i;
  1191.   }
  1192.   /* lprint_matrix(d,n); */
  1193.   /* lprint_matrix(ig,n); */
  1195.   lcopy_matrix(g,n,g_new);   /* restore g_new matrix */
  1196.   if(add>0)
  1197.     for(j=0;j<n;j++) {
  1198.       for(k=0;k<n;k++) {
  1199. g_new[j*n+k]+=add*ce[j]*ce[k];
  1200.       }
  1201.     }
  1203.   for(i=0;i<n;i++) {  /* fix linear dependent vectors */
  1204.     g0_new[i]=g0[i]+add*ce0[0]*ce[i];
  1205.   }
  1206.   if(m>0) ce0_new[0]=-ce0[0];
  1207.   for(i=0;i<n;i++) {  /* fix linear dependent vectors */
  1208.     if(lin_dependent[i]) {
  1209.       for(j=0;j<n;j++) {
  1210. if(!lin_dependent[j]) {
  1211.   g0_new[j]+=start[i]*g_new[i*n+j];
  1212. }
  1213.       }
  1214.       if(m>0) ce0_new[0]-=(start[i]*ce[i]);
  1215.     }
  1216.   }
  1217.   from=0;   /* remove linear dependent vectors */
  1218.   to=0;
  1219.   n_indep=0;
  1220.   for(i=0;i<n;i++) {
  1221.     if(!lin_dependent[i]) {
  1222.       g0_new[n_indep]=g0_new[i];
  1223.       ce_new[n_indep]=ce[i]; 
  1224.       low_new[n_indep]=low[i];
  1225.       up_new[n_indep]=up[i];
  1226.       primal[n_indep]=start[i];
  1227.       n_indep++;
  1228.     }
  1229.     for(j=0;j<n;j++) {
  1230.       if((!lin_dependent[i]) && (!lin_dependent[j])) {
  1231.         ig[to]=ig[from];
  1232.         g_new[to]=g_new[from];
  1233. to++;
  1234.       }
  1235.       from++;
  1236.     }
  1237.   }
  1239.   
  1240.   /* cannot optimize with only one variable */
  1241.   if((n_indep<=1) && (m>0) && (!changed)) { 
  1242.     for(i=n-1;i>=0;i--) {
  1243.       primal[i]=init[i];
  1244.     }
  1245.     return((int)ONLY_ONE_VARIABLE);
  1246.   }
  1248.   result=solve_dual(n_indep,m,precision,epsilon_crit,maxiter,g_new,g0_new,
  1249.     ce_new,ce0_new,low_new,up_new,primal,d,d0,ig,
  1250.     dual,dual_old,temp,goal);
  1251.   
  1252.   j=n_indep;
  1253.   for(i=n-1;i>=0;i--) {
  1254.     if(!lin_dependent[i]) {
  1255.       j--;
  1256.       primal[i]=primal[j];
  1257.     }
  1258.     else if((m==0) && (g[i*n+i]==0)) {
  1259.       /* if we use a biased hyperplane, each example with a zero diagonal */
  1260.       /* entry must have an alpha at the upper bound. Doing this */
  1261.       /* is essential for the HIDEO optimizer, since it cannot handle zero */
  1262.       /* diagonal entries in the hessian for the unbiased hyperplane case. */
  1263.       primal[i]=up[i];  
  1264.     }
  1265.     else {
  1266.       primal[i]=start[i];  /* leave as is */
  1267.     }
  1268.     temp[i]=primal[i];
  1269.   }
  1270.    
  1271.  
  1272.   return((int)result);
  1273. }
  1276. int CSVM::solve_dual(
  1277.      /* Solves the dual using the method of Hildreth and D'Espo. */
  1278.      /* Can only handle problems with zero or exactly one */
  1279.      /* equality constraints. */
  1281.      long   n,            /* number of variables */
  1282.      long   m,            /* number of linear equality constraints */
  1283.      double precision,    /* solve at least to this dual precision */
  1284.      double epsilon_crit, /* stop, if KT-Conditions approx fulfilled */
  1285.      long   maxiter,      /* stop after that many iterations */
  1286.      double *g,
  1287.      double *g0,          /* linear part of objective */
  1288.      double *ce,double *ce0,     /* linear equality constraints */
  1289.      double *low,double *up,     /* box constraints */
  1290.      double *primal,      /* variables (with initial values) */
  1291.      double *d,double *d0,double *ig,double *dual,double *dual_old,double *temp,       /* buffer  */
  1292.      long goal)
  1293. {
  1294.   long i,j,k,iter;
  1295.   double sum,w,maxviol,viol,temp1,temp2,isnantest;
  1296.   double model_b,dist;
  1297.   long retrain,maxfaktor,primal_optimal=0,at_bound,scalemaxiter;
  1298.   double epsilon_a=1E-15,epsilon_hideo;
  1299.   double eq; 
  1301.   if((m<0) || (m>1)) 
  1302.     printe("SOLVE DUAL: inappropriate number of eq-constrains!");
  1305.   for(i=0;i<2*(n+m);i++) {
  1306.     dual[i]=0;
  1307.     dual_old[i]=0;
  1308.   }
  1309.   for(i=0;i<n;i++) {   
  1310.     for(j=0;j<n;j++) {   /* dual hessian for box constraints */
  1311.       d[i*2*(n+m)+j]=ig[i*n+j];
  1312.       d[(i+n)*2*(n+m)+j]=-ig[i*n+j];
  1313.       d[i*2*(n+m)+j+n]=-ig[i*n+j];
  1314.       d[(i+n)*2*(n+m)+j+n]=ig[i*n+j];
  1315.     }
  1316.     if(m>0) {
  1317.       sum=0;              /* dual hessian for eq constraints */
  1318.       for(j=0;j<n;j++) {
  1319. sum+=(ce[j]*ig[i*n+j]);
  1320.       }
  1321.       d[i*2*(n+m)+2*n]=sum;
  1322.       d[i*2*(n+m)+2*n+1]=-sum;
  1323.       d[(n+i)*2*(n+m)+2*n]=-sum;
  1324.       d[(n+i)*2*(n+m)+2*n+1]=sum;
  1325.       d[(n+n)*2*(n+m)+i]=sum;
  1326.       d[(n+n+1)*2*(n+m)+i]=-sum;
  1327.       d[(n+n)*2*(n+m)+(n+i)]=-sum;
  1328.       d[(n+n+1)*2*(n+m)+(n+i)]=sum;
  1329.       
  1330.       sum=0;
  1331.       for(j=0;j<n;j++) {
  1332. for(k=0;k<n;k++) {
  1333.   sum+=(ce[k]*ce[j]*ig[j*n+k]);
  1334. }
  1335.       }
  1336.       d[(n+n)*2*(n+m)+2*n]=sum;
  1337.       d[(n+n)*2*(n+m)+2*n+1]=-sum;
  1338.       d[(n+n+1)*2*(n+m)+2*n]=-sum;
  1339.       d[(n+n+1)*2*(n+m)+2*n+1]=sum;
  1340.     } 
  1341.   }
  1343.   for(i=0;i<n;i++) {   /* dual linear component for the box constraints */
  1344.     w=0;
  1345.     for(j=0;j<n;j++) {
  1346.       w+=(ig[i*n+j]*g0[j]); 
  1347.     }
  1348.     d0[i]=up[i]+w;
  1349.     d0[i+n]=-low[i]-w;
  1350.   }
  1352.   if(m>0) {  
  1353.     sum=0;             /* dual linear component for eq constraints */
  1354.     for(j=0;j<n;j++) {
  1355.       for(k=0;k<n;k++) {
  1356. sum+=(ce[k]*ig[k*n+j]*g0[j]); 
  1357.       }
  1358.     }
  1359.     d0[2*n]=ce0[0]+sum;
  1360.     d0[2*n+1]=-ce0[0]-sum;
  1361.   }
  1363.   maxviol=999999;
  1364.   iter=0;
  1365.   retrain=1;
  1366.   maxfaktor=1;
  1367.   scalemaxiter=maxiter/5;
  1368.   while((retrain) && (maxviol > 0) && (iter < (scalemaxiter*maxfaktor))) {
  1369.     iter++;
  1370.     
  1371.     while((maxviol > precision) && (iter < (scalemaxiter*maxfaktor))) {
  1372.       iter++;
  1373.       maxviol=0;
  1374.       for(i=0;i<2*(n+m);i++) {
  1375. sum=d0[i];
  1376. for(j=0;j<2*(n+m);j++) {
  1377.   sum+=d[i*2*(n+m)+j]*dual_old[j];
  1378. }
  1379. sum-=d[i*2*(n+m)+i]*dual_old[i];
  1380. dual[i]=-sum/d[i*2*(n+m)+i];
  1381. if(dual[i]<0) dual[i]=0;
  1383. viol=fabs(dual[i]-dual_old[i]);
  1384. if(viol>maxviol) 
  1385.   maxviol=viol;
  1386. dual_old[i]=dual[i];
  1387.       }
  1388.       /*
  1389.       sprintf(temstr,"%d) maxviol=%20f precision=%fn",iter,maxviol,precision); 
  1390.       */
  1391.     }
  1392.   
  1393.     if(m>0) {
  1394.       for(i=0;i<n;i++) {
  1395. temp[i]=dual[i]-dual[i+n]+ce[i]*(dual[n+n]-dual[n+n+1])+g0[i];
  1396.       }
  1397.     } 
  1398.     else {
  1399.       for(i=0;i<n;i++) {
  1400. temp[i]=dual[i]-dual[i+n]+g0[i];
  1401.       }
  1402.     }
  1403.     for(i=0;i<n;i++) {
  1404.       primal[i]=0;             /* calc value of primal variables */
  1405.       for(j=0;j<n;j++) {
  1406. primal[i]+=ig[i*n+j]*temp[j];
  1407.       }
  1408.       primal[i]*=-1.0;
  1409.       if(primal[i]<=(low[i])) {  /* clip conservatively */
  1410. primal[i]=low[i];
  1411.       }
  1412.       else if(primal[i]>=(up[i])) {
  1413. primal[i]=up[i];
  1414.       }
  1415.     }
  1417.     if(m>0) 
  1418.       model_b=dual[n+n+1]-dual[n+n];
  1419.     else
  1420.       model_b=0;
  1422.     epsilon_hideo=EPSILON_HIDEO;
  1423.     for(i=0;i<n;i++) {           /* check precision of alphas */
  1424.       isnantest+=primal[j];
  1425.       dist=-model_b*ce[i]; 
  1426.       dist+=(g0[i]+1.0);
  1427.       for(j=0;j<i;j++) {
  1428. dist+=(primal[j]*g[j*n+i]);
  1429.       }
  1430.       for(j=i;j<n;j++) {
  1431. dist+=(primal[j]*g[i*n+j]);
  1432.       }
  1433.       if((primal[i]<(up[i]-epsilon_hideo)) && (dist < (1.0-epsilon_crit))) {
  1434. epsilon_hideo=(up[i]-primal[i])*2.0;
  1435.       }
  1436.       else if((primal[i]>(low[i]+epsilon_hideo)) &&(dist>(1.0+epsilon_crit))) {
  1437. epsilon_hideo=(primal[i]-low[i])*2.0;
  1438.       }
  1439.     }
  1440.      /*sprintf(temstr,"nEPSILON_HIDEO=%.30fn",epsilon_hideo); */
  1443.     for(i=0;i<n;i++) {           /* clip alphas to bounds */
  1444.       if(primal[i]<=(low[i]+epsilon_hideo)) {
  1445. primal[i]=low[i];
  1446.       }
  1447.       else if(primal[i]>=(up[i]-epsilon_hideo)) {
  1448. primal[i]=up[i];
  1449.       }
  1450.     }
  1452.     retrain=0;
  1453.     primal_optimal=1;
  1454.     at_bound=0;
  1455.     for(i=0;(i<n);i++) {  /* check primal KT-Conditions */
  1456.       dist=-model_b*ce[i]; 
  1457.       dist+=(g0[i]+1.0);
  1458.       for(j=0;j<i;j++) {
  1459. dist+=(primal[j]*g[j*n+i]);
  1460.       }
  1461.       for(j=i;j<n;j++) {
  1462. dist+=(primal[j]*g[i*n+j]);
  1463.       }
  1464.       if((primal[i]<(up[i]-epsilon_a)) && (dist < (1.0-epsilon_crit))) {
  1465. retrain=1;
  1466. primal_optimal=0;
  1467.       }
  1468.       else if((primal[i]>(low[i]+epsilon_a)) && (dist > (1.0+epsilon_crit))) {
  1469. retrain=1;
  1470. primal_optimal=0;
  1471.       }
  1472.       if((primal[i]<=(low[i]+epsilon_a)) || (primal[i]>=(up[i]-epsilon_a))) {
  1473. at_bound++;
  1474.       }
  1475.       /*    sprintf(temstr,"HIDEOtemp: a[%ld]=%.30f, dist=%.6f, b=%f, at_bound=%ldn",i,primal[i],dist,model_b,at_bound);  */
  1476.     }
  1477.     if(m>0) {
  1478.       eq=-ce0[0];               /* check precision of eq-constraint */
  1479.       for(i=0;i<n;i++) { 
  1480. eq+=(ce[i]*primal[i]);
  1481.       }
  1482.       if((EPSILON_EQ < fabs(eq)) 
  1483.  /*
  1484.  && !((goal==PRIMAL_OPTIMAL) 
  1485.        && (at_bound==n)) */
  1486.  ) {
  1487. retrain=1;
  1488. primal_optimal=0;
  1489.       }
  1490.       /* sprintf(temstr,"n eq=%.30f ce0=%f at-bound=%ldn",eq,ce0[0],at_bound);  */
  1491.     }
  1493.     if(retrain) {
  1494.       precision/=10;
  1495.       if(((goal == PRIMAL_OPTIMAL) && (maxfaktor < 50000))
  1496.  || (maxfaktor < 5)) {
  1497. maxfaktor++;
  1498.       }
  1499.     }
  1500.   }
  1502.   if(!primal_optimal) {
  1503.     for(i=0;i<n;i++) {
  1504.       primal[i]=0;             /* calc value of primal variables */
  1505.       for(j=0;j<n;j++) {
  1506. primal[i]+=ig[i*n+j]*temp[j];
  1507.       }
  1508.       primal[i]*=-1.0;
  1509.       if(primal[i]<=(low[i]+epsilon_a)) {  /* clip conservatively */
  1510. primal[i]=low[i];
  1511.       }
  1512.       else if(primal[i]>=(up[i]-epsilon_a)) {
  1513. primal[i]=up[i];
  1514.       }
  1515.     }
  1516.   }
  1518.   isnantest=0;
  1519.   for(i=0;i<n;i++) {           /* check for isnan */
  1520.     isnantest+=primal[i];
  1521.   }
  1523.   if(m>0) {
  1524.     temp1=dual[n+n+1];   /* copy the dual variables for the eq */
  1525.     temp2=dual[n+n];     /* constraints to a handier location */
  1526.     for(i=n+n+1;i>=2;i--) {
  1527.       dual[i]=dual[i-2];
  1528.     }
  1529.     dual[0]=temp2;
  1530.     dual[1]=temp1;
  1531.     isnantest+=temp1+temp2;
  1532.   }
  1534.   if(isnan(isnantest)) {
  1535.     return((int)NAN_SOLUTION);
  1536.   }
  1537.   else if(primal_optimal) {
  1538.     return((int)PRIMAL_OPTIMAL);
  1539.   }
  1540.   else if(maxviol == 0.0) {
  1541.     return((int)DUAL_OPTIMAL);
  1542.   }
  1543.   else {
  1544.     return((int)MAXITER_EXCEEDED);
  1545.   }
  1546. }
  1549. void CSVM::linvert_matrix(
  1550. double *matrix,
  1551. long depth,
  1552. double *inverse,double lindep_sensitivity,
  1553. long *lin_dependent)  /* indicates the active parts of matrix on 
  1554.  input and output*/
  1555. {
  1556.   long i,j,k;
  1557.   double factor;
  1559.   for(i=0;i<depth;i++) {
  1560.     /*    lin_dependent[i]=0; */
  1561.     for(j=0;j<depth;j++) {
  1562.       inverse[i*depth+j]=0.0;
  1563.     }
  1564.     inverse[i*depth+i]=1.0;
  1565.   }
  1566.   for(i=0;i<depth;i++) {
  1567.     if(lin_dependent[i] || (fabs(matrix[i*depth+i])<lindep_sensitivity)) {
  1568.       lin_dependent[i]=1;
  1569.     }
  1570.     else {
  1571.       for(j=i+1;j<depth;j++) {
  1572. factor=matrix[j*depth+i]/matrix[i*depth+i];
  1573. for(k=i;k<depth;k++) {
  1574.   matrix[j*depth+k]-=(factor*matrix[i*depth+k]);
  1575. }
  1576. for(k=0;k<depth;k++) {
  1577.   inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
  1578. }
  1579.       }
  1580.     }
  1581.   }
  1582.   for(i=depth-1;i>=0;i--) {
  1583.     if(!lin_dependent[i]) {
  1584.       factor=1/matrix[i*depth+i];
  1585.       for(k=0;k<depth;k++) {
  1586. inverse[i*depth+k]*=factor;
  1587.       }
  1588.       matrix[i*depth+i]=1;
  1589.       for(j=i-1;j>=0;j--) {
  1590. factor=matrix[j*depth+i];
  1591. matrix[j*depth+i]=0;
  1592. for(k=0;k<depth;k++) {
  1593.   inverse[j*depth+k]-=(factor*inverse[i*depth+k]);
  1594. }
  1595.       }
  1596.     }
  1597.   }
  1598. }
  1600. void CSVM::ladd_matrix(
  1601. double *matrix,
  1602. long depth,
  1603. double scalar)
  1604. {
  1605.   long i,j;
  1606.   for(i=0;i<depth;i++) {
  1607.     for(j=0;j<depth;j++) {
  1608.       matrix[i*depth+j]+=scalar;
  1609.     }
  1610.   }
  1611. }
  1613. void CSVM::lcopy_matrix(
  1614. double *matrix,
  1615. long depth,
  1616. double *matrix2)
  1617. {
  1618.   long i;
  1619.   
  1620.   for(i=0;i<(depth)*(depth);i++) {
  1621.     matrix2[i]=matrix[i];
  1622.   }
  1623. }
  1625. void CSVM::lswitch_rows_matrix(
  1626. double *matrix,
  1627. long depth,long r1,long r2)
  1628. {
  1629.   long i;
  1630.   double temp;
  1632.   for(i=0;i<depth;i++) {
  1633.     temp=matrix[r1*depth+i];
  1634.     matrix[r1*depth+i]=matrix[r2*depth+i];
  1635.     matrix[r2*depth+i]=temp;
  1636.   }
  1637. }
  1639. void CSVM::lswitchrk_matrix(
  1640. double *matrix,
  1641. long depth,long rk1,long rk2)
  1642. {
  1643.   long i;
  1644.   double temp;
  1646.   for(i=0;i<depth;i++) {
  1647.     temp=matrix[rk1*depth+i];
  1648.     matrix[rk1*depth+i]=matrix[rk2*depth+i];
  1649.     matrix[rk2*depth+i]=temp;
  1650.   }
  1651.   for(i=0;i<depth;i++) {
  1652.     temp=matrix[i*depth+rk1];
  1653.     matrix[i*depth+rk1]=matrix[i*depth+rk2];
  1654.     matrix[i*depth+rk2]=temp;
  1655.   }
  1656. }
  1658. double CSVM::calculate_qp_objective(
  1659. long opt_n,
  1660. double *opt_g,double *opt_g0,double *alpha)
  1661. {
  1662.   double obj;
  1663.   long i,j;
  1664.   obj=0;  /* calculate objective  */
  1665.   for(i=0;i<opt_n;i++) {
  1666.     obj+=(opt_g0[i]*alpha[i]);
  1667.     obj+=(0.5*alpha[i]*alpha[i]*opt_g[i*opt_n+i]);
  1668.     for(j=0;j<i;j++) {
  1669.       obj+=(alpha[j]*alpha[i]*opt_g[j*opt_n+i]);
  1670.     }
  1671.   }
  1672.   return(obj);
  1673. }
  1674. /********************************svm_hideo****************************/
  1677. /********************************svm_learn****************************/
  1678. /* Learns an SVM model based on the training data in docs/label. The resulting
  1679. model is returned in the structure model. */
  1681. void CSVM::svm_learn(
  1682.                DOC *docs,                
  1683.                long *label,               
  1684.                long totdoc,               
  1685.                long totwords,             
  1686.                LEARN_PARM *learn_parm,    
  1687.                KERNEL_PARM *kernel_parm, 
  1688.                KERNEL_CACHE *kernel_cache,
  1689.                MODEL *model             
  1690.                )
  1691. {
  1692.     long *inconsistent,i;
  1693.     long inconsistentnum;
  1694.     long misclassified,upsupvecnum;
  1695.     double loss,model_length,example_length;
  1696.     double maxdiff,*lin,*a;
  1697.     long runtime_start,runtime_end;
  1698.     long iterations;
  1699.     long *unlabeled,transduction;
  1700.     long heldout;
  1701.     long loo_count=0,loo_count_pos=0,loo_count_neg=0,trainpos=0,trainneg=0;
  1702.     long loocomputed=0,runtime_start_loo=0,runtime_start_xa=0;
  1703.     double heldout_c=0,r_delta_sq=0,r_delta,r_delta_avg;
  1704.     
  1705.     double *xi_fullset; /* buffer for storing xi on full sample in loo */
  1706.     double *a_fullset;  /* buffer for storing alpha on full sample in loo */
  1707.     TIMING timing_profile;
  1708.     SHRINK_STATE shrink_state;
  1709.     
  1710.     runtime_start=get_runtime();
  1712.     timing_profile.time_kernel=0;
  1713.     timing_profile.time_opti=0;
  1714.     timing_profile.time_shrink=0;
  1715.     timing_profile.time_update=0;
  1716.     timing_profile.time_model=0;
  1717.     timing_profile.time_check=0;
  1718.     timing_profile.time_select=0;
  1719.     
  1720. com_result.kernel_cache_statistic=0;
  1721.     
  1722.     learn_parm->totwords=totwords;
  1723.     
  1724.     /* make sure -n value is reasonable */
  1725.     if((learn_parm->svm_newvarsinqp < 2) || (learn_parm->svm_newvarsinqp > learn_parm->svm_maxqpsize))
  1726.     {
  1727.         learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;
  1728.     }
  1729.     
  1730.     init_shrink_state(&shrink_state,totdoc,(long)10000);
  1731.     
  1732.     inconsistent = (long *)my_malloc(sizeof(long)*totdoc);
  1733.     unlabeled = (long *)my_malloc(sizeof(long)*totdoc);
  1734.     a = (double *)my_malloc(sizeof(double)*totdoc);
  1735.     a_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  1736.     xi_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  1737.     lin = (double *)my_malloc(sizeof(double)*totdoc);
  1738.     learn_parm->svm_cost = (double *)my_malloc(sizeof(double)*totdoc);
  1739.     model->supvec = (DOC **)my_malloc(sizeof(DOC *)*(totdoc+2));
  1740.     model->alpha = (double *)my_malloc(sizeof(double)*(totdoc+2));
  1741.     model->index = (long *)my_malloc(sizeof(long)*(totdoc+2));
  1743.     
  1744.     model->at_upper_bound=0;
  1745.     model->b=0;        
  1746.     model->supvec[0]=0;  /* element 0 reserved and empty for now */
  1747.     model->alpha[0]=0;
  1748.     model->lin_weights=NULL;
  1749.     model->totwords=totwords;
  1750.     model->totdoc=totdoc;
  1751.     model->kernel_parm=(*kernel_parm);
  1752.     model->sv_num=1;
  1753.     model->loo_error=-1;
  1754.     model->loo_recall=-1;
  1755.     model->loo_precision=-1;
  1756.     model->xa_error=-1;
  1757.     model->xa_recall=-1;
  1758.     model->xa_precision=-1;
  1759.     inconsistentnum=0;
  1760.     transduction=0;
  1761.     
  1762.     r_delta=estimate_r_delta(docs,totdoc,kernel_parm);
  1763.     r_delta_sq=r_delta*r_delta;
  1764.     
  1765.     r_delta_avg=estimate_r_delta_average(docs,totdoc,kernel_parm);
  1766.     if(learn_parm->svm_c == 0.0)
  1767.     {  /* default value for C */
  1768.         learn_parm->svm_c=1.0/(r_delta_avg*r_delta_avg);
  1769. if (com_pro.show_compute_1)
  1770. {
  1771.         sprintf(temstr,"Setting default regularization parameter C=%.4fn",learn_parm->svm_c);
  1772.         printm(temstr);
  1773. }
  1774.     }
  1775.     
  1776.     for(i=0;i<totdoc;i++) 
  1777.     {    /* various inits */
  1778.         inconsistent[i]=0;
  1779.         a[i]=0;
  1780.         lin[i]=0;
  1781.         unlabeled[i]=0;
  1782.         if(label[i] == 0) 
  1783.         {
  1784.             unlabeled[i]=1;
  1785.             transduction=1;
  1786.         }
  1787.         if(label[i] > 0)
  1788.         {
  1789.             learn_parm->svm_cost[i]=learn_parm->svm_c*learn_parm->svm_costratio*
  1790.                 fabs((double)label[i]);
  1791.             label[i]=1;
  1792.             trainpos++;
  1793.         }
  1794.         else if(label[i] < 0) 
  1795.         {
  1796.             learn_parm->svm_cost[i]=learn_parm->svm_c*fabs((double)label[i]);
  1797.             label[i]=-1;
  1798.             trainneg++;
  1799.         }
  1800.         else
  1801.         {
  1802.             learn_parm->svm_cost[i]=0;
  1803.         }
  1804.     }
  1805.     
  1806.     /* caching makes no sense for linear kernel */
  1807.     if(kernel_parm->kernel_type == LINEAR)
  1808.     {
  1809.         kernel_cache = NULL;   
  1810.     } 
  1811.     
  1812.     if(transduction) 
  1813.     {
  1814.         learn_parm->svm_iter_to_shrink=99999999;
  1815.         sprintf(temstr,"nDeactivating Shrinking due to an incompatibility with the transductive nlearner in the current version.nn");
  1816.         printm(temstr);
  1817.     }
  1818.     
  1819.     if(transduction && learn_parm->compute_loo) 
  1820.     {
  1821.         learn_parm->compute_loo=0;
  1822.         sprintf(temstr,"nCannot compute leave-one-out estimates for transductive learner.nn");
  1823.         printm(temstr);
  1824.     }
  1825.     
  1826.     if(learn_parm->remove_inconsistent && learn_parm->compute_loo) 
  1827.     {
  1828.         learn_parm->compute_loo=0;
  1829.         sprintf(temstr,"nCannot compute leave-one-out estimates when removing inconsistent examples.nn");
  1830.         printm(temstr);
  1831.     }    
  1832.     
  1833.     if((trainpos == 1) || (trainneg == 1)) 
  1834.     {
  1835.         learn_parm->compute_loo=0;
  1836.         sprintf(temstr,"nCannot compute leave-one-out with only one example in one class.nn");
  1837.         printm(temstr);
  1838.     }    
  1839.     
  1840.     if (com_pro.show_action)
  1841. {
  1842. sprintf(temstr,"Optimizing..."); 
  1843. printm(temstr);
  1844. }
  1845.     
  1846.     /* train the svm */
  1847.     iterations=optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
  1848.         kernel_parm,kernel_cache,&shrink_state,model,inconsistent,unlabeled,a,lin,&timing_profile,  &maxdiff,(long)-1,(long)1);
  1849.     if (com_pro.show_action)
  1850. {
  1851. sprintf(temstr,"done. (%ld iterations) ",iterations);
  1852. printm(temstr);
  1853. }
  1854.     
  1855.     misclassified=0;
  1856.     for(i=0;(i<totdoc);i++)
  1857.     { /* get final statistic */
  1858.         if((lin[i]-model->b)*(double)label[i] <= 0.0) 
  1859.             misclassified++;
  1860.     }
  1861. if (com_pro.show_action)
  1862. {
  1863. printm("optimization finished");
  1864. }
  1865. if (com_pro.show_trainresult)
  1866. {
  1867. sprintf(temstr," (%ld misclassified, maxdiff=%.5f).n", misclassified,maxdiff); 
  1868. printm(temstr);
  1869. }
  1870. com_result.train_misclassify=misclassified;
  1871. com_result.max_difference=maxdiff;
  1872.              
  1873.     runtime_end=get_runtime();
  1874.              
  1875.     if (learn_parm->remove_inconsistent)
  1876.     {     
  1877.         inconsistentnum=0;
  1878.         for(i=0;i<totdoc;i++) 
  1879.             if(inconsistent[i]) 
  1880.                inconsistentnum++;
  1881.         sprintf(temstr,"Number of SV: %ld (plus %ld inconsistent examples)n", model->sv_num-1,inconsistentnum);
  1882.         printm(temstr);
  1883.     }
  1884.     
  1885.     else
  1886.     {
  1887.      upsupvecnum=0;
  1888.      for(i=1;i<model->sv_num;i++) 
  1889.      {
  1890.          if(fabs(model->alpha[i]) >= (learn_parm->svm_cost[(model->supvec[i])->docnum]-learn_parm->epsilon_a)) 
  1891.              upsupvecnum++;
  1892.      }
  1893.  if (com_pro.show_trainresult)
  1894.  {
  1895.  sprintf(temstr,"Number of SV: %ld (including %ld at upper bound)n", model->sv_num-1,upsupvecnum);
  1896.      printm(temstr);
  1897.  }
  1898.     }
  1900. if( (!learn_parm->skip_final_opt_check)) 
  1901. {
  1902. loss=0;
  1903. model_length=0; 
  1904. for(i=0;i<totdoc;i++)
  1905. {
  1906. if((lin[i]-model->b)*(double)label[i] < 1.0-learn_parm->epsilon_crit)
  1907. loss+=1.0-(lin[i]-model->b)*(double)label[i];
  1908. model_length+=a[i]*label[i]*lin[i];
  1909. }
  1910. model_length=sqrt(model_length);
  1911. sprintf(temstr,"L1 loss: loss=%.5fn",loss);   printm(temstr);
  1912. sprintf(temstr,"Norm of weight vector: |w|=%.5fn",model_length);printm(temstr);
  1913. example_length=estimate_sphere(model,kernel_parm); 
  1914. sprintf(temstr,"Norm of longest example vector: |x|=%.5fn",  length_of_longest_document_vector(docs,totdoc,kernel_parm));
  1915. printm(temstr);
  1916. sprintf(temstr,"Estimated VCdim of classifier: VCdim<=%.5fn",       estimate_margin_vcdim(model,model_length,example_length,    kernel_parm));
  1917. printm(temstr);
  1918. if((!learn_parm->remove_inconsistent) && (!transduction)) 
  1919. {
  1920. runtime_start_xa=get_runtime();
  1921.                      sprintf(temstr,"Computing XiAlpha-estimates..."); 
  1922.                      printm(temstr);
  1923.                      compute_xa_estimates(model,label,unlabeled,totdoc,docs,lin,a,
  1924.                          kernel_parm,learn_parm,&(model->xa_error),
  1925.                          &(model->xa_recall),&(model->xa_precision));
  1926.                      
  1927.                      
  1928.                      sprintf(temstr,"Runtime for XiAlpha-estimates in cpu-seconds: %.2fn",
  1929.                          (get_runtime()-runtime_start_xa)/100.0);
  1930.                      printm(temstr);
  1931.                      
  1932.                      fprintf(stdout,"XiAlpha-estimate of the error: error<=%.2f%% (rho=%.2f,depth=%ld)n",
  1933.                          model->xa_error,learn_parm->rho,learn_parm->xa_depth);
  1934.                      fprintf(stdout,"XiAlpha-estimate of the recall: recall=>%.2f%% (rho=%.2f,depth=%ld)n",
  1935.                          model->xa_recall,learn_parm->rho,learn_parm->xa_depth);
  1936.                      fprintf(stdout,"XiAlpha-estimate of the precision: precision=>%.2f%% (rho=%.2f,depth=%ld)n",
  1937.                          model->xa_precision,learn_parm->rho,learn_parm->xa_depth);
  1938.                  }
  1939.                  else if(!learn_parm->remove_inconsistent)
  1940.                  {
  1941.                      estimate_transduction_quality(model,label,unlabeled,totdoc,docs,lin);
  1942.                  }
  1943.              }
  1944. if (com_pro.show_trainresult)
  1945. {
  1946. sprintf(temstr,"Number of kernel evaluations: %ldn",com_result.kernel_cache_statistic);
  1947.         printm(temstr);
  1948. }
  1949.              /* leave-one-out testing starts now */
  1950.              if(learn_parm->compute_loo)
  1951.              {
  1952.                  /* save results of training on full dataset for leave-one-out */
  1953.                  runtime_start_loo=get_runtime();
  1954.                  for(i=0;i<totdoc;i++) 
  1955.                  {
  1956.                      xi_fullset[i]=1.0-((lin[i]-model->b)*(double)label[i]);
  1957.                      a_fullset[i]=a[i];
  1958.                  }
  1959.                  sprintf(temstr,"Computing leave-one-out");
  1960.                  printm(temstr);
  1961.                  
  1962.                  /* repeat this loop for every held-out example */
  1963.                  for(heldout=0;(heldout<totdoc);heldout++)
  1964.                  {
  1965.                      if(learn_parm->rho*a_fullset[heldout]*r_delta_sq+xi_fullset[heldout]
  1966.                          < 1.0) 
  1967.                      { 
  1968.                          /* guaranteed to not produce a leave-one-out error */
  1969.                          sprintf(temstr,"+"); 
  1970.                          printm(temstr);
  1971.                      }
  1972.                      else if(xi_fullset[heldout] > 1.0) 
  1973.                      {
  1974.                          /* guaranteed to produce a leave-one-out error */
  1975.                          loo_count++;
  1976.                          if(label[heldout] > 0) loo_count_pos++; else loo_count_neg++;
  1977.                          sprintf(temstr,"-");  printm(temstr);
  1978.                      }
  1979.                      else
  1980.                      {
  1981.                          loocomputed++;
  1982.                          heldout_c=learn_parm->svm_cost[heldout]; /* set upper bound to zero */
  1983.                          learn_parm->svm_cost[heldout]=0;
  1984.                          /* make sure heldout example is not currently  */
  1985.                          /* shrunk away. Assumes that lin is up to date! */
  1986.                          shrink_state.active[heldout]=1;  
  1987.                          
  1988.                          
  1989.                          optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
  1990.                              kernel_parm,
  1991.                              kernel_cache,&shrink_state,model,inconsistent,unlabeled,
  1992.                              a,lin,&timing_profile,
  1993.                              &maxdiff,heldout,(long)2);
  1994.                          
  1995.                          /* printf("%fn",(lin[heldout]-model->b)*(double)label[heldout]); */
  1996.                          
  1997.                          if(((lin[heldout]-model->b)*(double)label[heldout]) < 0.0)
  1998.                          { 
  1999.                              loo_count++;                           /* there was a loo-error */
  2000.                              if(label[heldout] > 0) loo_count_pos++; else loo_count_neg++;
  2001.                          }
  2002.                          else
  2003.                          {
  2004.                              
  2005.                          }
  2006.                          /* now we need to restore the original data set*/
  2007.                          learn_parm->svm_cost[heldout]=heldout_c; /* restore upper bound */
  2008.                      }
  2009.                  } /* end of leave-one-out loop */
  2010.                  
  2011.                  
  2012.                  sprintf(temstr,"nRetrain on full problem");  printm(temstr);
  2013.                  optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
  2014.                      kernel_parm,
  2015.                      kernel_cache,&shrink_state,model,inconsistent,unlabeled,
  2016.                      a,lin,&timing_profile,
  2017.                      &maxdiff,(long)-1,(long)1);
  2018.                  
  2019.                  
  2020.                  /* after all leave-one-out computed */
  2021.                  model->loo_error=100.0*loo_count/(double)totdoc;
  2022.                  model->loo_recall=(1.0-(double)loo_count_pos/(double)trainpos)*100.0;
  2023.                  model->loo_precision=(trainpos-loo_count_pos)/
  2024.                      (double)(trainpos-loo_count_pos+loo_count_neg)*100.0;
  2025.                  fprintf(stdout,"Leave-one-out estimate of the error: error=%.2f%%n",
  2026.                      model->loo_error);
  2027.                  fprintf(stdout,"Leave-one-out estimate of the recall: recall=%.2f%%n",
  2028.                      model->loo_recall);
  2029.                  fprintf(stdout,"Leave-one-out estimate of the precision: precision=%.2f%%n",
  2030.                      model->loo_precision);
  2031.                  fprintf(stdout,"Actual leave-one-outs computed:  %ld (rho=%.2f)n",
  2032.                      loocomputed,learn_parm->rho);
  2033.                  sprintf(temstr,"Runtime for leave-one-out in cpu-seconds: %.2fn",
  2034.                      (double)(get_runtime()-runtime_start_loo)/100.0);
  2035.                  printm(temstr);
  2036.              }
  2037.              
  2038.              // if(learn_parm->alphafile[0])
  2039.              //     write_alphas(learn_parm->alphafile,a,label,totdoc);
  2040.              
  2041.              shrink_state_cleanup(&shrink_state);
  2042.              
  2043.              free(inconsistent);
  2044.              free(unlabeled);
  2045.              free(a);
  2046.              free(a_fullset);
  2047.              free(xi_fullset);
  2048.              free(lin);
  2049.              free(learn_parm->svm_cost);
  2050. }
  2054. long CSVM::optimize_to_convergence(
  2055.                              DOC *docs,
  2056.                              long *label,
  2057.                              long totdoc,
  2058.                              long totwords,           
  2059.                              LEARN_PARM *learn_parm,    
  2060.                              KERNEL_PARM *kernel_parm,  
  2061.                              KERNEL_CACHE *kernel_cache,
  2062.                              SHRINK_STATE *shrink_state,
  2063.                              MODEL *model,              
  2064.                              long *inconsistent,
  2065.                              long *unlabeled,
  2066.                              double *a,
  2067.                              double *lin,
  2068.                              TIMING *timing_profile,
  2069.                              double *maxdiff,
  2070.                              long heldout,
  2071.                              long retrain)
  2072. {
  2073.     long *chosen,*key,i,j,jj,*last_suboptimal_at,noshrink;
  2074.     long inconsistentnum,choosenum,already_chosen=0,iteration;
  2075.     long misclassified,supvecnum=0,*active2dnum,inactivenum;
  2076.     long *working2dnum,*selexam;
  2077.     long activenum;
  2078.     double criterion,eq;
  2079.     double *a_old;
  2080.     long t0=0,t1=0,t2=0,t3=0,t4=0,t5=0,t6=0; /* timing */
  2081.     long transductcycle;
  2082.     long transduction;
  2083.     double epsilon_crit_org; 
  2084.     
  2085.     double *selcrit;  /* buffer for sorting */        
  2086.     CFLOAT *aicache;  /* buffer to keep one row of hessian */
  2087.     double *weights;  /* buffer for weight vector in linear case */
  2088.     QP qp;            /* buffer for one quadratic program */
  2089.     
  2090.     epsilon_crit_org=learn_parm->epsilon_crit; /* save org */
  2091.     if(kernel_parm->kernel_type == LINEAR)
  2092.     {
  2093.         learn_parm->epsilon_crit=2.0;
  2094.         kernel_cache=NULL;   /* caching makes no sense for linear kernel */
  2095.     } 
  2096.     learn_parm->epsilon_shrink=2;
  2097.     (*maxdiff)=1;
  2098.     
  2099.     learn_parm->totwords=totwords;
  2100.     
  2101.     chosen = (long *)my_malloc(sizeof(long)*totdoc);
  2102.     last_suboptimal_at = (long *)my_malloc(sizeof(long)*totdoc);
  2103.     key = (long *)my_malloc(sizeof(long)*(totdoc+11)); 
  2104.     selcrit = (double *)my_malloc(sizeof(double)*totdoc);
  2105.     selexam = (long *)my_malloc(sizeof(long)*totdoc);
  2106.     a_old = (double *)my_malloc(sizeof(double)*totdoc);
  2107.     aicache = (CFLOAT *)my_malloc(sizeof(CFLOAT)*totdoc);
  2108.     working2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  2109.     active2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  2110.     qp.opt_ce = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  2111.     qp.opt_ce0 = (double *)my_malloc(sizeof(double));
  2112.     qp.opt_g = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize
  2113.         *learn_parm->svm_maxqpsize);
  2114.     qp.opt_g0 = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  2115.     qp.opt_xinit = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  2116.     qp.opt_low=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  2117.     qp.opt_up=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  2118.     weights=(double *)my_malloc(sizeof(double)*(totwords+1));
  2119.     
  2120.     choosenum=0;
  2121.     inconsistentnum=0;
  2122.     transductcycle=0;
  2123.     transduction=0;
  2124.     if(!retrain) retrain=1;
  2125.     iteration=1;
  2126.     
  2127.     if(kernel_cache)
  2128.     {
  2129.         kernel_cache->time=iteration;  /* for lru cache */
  2130.         kernel_cache_reset_lru(kernel_cache);
  2131.     }
  2132.     
  2133.     for(i=0;i<totdoc;i++)
  2134.     {    /* various inits */
  2135.         chosen[i]=0;
  2136.         a_old[i]=a[i];
  2137.         last_suboptimal_at[i]=1;
  2138.         if(inconsistent[i]) 
  2139.             inconsistentnum++;
  2140.         if(unlabeled[i]) 
  2141.         {
  2142.             transduction=1;
  2143.         }
  2144.     }
  2145.     activenum=compute_index(shrink_state->active,totdoc,active2dnum);
  2146.     inactivenum=totdoc-activenum;
  2147.     clear_index(working2dnum);
  2148.     
  2149.     /* repeat this loop until we have convergence */
  2150.     for(;retrain;iteration++) 
  2151.     {
  2152.         
  2153.         if(kernel_cache)
  2154.             kernel_cache->time=iteration;  /* for lru cache */
  2155.         i=0;
  2156.         for(jj=0;(j=working2dnum[jj])>=0;jj++)
  2157.         { /* clear working set */
  2158.             if((chosen[j]>=(learn_parm->svm_maxqpsize/
  2159.                 minl(learn_parm->svm_maxqpsize,
  2160.                 learn_parm->svm_newvarsinqp))) 
  2161.                 || (inconsistent[j])
  2162.                 || (j == heldout)) 
  2163.             {
  2164.                 chosen[j]=0; 
  2165.                 choosenum--; 
  2166.             }
  2167.             else 
  2168.             {
  2169.                 chosen[j]++;
  2170.                 working2dnum[i++]=j;
  2171.             }
  2172.         }
  2173.         working2dnum[i]=-1;
  2174.         
  2175.         if(retrain == 2)
  2176.         {
  2177.             choosenum=0;
  2178.             for(jj=0;(j=working2dnum[jj])>=0;jj++) 
  2179.             { /* fully clear working set */
  2180.                 chosen[j]=0; 
  2181.             }
  2182.             clear_index(working2dnum);
  2183.             for(i=0;i<totdoc;i++) 
  2184.             { /* set inconsistent examples to zero (-i 1) */
  2185.                 if((inconsistent[i] || (heldout==i)) && (a[i] != 0.0)) 
  2186.                 {
  2187.                     chosen[i]=99999;
  2188.                     choosenum++;
  2189.                     a[i]=0;
  2190.                 }
  2191.             }
  2192.             if(learn_parm->biased_hyperplane)
  2193.             {
  2194.                 eq=0;
  2195.                 for(i=0;i<totdoc;i++) 
  2196.                 { /* make sure we fulfill equality constraint */
  2197.                     eq+=a[i]*label[i];
  2198.                 }
  2199.                 for(i=0;(i<totdoc) && (fabs(eq) > learn_parm->epsilon_a);i++) 
  2200.                 {
  2201.                     if((eq*label[i] > 0) && (a[i] > 0)) 
  2202.                     {
  2203.                         chosen[i]=88888;
  2204.                         choosenum++;
  2205.                         if((eq*label[i]) > a[i]) 
  2206.                         {
  2207.                             eq-=(a[i]*label[i]);
  2208.                             a[i]=0;
  2209.                         }
  2210.                         else 
  2211.                         {
  2212.                             a[i]-=(eq*label[i]);
  2213.                             eq=0;
  2214.                         }
  2215.                     }
  2216.                 }
  2217.             }
  2218.             compute_index(chosen,totdoc,working2dnum);
  2219.         }
  2220.         else
  2221.         {      /* select working set according to steepest gradient */ 
  2222.             if((minl(learn_parm->svm_newvarsinqp,learn_parm->svm_maxqpsize)>=4) 
  2223.                 && (kernel_parm->kernel_type != LINEAR)) 
  2224.             {
  2225.                 /* select part of the working set from cache */
  2226.                 already_chosen=select_next_qp_subproblem_grad_cache(
  2227.                     label,unlabeled,a,lin,totdoc,
  2228.                     minl((long)(learn_parm->svm_maxqpsize-choosenum),
  2229.                     (long)(learn_parm->svm_newvarsinqp/2)),
  2230.                     learn_parm,inconsistent,active2dnum,
  2231.                     working2dnum,selcrit,selexam,kernel_cache,
  2232.                     key,chosen);
  2233.                 choosenum+=already_chosen;
  2234.             }
  2235.             choosenum+=select_next_qp_subproblem_grad(label,unlabeled,a,lin,totdoc,
  2236.                 minl((long)(learn_parm->svm_maxqpsize-choosenum),
  2237.                 (long)(learn_parm->svm_newvarsinqp-already_chosen)),
  2238.                 learn_parm,inconsistent,active2dnum,
  2239.                 working2dnum,selcrit,selexam,kernel_cache,key,
  2240.                 chosen);
  2241.         }
  2242.         
  2243.         //      sprintf(temstr," %ld vectors chosenn",choosenum);  printm(temstr);
  2244.         
  2245.         
  2246.         t1=get_runtime();
  2247.         
  2248.         if(kernel_cache) 
  2249.             cache_multiple_kernel_rows(kernel_cache,docs,working2dnum,choosenum,kernel_parm); 
  2250.         
  2251.         
  2252.         t2=get_runtime();
  2253.         if(retrain != 2)
  2254.         {
  2255.             optimize_svm(docs,label,unlabeled,chosen,active2dnum,model,totdoc,
  2256.                 working2dnum,choosenum,a,lin,learn_parm,aicache,
  2257.                 kernel_parm,&qp,&epsilon_crit_org);
  2258.         }
  2259.         
  2260.         
  2261.         t3=get_runtime();
  2262.         update_linear_component(docs,label,active2dnum,a,a_old,working2dnum,totdoc,
  2263.             totwords,kernel_parm,kernel_cache,lin,aicache,weights);
  2264.         
  2265.         
  2266.         
  2267.         t4=get_runtime();
  2268.         supvecnum=calculate_svm_model(docs,label,unlabeled,lin,a,a_old,learn_parm,
  2269.             working2dnum,model);
  2270.         
  2271.         
  2272.         t5=get_runtime();
  2273.         
  2274.         /* The following computation of the objective function works only */
  2275.         /* relative to the active variables */
  2276.         
  2277.         criterion=compute_objective_function(a,lin,label,active2dnum);
  2278.         //  sprintf(temstr,"Objective function (over active variables): %.16fn",criterion);
  2279.         //  printm(temstr);
  2280.         
  2281.         for(jj=0;(i=working2dnum[jj])>=0;jj++)