view toolboxes/SVM-light/src/svm_learn.c @ 0:e9a9cd732c1e tip

first hg version after svn

/***********************************************************************/
/*                                                                     */
/*   svm_learn.c                                                       */
/*                                                                     */
/*   Learning module of Support Vector Machine.                        */
/*                                                                     */
/*   Author: Thorsten Joachims                                         */
/*   Date: 02.07.02                                                    */
/*                                                                     */
/*   Copyright (c) 2002  Thorsten Joachims - All rights reserved       */
/*                                                                     */
/*   This software is available for non-commercial use only. It must   */
/*   not be modified and distributed without prior permission of the   */
/*   author. The author is not responsible for implications from the   */
/*   use of this software.                                             */
/*                                                                     */
/***********************************************************************/


# include "svm_common.h"
# include "svm_learn.h"


/* interface to QP-solver */
double *optimize_qp(QP *, double *, long, double *, LEARN_PARM *);

/*---------------------------------------------------------------------------*/

/* Learns an SVM classification model based on the training data in
   docs/label. The resulting model is returned in the structure
   model. */

void svm_learn_classification(DOC **docs, double *class, long int
			      totdoc, long int totwords, 
			      LEARN_PARM *learn_parm, 
			      KERNEL_PARM *kernel_parm, 
			      KERNEL_CACHE *kernel_cache, 
			      MODEL *model,
			      double *alpha)
     /* docs:        Training vectors (x-part) */
     /* class:       Training labels (y-part, zero if test example for
                     transduction) */
     /* totdoc:      Number of examples in docs/label */
     /* totwords:    Number of features (i.e. highest feature index) */
     /* learn_parm:  Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache:Initialized Cache of size totdoc, if using a kernel. 
                     NULL if linear.*/
     /* model:       Returns learning result (assumed empty before called) */
     /* alpha:       Start values for the alpha variables or NULL
	             pointer. The new alpha values are returned after 
		     optimization if not NULL. Array must be of size totdoc. */
{
  long *inconsistent,i,*label;
  long inconsistentnum;
  long misclassified,upsupvecnum;
  double loss,model_length,example_length;
  double maxdiff,*lin,*a,*c;
  long runtime_start,runtime_end;
  long iterations;
  long *unlabeled,transduction;
  long heldout;
  long loo_count=0,loo_count_pos=0,loo_count_neg=0,trainpos=0,trainneg=0;
  long loocomputed=0,runtime_start_loo=0,runtime_start_xa=0;
  double heldout_c=0,r_delta_sq=0,r_delta,r_delta_avg;
  long *index,*index2dnum;
  double *weights;
  CFLOAT *aicache;  /* buffer to keep one row of hessian */

  double *xi_fullset; /* buffer for storing xi on full sample in loo */
  double *a_fullset;  /* buffer for storing alpha on full sample in loo */
  TIMING timing_profile;
  SHRINK_STATE shrink_state;

  runtime_start=get_runtime();
  timing_profile.time_kernel=0;
  timing_profile.time_opti=0;
  timing_profile.time_shrink=0;
  timing_profile.time_update=0;
  timing_profile.time_model=0;
  timing_profile.time_check=0;
  timing_profile.time_select=0;
  kernel_cache_statistic=0;

  learn_parm->totwords=totwords;

  /* make sure -n value is reasonable */
  if((learn_parm->svm_newvarsinqp < 2) 
     || (learn_parm->svm_newvarsinqp > learn_parm->svm_maxqpsize)) {
    learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;
  }

  init_shrink_state(&shrink_state,totdoc,(long)MAXSHRINK);

  label = (long *)my_malloc(sizeof(long)*totdoc);
  inconsistent = (long *)my_malloc(sizeof(long)*totdoc);
  unlabeled = (long *)my_malloc(sizeof(long)*totdoc);
  c = (double *)my_malloc(sizeof(double)*totdoc);
  a = (double *)my_malloc(sizeof(double)*totdoc);
  a_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  xi_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  lin = (double *)my_malloc(sizeof(double)*totdoc);
  learn_parm->svm_cost = (double *)my_malloc(sizeof(double)*totdoc);
  model->supvec = (DOC **)my_malloc(sizeof(DOC *)*(totdoc+2));
  model->alpha = (double *)my_malloc(sizeof(double)*(totdoc+2));
  model->index = (long *)my_malloc(sizeof(long)*(totdoc+2));

  model->at_upper_bound=0;
  model->b=0;	       
  model->supvec[0]=0;  /* element 0 reserved and empty for now */
  model->alpha[0]=0;
  model->lin_weights=NULL;
  model->totwords=totwords;
  model->totdoc=totdoc;
  model->kernel_parm=(*kernel_parm);
  model->sv_num=1;
  model->loo_error=-1;
  model->loo_recall=-1;
  model->loo_precision=-1;
  model->xa_error=-1;
  model->xa_recall=-1;
  model->xa_precision=-1;
  inconsistentnum=0;
  transduction=0;

  r_delta=estimate_r_delta(docs,totdoc,kernel_parm);
  r_delta_sq=r_delta*r_delta;

  r_delta_avg=estimate_r_delta_average(docs,totdoc,kernel_parm);
  if(learn_parm->svm_c == 0.0) {  /* default value for C */
    learn_parm->svm_c=1.0/(r_delta_avg*r_delta_avg);
    if(verbosity>=1) 
      printf("Setting default regularization parameter C=%.4f\n",
	     learn_parm->svm_c);
  }

  learn_parm->eps=-1.0;      /* equivalent regression epsilon for
				classification */

  for(i=0;i<totdoc;i++) {    /* various inits */
    docs[i]->docnum=i;
    inconsistent[i]=0;
    a[i]=0;
    lin[i]=0;
    c[i]=0.0;
    unlabeled[i]=0;
    if(class[i] == 0) {
      unlabeled[i]=1;
      label[i]=0;
      transduction=1;
    }
    if(class[i] > 0) {
      learn_parm->svm_cost[i]=learn_parm->svm_c*learn_parm->svm_costratio*
	docs[i]->costfactor;
      label[i]=1;
      trainpos++;
    }
    else if(class[i] < 0) {
      learn_parm->svm_cost[i]=learn_parm->svm_c*docs[i]->costfactor;
      label[i]=-1;
      trainneg++;
    }
    else {
      learn_parm->svm_cost[i]=0;
    }
  }
  if(verbosity>=2) {
    printf("%ld positive, %ld negative, and %ld unlabeled examples.\n",trainpos,trainneg,totdoc-trainpos-trainneg); fflush(stdout);
  }

  /* caching makes no sense for linear kernel */
  if(kernel_parm->kernel_type == LINEAR) {
    kernel_cache = NULL;   
  } 

  /* compute starting state for initial alpha values */
  if(alpha) {
    if(verbosity>=1) {
      printf("Computing starting state..."); fflush(stdout);
    }
    index = (long *)my_malloc(sizeof(long)*totdoc);
    index2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
    weights=(double *)my_malloc(sizeof(double)*(totwords+1));
    aicache = (CFLOAT *)my_malloc(sizeof(CFLOAT)*totdoc);
    for(i=0;i<totdoc;i++) {    /* create full index and clip alphas */
      index[i]=1;
      alpha[i]=fabs(alpha[i]);
      if(alpha[i]<0) alpha[i]=0;
      if(alpha[i]>learn_parm->svm_cost[i]) alpha[i]=learn_parm->svm_cost[i];
    }
    if(kernel_parm->kernel_type != LINEAR) {
      for(i=0;i<totdoc;i++)     /* fill kernel cache with unbounded SV */
	if((alpha[i]>0) && (alpha[i]<learn_parm->svm_cost[i]) 
	   && (kernel_cache_space_available(kernel_cache))) 
	  cache_kernel_row(kernel_cache,docs,i,kernel_parm);
      for(i=0;i<totdoc;i++)     /* fill rest of kernel cache with bounded SV */
	if((alpha[i]==learn_parm->svm_cost[i]) 
	   && (kernel_cache_space_available(kernel_cache))) 
	  cache_kernel_row(kernel_cache,docs,i,kernel_parm);
    }
    (void)compute_index(index,totdoc,index2dnum);
    update_linear_component(docs,label,index2dnum,alpha,a,index2dnum,totdoc,
			    totwords,kernel_parm,kernel_cache,lin,aicache,
			    weights);
    (void)calculate_svm_model(docs,label,unlabeled,lin,alpha,a,c,
			      learn_parm,index2dnum,index2dnum,model);
    for(i=0;i<totdoc;i++) {    /* copy initial alphas */
      a[i]=alpha[i];
    }
    free(index);
    free(index2dnum);
    free(weights);
    free(aicache);
    if(verbosity>=1) {
      printf("done.\n");  fflush(stdout);
    }   
  } 

  if(transduction) {
    learn_parm->svm_iter_to_shrink=99999999;
    if(verbosity >= 1)
      printf("\nDeactivating Shrinking due to an incompatibility with the transductive \nlearner in the current version.\n\n");
  }

  if(transduction && learn_parm->compute_loo) {
    learn_parm->compute_loo=0;
    if(verbosity >= 1)
      printf("\nCannot compute leave-one-out estimates for transductive learner.\n\n");
  }    

  if(learn_parm->remove_inconsistent && learn_parm->compute_loo) {
    learn_parm->compute_loo=0;
    printf("\nCannot compute leave-one-out estimates when removing inconsistent examples.\n\n");
  }    

  if(learn_parm->compute_loo && ((trainpos == 1) || (trainneg == 1))) {
    learn_parm->compute_loo=0;
    printf("\nCannot compute leave-one-out with only one example in one class.\n\n");
  }    


  if(verbosity==1) {
    printf("Optimizing"); fflush(stdout);
  }

  /* train the svm */
  iterations=optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
				     kernel_parm,kernel_cache,&shrink_state,model,
				     inconsistent,unlabeled,a,lin,
				     c,&timing_profile,
				     &maxdiff,(long)-1,
				     (long)1);
  
  if(verbosity>=1) {
    if(verbosity==1) printf("done. (%ld iterations)\n",iterations);

    misclassified=0;
    for(i=0;(i<totdoc);i++) { /* get final statistic */
      if((lin[i]-model->b)*(double)label[i] <= 0.0) 
	misclassified++;
    }

    printf("Optimization finished (%ld misclassified, maxdiff=%.5f).\n",
	   misclassified,maxdiff); 

    runtime_end=get_runtime();
    if(verbosity>=2) {
      printf("Runtime in cpu-seconds: %.2f (%.2f%% for kernel/%.2f%% for optimizer/%.2f%% for final/%.2f%% for update/%.2f%% for model/%.2f%% for check/%.2f%% for select)\n",
        ((float)runtime_end-(float)runtime_start)/100.0,
        (100.0*timing_profile.time_kernel)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_opti)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_shrink)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_update)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_model)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_check)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_select)/(float)(runtime_end-runtime_start));
    }
    else {
      printf("Runtime in cpu-seconds: %.2f\n",
	     (runtime_end-runtime_start)/100.0);
    }

    if(learn_parm->remove_inconsistent) {	  
      inconsistentnum=0;
      for(i=0;i<totdoc;i++) 
	if(inconsistent[i]) 
	  inconsistentnum++;
      printf("Number of SV: %ld (plus %ld inconsistent examples)\n",
	     model->sv_num-1,inconsistentnum);
    }
    else {
      upsupvecnum=0;
      for(i=1;i<model->sv_num;i++) {
	if(fabs(model->alpha[i]) >= 
	   (learn_parm->svm_cost[(model->supvec[i])->docnum]-
	    learn_parm->epsilon_a)) 
	  upsupvecnum++;
      }
      printf("Number of SV: %ld (including %ld at upper bound)\n",
	     model->sv_num-1,upsupvecnum);
    }
    
    if((verbosity>=1) && (!learn_parm->skip_final_opt_check)) {
      loss=0;
      model_length=0; 
      for(i=0;i<totdoc;i++) {
	if((lin[i]-model->b)*(double)label[i] < 1.0-learn_parm->epsilon_crit)
	  loss+=1.0-(lin[i]-model->b)*(double)label[i];
	model_length+=a[i]*label[i]*lin[i];
      }
      model_length=sqrt(model_length);
      fprintf(stdout,"L1 loss: loss=%.5f\n",loss);
      fprintf(stdout,"Norm of weight vector: |w|=%.5f\n",model_length);
      example_length=estimate_sphere(model,kernel_parm); 
      fprintf(stdout,"Norm of longest example vector: |x|=%.5f\n",
	      length_of_longest_document_vector(docs,totdoc,kernel_parm));
      fprintf(stdout,"Estimated VCdim of classifier: VCdim<=%.5f\n",
	      estimate_margin_vcdim(model,model_length,example_length,
				    kernel_parm));
      if((!learn_parm->remove_inconsistent) && (!transduction)) {
	runtime_start_xa=get_runtime();
	if(verbosity>=1) {
	  printf("Computing XiAlpha-estimates..."); fflush(stdout);
	}
	compute_xa_estimates(model,label,unlabeled,totdoc,docs,lin,a,
			     kernel_parm,learn_parm,&(model->xa_error),
			     &(model->xa_recall),&(model->xa_precision));
	if(verbosity>=1) {
	  printf("done\n");
	}
	printf("Runtime for XiAlpha-estimates in cpu-seconds: %.2f\n",
	       (get_runtime()-runtime_start_xa)/100.0);
	
	fprintf(stdout,"XiAlpha-estimate of the error: error<=%.2f%% (rho=%.2f,depth=%ld)\n",
		model->xa_error,learn_parm->rho,learn_parm->xa_depth);
	fprintf(stdout,"XiAlpha-estimate of the recall: recall=>%.2f%% (rho=%.2f,depth=%ld)\n",
		model->xa_recall,learn_parm->rho,learn_parm->xa_depth);
	fprintf(stdout,"XiAlpha-estimate of the precision: precision=>%.2f%% (rho=%.2f,depth=%ld)\n",
		model->xa_precision,learn_parm->rho,learn_parm->xa_depth);
      }
      else if(!learn_parm->remove_inconsistent) {
	estimate_transduction_quality(model,label,unlabeled,totdoc,docs,lin);
      }
    }
    if(verbosity>=1) {
      printf("Number of kernel evaluations: %ld\n",kernel_cache_statistic);
    }
  }


  /* leave-one-out testing starts now */
  if(learn_parm->compute_loo) {
    /* save results of training on full dataset for leave-one-out */
    runtime_start_loo=get_runtime();
    for(i=0;i<totdoc;i++) {
      xi_fullset[i]=1.0-((lin[i]-model->b)*(double)label[i]);
      if(xi_fullset[i]<0) xi_fullset[i]=0;
      a_fullset[i]=a[i];
    }
    if(verbosity>=1) {
      printf("Computing leave-one-out");
    }
    
    /* repeat this loop for every held-out example */
    for(heldout=0;(heldout<totdoc);heldout++) {
      if(learn_parm->rho*a_fullset[heldout]*r_delta_sq+xi_fullset[heldout]
	 < 1.0) { 
	/* guaranteed to not produce a leave-one-out error */
	if(verbosity==1) {
	  printf("+"); fflush(stdout); 
	}
      }
      else if(xi_fullset[heldout] > 1.0) {
	/* guaranteed to produce a leave-one-out error */
	loo_count++;
	if(label[heldout] > 0)  loo_count_pos++; else loo_count_neg++;
	if(verbosity==1) {
	  printf("-"); fflush(stdout); 
	}
      }
      else {
	loocomputed++;
	heldout_c=learn_parm->svm_cost[heldout]; /* set upper bound to zero */
	learn_parm->svm_cost[heldout]=0;
	/* make sure heldout example is not currently  */
	/* shrunk away. Assumes that lin is up to date! */
	shrink_state.active[heldout]=1;  
	if(verbosity>=2) 
	  printf("\nLeave-One-Out test on example %ld\n",heldout);
	if(verbosity>=1) {
	  printf("(?[%ld]",heldout); fflush(stdout); 
	}
	
	optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
				kernel_parm,
				kernel_cache,&shrink_state,model,inconsistent,unlabeled,
				a,lin,c,&timing_profile,
				&maxdiff,heldout,(long)2);

	/* printf("%.20f\n",(lin[heldout]-model->b)*(double)label[heldout]); */

	if(((lin[heldout]-model->b)*(double)label[heldout]) <= 0.0) { 
	  loo_count++;                            /* there was a loo-error */
	  if(label[heldout] > 0)  loo_count_pos++; else loo_count_neg++;
	  if(verbosity>=1) {
	    printf("-)"); fflush(stdout); 
	  }
	}
	else {
	  if(verbosity>=1) {
	    printf("+)"); fflush(stdout); 
	  }
	}
	/* now we need to restore the original data set*/
	learn_parm->svm_cost[heldout]=heldout_c; /* restore upper bound */
      }
    } /* end of leave-one-out loop */


    if(verbosity>=1) {
      printf("\nRetrain on full problem"); fflush(stdout); 
    }
    optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
			    kernel_parm,
			    kernel_cache,&shrink_state,model,inconsistent,unlabeled,
			    a,lin,c,&timing_profile,
			    &maxdiff,(long)-1,(long)1);
    if(verbosity >= 1) 
      printf("done.\n");
    
    
    /* after all leave-one-out computed */
    model->loo_error=100.0*loo_count/(double)totdoc;
    model->loo_recall=(1.0-(double)loo_count_pos/(double)trainpos)*100.0;
    model->loo_precision=(trainpos-loo_count_pos)/
      (double)(trainpos-loo_count_pos+loo_count_neg)*100.0;
    if(verbosity >= 1) {
      fprintf(stdout,"Leave-one-out estimate of the error: error=%.2f%%\n",
	      model->loo_error);
      fprintf(stdout,"Leave-one-out estimate of the recall: recall=%.2f%%\n",
	      model->loo_recall);
      fprintf(stdout,"Leave-one-out estimate of the precision: precision=%.2f%%\n",
	      model->loo_precision);
      fprintf(stdout,"Actual leave-one-outs computed:  %ld (rho=%.2f)\n",
	      loocomputed,learn_parm->rho);
      printf("Runtime for leave-one-out in cpu-seconds: %.2f\n",
	     (double)(get_runtime()-runtime_start_loo)/100.0);
    }
  }
    
  if(learn_parm->alphafile[0])
    write_alphas(learn_parm->alphafile,a,label,totdoc);
  
  shrink_state_cleanup(&shrink_state);
  free(label);
  free(inconsistent);
  free(unlabeled);
  free(c);
  free(a);
  free(a_fullset);
  free(xi_fullset);
  free(lin);
  free(learn_parm->svm_cost);
}


/* Learns an SVM regression model based on the training data in
   docs/label. The resulting model is returned in the structure
   model. */

void svm_learn_regression(DOC **docs, double *value, long int totdoc, 
			  long int totwords, LEARN_PARM *learn_parm, 
			  KERNEL_PARM *kernel_parm, 
			  KERNEL_CACHE **kernel_cache, MODEL *model)
     /* docs:        Training vectors (x-part) */
     /* class:       Training value (y-part) */
     /* totdoc:      Number of examples in docs/label */
     /* totwords:    Number of features (i.e. highest feature index) */
     /* learn_parm:  Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache:Initialized Cache, if using a kernel. NULL if
                     linear. Note that it will be free'd and reassigned */
     /* model:       Returns learning result (assumed empty before called) */
{
  long *inconsistent,i,j;
  long inconsistentnum;
  long upsupvecnum;
  double loss,model_length,example_length;
  double maxdiff,*lin,*a,*c;
  long runtime_start,runtime_end;
  long iterations,kernel_cache_size;
  long *unlabeled;
  double r_delta_sq=0,r_delta,r_delta_avg;
  double *xi_fullset; /* buffer for storing xi on full sample in loo */
  double *a_fullset;  /* buffer for storing alpha on full sample in loo */
  TIMING timing_profile;
  SHRINK_STATE shrink_state;
  DOC **docs_org;
  long *label;

  /* set up regression problem in standard form */
  docs_org=docs;
  docs = (DOC **)my_malloc(sizeof(DOC)*2*totdoc);
  label = (long *)my_malloc(sizeof(long)*2*totdoc);
  c = (double *)my_malloc(sizeof(double)*2*totdoc);
  for(i=0;i<totdoc;i++) {   
    j=2*totdoc-1-i;
    docs[i]=create_example(i,0,0,docs_org[i]->costfactor,docs_org[i]->fvec);
    label[i]=+1;
    c[i]=value[i];
    docs[j]=create_example(j,0,0,docs_org[i]->costfactor,docs_org[i]->fvec);
    label[j]=-1;
    c[j]=value[i];
  }
  totdoc*=2;

  /* need to get a bigger kernel cache */
  if(*kernel_cache) {
    kernel_cache_size=(*kernel_cache)->buffsize*sizeof(CFLOAT)/(1024*1024);
    kernel_cache_cleanup(*kernel_cache);
    (*kernel_cache)=kernel_cache_init(totdoc,kernel_cache_size);
  }

  runtime_start=get_runtime();
  timing_profile.time_kernel=0;
  timing_profile.time_opti=0;
  timing_profile.time_shrink=0;
  timing_profile.time_update=0;
  timing_profile.time_model=0;
  timing_profile.time_check=0;
  timing_profile.time_select=0;
  kernel_cache_statistic=0;

  learn_parm->totwords=totwords;

  /* make sure -n value is reasonable */
  if((learn_parm->svm_newvarsinqp < 2) 
     || (learn_parm->svm_newvarsinqp > learn_parm->svm_maxqpsize)) {
    learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;
  }

  init_shrink_state(&shrink_state,totdoc,(long)MAXSHRINK);

  inconsistent = (long *)my_malloc(sizeof(long)*totdoc);
  unlabeled = (long *)my_malloc(sizeof(long)*totdoc);
  a = (double *)my_malloc(sizeof(double)*totdoc);
  a_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  xi_fullset = (double *)my_malloc(sizeof(double)*totdoc);
  lin = (double *)my_malloc(sizeof(double)*totdoc);
  learn_parm->svm_cost = (double *)my_malloc(sizeof(double)*totdoc);
  model->supvec = (DOC **)my_malloc(sizeof(DOC *)*(totdoc+2));
  model->alpha = (double *)my_malloc(sizeof(double)*(totdoc+2));
  model->index = (long *)my_malloc(sizeof(long)*(totdoc+2));

  model->at_upper_bound=0;
  model->b=0;	       
  model->supvec[0]=0;  /* element 0 reserved and empty for now */
  model->alpha[0]=0;
  model->lin_weights=NULL;
  model->totwords=totwords;
  model->totdoc=totdoc;
  model->kernel_parm=(*kernel_parm);
  model->sv_num=1;
  model->loo_error=-1;
  model->loo_recall=-1;
  model->loo_precision=-1;
  model->xa_error=-1;
  model->xa_recall=-1;
  model->xa_precision=-1;
  inconsistentnum=0;

  r_delta=estimate_r_delta(docs,totdoc,kernel_parm);
  r_delta_sq=r_delta*r_delta;

  r_delta_avg=estimate_r_delta_average(docs,totdoc,kernel_parm);
  if(learn_parm->svm_c == 0.0) {  /* default value for C */
    learn_parm->svm_c=1.0/(r_delta_avg*r_delta_avg);
    if(verbosity>=1) 
      printf("Setting default regularization parameter C=%.4f\n",
	     learn_parm->svm_c);
  }

  for(i=0;i<totdoc;i++) {    /* various inits */
    inconsistent[i]=0;
    a[i]=0;
    lin[i]=0;
    unlabeled[i]=0;
    if(label[i] > 0) {
      learn_parm->svm_cost[i]=learn_parm->svm_c*learn_parm->svm_costratio*
	docs[i]->costfactor;
    }
    else if(label[i] < 0) {
      learn_parm->svm_cost[i]=learn_parm->svm_c*docs[i]->costfactor;
    }
  }

  /* caching makes no sense for linear kernel */
  if((kernel_parm->kernel_type == LINEAR) && (*kernel_cache)) {
    printf("WARNING: Using a kernel cache for linear case will slow optimization down!\n");
  } 

  if(verbosity==1) {
    printf("Optimizing"); fflush(stdout);
  }

  /* train the svm */
  iterations=optimize_to_convergence(docs,label,totdoc,totwords,learn_parm,
				     kernel_parm,*kernel_cache,&shrink_state,
				     model,inconsistent,unlabeled,a,lin,c,
				     &timing_profile,&maxdiff,(long)-1,
				     (long)1);
  
  if(verbosity>=1) {
    if(verbosity==1) printf("done. (%ld iterations)\n",iterations);

    printf("Optimization finished (maxdiff=%.5f).\n",maxdiff); 

    runtime_end=get_runtime();
    if(verbosity>=2) {
      printf("Runtime in cpu-seconds: %.2f (%.2f%% for kernel/%.2f%% for optimizer/%.2f%% for final/%.2f%% for update/%.2f%% for model/%.2f%% for check/%.2f%% for select)\n",
        ((float)runtime_end-(float)runtime_start)/100.0,
        (100.0*timing_profile.time_kernel)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_opti)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_shrink)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_update)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_model)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_check)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_select)/(float)(runtime_end-runtime_start));
    }
    else {
      printf("Runtime in cpu-seconds: %.2f\n",
	     (runtime_end-runtime_start)/100.0);
    }

    if(learn_parm->remove_inconsistent) {	  
      inconsistentnum=0;
      for(i=0;i<totdoc;i++) 
	if(inconsistent[i]) 
	  inconsistentnum++;
      printf("Number of SV: %ld (plus %ld inconsistent examples)\n",
	     model->sv_num-1,inconsistentnum);
    }
    else {
      upsupvecnum=0;
      for(i=1;i<model->sv_num;i++) {
	if(fabs(model->alpha[i]) >= 
	   (learn_parm->svm_cost[(model->supvec[i])->docnum]-
	    learn_parm->epsilon_a)) 
	  upsupvecnum++;
      }
      printf("Number of SV: %ld (including %ld at upper bound)\n",
	     model->sv_num-1,upsupvecnum);
    }
    
    if((verbosity>=1) && (!learn_parm->skip_final_opt_check)) {
      loss=0;
      model_length=0; 
      for(i=0;i<totdoc;i++) {
	if((lin[i]-model->b)*(double)label[i] < (-learn_parm->eps+(double)label[i]*c[i])-learn_parm->epsilon_crit)
	  loss+=-learn_parm->eps+(double)label[i]*c[i]-(lin[i]-model->b)*(double)label[i];
	model_length+=a[i]*label[i]*lin[i];
      }
      model_length=sqrt(model_length);
      fprintf(stdout,"L1 loss: loss=%.5f\n",loss);
      fprintf(stdout,"Norm of weight vector: |w|=%.5f\n",model_length);
      example_length=estimate_sphere(model,kernel_parm); 
      fprintf(stdout,"Norm of longest example vector: |x|=%.5f\n",
	      length_of_longest_document_vector(docs,totdoc,kernel_parm));
    }
    if(verbosity>=1) {
      printf("Number of kernel evaluations: %ld\n",kernel_cache_statistic);
    }
  }
    
  if(learn_parm->alphafile[0])
    write_alphas(learn_parm->alphafile,a,label,totdoc);

  /* this makes sure the model we return does not contain pointers to the 
     temporary documents */
  for(i=1;i<model->sv_num;i++) { 
    j=model->supvec[i]->docnum;
    if(j >= (totdoc/2)) {
      j=totdoc-j-1;
    }
    model->supvec[i]=docs_org[j];
  }
  
  shrink_state_cleanup(&shrink_state);
  for(i=0;i<totdoc;i++)
    free_example(docs[i],0);
  free(docs);
  free(label);
  free(inconsistent);
  free(unlabeled);
  free(c);
  free(a);
  free(a_fullset);
  free(xi_fullset);
  free(lin);
  free(learn_parm->svm_cost);
}

void svm_learn_ranking(DOC **docs, double *rankvalue, long int totdoc, 
		       long int totwords, LEARN_PARM *learn_parm, 
		       KERNEL_PARM *kernel_parm, KERNEL_CACHE **kernel_cache, 
		       MODEL *model)
     /* docs:        Training vectors (x-part) */
     /* rankvalue:   Training target values that determine the ranking */
     /* totdoc:      Number of examples in docs/label */
     /* totwords:    Number of features (i.e. highest feature index) */
     /* learn_parm:  Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache:Initialized pointer to Cache of size 1*totdoc, if 
	             using a kernel. NULL if linear. NOTE: Cache is 
                     getting reinitialized in this function */
     /* model:       Returns learning result (assumed empty before called) */
{
  DOC **docdiff;
  long i,j,k,totpair,kernel_cache_size;
  double *target,*alpha,cost;
  long *greater,*lesser;
  MODEL *pairmodel;
  SVECTOR *flow,*fhigh;

  totpair=0;
  for(i=0;i<totdoc;i++) {
    for(j=i+1;j<totdoc;j++) {
      if((docs[i]->queryid==docs[j]->queryid) && (rankvalue[i] != rankvalue[j])) {
	totpair++;
      }
    }
  }

  printf("Constructing %ld rank constraints...",totpair); fflush(stdout);
  docdiff=(DOC **)my_malloc(sizeof(DOC)*totpair);
  target=(double *)my_malloc(sizeof(double)*totpair); 
  greater=(long *)my_malloc(sizeof(long)*totpair); 
  lesser=(long *)my_malloc(sizeof(long)*totpair); 

  k=0;
  for(i=0;i<totdoc;i++) {
    for(j=i+1;j<totdoc;j++) {
      if(docs[i]->queryid == docs[j]->queryid) {
	cost=(docs[i]->costfactor+docs[j]->costfactor)/2.0;
	if(rankvalue[i] > rankvalue[j]) {
	  if(kernel_parm->kernel_type == LINEAR)
	    docdiff[k]=create_example(k,0,0,cost,
				      sub_ss(docs[i]->fvec,docs[j]->fvec));
	  else {
	    flow=copy_svector(docs[j]->fvec);
	    flow->factor=-1.0;
	    flow->next=NULL;
	    fhigh=copy_svector(docs[i]->fvec);
	    fhigh->factor=1.0;
	    fhigh->next=flow;
	    docdiff[k]=create_example(k,0,0,cost,fhigh);
	  }
	  target[k]=1;
	  greater[k]=i;
	  lesser[k]=j;
	  k++;
	}
	else if(rankvalue[i] < rankvalue[j]) {
	  if(kernel_parm->kernel_type == LINEAR)
	    docdiff[k]=create_example(k,0,0,cost,
				      sub_ss(docs[i]->fvec,docs[j]->fvec));
	  else {
	    flow=copy_svector(docs[j]->fvec);
	    flow->factor=-1.0;
	    flow->next=NULL;
	    fhigh=copy_svector(docs[i]->fvec);
	    fhigh->factor=1.0;
	    fhigh->next=flow;
	    docdiff[k]=create_example(k,0,0,cost,fhigh);
	  }
	  target[k]=-1;
	  greater[k]=i;
	  lesser[k]=j;
	  k++;
	}
      }
    }
  }
  printf("done.\n"); fflush(stdout);

  /* need to get a bigger kernel cache */
  if(*kernel_cache) {
    kernel_cache_size=(*kernel_cache)->buffsize*sizeof(CFLOAT)/(1024*1024);
    kernel_cache_cleanup(*kernel_cache);
    (*kernel_cache)=kernel_cache_init(totpair,kernel_cache_size);
  }

  /* must use unbiased hyperplane on difference vectors */
  learn_parm->biased_hyperplane=0;
  pairmodel=(MODEL *)my_malloc(sizeof(MODEL));
  svm_learn_classification(docdiff,target,totpair,totwords,learn_parm,
			   kernel_parm,(*kernel_cache),pairmodel,NULL);

  /* Transfer the result into a more compact model. If you would like
     to output the original model on pairs of documents, see below. */
  alpha=(double *)my_malloc(sizeof(double)*totdoc); 
  for(i=0;i<totdoc;i++) {
    alpha[i]=0;
  }
  for(i=1;i<pairmodel->sv_num;i++) {
    alpha[lesser[(pairmodel->supvec[i])->docnum]]-=pairmodel->alpha[i];
    alpha[greater[(pairmodel->supvec[i])->docnum]]+=pairmodel->alpha[i];
  }
  model->supvec = (DOC **)my_malloc(sizeof(DOC *)*(totdoc+2));
  model->alpha = (double *)my_malloc(sizeof(double)*(totdoc+2));
  model->index = (long *)my_malloc(sizeof(long)*(totdoc+2));
  model->supvec[0]=0;  /* element 0 reserved and empty for now */
  model->alpha[0]=0;
  model->sv_num=1;
  for(i=0;i<totdoc;i++) {
    if(alpha[i]) {
      model->supvec[model->sv_num]=docs[i];
      model->alpha[model->sv_num]=alpha[i];
      model->index[i]=model->sv_num;
      model->sv_num++;
    }
    else {
      model->index[i]=-1;
    }
  }
  model->at_upper_bound=0;
  model->b=0;	       
  model->lin_weights=NULL;
  model->totwords=totwords;
  model->totdoc=totdoc;
  model->kernel_parm=(*kernel_parm);
  model->loo_error=-1;
  model->loo_recall=-1;
  model->loo_precision=-1;
  model->xa_error=-1;
  model->xa_recall=-1;
  model->xa_precision=-1;

  free(alpha);
  free(greater);
  free(lesser);
  free(target);

  /* If you would like to output the original model on pairs of
     document, replace the following lines with '(*model)=(*pairmodel);' */
  for(i=0;i<totpair;i++)
    free_example(docdiff[i],1);
  free(docdiff);
  free_model(pairmodel,0);
}


/* The following solves a freely defined and given set of
   inequalities. The optimization problem is of the following form:

   min 0.5 w*w + C sum_i C_i \xi_i
   s.t. x_i * w > rhs_i - \xi_i

   This corresponds to the -z o option. */

void svm_learn_optimization(DOC **docs, double *rhs, long int
			    totdoc, long int totwords, 
			    LEARN_PARM *learn_parm, 
			    KERNEL_PARM *kernel_parm, 
			    KERNEL_CACHE *kernel_cache, MODEL *model,
			    double *alpha)
     /* docs:        Left-hand side of inequalities (x-part) */
     /* rhs:         Right-hand side of inequalities */
     /* totdoc:      Number of examples in docs/label */
     /* totwords:    Number of features (i.e. highest feature index) */
     /* learn_parm:  Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache:Initialized Cache of size 1*totdoc, if using a kernel. 
                     NULL if linear.*/
     /* model:       Returns solution as SV expansion (assumed empty before called) */
     /* alpha:       Start values for the alpha variables or NULL
	             pointer. The new alpha values are returned after 
		     optimization if not NULL. Array must be of size totdoc. */
{
  long i,*label;
  long misclassified,upsupvecnum;
  double loss,model_length,example_length;
  double maxdiff,*lin,*a,*c;
  long runtime_start,runtime_end;
  long iterations,maxslackid,svsetnum;
  long *unlabeled,*inconsistent;
  double r_delta_sq=0,r_delta,r_delta_avg;
  long *index,*index2dnum;
  double *weights,*slack,*alphaslack;
  CFLOAT *aicache;  /* buffer to keep one row of hessian */

  TIMING timing_profile;
  SHRINK_STATE shrink_state;

  runtime_start=get_runtime();
  timing_profile.time_kernel=0;
  timing_profile.time_opti=0;
  timing_profile.time_shrink=0;
  timing_profile.time_update=0;
  timing_profile.time_model=0;
  timing_profile.time_check=0;
  timing_profile.time_select=0;
  kernel_cache_statistic=0;

  learn_parm->totwords=totwords;

  /* make sure -n value is reasonable */
  if((learn_parm->svm_newvarsinqp < 2) 
     || (learn_parm->svm_newvarsinqp > learn_parm->svm_maxqpsize)) {
    learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;
  }

  init_shrink_state(&shrink_state,totdoc,(long)MAXSHRINK);

  label = (long *)my_malloc(sizeof(long)*totdoc);
  unlabeled = (long *)my_malloc(sizeof(long)*totdoc);
  inconsistent = (long *)my_malloc(sizeof(long)*totdoc);
  c = (double *)my_malloc(sizeof(double)*totdoc);
  a = (double *)my_malloc(sizeof(double)*totdoc);
  lin = (double *)my_malloc(sizeof(double)*totdoc);
  learn_parm->svm_cost = (double *)my_malloc(sizeof(double)*totdoc);
  model->supvec = (DOC **)my_malloc(sizeof(DOC *)*(totdoc+2));
  model->alpha = (double *)my_malloc(sizeof(double)*(totdoc+2));
  model->index = (long *)my_malloc(sizeof(long)*(totdoc+2));

  model->at_upper_bound=0;
  model->b=0;	       
  model->supvec[0]=0;  /* element 0 reserved and empty for now */
  model->alpha[0]=0;
  model->lin_weights=NULL;
  model->totwords=totwords;
  model->totdoc=totdoc;
  model->kernel_parm=(*kernel_parm);
  model->sv_num=1;
  model->loo_error=-1;
  model->loo_recall=-1;
  model->loo_precision=-1;
  model->xa_error=-1;
  model->xa_recall=-1;
  model->xa_precision=-1;

  r_delta=estimate_r_delta(docs,totdoc,kernel_parm);
  r_delta_sq=r_delta*r_delta;

  r_delta_avg=estimate_r_delta_average(docs,totdoc,kernel_parm);
  if(learn_parm->svm_c == 0.0) {  /* default value for C */
    learn_parm->svm_c=1.0/(r_delta_avg*r_delta_avg);
    if(verbosity>=1) 
      printf("Setting default regularization parameter C=%.4f\n",
	     learn_parm->svm_c);
  }

  learn_parm->biased_hyperplane=0; /* learn an unbiased hyperplane */

  learn_parm->eps=0.0;      /* No margin, unless explicitly handcoded
                               in the right-hand side in the training
                               set.  */

  for(i=0;i<totdoc;i++) {    /* various inits */
    docs[i]->docnum=i;
    a[i]=0;
    lin[i]=0;
    c[i]=rhs[i];       /* set right-hand side */
    unlabeled[i]=0;
    inconsistent[i]=0;
    learn_parm->svm_cost[i]=learn_parm->svm_c*learn_parm->svm_costratio*
      docs[i]->costfactor;
    label[i]=1;
  }
  if(learn_parm->sharedslack) /* if shared slacks are used, they must */
    for(i=0;i<totdoc;i++)     /*  be used on every constraint */
      if(!docs[i]->slackid) {
	perror("Error: Missing shared slacks definitions in some of the examples.");
	exit(0);
      }
      
  /* compute starting state for initial alpha values */
  if(alpha) {
    if(verbosity>=1) {
      printf("Computing starting state..."); fflush(stdout);
    }
    index = (long *)my_malloc(sizeof(long)*totdoc);
    index2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
    weights=(double *)my_malloc(sizeof(double)*(totwords+1));
    aicache = (CFLOAT *)my_malloc(sizeof(CFLOAT)*totdoc);
    for(i=0;i<totdoc;i++) {    /* create full index and clip alphas */
      index[i]=1;
      alpha[i]=fabs(alpha[i]);
      if(alpha[i]<0) alpha[i]=0;
      if(alpha[i]>learn_parm->svm_cost[i]) alpha[i]=learn_parm->svm_cost[i];
    }
    if(kernel_parm->kernel_type != LINEAR) {
      for(i=0;i<totdoc;i++)     /* fill kernel cache with unbounded SV */
	if((alpha[i]>0) && (alpha[i]<learn_parm->svm_cost[i]) 
	   && (kernel_cache_space_available(kernel_cache))) 
	  cache_kernel_row(kernel_cache,docs,i,kernel_parm);
      for(i=0;i<totdoc;i++)     /* fill rest of kernel cache with bounded SV */
	if((alpha[i]==learn_parm->svm_cost[i]) 
	   && (kernel_cache_space_available(kernel_cache))) 
	  cache_kernel_row(kernel_cache,docs,i,kernel_parm);
    }
    (void)compute_index(index,totdoc,index2dnum);
    update_linear_component(docs,label,index2dnum,alpha,a,index2dnum,totdoc,
			    totwords,kernel_parm,kernel_cache,lin,aicache,
			    weights);
    (void)calculate_svm_model(docs,label,unlabeled,lin,alpha,a,c,
			      learn_parm,index2dnum,index2dnum,model);
    for(i=0;i<totdoc;i++) {    /* copy initial alphas */
      a[i]=alpha[i];
    }
    free(index);
    free(index2dnum);
    free(weights);
    free(aicache);
    if(verbosity>=1) {
      printf("done.\n");  fflush(stdout);
    }   
  } 

  /* removing inconsistent does not work for general optimization problem */
  if(learn_parm->remove_inconsistent) {	  
    learn_parm->remove_inconsistent = 0;
    printf("'remove inconsistent' not available in this mode. Switching option off!"); fflush(stdout);
  }

  /* caching makes no sense for linear kernel */
  if(kernel_parm->kernel_type == LINEAR) {
    kernel_cache = NULL;   
  } 

  if(verbosity==1) {
    printf("Optimizing"); fflush(stdout);
  }

  /* train the svm */
  if(learn_parm->sharedslack)
    iterations=optimize_to_convergence_sharedslack(docs,label,totdoc,
				     totwords,learn_parm,kernel_parm,
				     kernel_cache,&shrink_state,model,
				     a,lin,c,&timing_profile,
				     &maxdiff);
  else
    iterations=optimize_to_convergence(docs,label,totdoc,
				     totwords,learn_parm,kernel_parm,
				     kernel_cache,&shrink_state,model,
				     inconsistent,unlabeled,
				     a,lin,c,&timing_profile,
				     &maxdiff,(long)-1,(long)1);
  
  if(verbosity>=1) {
    if(verbosity==1) printf("done. (%ld iterations)\n",iterations);

    misclassified=0;
    for(i=0;(i<totdoc);i++) { /* get final statistic */
      if((lin[i]-model->b)*(double)label[i] <= 0.0) 
	misclassified++;
    }

    printf("Optimization finished (maxdiff=%.5f).\n",maxdiff); 

    runtime_end=get_runtime();
    if(verbosity>=2) {
      printf("Runtime in cpu-seconds: %.2f (%.2f%% for kernel/%.2f%% for optimizer/%.2f%% for final/%.2f%% for update/%.2f%% for model/%.2f%% for check/%.2f%% for select)\n",
        ((float)runtime_end-(float)runtime_start)/100.0,
        (100.0*timing_profile.time_kernel)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_opti)/(float)(runtime_end-runtime_start),
	(100.0*timing_profile.time_shrink)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_update)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_model)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_check)/(float)(runtime_end-runtime_start),
        (100.0*timing_profile.time_select)/(float)(runtime_end-runtime_start));
    }
    else {
      printf("Runtime in cpu-seconds: %.2f\n",
	     (runtime_end-runtime_start)/100.0);
    }
  }
  if((verbosity>=1) && (!learn_parm->skip_final_opt_check)) {
    loss=0;
    model_length=0; 
    for(i=0;i<totdoc;i++) {
      if((lin[i]-model->b)*(double)label[i] < c[i]-learn_parm->epsilon_crit)
	loss+=c[i]-(lin[i]-model->b)*(double)label[i];
      model_length+=a[i]*label[i]*lin[i];
    }
    model_length=sqrt(model_length);
    fprintf(stdout,"Norm of weight vector: |w|=%.5f\n",model_length);
  }
  
  if(learn_parm->sharedslack) {
    index = (long *)my_malloc(sizeof(long)*totdoc);
    index2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
    maxslackid=0;
    for(i=0;i<totdoc;i++) {    /* create full index */
      index[i]=1;
      if(maxslackid<docs[i]->slackid)
	maxslackid=docs[i]->slackid;
    }
    (void)compute_index(index,totdoc,index2dnum);
    slack=(double *)my_malloc(sizeof(double)*(maxslackid+1));
    alphaslack=(double *)my_malloc(sizeof(double)*(maxslackid+1));
    for(i=0;i<=maxslackid;i++) {    /* init shared slacks */
      slack[i]=0;
      alphaslack[i]=0;
    }
    compute_shared_slacks(docs,label,a,lin,c,index2dnum,learn_parm,
			  slack,alphaslack);
    loss=0;
    model->at_upper_bound=0;
    svsetnum=0;
    for(i=0;i<=maxslackid;i++) {    /* create full index */
      loss+=slack[i];
      if(alphaslack[i] > (learn_parm->svm_c - learn_parm->epsilon_a)) 
	model->at_upper_bound++;
      if(alphaslack[i] > learn_parm->epsilon_a)
	svsetnum++;
    }
    free(index);
    free(index2dnum);
    free(slack);
    free(alphaslack);
  }
  
  if((verbosity>=1) && (!learn_parm->skip_final_opt_check)) {
    if(learn_parm->sharedslack) {
      printf("Number of SV: %ld\n",
	     model->sv_num-1);
      printf("Number of non-zero slack variables: %ld (out of %ld)\n",
	     model->at_upper_bound,svsetnum);
      fprintf(stdout,"L1 loss: loss=%.5f\n",loss);
    }
    else {
      upsupvecnum=0;
      for(i=1;i<model->sv_num;i++) {
	if(fabs(model->alpha[i]) >= 
	   (learn_parm->svm_cost[(model->supvec[i])->docnum]-
	    learn_parm->epsilon_a)) 
	  upsupvecnum++;
      }
      printf("Number of SV: %ld (including %ld at upper bound)\n",
	     model->sv_num-1,upsupvecnum);
      fprintf(stdout,"L1 loss: loss=%.5f\n",loss);
    }
    example_length=estimate_sphere(model,kernel_parm); 
    fprintf(stdout,"Norm of longest example vector: |x|=%.5f\n",
	    length_of_longest_document_vector(docs,totdoc,kernel_parm));
  }
  if(verbosity>=1) {
    printf("Number of kernel evaluations: %ld\n",kernel_cache_statistic);
  }
    
  if(alpha) {
    for(i=0;i<totdoc;i++) {    /* copy final alphas */
      alpha[i]=a[i];
    }
  }
 
  if(learn_parm->alphafile[0])
    write_alphas(learn_parm->alphafile,a,label,totdoc);
  
  shrink_state_cleanup(&shrink_state);
  free(label);
  free(unlabeled);
  free(inconsistent);
  free(c);
  free(a);
  free(lin);
  free(learn_parm->svm_cost);
}


long optimize_to_convergence(DOC **docs, long int *label, long int totdoc, 
			     long int totwords, LEARN_PARM *learn_parm, 
			     KERNEL_PARM *kernel_parm, 
			     KERNEL_CACHE *kernel_cache, 
			     SHRINK_STATE *shrink_state, MODEL *model, 
			     long int *inconsistent, long int *unlabeled, 
			     double *a, double *lin, double *c, 
			     TIMING *timing_profile, double *maxdiff, 
			     long int heldout, long int retrain)
     /* docs: Training vectors (x-part) */
     /* label: Training labels/value (y-part, zero if test example for
			      transduction) */
     /* totdoc: Number of examples in docs/label */
     /* totwords: Number of features (i.e. highest feature index) */
     /* laern_parm: Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache: Initialized/partly filled Cache, if using a kernel. 
                      NULL if linear. */
     /* shrink_state: State of active variables */
     /* model: Returns learning result */
     /* inconsistent: examples thrown out as inconstistent */
     /* unlabeled: test examples for transduction */
     /* a: alphas */
     /* lin: linear component of gradient */
     /* c: right hand side of inequalities (margin) */
     /* maxdiff: returns maximum violation of KT-conditions */
     /* heldout: marks held-out example for leave-one-out (or -1) */
     /* retrain: selects training mode (1=regular / 2=holdout) */
{
  long *chosen,*key,i,j,jj,*last_suboptimal_at,noshrink;
  long inconsistentnum,choosenum,already_chosen=0,iteration;
  long misclassified,supvecnum=0,*active2dnum,inactivenum;
  long *working2dnum,*selexam;
  long activenum;
  double criterion,eq;
  double *a_old;
  long t0=0,t1=0,t2=0,t3=0,t4=0,t5=0,t6=0; /* timing */
  long transductcycle;
  long transduction;
  double epsilon_crit_org; 
  double bestmaxdiff;
  long   bestmaxdiffiter,terminate;

  double *selcrit;  /* buffer for sorting */        
  CFLOAT *aicache;  /* buffer to keep one row of hessian */
  double *weights;  /* buffer for weight vector in linear case */
  QP qp;            /* buffer for one quadratic program */

  epsilon_crit_org=learn_parm->epsilon_crit; /* save org */
  if(kernel_parm->kernel_type == LINEAR) {
    learn_parm->epsilon_crit=2.0;
    kernel_cache=NULL;   /* caching makes no sense for linear kernel */
  } 
  learn_parm->epsilon_shrink=2;
  (*maxdiff)=1;

  learn_parm->totwords=totwords;

  chosen = (long *)my_malloc(sizeof(long)*totdoc);
  last_suboptimal_at = (long *)my_malloc(sizeof(long)*totdoc);
  key = (long *)my_malloc(sizeof(long)*(totdoc+11)); 
  selcrit = (double *)my_malloc(sizeof(double)*totdoc);
  selexam = (long *)my_malloc(sizeof(long)*totdoc);
  a_old = (double *)my_malloc(sizeof(double)*totdoc);
  aicache = (CFLOAT *)my_malloc(sizeof(CFLOAT)*totdoc);
  working2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  active2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  qp.opt_ce = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_ce0 = (double *)my_malloc(sizeof(double));
  qp.opt_g = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize
				 *learn_parm->svm_maxqpsize);
  qp.opt_g0 = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_xinit = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_low=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_up=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  weights=(double *)my_malloc(sizeof(double)*(totwords+1));

  choosenum=0;
  inconsistentnum=0;
  transductcycle=0;
  transduction=0;
  if(!retrain) retrain=1;
  iteration=1;
  bestmaxdiffiter=1;
  bestmaxdiff=999999999;
  terminate=0;

  if(kernel_cache) {
    kernel_cache->time=iteration;  /* for lru cache */
    kernel_cache_reset_lru(kernel_cache);
  }

  for(i=0;i<totdoc;i++) {    /* various inits */
    chosen[i]=0;
    a_old[i]=a[i];
    last_suboptimal_at[i]=1;
    if(inconsistent[i]) 
      inconsistentnum++;
    if(unlabeled[i]) {
      transduction=1;
    }
  }
  activenum=compute_index(shrink_state->active,totdoc,active2dnum);
  inactivenum=totdoc-activenum;
  clear_index(working2dnum);

                            /* repeat this loop until we have convergence */
  for(;retrain && (!terminate);iteration++) {

    if(kernel_cache)
      kernel_cache->time=iteration;  /* for lru cache */
    if(verbosity>=2) {
      printf(
	"Iteration %ld: ",iteration); fflush(stdout);
    }
    else if(verbosity==1) {
      printf("."); fflush(stdout);
    }

    if(verbosity>=2) t0=get_runtime();
    if(verbosity>=3) {
      printf("\nSelecting working set... "); fflush(stdout); 
    }

    if(learn_parm->svm_newvarsinqp>learn_parm->svm_maxqpsize) 
      learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;

    i=0;
    for(jj=0;(j=working2dnum[jj])>=0;jj++) { /* clear working set */
      if((chosen[j]>=(learn_parm->svm_maxqpsize/
		      minl(learn_parm->svm_maxqpsize,
			   learn_parm->svm_newvarsinqp))) 
	 || (inconsistent[j])
	 || (j == heldout)) {
	chosen[j]=0; 
	choosenum--; 
      }
      else {
	chosen[j]++;
	working2dnum[i++]=j;
      }
    }
    working2dnum[i]=-1;

    if(retrain == 2) {
      choosenum=0;
      for(jj=0;(j=working2dnum[jj])>=0;jj++) { /* fully clear working set */
	chosen[j]=0; 
      }
      clear_index(working2dnum);
      for(i=0;i<totdoc;i++) { /* set inconsistent examples to zero (-i 1) */
	if((inconsistent[i] || (heldout==i)) && (a[i] != 0.0)) {
	  chosen[i]=99999;
	  choosenum++;
	  a[i]=0;
	}
      }
      if(learn_parm->biased_hyperplane) {
	eq=0;
	for(i=0;i<totdoc;i++) { /* make sure we fulfill equality constraint */
	  eq+=a[i]*label[i];
	}
	for(i=0;(i<totdoc) && (fabs(eq) > learn_parm->epsilon_a);i++) {
	  if((eq*label[i] > 0) && (a[i] > 0)) {
	    chosen[i]=88888;
	    choosenum++;
	    if((eq*label[i]) > a[i]) {
	      eq-=(a[i]*label[i]);
	      a[i]=0;
	    }
	    else {
	      a[i]-=(eq*label[i]);
	      eq=0;
	    }
	  }
	}
      }
      compute_index(chosen,totdoc,working2dnum);
    }
    else {      /* select working set according to steepest gradient */
      if(iteration % 101) {
        already_chosen=0;
	if((minl(learn_parm->svm_newvarsinqp,
		 learn_parm->svm_maxqpsize-choosenum)>=4) 
	   && (kernel_parm->kernel_type != LINEAR)) {
	  /* select part of the working set from cache */
	  already_chosen=select_next_qp_subproblem_grad(
			      label,unlabeled,a,lin,c,totdoc,
			      (long)(minl(learn_parm->svm_maxqpsize-choosenum,
					  learn_parm->svm_newvarsinqp)
				     /2),
			      learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,1,
			      key,chosen);
	  choosenum+=already_chosen;
	}
	choosenum+=select_next_qp_subproblem_grad(
                              label,unlabeled,a,lin,c,totdoc,
                              minl(learn_parm->svm_maxqpsize-choosenum,
				   learn_parm->svm_newvarsinqp-already_chosen),
                              learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,0,key,
			      chosen);
      }
      else { /* once in a while, select a somewhat random working set
		to get unlocked of infinite loops due to numerical
		inaccuracies in the core qp-solver */
	choosenum+=select_next_qp_subproblem_rand(
                              label,unlabeled,a,lin,c,totdoc,
                              minl(learn_parm->svm_maxqpsize-choosenum,
				   learn_parm->svm_newvarsinqp),
                              learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,key,
			      chosen,iteration);
      }
    }

    if(verbosity>=2) {
      printf(" %ld vectors chosen\n",choosenum); fflush(stdout); 
    }

    if(verbosity>=2) t1=get_runtime();

    if(kernel_cache) 
      cache_multiple_kernel_rows(kernel_cache,docs,working2dnum,
				 choosenum,kernel_parm); 
    
    if(verbosity>=2) t2=get_runtime();
    if(retrain != 2) {
      optimize_svm(docs,label,unlabeled,inconsistent,0.0,chosen,active2dnum,
		   model,totdoc,working2dnum,choosenum,a,lin,c,learn_parm,
		   aicache,kernel_parm,&qp,&epsilon_crit_org);
    }

    if(verbosity>=2) t3=get_runtime();
    update_linear_component(docs,label,active2dnum,a,a_old,working2dnum,totdoc,
			    totwords,kernel_parm,kernel_cache,lin,aicache,
			    weights);

    if(verbosity>=2) t4=get_runtime();
    supvecnum=calculate_svm_model(docs,label,unlabeled,lin,a,a_old,c,
		                  learn_parm,working2dnum,active2dnum,model);

    if(verbosity>=2) t5=get_runtime();

    /* The following computation of the objective function works only */
    /* relative to the active variables */
    if(verbosity>=3) {
      criterion=compute_objective_function(a,lin,c,learn_parm->eps,label,
		                           active2dnum);
      printf("Objective function (over active variables): %.16f\n",criterion);
      fflush(stdout); 
    }

    for(jj=0;(i=working2dnum[jj])>=0;jj++) {
      a_old[i]=a[i];
    }

    if(retrain == 2) {  /* reset inconsistent unlabeled examples */
      for(i=0;(i<totdoc);i++) {
	if(inconsistent[i] && unlabeled[i]) {
	  inconsistent[i]=0;
	  label[i]=0;
	}
      }
    }

    retrain=check_optimality(model,label,unlabeled,a,lin,c,totdoc,learn_parm,
			     maxdiff,epsilon_crit_org,&misclassified,
			     inconsistent,active2dnum,last_suboptimal_at,
			     iteration,kernel_parm);

    if(verbosity>=2) {
      t6=get_runtime();
      timing_profile->time_select+=t1-t0;
      timing_profile->time_kernel+=t2-t1;
      timing_profile->time_opti+=t3-t2;
      timing_profile->time_update+=t4-t3;
      timing_profile->time_model+=t5-t4;
      timing_profile->time_check+=t6-t5;
    }

    /* checking whether optimizer got stuck */
    if((*maxdiff) < bestmaxdiff) {
      bestmaxdiff=(*maxdiff);
      bestmaxdiffiter=iteration;
    }
    if(iteration > (bestmaxdiffiter+learn_parm->maxiter)) { 
      /* long time no progress? */
      terminate=1;
      retrain=0;
      if(verbosity>=1) 
	printf("\nWARNING: Relaxing KT-Conditions due to slow progress! Terminating!\n");
    }

    noshrink=0;
    if((!retrain) && (inactivenum>0) 
       && ((!learn_parm->skip_final_opt_check) 
	   || (kernel_parm->kernel_type == LINEAR))) { 
      if(((verbosity>=1) && (kernel_parm->kernel_type != LINEAR)) 
	 || (verbosity>=2)) {
	if(verbosity==1) {
	  printf("\n");
	}
	printf(" Checking optimality of inactive variables..."); 
	fflush(stdout);
      }
      t1=get_runtime();
      reactivate_inactive_examples(label,unlabeled,a,shrink_state,lin,c,totdoc,
				   totwords,iteration,learn_parm,inconsistent,
				   docs,kernel_parm,kernel_cache,model,aicache,
				   weights,maxdiff);
      /* Update to new active variables. */
      activenum=compute_index(shrink_state->active,totdoc,active2dnum);
      inactivenum=totdoc-activenum;
      /* reset watchdog */
      bestmaxdiff=(*maxdiff);
      bestmaxdiffiter=iteration;
      /* termination criterion */
      noshrink=1;
      retrain=0;
      if((*maxdiff) > learn_parm->epsilon_crit) 
	retrain=1;
      timing_profile->time_shrink+=get_runtime()-t1;
      if(((verbosity>=1) && (kernel_parm->kernel_type != LINEAR)) 
	 || (verbosity>=2)) {
	printf("done.\n");  fflush(stdout);
        printf(" Number of inactive variables = %ld\n",inactivenum);
      }		  
    }

    if((!retrain) && (learn_parm->epsilon_crit>(*maxdiff))) 
      learn_parm->epsilon_crit=(*maxdiff);
    if((!retrain) && (learn_parm->epsilon_crit>epsilon_crit_org)) {
      learn_parm->epsilon_crit/=2.0;
      retrain=1;
      noshrink=1;
    }
    if(learn_parm->epsilon_crit<epsilon_crit_org) 
      learn_parm->epsilon_crit=epsilon_crit_org;
    
    if(verbosity>=2) {
      printf(" => (%ld SV (incl. %ld SV at u-bound), max violation=%.5f)\n",
	     supvecnum,model->at_upper_bound,(*maxdiff)); 
      fflush(stdout);
    }
    if(verbosity>=3) {
      printf("\n");
    }

    if((!retrain) && (transduction)) {
      for(i=0;(i<totdoc);i++) {
	shrink_state->active[i]=1;
      }
      activenum=compute_index(shrink_state->active,totdoc,active2dnum);
      inactivenum=0;
      if(verbosity==1) printf("done\n");
      retrain=incorporate_unlabeled_examples(model,label,inconsistent,
					     unlabeled,a,lin,totdoc,
					     selcrit,selexam,key,
					     transductcycle,kernel_parm,
					     learn_parm);
      epsilon_crit_org=learn_parm->epsilon_crit;
      if(kernel_parm->kernel_type == LINEAR)
	learn_parm->epsilon_crit=1; 
      transductcycle++;
      /* reset watchdog */
      bestmaxdiff=(*maxdiff);
      bestmaxdiffiter=iteration;
    } 
    else if(((iteration % 10) == 0) && (!noshrink)) {
      activenum=shrink_problem(docs,learn_parm,shrink_state,kernel_parm,
			       active2dnum,last_suboptimal_at,iteration,totdoc,
			       maxl((long)(activenum/10),
				    maxl((long)(totdoc/500),100)),
			       a,inconsistent);
      inactivenum=totdoc-activenum;
      if((kernel_cache)
	 && (supvecnum>kernel_cache->max_elems)
	 && ((kernel_cache->activenum-activenum)>maxl((long)(activenum/10),500))) {
	kernel_cache_shrink(kernel_cache,totdoc,
			    minl((kernel_cache->activenum-activenum),
				 (kernel_cache->activenum-supvecnum)),
			    shrink_state->active); 
      }
    }

    if((!retrain) && learn_parm->remove_inconsistent) {
      if(verbosity>=1) {
	printf(" Moving training errors to inconsistent examples...");
	fflush(stdout);
      }
      if(learn_parm->remove_inconsistent == 1) {
	retrain=identify_inconsistent(a,label,unlabeled,totdoc,learn_parm,
				      &inconsistentnum,inconsistent); 
      }
      else if(learn_parm->remove_inconsistent == 2) {
	retrain=identify_misclassified(lin,label,unlabeled,totdoc,
				       model,&inconsistentnum,inconsistent); 
      }
      else if(learn_parm->remove_inconsistent == 3) {
	retrain=identify_one_misclassified(lin,label,unlabeled,totdoc,
				   model,&inconsistentnum,inconsistent);
      }
      if(retrain) {
	if(kernel_parm->kernel_type == LINEAR) { /* reinit shrinking */
	  learn_parm->epsilon_crit=2.0;
	} 
      }
      if(verbosity>=1) {
	printf("done.\n");
	if(retrain) {
	  printf(" Now %ld inconsistent examples.\n",inconsistentnum);
	}
      }
    }
  } /* end of loop */

  free(chosen);
  free(last_suboptimal_at);
  free(key);
  free(selcrit);
  free(selexam);
  free(a_old);
  free(aicache);
  free(working2dnum);
  free(active2dnum);
  free(qp.opt_ce);
  free(qp.opt_ce0);
  free(qp.opt_g);
  free(qp.opt_g0);
  free(qp.opt_xinit);
  free(qp.opt_low);
  free(qp.opt_up);
  free(weights);

  learn_parm->epsilon_crit=epsilon_crit_org; /* restore org */
  model->maxdiff=(*maxdiff);

  return(iteration);
}

long optimize_to_convergence_sharedslack(DOC **docs, long int *label, 
			     long int totdoc, 
			     long int totwords, LEARN_PARM *learn_parm, 
			     KERNEL_PARM *kernel_parm, 
			     KERNEL_CACHE *kernel_cache, 
			     SHRINK_STATE *shrink_state, MODEL *model, 
			     double *a, double *lin, double *c, 
			     TIMING *timing_profile, double *maxdiff)
     /* docs: Training vectors (x-part) */
     /* label: Training labels/value (y-part, zero if test example for
			      transduction) */
     /* totdoc: Number of examples in docs/label */
     /* totwords: Number of features (i.e. highest feature index) */
     /* learn_parm: Learning paramenters */
     /* kernel_parm: Kernel paramenters */
     /* kernel_cache: Initialized/partly filled Cache, if using a kernel. 
                      NULL if linear. */
     /* shrink_state: State of active variables */
     /* model: Returns learning result */
     /* a: alphas */
     /* lin: linear component of gradient */
     /* c: right hand side of inequalities (margin) */
     /* maxdiff: returns maximum violation of KT-conditions */
{
  long *chosen,*key,i,j,jj,*last_suboptimal_at,noshrink,*unlabeled;
  long *inconsistent,choosenum,already_chosen=0,iteration;
  long misclassified,supvecnum=0,*active2dnum,inactivenum;
  long *working2dnum,*selexam,*ignore;
  long activenum,retrain,maxslackid,slackset,jointstep;
  double criterion,eq_target;
  double *a_old,*alphaslack;
  long t0=0,t1=0,t2=0,t3=0,t4=0,t5=0,t6=0; /* timing */
  double epsilon_crit_org,maxsharedviol; 
  double bestmaxdiff;
  long   bestmaxdiffiter,terminate;

  double *selcrit;  /* buffer for sorting */        
  CFLOAT *aicache;  /* buffer to keep one row of hessian */
  double *weights;  /* buffer for weight vector in linear case */
  QP qp;            /* buffer for one quadratic program */
  double *slack;    /* vector of slack variables for optimization with
		       shared slacks */

  epsilon_crit_org=learn_parm->epsilon_crit; /* save org */
  if(kernel_parm->kernel_type == LINEAR) {
    learn_parm->epsilon_crit=2.0;
    kernel_cache=NULL;   /* caching makes no sense for linear kernel */
  } 
  learn_parm->epsilon_shrink=2;
  (*maxdiff)=1;

  learn_parm->totwords=totwords;

  chosen = (long *)my_malloc(sizeof(long)*totdoc);
  unlabeled = (long *)my_malloc(sizeof(long)*totdoc);
  inconsistent = (long *)my_malloc(sizeof(long)*totdoc);
  ignore = (long *)my_malloc(sizeof(long)*totdoc);
  key = (long *)my_malloc(sizeof(long)*(totdoc+11)); 
  selcrit = (double *)my_malloc(sizeof(double)*totdoc);
  selexam = (long *)my_malloc(sizeof(long)*totdoc);
  a_old = (double *)my_malloc(sizeof(double)*totdoc);
  aicache = (CFLOAT *)my_malloc(sizeof(CFLOAT)*totdoc);
  working2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  active2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
  qp.opt_ce = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_ce0 = (double *)my_malloc(sizeof(double));
  qp.opt_g = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize
				 *learn_parm->svm_maxqpsize);
  qp.opt_g0 = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_xinit = (double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_low=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  qp.opt_up=(double *)my_malloc(sizeof(double)*learn_parm->svm_maxqpsize);
  weights=(double *)my_malloc(sizeof(double)*(totwords+1));
  maxslackid=0;
  for(i=0;i<totdoc;i++) {    /* determine size of slack array */
    if(maxslackid<docs[i]->slackid)
      maxslackid=docs[i]->slackid;
  }
  slack=(double *)my_malloc(sizeof(double)*(maxslackid+1));
  alphaslack=(double *)my_malloc(sizeof(double)*(maxslackid+1));
  last_suboptimal_at = (long *)my_malloc(sizeof(long)*(maxslackid+1));
  for(i=0;i<=maxslackid;i++) {    /* init shared slacks */
    slack[i]=0;
    alphaslack[i]=0;
    last_suboptimal_at[i]=1;
  }

  choosenum=0;
  retrain=1;
  iteration=1;
  bestmaxdiffiter=1;
  bestmaxdiff=999999999;
  terminate=0;

  if(kernel_cache) {
    kernel_cache->time=iteration;  /* for lru cache */
    kernel_cache_reset_lru(kernel_cache);
  }

  for(i=0;i<totdoc;i++) {    /* various inits */
    chosen[i]=0;
    unlabeled[i]=0;
    inconsistent[i]=0;
    ignore[i]=0;
    a_old[i]=a[i];
  }
  activenum=compute_index(shrink_state->active,totdoc,active2dnum);
  inactivenum=totdoc-activenum;
  clear_index(working2dnum);

  /* call to init slack and alphaslack */
  compute_shared_slacks(docs,label,a,lin,c,active2dnum,learn_parm,
			slack,alphaslack);

                            /* repeat this loop until we have convergence */
  for(;retrain && (!terminate);iteration++) {

    if(kernel_cache)
      kernel_cache->time=iteration;  /* for lru cache */
    if(verbosity>=2) {
      printf(
	"Iteration %ld: ",iteration); fflush(stdout);
    }
    else if(verbosity==1) {
      printf("."); fflush(stdout);
    }

    if(verbosity>=2) t0=get_runtime();
    if(verbosity>=3) {
      printf("\nSelecting working set... "); fflush(stdout); 
    }

    if(learn_parm->svm_newvarsinqp>learn_parm->svm_maxqpsize) 
      learn_parm->svm_newvarsinqp=learn_parm->svm_maxqpsize;

    /* select working set according to steepest gradient */
    jointstep=0;
    eq_target=0;
    if(iteration % 101) {
      slackset=select_next_qp_slackset(docs,label,a,lin,slack,alphaslack,c,
				       learn_parm,active2dnum,&maxsharedviol);
      if((iteration % 2) 
	 || (!slackset) || (maxsharedviol<learn_parm->epsilon_crit)){
	/* do a step with examples from different slack sets */
	if(verbosity >= 2) {
	  printf("(i-step)"); fflush(stdout);
	}
	i=0;
	for(jj=0;(j=working2dnum[jj])>=0;jj++) { /* clear old part of working set */
	  if((chosen[j]>=(learn_parm->svm_maxqpsize/
			  minl(learn_parm->svm_maxqpsize,
			       learn_parm->svm_newvarsinqp)))) {
	    chosen[j]=0; 
	    choosenum--; 
	  }
	  else {
	    chosen[j]++;
	    working2dnum[i++]=j;
	  }
	}
	working2dnum[i]=-1;
	
	already_chosen=0;
	if((minl(learn_parm->svm_newvarsinqp,
		 learn_parm->svm_maxqpsize-choosenum)>=4) 
	   && (kernel_parm->kernel_type != LINEAR)) {
	  /* select part of the working set from cache */
	  already_chosen=select_next_qp_subproblem_grad(
			      label,unlabeled,a,lin,c,totdoc,
			      (long)(minl(learn_parm->svm_maxqpsize-choosenum,
					  learn_parm->svm_newvarsinqp)
				     /2),
			      learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,
			      (long)1,key,chosen);
	  choosenum+=already_chosen;
	}
	choosenum+=select_next_qp_subproblem_grad(
                              label,unlabeled,a,lin,c,totdoc,
                              minl(learn_parm->svm_maxqpsize-choosenum,
				   learn_parm->svm_newvarsinqp-already_chosen),
                              learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,
			      (long)0,key,chosen);
      }
      else { /* do a step with all examples from same slack set */
	if(verbosity >= 2) {
	  printf("(j-step on %ld)",slackset); fflush(stdout);
	}
	jointstep=1;
	for(jj=0;(j=working2dnum[jj])>=0;jj++) { /* clear working set */
	    chosen[j]=0; 
	}
	working2dnum[0]=-1;
	eq_target=alphaslack[slackset];
	for(j=0;j<totdoc;j++) {                  /* mask all but slackset */
	  /* for(jj=0;(j=active2dnum[jj])>=0;jj++) { */
	  if(docs[j]->slackid != slackset)
	    ignore[j]=1; 
	  else {
	    ignore[j]=0; 
	    learn_parm->svm_cost[j]=learn_parm->svm_c;
	    /* printf("Inslackset(%ld,%ld)",j,shrink_state->active[j]); */
	  }
	}
	learn_parm->biased_hyperplane=1;
	choosenum=select_next_qp_subproblem_grad(
                              label,unlabeled,a,lin,c,totdoc,
                              learn_parm->svm_maxqpsize,
                              learn_parm,ignore,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,
			      (long)0,key,chosen);
	learn_parm->biased_hyperplane=0;
      }
    }
    else { /* once in a while, select a somewhat random working set
	      to get unlocked of infinite loops due to numerical
	      inaccuracies in the core qp-solver */
      choosenum+=select_next_qp_subproblem_rand(
                              label,unlabeled,a,lin,c,totdoc,
                              minl(learn_parm->svm_maxqpsize-choosenum,
				   learn_parm->svm_newvarsinqp),
                              learn_parm,inconsistent,active2dnum,
			      working2dnum,selcrit,selexam,kernel_cache,key,
			      chosen,iteration);
    }

    if(verbosity>=2) {
      printf(" %ld vectors chosen\n",choosenum); fflush(stdout); 
    }

    if(verbosity>=2) t1=get_runtime();

    if(kernel_cache) 
      cache_multiple_kernel_rows(kernel_cache,docs,working2dnum,
				 choosenum,kernel_parm); 

    if(verbosity>=2) t2=get_runtime();
    if(jointstep) learn_parm->biased_hyperplane=1;
    optimize_svm(docs,label,unlabeled,ignore,eq_target,chosen,active2dnum,
		 model,totdoc,working2dnum,choosenum,a,lin,c,learn_parm,
		 aicache,kernel_parm,&qp,&epsilon_crit_org);
    learn_parm->biased_hyperplane=0;

    for(jj=0;(i=working2dnum[jj])>=0;jj++)   /* recompute sums of alphas */
      alphaslack[docs[i]->slackid]+=(a[i]-a_old[i]);
    for(jj=0;(i=working2dnum[jj])>=0;jj++) { /* reduce alpha to fulfill
						constraints */
      if(alphaslack[docs[i]->slackid] > learn_parm->svm_c) {
	if(a[i] < (alphaslack[docs[i]->slackid]-learn_parm->svm_c)) {
	  alphaslack[docs[i]->slackid]-=a[i];
	  a[i]=0;
	}
	else {
	  a[i]-=(alphaslack[docs[i]->slackid]-learn_parm->svm_c);
	  alphaslack[docs[i]->slackid]=learn_parm->svm_c;
	}
      }
    }
    for(jj=0;(i=active2dnum[jj])>=0;jj++) 
      learn_parm->svm_cost[i]=a[i]+(learn_parm->svm_c
				    -alphaslack[docs[i]->slackid]);

    if(verbosity>=2) t3=get_runtime();
    update_linear_component(docs,label,active2dnum,a,a_old,working2dnum,totdoc,
			    totwords,kernel_parm,kernel_cache,lin,aicache,
			    weights);
    compute_shared_slacks(docs,label,a,lin,c,active2dnum,learn_parm,
			  slack,alphaslack);

    if(verbosity>=2) t4=get_runtime();
    supvecnum=calculate_svm_model(docs,label,unlabeled,lin,a,a_old,c,
		                  learn_parm,working2dnum,active2dnum,model);

    if(verbosity>=2) t5=get_runtime();

    /* The following computation of the objective function works only */
    /* relative to the active variables */
    if(verbosity>=3) {
      criterion=compute_objective_function(a,lin,c,learn_parm->eps,label,
		                           active2dnum);
      printf("Objective function (over active variables): %.16f\n",criterion);
      fflush(stdout); 
    }

    for(jj=0;(i=working2dnum[jj])>=0;jj++) {
      a_old[i]=a[i];
    }

    retrain=check_optimality_sharedslack(docs,model,label,a,lin,c,
                             slack,alphaslack,totdoc,learn_parm,
			     maxdiff,epsilon_crit_org,&misclassified,
			     active2dnum,last_suboptimal_at,
			     iteration,kernel_parm);

    if(verbosity>=2) {
      t6=get_runtime();
      timing_profile->time_select+=t1-t0;
      timing_profile->time_kernel+=t2-t1;
      timing_profile->time_opti+=t3-t2;
      timing_profile->time_update+=t4-t3;
      timing_profile->time_model+=t5-t4;
      timing_profile->time_check+=t6-t5;
    }

    /* checking whether optimizer got stuck */
    if((*maxdiff) < bestmaxdiff) {
      bestmaxdiff=(*maxdiff);
      bestmaxdiffiter=iteration;
    }
    if(iteration > (bestmaxdiffiter+learn_parm->maxiter)) { 
      /* long time no progress? */
      terminate=1;
      retrain=0;
      if(verbosity>=1) 
	printf("\nWARNING: Relaxing KT-Conditions due to slow progress! Terminating!\n");
    }

    noshrink=0; 

    if((!retrain) && (inactivenum>0) 
       && ((!learn_parm->skip_final_opt_check) 
	   || (kernel_parm->kernel_type == LINEAR))) { 
      if(((verbosity>=1) && (kernel_parm->kernel_type != LINEAR)) 
	 || (verbosity>=2)) {
	if(verbosity==1) {
	  printf("\n");
	}
	printf(" Checking optimality of inactive variables..."); 
	fflush(stdout);
      }
      t1=get_runtime();
      reactivate_inactive_examples(label,unlabeled,a,shrink_state,lin,c,totdoc,
				   totwords,iteration,learn_parm,inconsistent,
				   docs,kernel_parm,kernel_cache,model,aicache,
				   weights,maxdiff);
      /* Update to new active variables. */
      activenum=compute_index(shrink_state->active,totdoc,active2dnum);
      inactivenum=totdoc-activenum;
      /* check optimality, since check in reactivate does not work for
	 sharedslacks */
      retrain=check_optimality_sharedslack(docs,model,label,a,lin,c,
			     slack,alphaslack,totdoc,learn_parm,
			     maxdiff,epsilon_crit_org,&misclassified,
			     active2dnum,last_suboptimal_at,
			     iteration,kernel_parm);

      /* reset watchdog */
      bestmaxdiff=(*maxdiff);
      bestmaxdiffiter=iteration;
      /* termination criterion */
      noshrink=1;
      retrain=0;
      if((*maxdiff) > learn_parm->epsilon_crit) 
	retrain=1;
      timing_profile->time_shrink+=get_runtime()-t1;
      if(((verbosity>=1) && (kernel_parm->kernel_type != LINEAR)) 
	 || (verbosity>=2)) {
	printf("done.\n");  fflush(stdout);
        printf(" Number of inactive variables = %ld\n",inactivenum);
      }		  
    }

    if((!retrain) && (learn_parm->epsilon_crit>(*maxdiff))) 
      learn_parm->epsilon_crit=(*maxdiff);
    if((!retrain) && (learn_parm->epsilon_crit>epsilon_crit_org)) {
      learn_parm->epsilon_crit/=2.0;
      retrain=1;
      noshrink=1;
    }
    if(learn_parm->epsilon_crit<epsilon_crit_org) 
      learn_parm->epsilon_crit=epsilon_crit_org;
    
    if(verbosity>=2) {
      printf(" => (%ld SV (incl. %ld SV at u-bound), max violation=%.5f)\n",
	     supvecnum,model->at_upper_bound,(*maxdiff)); 
      fflush(stdout);
    }
    if(verbosity>=3) {
      printf("\n");
    }

    if(((iteration % 10) == 0) && (!noshrink)) {
      activenum=shrink_problem(docs,learn_parm,shrink_state,
			       kernel_parm,active2dnum,
			       last_suboptimal_at,iteration,totdoc,
			       maxl((long)(activenum/10),
				    maxl((long)(totdoc/500),100)),
			       a,inconsistent);
      inactivenum=totdoc-activenum;
      if((kernel_cache)
	 && (supvecnum>kernel_cache->max_elems)
	 && ((kernel_cache->activenum-activenum)>maxl((long)(activenum/10),500))) {
	kernel_cache_shrink(kernel_cache,totdoc,
			    minl((kernel_cache->activenum-activenum),
				 (kernel_cache->activenum-supvecnum)),
			    shrink_state->active); 
      }
    }

  } /* end of loop */


  free(alphaslack);
  free(slack);
  free(chosen);
  free(unlabeled);
  free(inconsistent);
  free(ignore);
  free(last_suboptimal_at);
  free(key);
  free(selcrit);
  free(selexam);
  free(a_old);
  free(aicache);
  free(working2dnum);
  free(active2dnum);
  free(qp.opt_ce);
  free(qp.opt_ce0);
  free(qp.opt_g);
  free(qp.opt_g0);
  free(qp.opt_xinit);
  free(qp.opt_low);
  free(qp.opt_up);
  free(weights);

  learn_parm->epsilon_crit=epsilon_crit_org; /* restore org */
  model->maxdiff=(*maxdiff);

  return(iteration);
}


double compute_objective_function(double *a, double *lin, double *c, 
				  double eps, long int *label, 
				  long int *active2dnum)
     /* Return value of objective function. */
     /* Works only relative to the active variables! */
{
  long i,ii;
  double criterion;
  /* calculate value of objective function */
  criterion=0;
  for(ii=0;active2dnum[ii]>=0;ii++) {
    i=active2dnum[ii];
    criterion=criterion+(eps-(double)label[i]*c[i])*a[i]+0.5*a[i]*label[i]*lin[i];
  } 
  return(criterion);
}

void clear_index(long int *index)
              /* initializes and empties index */
{
  index[0]=-1;
} 

void add_to_index(long int *index, long int elem)
     /* initializes and empties index */
{
  register long i;
  for(i=0;index[i] != -1;i++);
  index[i]=elem;
  index[i+1]=-1;
}

long compute_index(long int *binfeature, long int range, long int *index)
     /* create an inverted index of binfeature */
{               
  register long i,ii;

  ii=0;
  for(i=0;i<range;i++) {
    if(binfeature[i]) {
      index[ii]=i;
      ii++;
    }
  }
  for(i=0;i<4;i++) {
    index[ii+i]=-1;
  }
  return(ii);
}


void optimize_svm(DOC **docs, long int *label, long int *unlabeled, 
		  long int *exclude_from_eq_const, double eq_target,
		  long int *chosen, long int *active2dnum, MODEL *model, 
		  long int totdoc, long int *working2dnum, long int varnum, 
		  double *a, double *lin, double *c, LEARN_PARM *learn_parm, 
		  CFLOAT *aicache, KERNEL_PARM *kernel_parm, QP *qp, 
		  double *epsilon_crit_target)
     /* Do optimization on the working set. */
{
    long i;
    double *a_v;

    compute_matrices_for_optimization(docs,label,unlabeled,
				      exclude_from_eq_const,eq_target,chosen,
				      active2dnum,working2dnum,model,a,lin,c,
				      varnum,totdoc,learn_parm,aicache,
				      kernel_parm,qp);

    if(verbosity>=3) {
      printf("Running optimizer..."); fflush(stdout);
    }
    /* call the qp-subsolver */
    a_v=optimize_qp(qp,epsilon_crit_target,
		    learn_parm->svm_maxqpsize,
		    &(model->b),   /* in case the optimizer gives us */
                                   /* the threshold for free. otherwise */
                                   /* b is calculated in calculate_model. */
		    learn_parm);
    if(verbosity>=3) {         
      printf("done\n");
    }

    for(i=0;i<varnum;i++) {
      a[working2dnum[i]]=a_v[i];
      /*
      if(a_v[i]<=(0+learn_parm->epsilon_a)) {
	a[working2dnum[i]]=0;
      }
      else if(a_v[i]>=(learn_parm->svm_cost[working2dnum[i]]-learn_parm->epsilon_a)) {
	a[working2dnum[i]]=learn_parm->svm_cost[working2dnum[i]];
      }
      */
    }
}

void compute_matrices_for_optimization(DOC **docs, long int *label, 
          long int *unlabeled, long *exclude_from_eq_const, double eq_target,
	  long int *chosen, long int *active2dnum, 
          long int *key, MODEL *model, double *a, double *lin, double *c, 
	  long int varnum, long int totdoc, LEARN_PARM *learn_parm, 
          CFLOAT *aicache, KERNEL_PARM *kernel_parm, QP *qp)
{
  register long ki,kj,i,j;
  register double kernel_temp;

  if(verbosity>=3) {
    fprintf(stdout,"Computing qp-matrices (type %ld kernel [degree %ld, rbf_gamma %f, coef_lin %f, coef_const %f])...",kernel_parm->kernel_type,kernel_parm->poly_degree,kernel_parm->rbf_gamma,kernel_parm->coef_lin,kernel_parm->coef_const); 
    fflush(stdout);
  }

  qp->opt_n=varnum;
  qp->opt_ce0[0]=-eq_target; /* compute the constant for equality constraint */
  for(j=1;j<model->sv_num;j++) { /* start at 1 */
    if((!chosen[(model->supvec[j])->docnum])
       && (!exclude_from_eq_const[(model->supvec[j])->docnum])) {
      qp->opt_ce0[0]+=model->alpha[j];
    }
  } 
  if(learn_parm->biased_hyperplane) 
    qp->opt_m=1;
  else 
    qp->opt_m=0;  /* eq-constraint will be ignored */

  /* init linear part of objective function */
  for(i=0;i<varnum;i++) {
    qp->opt_g0[i]=lin[key[i]];
  }

  for(i=0;i<varnum;i++) {
    ki=key[i];

    /* Compute the matrix for equality constraints */
    qp->opt_ce[i]=label[ki];
    qp->opt_low[i]=0;
    qp->opt_up[i]=learn_parm->svm_cost[ki];

    kernel_temp=(double)kernel(kernel_parm,docs[ki],docs[ki]); 
    /* compute linear part of objective function */
    qp->opt_g0[i]-=(kernel_temp*a[ki]*(double)label[ki]); 
    /* compute quadratic part of objective function */
    qp->opt_g[varnum*i+i]=kernel_temp;
    for(j=i+1;j<varnum;j++) {
      kj=key[j];
      kernel_temp=(double)kernel(kernel_parm,docs[ki],docs[kj]);
      /* compute linear part of objective function */
      qp->opt_g0[i]-=(kernel_temp*a[kj]*(double)label[kj]);
      qp->opt_g0[j]-=(kernel_temp*a[ki]*(double)label[ki]); 
      /* compute quadratic part of objective function */
      qp->opt_g[varnum*i+j]=(double)label[ki]*(double)label[kj]*kernel_temp;
      qp->opt_g[varnum*j+i]=(double)label[ki]*(double)label[kj]*kernel_temp;
    }

    if(verbosity>=3) {
      if(i % 20 == 0) {
	fprintf(stdout,"%ld..",i); fflush(stdout);
      }
    }
  }

  for(i=0;i<varnum;i++) {
    /* assure starting at feasible point */
    qp->opt_xinit[i]=a[key[i]];
    /* set linear part of objective function */
    qp->opt_g0[i]=(learn_parm->eps-(double)label[key[i]]*c[key[i]])+qp->opt_g0[i]*(double)label[key[i]];    
  }

  if(verbosity>=3) {
    fprintf(stdout,"done\n");
  }
}

long calculate_svm_model(DOC **docs, long int *label, long int *unlabeled, 
			 double *lin, double *a, double *a_old, double *c, 
			 LEARN_PARM *learn_parm, long int *working2dnum, 
			 long int *active2dnum, MODEL *model)
     /* Compute decision function based on current values */
     /* of alpha. */
{
  long i,ii,pos,b_calculated=0,first_low,first_high;
  double ex_c,b_temp,b_low,b_high;

  if(verbosity>=3) {
    printf("Calculating model..."); fflush(stdout);
  }

  if(!learn_parm->biased_hyperplane) {
    model->b=0;
    b_calculated=1;
  }

  for(ii=0;(i=working2dnum[ii])>=0;ii++) {
    if((a_old[i]>0) && (a[i]==0)) { /* remove from model */
      pos=model->index[i]; 
      model->index[i]=-1;
      (model->sv_num)--;
      model->supvec[pos]=model->supvec[model->sv_num];
      model->alpha[pos]=model->alpha[model->sv_num];
      model->index[(model->supvec[pos])->docnum]=pos;
    }
    else if((a_old[i]==0) && (a[i]>0)) { /* add to model */
      model->supvec[model->sv_num]=docs[i];
      model->alpha[model->sv_num]=a[i]*(double)label[i];
      model->index[i]=model->sv_num;
      (model->sv_num)++;
    }
    else if(a_old[i]==a[i]) { /* nothing to do */
    }
    else {  /* just update alpha */
      model->alpha[model->index[i]]=a[i]*(double)label[i];
    }
      
    ex_c=learn_parm->svm_cost[i]-learn_parm->epsilon_a;
    if((a_old[i]>=ex_c) && (a[i]<ex_c)) { 
      (model->at_upper_bound)--;
    }
    else if((a_old[i]<ex_c) && (a[i]>=ex_c)) { 
      (model->at_upper_bound)++;
    }

    if((!b_calculated) 
       && (a[i]>learn_parm->epsilon_a) && (a[i]<ex_c)) {   /* calculate b */
     	model->b=((double)label[i]*learn_parm->eps-c[i]+lin[i]); 
	/* model->b=(-(double)label[i]+lin[i]); */
	b_calculated=1;
    }
  }      

  /* No alpha in the working set not at bounds, so b was not
     calculated in the usual way. The following handles this special
     case. */
  if(learn_parm->biased_hyperplane 
     && (!b_calculated)
     && (model->sv_num-1 == model->at_upper_bound)) { 
    first_low=1;
    first_high=1;
    b_low=0;
    b_high=0;
    for(ii=0;(i=active2dnum[ii])>=0;ii++) {
      ex_c=learn_parm->svm_cost[i]-learn_parm->epsilon_a;
      if(a[i]<ex_c) { 
	if(label[i]>0)  {
	  b_temp=-(learn_parm->eps-c[i]+lin[i]);
	  if((b_temp>b_low) || (first_low)) {
	    b_low=b_temp;
	    first_low=0;
	  }
	}
	else {
	  b_temp=-(-learn_parm->eps-c[i]+lin[i]);
	  if((b_temp<b_high) || (first_high)) {
	    b_high=b_temp;
	    first_high=0;
	  }
	}
      }
      else {
	if(label[i]<0)  {
	  b_temp=-(-learn_parm->eps-c[i]+lin[i]);
	  if((b_temp>b_low) || (first_low)) {
	    b_low=b_temp;
	    first_low=0;
	  }
	}
	else {
	  b_temp=-(learn_parm->eps-c[i]+lin[i]);
	  if((b_temp<b_high) || (first_high)) {
	    b_high=b_temp;
	    first_high=0;
	  }
	}
      }
    }
    if(first_high) {
      model->b=-b_low;
    }
    else if(first_low) {
      model->b=-b_high;
    }
    else {
      model->b=-(b_high+b_low)/2.0;  /* select b as the middle of range */
      /* printf("\nb_low=%f, b_high=%f,b=%f\n",b_low,b_high,model->b); */
    }
  }

  if(verbosity>=3) {
    printf("done\n"); fflush(stdout);
  }

  return(model->sv_num-1); /* have to substract one, since element 0 is empty*/
}

long check_optimality(MODEL *model, long int *label, long int *unlabeled, 
		      double *a, double *lin, double *c, long int totdoc, 
		      LEARN_PARM *learn_parm, double *maxdiff, 
		      double epsilon_crit_org, long int *misclassified, 
		      long int *inconsistent, long int *active2dnum,
		      long int *last_suboptimal_at, 
		      long int iteration, KERNEL_PARM *kernel_parm)
     /* Check KT-conditions */
{
  long i,ii,retrain;
  double dist,ex_c,target;

  if(kernel_parm->kernel_type == LINEAR) {  /* be optimistic */
    learn_parm->epsilon_shrink=-learn_parm->epsilon_crit+epsilon_crit_org;  
  }
  else {  /* be conservative */
    learn_parm->epsilon_shrink=learn_parm->epsilon_shrink*0.7+(*maxdiff)*0.3; 
  }
  retrain=0;
  (*maxdiff)=0;
  (*misclassified)=0;
  for(ii=0;(i=active2dnum[ii])>=0;ii++) {
    if((!inconsistent[i]) && label[i]) {
      dist=(lin[i]-model->b)*(double)label[i];/* 'distance' from
						 hyperplane*/
      target=-(learn_parm->eps-(double)label[i]*c[i]);
      ex_c=learn_parm->svm_cost[i]-learn_parm->epsilon_a;
      if(dist <= 0) {       
	(*misclassified)++;  /* does not work due to deactivation of var */
      }
      if((a[i]>learn_parm->epsilon_a) && (dist > target)) {
	if((dist-target)>(*maxdiff))  /* largest violation */
	  (*maxdiff)=dist-target;
      }
      else if((a[i]<ex_c) && (dist < target)) {
	if((target-dist)>(*maxdiff))  /* largest violation */
	  (*maxdiff)=target-dist;
      }
      /* Count how long a variable was at lower/upper bound (and optimal).*/
      /* Variables, which were at the bound and optimal for a long */
      /* time are unlikely to become support vectors. In case our */
      /* cache is filled up, those variables are excluded to save */
      /* kernel evaluations. (See chapter 'Shrinking').*/ 
      if((a[i]>(learn_parm->epsilon_a)) 
	 && (a[i]<ex_c)) { 
	last_suboptimal_at[i]=iteration;         /* not at bound */
      }
      else if((a[i]<=(learn_parm->epsilon_a)) 
	      && (dist < (target+learn_parm->epsilon_shrink))) {
	last_suboptimal_at[i]=iteration;         /* not likely optimal */
      }
      else if((a[i]>=ex_c)
	      && (dist > (target-learn_parm->epsilon_shrink)))  { 
	last_suboptimal_at[i]=iteration;         /* not likely optimal */
      }
    }   
  }
  /* termination criterion */
  if((!retrain) && ((*maxdiff) > learn_parm->epsilon_crit)) {  
    retrain=1;
  }
  return(retrain);
}

long check_optimality_sharedslack(DOC **docs, MODEL *model, long int *label,
		      double *a, double *lin, double *c, double *slack,
		      double *alphaslack,
		      long int totdoc, 
		      LEARN_PARM *learn_parm, double *maxdiff, 
		      double epsilon_crit_org, long int *misclassified, 
		      long int *active2dnum,
		      long int *last_suboptimal_at, 
		      long int iteration, KERNEL_PARM *kernel_parm)
     /* Check KT-conditions */
{
  long i,ii,retrain;
  double dist,ex_c=0,target;

  if(kernel_parm->kernel_type == LINEAR) {  /* be optimistic */
    learn_parm->epsilon_shrink=-learn_parm->epsilon_crit+epsilon_crit_org;  
  }
  else {  /* be conservative */
    learn_parm->epsilon_shrink=learn_parm->epsilon_shrink*0.7+(*maxdiff)*0.3; 
  }

  retrain=0;
  (*maxdiff)=0;
  (*misclassified)=0;
  for(ii=0;(i=active2dnum[ii])>=0;ii++) {
    /* 'distance' from hyperplane*/
    dist=(lin[i]-model->b)*(double)label[i]+slack[docs[i]->slackid];
    target=-(learn_parm->eps-(double)label[i]*c[i]);
    ex_c=learn_parm->svm_c-learn_parm->epsilon_a;
    if((a[i]>learn_parm->epsilon_a) && (dist > target)) {
      if((dist-target)>(*maxdiff)) {  /* largest violation */
	(*maxdiff)=dist-target;
	if(verbosity>=5) printf("sid %ld: dist=%.2f, target=%.2f, slack=%.2f, a=%f, alphaslack=%f\n",docs[i]->slackid,dist,target,slack[docs[i]->slackid],a[i],alphaslack[docs[i]->slackid]);
	if(verbosity>=5) printf(" (single %f)\n",(*maxdiff));
      }
    }
    if((alphaslack[docs[i]->slackid]<ex_c) && (slack[docs[i]->slackid]>0)) {
      if((slack[docs[i]->slackid])>(*maxdiff)) { /* largest violation */
	(*maxdiff)=slack[docs[i]->slackid];
	if(verbosity>=5) printf("sid %ld: dist=%.2f, target=%.2f, slack=%.2f, a=%f, alphaslack=%f\n",docs[i]->slackid,dist,target,slack[docs[i]->slackid],a[i],alphaslack[docs[i]->slackid]);
	if(verbosity>=5) printf(" (joint %f)\n",(*maxdiff));
      }
    }
    /* Count how long a variable was at lower/upper bound (and optimal).*/
    /* Variables, which were at the bound and optimal for a long */
    /* time are unlikely to become support vectors. In case our */
    /* cache is filled up, those variables are excluded to save */
    /* kernel evaluations. (See chapter 'Shrinking').*/ 
    if((a[i]>(learn_parm->epsilon_a)) 
       && (a[i]<ex_c)) { 
      last_suboptimal_at[docs[i]->slackid]=iteration;  /* not at bound */
    }
    else if((a[i]<=(learn_parm->epsilon_a)) 
	    && (dist < (target+learn_parm->epsilon_shrink))) {
      last_suboptimal_at[docs[i]->slackid]=iteration;  /* not likely optimal */
    }
    else if((a[i]>=ex_c)
	    && (slack[docs[i]->slackid] < learn_parm->epsilon_shrink))  { 
      last_suboptimal_at[docs[i]->slackid]=iteration;  /* not likely optimal */
    }
  }   
  /* termination criterion */
  if((!retrain) && ((*maxdiff) > learn_parm->epsilon_crit)) {  
    retrain=1;
  }
  return(retrain);
}

void compute_shared_slacks(DOC **docs, long int *label, 
			   double *a, double *lin, 
			   double *c, long int *active2dnum,
			   LEARN_PARM *learn_parm, 
			   double *slack, double *alphaslack)
     /* compute the value of shared slacks and the joint alphas */
{
  long jj,i;
  double dist,target;

  for(jj=0;(i=active2dnum[jj])>=0;jj++) { /* clear slack variables */
    slack[docs[i]->slackid]=0.0;
    alphaslack[docs[i]->slackid]=0.0;
  }
  for(jj=0;(i=active2dnum[jj])>=0;jj++) { /* recompute slack variables */
    dist=(lin[i])*(double)label[i];
    target=-(learn_parm->eps-(double)label[i]*c[i]);
    if((target-dist) > slack[docs[i]->slackid])
      slack[docs[i]->slackid]=target-dist;
    alphaslack[docs[i]->slackid]+=a[i];
  }
}


long identify_inconsistent(double *a, long int *label, 
			   long int *unlabeled, long int totdoc, 
			   LEARN_PARM *learn_parm, 
			   long int *inconsistentnum, long int *inconsistent)
{
  long i,retrain;

  /* Throw out examples with multipliers at upper bound. This */
  /* corresponds to the -i 1 option. */
  /* ATTENTION: this is just a heuristic for finding a close */
  /*            to minimum number of examples to exclude to */
  /*            make the problem separable with desired margin */
  retrain=0;
  for(i=0;i<totdoc;i++) {
    if((!inconsistent[i]) && (!unlabeled[i]) 
       && (a[i]>=(learn_parm->svm_cost[i]-learn_parm->epsilon_a))) { 
	(*inconsistentnum)++;
	inconsistent[i]=1;  /* never choose again */
	retrain=2;          /* start over */
	if(verbosity>=3) {
	  printf("inconsistent(%ld)..",i); fflush(stdout);
	}
    }
  }
  return(retrain);
}

long identify_misclassified(double *lin, long int *label, 
			    long int *unlabeled, long int totdoc, 
			    MODEL *model, long int *inconsistentnum, 
			    long int *inconsistent)
{
  long i,retrain;
  double dist;

  /* Throw out misclassified examples. This */
  /* corresponds to the -i 2 option. */
  /* ATTENTION: this is just a heuristic for finding a close */
  /*            to minimum number of examples to exclude to */
  /*            make the problem separable with desired margin */
  retrain=0;
  for(i=0;i<totdoc;i++) {
    dist=(lin[i]-model->b)*(double)label[i]; /* 'distance' from hyperplane*/  
    if((!inconsistent[i]) && (!unlabeled[i]) && (dist <= 0)) { 
	(*inconsistentnum)++;
	inconsistent[i]=1;  /* never choose again */
	retrain=2;          /* start over */
	if(verbosity>=3) {
	  printf("inconsistent(%ld)..",i); fflush(stdout);
	}
    }
  }
  return(retrain);
}

long identify_one_misclassified(double *lin, long int *label, 
				long int *unlabeled, 
				long int totdoc, MODEL *model, 
				long int *inconsistentnum, 
				long int *inconsistent)
{
  long i,retrain,maxex=-1;
  double dist,maxdist=0;

  /* Throw out the 'most misclassified' example. This */
  /* corresponds to the -i 3 option. */
  /* ATTENTION: this is just a heuristic for finding a close */
  /*            to minimum number of examples to exclude to */
  /*            make the problem separable with desired margin */
  retrain=0;
  for(i=0;i<totdoc;i++) {
    if((!inconsistent[i]) && (!unlabeled[i])) {
      dist=(lin[i]-model->b)*(double)label[i];/* 'distance' from hyperplane*/  
      if(dist<maxdist) {
	maxdist=dist;
	maxex=i;
      }
    }
  }
  if(maxex>=0) {
    (*inconsistentnum)++;
    inconsistent[maxex]=1;  /* never choose again */
    retrain=2;          /* start over */
    if(verbosity>=3) {
      printf("inconsistent(%ld)..",i); fflush(stdout);
    }
  }
  return(retrain);
}

void update_linear_component(DOC **docs, long int *label, 
			     long int *active2dnum, double *a, 
			     double *a_old, long int *working2dnum, 
			     long int totdoc, long int totwords, 
			     KERNEL_PARM *kernel_parm, 
			     KERNEL_CACHE *kernel_cache, 
			     double *lin, CFLOAT *aicache, double *weights)
     /* keep track of the linear component */
     /* lin of the gradient etc. by updating */
     /* based on the change of the variables */
     /* in the current working set */
{
  register long i,ii,j,jj;
  register double tec;
  SVECTOR *f;

  if(kernel_parm->kernel_type==0) { /* special linear case */
    clear_vector_n(weights,totwords);
    for(ii=0;(i=working2dnum[ii])>=0;ii++) {
      if(a[i] != a_old[i]) {
	for(f=docs[i]->fvec;f;f=f->next)  
	  add_vector_ns(weights,f,
			f->factor*((a[i]-a_old[i])*(double)label[i]));
      }
    }
    for(jj=0;(j=active2dnum[jj])>=0;jj++) {
      for(f=docs[j]->fvec;f;f=f->next)  
	lin[j]+=f->factor*sprod_ns(weights,f);
    }
  }
  else {                            /* general case */
    for(jj=0;(i=working2dnum[jj])>=0;jj++) {
      if(a[i] != a_old[i]) {
	get_kernel_row(kernel_cache,docs,i,totdoc,active2dnum,aicache,
		       kernel_parm);
	for(ii=0;(j=active2dnum[ii])>=0;ii++) {
	  tec=aicache[j];
	  lin[j]+=(((a[i]*tec)-(a_old[i]*tec))*(double)label[i]);
	}
      }
    }
  }
}


long incorporate_unlabeled_examples(MODEL *model, long int *label, 
				    long int *inconsistent, 
				    long int *unlabeled, 
				    double *a, double *lin, 
				    long int totdoc, double *selcrit, 
				    long int *select, long int *key, 
				    long int transductcycle, 
				    KERNEL_PARM *kernel_parm, 
				    LEARN_PARM *learn_parm)
{
  long i,j,k,j1,j2,j3,j4,unsupaddnum1=0,unsupaddnum2=0;
  long pos,neg,upos,uneg,orgpos,orgneg,nolabel,newpos,newneg,allunlab;
  double dist,model_length,posratio,negratio;
  long check_every=2;
  double loss;
  static double switchsens=0.0,switchsensorg=0.0;
  double umin,umax,sumalpha;
  long imin=0,imax=0;
  static long switchnum=0;

  switchsens/=1.2;

  /* assumes that lin[] is up to date -> no inactive vars */

  orgpos=0;
  orgneg=0;
  newpos=0;
  newneg=0;
  nolabel=0;
  allunlab=0;
  for(i=0;i<totdoc;i++) {
    if(!unlabeled[i]) {
      if(label[i] > 0) {
	orgpos++;
      }
      else {
	orgneg++;
      }
    }
    else {
      allunlab++;
      if(unlabeled[i]) {
	if(label[i] > 0) {
	  newpos++;
	}
	else if(label[i] < 0) {
	  newneg++;
	}
      }
    }
    if(label[i]==0) {
      nolabel++;
    }
  }

  if(learn_parm->transduction_posratio >= 0) {
    posratio=learn_parm->transduction_posratio;
  }
  else {
    posratio=(double)orgpos/(double)(orgpos+orgneg); /* use ratio of pos/neg */
  }                                                  /* in training data */
  negratio=1.0-posratio;

  learn_parm->svm_costratio=1.0;                     /* global */
  if(posratio>0) {
    learn_parm->svm_costratio_unlab=negratio/posratio;
  }
  else {
    learn_parm->svm_costratio_unlab=1.0;
  }
  
  pos=0;
  neg=0;
  upos=0;
  uneg=0;
  for(i=0;i<totdoc;i++) {
    dist=(lin[i]-model->b);  /* 'distance' from hyperplane*/
    if(dist>0) {
      pos++;
    }
    else {
      neg++;
    }
    if(unlabeled[i]) {
      if(dist>0) {
	upos++;
      }
      else {
	uneg++;
      }
    }
    if((!unlabeled[i]) && (a[i]>(learn_parm->svm_cost[i]-learn_parm->epsilon_a))) {
      /*      printf("Ubounded %ld (class %ld, unlabeled %ld)\n",i,label[i],unlabeled[i]); */
    }
  }
  if(verbosity>=2) {
    printf("POS=%ld, ORGPOS=%ld, ORGNEG=%ld\n",pos,orgpos,orgneg);
    printf("POS=%ld, NEWPOS=%ld, NEWNEG=%ld\n",pos,newpos,newneg);
    printf("pos ratio = %f (%f).\n",(double)(upos)/(double)(allunlab),posratio);
    fflush(stdout);
  }

  if(transductcycle == 0) {
    j1=0; 
    j2=0;
    j4=0;
    for(i=0;i<totdoc;i++) {
      dist=(lin[i]-model->b);  /* 'distance' from hyperplane*/
      if((label[i]==0) && (unlabeled[i])) {
	selcrit[j4]=dist;
	key[j4]=i;
	j4++;
      }
    }
    unsupaddnum1=0;	
    unsupaddnum2=0;	
    select_top_n(selcrit,j4,select,(long)(allunlab*posratio+0.5));
    for(k=0;(k<(long)(allunlab*posratio+0.5));k++) {
      i=key[select[k]];
      label[i]=1;
      unsupaddnum1++;	
      j1++;
    }
    for(i=0;i<totdoc;i++) {
      if((label[i]==0) && (unlabeled[i])) {
	label[i]=-1;
	j2++;
	unsupaddnum2++;
      }
    }
    for(i=0;i<totdoc;i++) {  /* set upper bounds on vars */
      if(unlabeled[i]) {
	if(label[i] == 1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_costratio_unlab*learn_parm->svm_unlabbound;
	}
	else if(label[i] == -1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_unlabbound;
	}
      }
    }
    if(verbosity>=1) {
      /* printf("costratio %f, costratio_unlab %f, unlabbound %f\n",
	 learn_parm->svm_costratio,learn_parm->svm_costratio_unlab,
	 learn_parm->svm_unlabbound); */
      printf("Classifying unlabeled data as %ld POS / %ld NEG.\n",
	     unsupaddnum1,unsupaddnum2); 
      fflush(stdout);
    }
    if(verbosity >= 1) 
      printf("Retraining.");
    if(verbosity >= 2) printf("\n");
    return((long)3);
  }
  if((transductcycle % check_every) == 0) {
    if(verbosity >= 1) 
      printf("Retraining.");
    if(verbosity >= 2) printf("\n");
    j1=0;
    j2=0;
    unsupaddnum1=0;
    unsupaddnum2=0;
    for(i=0;i<totdoc;i++) {
      if((unlabeled[i] == 2)) {
	unlabeled[i]=1;
	label[i]=1;
	j1++;
	unsupaddnum1++;
      }
      else if((unlabeled[i] == 3)) {
	unlabeled[i]=1;
	label[i]=-1;
	j2++;
	unsupaddnum2++;
      }
    }
    for(i=0;i<totdoc;i++) {  /* set upper bounds on vars */
      if(unlabeled[i]) {
	if(label[i] == 1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_costratio_unlab*learn_parm->svm_unlabbound;
	}
	else if(label[i] == -1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_unlabbound;
	}
      }
    }

    if(verbosity>=2) {
      /* printf("costratio %f, costratio_unlab %f, unlabbound %f\n",
	     learn_parm->svm_costratio,learn_parm->svm_costratio_unlab,
	     learn_parm->svm_unlabbound); */
      printf("%ld positive -> Added %ld POS / %ld NEG unlabeled examples.\n",
	     upos,unsupaddnum1,unsupaddnum2); 
      fflush(stdout);
    }

    if(learn_parm->svm_unlabbound == 1) {
      learn_parm->epsilon_crit=0.001; /* do the last run right */
    }
    else {
      learn_parm->epsilon_crit=0.01; /* otherwise, no need to be so picky */
    }

    return((long)3);
  }
  else if(((transductcycle % check_every) < check_every)) { 
    model_length=0;
    sumalpha=0;
    loss=0;
    for(i=0;i<totdoc;i++) {
      model_length+=a[i]*label[i]*lin[i];
      sumalpha+=a[i];
      dist=(lin[i]-model->b);  /* 'distance' from hyperplane*/
      if((label[i]*dist)<(1.0-learn_parm->epsilon_crit)) {
	loss+=(1.0-(label[i]*dist))*learn_parm->svm_cost[i]; 
      }
    }
    model_length=sqrt(model_length); 
    if(verbosity>=2) {
      printf("Model-length = %f (%f), loss = %f, objective = %f\n",
	     model_length,sumalpha,loss,loss+0.5*model_length*model_length);
      fflush(stdout);
    }
    j1=0;
    j2=0;
    j3=0;
    j4=0;
    unsupaddnum1=0;	
    unsupaddnum2=0;	
    umin=99999;
    umax=-99999;
    j4=1;
    while(j4) {
      umin=99999;
      umax=-99999;
      for(i=0;(i<totdoc);i++) { 
	dist=(lin[i]-model->b);  
	if((label[i]>0) && (unlabeled[i]) && (!inconsistent[i]) 
	   && (dist<umin)) {
	  umin=dist;
	  imin=i;
	}
	if((label[i]<0) && (unlabeled[i])  && (!inconsistent[i]) 
	   && (dist>umax)) {
	  umax=dist;
	  imax=i;
	}
      }
      if((umin < (umax+switchsens-1E-4))) {
	j1++;
	j2++;
	unsupaddnum1++;	
	unlabeled[imin]=3;
	inconsistent[imin]=1;
	unsupaddnum2++;	
	unlabeled[imax]=2;
	inconsistent[imax]=1;
      }
      else
	j4=0;
      j4=0;
    }
    for(j=0;(j<totdoc);j++) {
      if(unlabeled[j] && (!inconsistent[j])) {
	if(label[j]>0) {
	  unlabeled[j]=2;
	}
	else if(label[j]<0) {
	  unlabeled[j]=3;
	}
	/* inconsistent[j]=1; */
	j3++;
      }
    }
    switchnum+=unsupaddnum1+unsupaddnum2;

    /* stop and print out current margin
       printf("switchnum %ld %ld\n",switchnum,kernel_parm->poly_degree);
       if(switchnum == 2*kernel_parm->poly_degree) {
       learn_parm->svm_unlabbound=1;
       }
       */

    if((!unsupaddnum1) && (!unsupaddnum2)) {
      if((learn_parm->svm_unlabbound>=1) && ((newpos+newneg) == allunlab)) {
	for(j=0;(j<totdoc);j++) {
	  inconsistent[j]=0;
	  if(unlabeled[j]) unlabeled[j]=1;
	}
	write_prediction(learn_parm->predfile,model,lin,a,unlabeled,label,
			 totdoc,learn_parm);  
	if(verbosity>=1)
	  printf("Number of switches: %ld\n",switchnum);
	return((long)0);
      }
      switchsens=switchsensorg;
      learn_parm->svm_unlabbound*=1.5;
      if(learn_parm->svm_unlabbound>1) {
	learn_parm->svm_unlabbound=1;
      }
      model->at_upper_bound=0; /* since upper bound increased */
      if(verbosity>=1) 
	printf("Increasing influence of unlabeled examples to %f%% .",
	       learn_parm->svm_unlabbound*100.0);
    }
    else if(verbosity>=1) {
      printf("%ld positive -> Switching labels of %ld POS / %ld NEG unlabeled examples.",
	     upos,unsupaddnum1,unsupaddnum2); 
      fflush(stdout);
    }

    if(verbosity >= 2) printf("\n");
    
    learn_parm->epsilon_crit=0.5; /* don't need to be so picky */

    for(i=0;i<totdoc;i++) {  /* set upper bounds on vars */
      if(unlabeled[i]) {
	if(label[i] == 1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_costratio_unlab*learn_parm->svm_unlabbound;
	}
	else if(label[i] == -1) {
	  learn_parm->svm_cost[i]=learn_parm->svm_c*
	    learn_parm->svm_unlabbound;
	}
      }
    }

    return((long)2);
  }

  return((long)0); 
}

/*************************** Working set selection ***************************/

long select_next_qp_subproblem_grad(long int *label, 
				    long int *unlabeled, 
				    double *a, double *lin, 
				    double *c, long int totdoc, 
				    long int qp_size, 
				    LEARN_PARM *learn_parm, 
				    long int *inconsistent, 
				    long int *active2dnum, 
				    long int *working2dnum, 
				    double *selcrit, 
				    long int *select, 
				    KERNEL_CACHE *kernel_cache, 
				    long int cache_only,
				    long int *key, long int *chosen)
     /* Use the feasible direction approach to select the next
      qp-subproblem (see chapter 'Selecting a good working set'). If
      'cache_only' is true, then the variables are selected only among
      those for which the kernel evaluations are cached. */
{
  long choosenum,i,j,k,activedoc,inum,valid;
  double s;

  for(inum=0;working2dnum[inum]>=0;inum++); /* find end of index */
  choosenum=0;
  activedoc=0;
  for(i=0;(j=active2dnum[i])>=0;i++) {
    s=-label[j];
    if(kernel_cache && cache_only) 
      valid=(kernel_cache->index[j]>=0);
    else
      valid=1;
    if(valid
       && (!((a[j]<=(0+learn_parm->epsilon_a)) && (s<0)))
       && (!((a[j]>=(learn_parm->svm_cost[j]-learn_parm->epsilon_a)) 
	     && (s>0)))
       && (!chosen[j]) 
       && (label[j])
       && (!inconsistent[j]))
      {
      selcrit[activedoc]=(double)label[j]*(learn_parm->eps-(double)label[j]*c[j]+(double)label[j]*lin[j]);
      /*      selcrit[activedoc]=(double)label[j]*(-1.0+(double)label[j]*lin[j]); */
      key[activedoc]=j;
      activedoc++;
    }
  }
  select_top_n(selcrit,activedoc,select,(long)(qp_size/2));
  for(k=0;(choosenum<(qp_size/2)) && (k<(qp_size/2)) && (k<activedoc);k++) {
    /* if(learn_parm->biased_hyperplane || (selcrit[select[k]] > 0)) { */
      i=key[select[k]];
      chosen[i]=1;
      working2dnum[inum+choosenum]=i;
      choosenum+=1;
      if(kernel_cache)
	kernel_cache_touch(kernel_cache,i); /* make sure it does not get
					       kicked out of cache */
      /* } */
  }

  activedoc=0;
  for(i=0;(j=active2dnum[i])>=0;i++) {
    s=label[j];
    if(kernel_cache && cache_only) 
      valid=(kernel_cache->index[j]>=0);
    else
      valid=1;
    if(valid
       && (!((a[j]<=(0+learn_parm->epsilon_a)) && (s<0)))
       && (!((a[j]>=(learn_parm->svm_cost[j]-learn_parm->epsilon_a)) 
	     && (s>0))) 
       && (!chosen[j]) 
       && (label[j])
       && (!inconsistent[j])) 
      {
      selcrit[activedoc]=-(double)label[j]*(learn_parm->eps-(double)label[j]*c[j]+(double)label[j]*lin[j]);
      /*  selcrit[activedoc]=-(double)(label[j]*(-1.0+(double)label[j]*lin[j])); */
      key[activedoc]=j;
      activedoc++;
    }
  }
  select_top_n(selcrit,activedoc,select,(long)(qp_size/2));
  for(k=0;(choosenum<qp_size) && (k<(qp_size/2)) && (k<activedoc);k++) {
    /* if(learn_parm->biased_hyperplane || (selcrit[select[k]] > 0)) { */
      i=key[select[k]];
      chosen[i]=1;
      working2dnum[inum+choosenum]=i;
      choosenum+=1;
      if(kernel_cache)
	kernel_cache_touch(kernel_cache,i); /* make sure it does not get
					       kicked out of cache */
      /* } */
  } 
  working2dnum[inum+choosenum]=-1; /* complete index */
  return(choosenum);
}

long select_next_qp_subproblem_rand(long int *label, 
				    long int *unlabeled, 
				    double *a, double *lin, 
				    double *c, long int totdoc, 
				    long int qp_size, 
				    LEARN_PARM *learn_parm, 
				    long int *inconsistent, 
				    long int *active2dnum, 
				    long int *working2dnum, 
				    double *selcrit, 
				    long int *select, 
				    KERNEL_CACHE *kernel_cache, 
				    long int *key, 
				    long int *chosen, 
				    long int iteration)
/* Use the feasible direction approach to select the next
   qp-subproblem (see section 'Selecting a good working set'). Chooses
   a feasible direction at (pseudo) random to help jump over numerical
   problem. */
{
  long choosenum,i,j,k,activedoc,inum;
  double s;

  for(inum=0;working2dnum[inum]>=0;inum++); /* find end of index */
  choosenum=0;
  activedoc=0;
  for(i=0;(j=active2dnum[i])>=0;i++) {
    s=-label[j];
    if((!((a[j]<=(0+learn_parm->epsilon_a)) && (s<0)))
       && (!((a[j]>=(learn_parm->svm_cost[j]-learn_parm->epsilon_a)) 
	     && (s>0)))
       && (!inconsistent[j]) 
       && (label[j])
       && (!chosen[j])) {
      selcrit[activedoc]=(j+iteration) % totdoc;
      key[activedoc]=j;
      activedoc++;
    }
  }
  select_top_n(selcrit,activedoc,select,(long)(qp_size/2));
  for(k=0;(choosenum<(qp_size/2)) && (k<(qp_size/2)) && (k<activedoc);k++) {
    i=key[select[k]];
    chosen[i]=1;
    working2dnum[inum+choosenum]=i;
    choosenum+=1;
    kernel_cache_touch(kernel_cache,i); /* make sure it does not get kicked */
                                        /* out of cache */
  }

  activedoc=0;
  for(i=0;(j=active2dnum[i])>=0;i++) {
    s=label[j];
    if((!((a[j]<=(0+learn_parm->epsilon_a)) && (s<0)))
       && (!((a[j]>=(learn_parm->svm_cost[j]-learn_parm->epsilon_a)) 
	     && (s>0))) 
       && (!inconsistent[j]) 
       && (label[j])
       && (!chosen[j])) {
      selcrit[activedoc]=(j+iteration) % totdoc;
      key[activedoc]=j;
      activedoc++;
    }
  }
  select_top_n(selcrit,activedoc,select,(long)(qp_size/2));
  for(k=0;(choosenum<qp_size) && (k<(qp_size/2)) && (k<activedoc);k++) {
    i=key[select[k]];
    chosen[i]=1;
    working2dnum[inum+choosenum]=i;
    choosenum+=1;
    kernel_cache_touch(kernel_cache,i); /* make sure it does not get kicked */
                                        /* out of cache */
  } 
  working2dnum[inum+choosenum]=-1; /* complete index */
  return(choosenum);
}

long select_next_qp_slackset(DOC **docs, long int *label, 
			     double *a, double *lin, 
			     double *slack, double *alphaslack, 
			     double *c,
			     LEARN_PARM *learn_parm, 
			     long int *active2dnum, double *maxviol)
     /* returns the slackset with the largest internal violation */
{
  long i,ii,maxdiffid;
  double dist,target,maxdiff,ex_c;

  maxdiff=0;
  maxdiffid=0;
  for(ii=0;(i=active2dnum[ii])>=0;ii++) {
    ex_c=learn_parm->svm_c-learn_parm->epsilon_a;
    if(alphaslack[docs[i]->slackid] >= ex_c) {
      dist=(lin[i])*(double)label[i]+slack[docs[i]->slackid]; /* distance */
      target=-(learn_parm->eps-(double)label[i]*c[i]); /* rhs of constraint */
      if((a[i]>learn_parm->epsilon_a) && (dist > target)) {
	if((dist-target)>maxdiff) { /* largest violation */
	  maxdiff=dist-target;
	  maxdiffid=docs[i]->slackid;
	}
      }
    }
  }
  (*maxviol)=maxdiff;
  return(maxdiffid);
}


void select_top_n(double *selcrit, long int range, long int *select, 
		  long int n)
{
  register long i,j;

  for(i=0;(i<n) && (i<range);i++) { /* Initialize with the first n elements */
    for(j=i;j>=0;j--) {
      if((j>0) && (selcrit[select[j-1]]<selcrit[i])){
	select[j]=select[j-1];
      }
      else {
	select[j]=i;
	j=-1;
      }
    }
  }
  if(n>0) {
    for(i=n;i<range;i++) {  
      if(selcrit[i]>selcrit[select[n-1]]) {
	for(j=n-1;j>=0;j--) {
	  if((j>0) && (selcrit[select[j-1]]<selcrit[i])) {
	    select[j]=select[j-1];
	  }
	  else {
	    select[j]=i;
	    j=-1;
	  }
	}
      }
    }
  }
}      
      

/******************************** Shrinking  *********************************/

void init_shrink_state(SHRINK_STATE *shrink_state, long int totdoc, 
		       long int maxhistory)
{
  long i;

  shrink_state->deactnum=0;
  shrink_state->active = (long *)my_malloc(sizeof(long)*totdoc);
  shrink_state->inactive_since = (long *)my_malloc(sizeof(long)*totdoc);
  shrink_state->a_history = (double **)my_malloc(sizeof(double *)*maxhistory);
  shrink_state->maxhistory=maxhistory;
  shrink_state->last_lin = (double *)my_malloc(sizeof(double)*totdoc);
  shrink_state->last_a = (double *)my_malloc(sizeof(double)*totdoc);

  for(i=0;i<totdoc;i++) { 
    shrink_state->active[i]=1;
    shrink_state->inactive_since[i]=0;
    shrink_state->last_a[i]=0;
    shrink_state->last_lin[i]=0;
  }
}

void shrink_state_cleanup(SHRINK_STATE *shrink_state)
{
  free(shrink_state->active);
  free(shrink_state->inactive_since);
  if(shrink_state->deactnum > 0) 
    free(shrink_state->a_history[shrink_state->deactnum-1]);
  free(shrink_state->a_history);
  free(shrink_state->last_a);
  free(shrink_state->last_lin);
}

long shrink_problem(DOC **docs,
		    LEARN_PARM *learn_parm, 
		    SHRINK_STATE *shrink_state, 
		    KERNEL_PARM *kernel_parm,
		    long int *active2dnum, 
		    long int *last_suboptimal_at, 
		    long int iteration, 
		    long int totdoc, 
		    long int minshrink, 
		    double *a, 
		    long int *inconsistent)
     /* Shrink some variables away.  Do the shrinking only if at least
        minshrink variables can be removed. */
{
  long i,ii,change,activenum,lastiter;
  double *a_old;
  
  activenum=0;
  change=0;
  for(ii=0;active2dnum[ii]>=0;ii++) {
    i=active2dnum[ii];
    activenum++;
    if(learn_parm->sharedslack)
      lastiter=last_suboptimal_at[docs[i]->slackid];
    else
      lastiter=last_suboptimal_at[i];
    if(((iteration-lastiter) > learn_parm->svm_iter_to_shrink) 
       || (inconsistent[i])) {
      change++;
    }
  }
  if((change>=minshrink) /* shrink only if sufficiently many candidates */
     && (shrink_state->deactnum<shrink_state->maxhistory)) { /* and enough memory */
    /* Shrink problem by removing those variables which are */
    /* optimal at a bound for a minimum number of iterations */
    if(verbosity>=2) {
      printf(" Shrinking..."); fflush(stdout);
    }
    if(kernel_parm->kernel_type != LINEAR) { /*  non-linear case save alphas */
      a_old=(double *)my_malloc(sizeof(double)*totdoc);
      shrink_state->a_history[shrink_state->deactnum]=a_old;
      for(i=0;i<totdoc;i++) {
	a_old[i]=a[i];
      }
    }
    for(ii=0;active2dnum[ii]>=0;ii++) {
      i=active2dnum[ii];
      if(learn_parm->sharedslack)
	lastiter=last_suboptimal_at[docs[i]->slackid];
      else
	lastiter=last_suboptimal_at[i];
      if(((iteration-lastiter) > learn_parm->svm_iter_to_shrink) 
	 || (inconsistent[i])) {
	shrink_state->active[i]=0;
	shrink_state->inactive_since[i]=shrink_state->deactnum;
      }
    }
    activenum=compute_index(shrink_state->active,totdoc,active2dnum);
    shrink_state->deactnum++;
    if(kernel_parm->kernel_type == LINEAR) { 
      shrink_state->deactnum=0;
    }
    if(verbosity>=2) {
      printf("done.\n"); fflush(stdout);
      printf(" Number of inactive variables = %ld\n",totdoc-activenum);
    }
  }
  return(activenum);
} 


void reactivate_inactive_examples(long int *label, 
				  long int *unlabeled, 
				  double *a, 
				  SHRINK_STATE *shrink_state, 
				  double *lin, 
				  double *c, 
				  long int totdoc, 
				  long int totwords, 
				  long int iteration, 
				  LEARN_PARM *learn_parm, 
				  long int *inconsistent, 
				  DOC **docs, 
				  KERNEL_PARM *kernel_parm, 
				  KERNEL_CACHE *kernel_cache, 
				  MODEL *model, 
				  CFLOAT *aicache, 
				  double *weights, 
				  double *maxdiff)
     /* Make all variables active again which had been removed by
        shrinking. */
     /* Computes lin for those variables from scratch. */
{
  register long i,j,ii,jj,t,*changed2dnum,*inactive2dnum;
  long *changed,*inactive;
  register double kernel_val,*a_old,dist;
  double ex_c,target;
  SVECTOR *f;

  if(kernel_parm->kernel_type == LINEAR) { /* special linear case */
    a_old=shrink_state->last_a;    
    clear_vector_n(weights,totwords);
    for(i=0;i<totdoc;i++) {
      if(a[i] != a_old[i]) {
	for(f=docs[i]->fvec;f;f=f->next)  
	  add_vector_ns(weights,f,
			f->factor*((a[i]-a_old[i])*(double)label[i]));
	a_old[i]=a[i];
      }
    }
    for(i=0;i<totdoc;i++) {
      if(!shrink_state->active[i]) {
	for(f=docs[i]->fvec;f;f=f->next)  
	  lin[i]=shrink_state->last_lin[i]+f->factor*sprod_ns(weights,f);
      }
      shrink_state->last_lin[i]=lin[i];
    }
  }
  else {
    changed=(long *)my_malloc(sizeof(long)*totdoc);
    changed2dnum=(long *)my_malloc(sizeof(long)*(totdoc+11));
    inactive=(long *)my_malloc(sizeof(long)*totdoc);
    inactive2dnum=(long *)my_malloc(sizeof(long)*(totdoc+11));
    for(t=shrink_state->deactnum-1;(t>=0) && shrink_state->a_history[t];t--) {
      if(verbosity>=2) {
	printf("%ld..",t); fflush(stdout);
      }
      a_old=shrink_state->a_history[t];    
      for(i=0;i<totdoc;i++) {
	inactive[i]=((!shrink_state->active[i]) 
		     && (shrink_state->inactive_since[i] == t));
	changed[i]= (a[i] != a_old[i]);
      }
      compute_index(inactive,totdoc,inactive2dnum);
      compute_index(changed,totdoc,changed2dnum);
      
      for(ii=0;(i=changed2dnum[ii])>=0;ii++) {
	get_kernel_row(kernel_cache,docs,i,totdoc,inactive2dnum,aicache,
		       kernel_parm);
	for(jj=0;(j=inactive2dnum[jj])>=0;jj++) {
	  kernel_val=aicache[j];
	  lin[j]+=(((a[i]*kernel_val)-(a_old[i]*kernel_val))*(double)label[i]);
	}
      }
    }
    free(changed);
    free(changed2dnum);
    free(inactive);
    free(inactive2dnum);
  }
  (*maxdiff)=0;
  for(i=0;i<totdoc;i++) {
    shrink_state->inactive_since[i]=shrink_state->deactnum-1;
    if(!inconsistent[i]) {
      dist=(lin[i]-model->b)*(double)label[i];
      target=-(learn_parm->eps-(double)label[i]*c[i]);
      ex_c=learn_parm->svm_cost[i]-learn_parm->epsilon_a;
      if((a[i]>learn_parm->epsilon_a) && (dist > target)) {
	if((dist-target)>(*maxdiff))  /* largest violation */
	  (*maxdiff)=dist-target;
      }
      else if((a[i]<ex_c) && (dist < target)) {
	if((target-dist)>(*maxdiff))  /* largest violation */
	  (*maxdiff)=target-dist;
      }
      if((a[i]>(0+learn_parm->epsilon_a)) 
	 && (a[i]<ex_c)) { 
	shrink_state->active[i]=1;                         /* not at bound */
      }
      else if((a[i]<=(0+learn_parm->epsilon_a)) && (dist < (target+learn_parm->epsilon_shrink))) {
	shrink_state->active[i]=1;
      }
      else if((a[i]>=ex_c)
	      && (dist > (target-learn_parm->epsilon_shrink))) {
	shrink_state->active[i]=1;
      }
      else if(learn_parm->sharedslack) { /* make all active when sharedslack */
	shrink_state->active[i]=1;
      }
    }
  }
  if(kernel_parm->kernel_type != LINEAR) { /* update history for non-linear */
    for(i=0;i<totdoc;i++) {
      (shrink_state->a_history[shrink_state->deactnum-1])[i]=a[i];
    }
    for(t=shrink_state->deactnum-2;(t>=0) && shrink_state->a_history[t];t--) {
      free(shrink_state->a_history[t]);
      shrink_state->a_history[t]=0;
    }
  }
}

/****************************** Cache handling *******************************/

void get_kernel_row(KERNEL_CACHE *kernel_cache, DOC **docs, 
		    long int docnum, long int totdoc, 
		    long int *active2dnum, CFLOAT *buffer, 
		    KERNEL_PARM *kernel_parm)
     /* Get's a row of the matrix of kernel values This matrix has the
      same form as the Hessian, just that the elements are not
      multiplied by */
     /* y_i * y_j * a_i * a_j */
     /* Takes the values from the cache if available. */
{
  register long i,j,start;
  DOC *ex;

  ex=docs[docnum];

  if(kernel_cache->index[docnum] != -1) { /* row is cached? */
    kernel_cache->lru[kernel_cache->index[docnum]]=kernel_cache->time; /* lru */
    start=kernel_cache->activenum*kernel_cache->index[docnum];
    for(i=0;(j=active2dnum[i])>=0;i++) {
      if(kernel_cache->totdoc2active[j] >= 0) { /* column is cached? */
	buffer[j]=kernel_cache->buffer[start+kernel_cache->totdoc2active[j]];
      }
      else {
	buffer[j]=(CFLOAT)kernel(kernel_parm,ex,docs[j]);
      }
    }
  }
  else {
    for(i=0;(j=active2dnum[i])>=0;i++) {
      buffer[j]=(CFLOAT)kernel(kernel_parm,ex,docs[j]);
    }
  }
}


void cache_kernel_row(KERNEL_CACHE *kernel_cache, DOC **docs, 
		      long int m, KERNEL_PARM *kernel_parm)
     /* Fills cache for the row m */
{
  register DOC *ex;
  register long j,k,l;
  register CFLOAT *cache;

  if(!kernel_cache_check(kernel_cache,m)) {  /* not cached yet*/
    cache = kernel_cache_clean_and_malloc(kernel_cache,m);
    if(cache) {
      l=kernel_cache->totdoc2active[m];
      ex=docs[m];
      for(j=0;j<kernel_cache->activenum;j++) {  /* fill cache */
	k=kernel_cache->active2totdoc[j];
	if((kernel_cache->index[k] != -1) && (l != -1) && (k != m)) {
	  cache[j]=kernel_cache->buffer[kernel_cache->activenum
				       *kernel_cache->index[k]+l];
	}
	else {
	  cache[j]=kernel(kernel_parm,ex,docs[k]);
	} 
      }
    }
    else {
      perror("Error: Kernel cache full! => increase cache size");
    }
  }
}

 
void cache_multiple_kernel_rows(KERNEL_CACHE *kernel_cache, DOC **docs, 
				long int *key, long int varnum, 
				KERNEL_PARM *kernel_parm)
     /* Fills cache for the rows in key */
{
  register long i;

  for(i=0;i<varnum;i++) {  /* fill up kernel cache */
    cache_kernel_row(kernel_cache,docs,key[i],kernel_parm);
  }
}

 
void kernel_cache_shrink(KERNEL_CACHE *kernel_cache, long int totdoc, 
			 long int numshrink, long int *after)
     /* Remove numshrink columns in the cache which correspond to
        examples marked 0 in after. */
{
  register long i,j,jj,from=0,to=0,scount;  
  long *keep;

  if(verbosity>=2) {
    printf(" Reorganizing cache..."); fflush(stdout);
  }

  keep=(long *)my_malloc(sizeof(long)*totdoc);
  for(j=0;j<totdoc;j++) {
    keep[j]=1;
  }
  scount=0;
  for(jj=0;(jj<kernel_cache->activenum) && (scount<numshrink);jj++) {
    j=kernel_cache->active2totdoc[jj];
    if(!after[j]) {
      scount++;
      keep[j]=0;
    }
  }

  for(i=0;i<kernel_cache->max_elems;i++) {
    for(jj=0;jj<kernel_cache->activenum;jj++) {
      j=kernel_cache->active2totdoc[jj];
      if(!keep[j]) {
	from++;
      }
      else {
	kernel_cache->buffer[to]=kernel_cache->buffer[from];
	to++;
	from++;
      }
    }
  }

  kernel_cache->activenum=0;
  for(j=0;j<totdoc;j++) {
    if((keep[j]) && (kernel_cache->totdoc2active[j] != -1)) {
      kernel_cache->active2totdoc[kernel_cache->activenum]=j;
      kernel_cache->totdoc2active[j]=kernel_cache->activenum;
      kernel_cache->activenum++;
    }
    else {
      kernel_cache->totdoc2active[j]=-1;
    }
  }

  kernel_cache->max_elems=(long)(kernel_cache->buffsize/kernel_cache->activenum);
  if(kernel_cache->max_elems>totdoc) {
    kernel_cache->max_elems=totdoc;
  }

  free(keep);

  if(verbosity>=2) {
    printf("done.\n"); fflush(stdout);
    printf(" Cache-size in rows = %ld\n",kernel_cache->max_elems);
  }
}

KERNEL_CACHE *kernel_cache_init(long int totdoc, long int buffsize)
{
  long i;
  KERNEL_CACHE *kernel_cache;

  kernel_cache=(KERNEL_CACHE *)my_malloc(sizeof(KERNEL_CACHE));
  kernel_cache->index = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->occu = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->lru = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->invindex = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->active2totdoc = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->totdoc2active = (long *)my_malloc(sizeof(long)*totdoc);
  kernel_cache->buffer = (CFLOAT *)my_malloc((size_t)(buffsize)*1024*1024);

  kernel_cache->buffsize=(long)(buffsize/sizeof(CFLOAT)*1024*1024);

  kernel_cache->max_elems=(long)(kernel_cache->buffsize/totdoc);
  if(kernel_cache->max_elems>totdoc) {
    kernel_cache->max_elems=totdoc;
  }

  if(verbosity>=2) {
    printf(" Cache-size in rows = %ld\n",kernel_cache->max_elems);
    printf(" Kernel evals so far: %ld\n",kernel_cache_statistic);    
  }

  kernel_cache->elems=0;   /* initialize cache */
  for(i=0;i<totdoc;i++) {
    kernel_cache->index[i]=-1;
    kernel_cache->lru[i]=0;
  }
  for(i=0;i<totdoc;i++) {
    kernel_cache->occu[i]=0;
    kernel_cache->invindex[i]=-1;
  }

  kernel_cache->activenum=totdoc;;
  for(i=0;i<totdoc;i++) {
      kernel_cache->active2totdoc[i]=i;
      kernel_cache->totdoc2active[i]=i;
  }

  kernel_cache->time=0;  

  return(kernel_cache);
} 

void kernel_cache_reset_lru(KERNEL_CACHE *kernel_cache)
{
  long maxlru=0,k;
  
  for(k=0;k<kernel_cache->max_elems;k++) {
    if(maxlru < kernel_cache->lru[k]) 
      maxlru=kernel_cache->lru[k];
  }
  for(k=0;k<kernel_cache->max_elems;k++) {
      kernel_cache->lru[k]-=maxlru;
  }
}

void kernel_cache_cleanup(KERNEL_CACHE *kernel_cache)
{
  free(kernel_cache->index);
  free(kernel_cache->occu);
  free(kernel_cache->lru);
  free(kernel_cache->invindex);
  free(kernel_cache->active2totdoc);
  free(kernel_cache->totdoc2active);
  free(kernel_cache->buffer);
  free(kernel_cache);
}

long kernel_cache_malloc(KERNEL_CACHE *kernel_cache)
{
  long i;

  if(kernel_cache_space_available(kernel_cache)) {
    for(i=0;i<kernel_cache->max_elems;i++) {
      if(!kernel_cache->occu[i]) {
	kernel_cache->occu[i]=1;
	kernel_cache->elems++;
	return(i);
      }
    }
  }
  return(-1);
}

void kernel_cache_free(KERNEL_CACHE *kernel_cache, long int i)
{
  kernel_cache->occu[i]=0;
  kernel_cache->elems--;
}

long kernel_cache_free_lru(KERNEL_CACHE *kernel_cache) 
     /* remove least recently used cache element */
{                                     
  register long k,least_elem=-1,least_time;

  least_time=kernel_cache->time+1;
  for(k=0;k<kernel_cache->max_elems;k++) {
    if(kernel_cache->invindex[k] != -1) {
      if(kernel_cache->lru[k]<least_time) {
	least_time=kernel_cache->lru[k];
	least_elem=k;
      }
    }
  }
  if(least_elem != -1) {
    kernel_cache_free(kernel_cache,least_elem);
    kernel_cache->index[kernel_cache->invindex[least_elem]]=-1;
    kernel_cache->invindex[least_elem]=-1;
    return(1);
  }
  return(0);
}


CFLOAT *kernel_cache_clean_and_malloc(KERNEL_CACHE *kernel_cache, 
				      long int docnum)
     /* Get a free cache entry. In case cache is full, the lru element
        is removed. */
{
  long result;
  if((result = kernel_cache_malloc(kernel_cache)) == -1) {
    if(kernel_cache_free_lru(kernel_cache)) {
      result = kernel_cache_malloc(kernel_cache);
    }
  }
  kernel_cache->index[docnum]=result;
  if(result == -1) {
    return(0);
  }
  kernel_cache->invindex[result]=docnum;
  kernel_cache->lru[kernel_cache->index[docnum]]=kernel_cache->time; /* lru */
  return((CFLOAT *)((long)kernel_cache->buffer
		    +(kernel_cache->activenum*sizeof(CFLOAT)*
		      kernel_cache->index[docnum])));
}

long kernel_cache_touch(KERNEL_CACHE *kernel_cache, long int docnum)
     /* Update lru time to avoid removal from cache. */
{
  if(kernel_cache && kernel_cache->index[docnum] != -1) {
    kernel_cache->lru[kernel_cache->index[docnum]]=kernel_cache->time; /* lru */
    return(1);
  }
  return(0);
}
  
long kernel_cache_check(KERNEL_CACHE *kernel_cache, long int docnum)
     /* Is that row cached? */
{
  return(kernel_cache->index[docnum] != -1);
}
  
long kernel_cache_space_available(KERNEL_CACHE *kernel_cache)
     /* Is there room for one more row? */
{
  return(kernel_cache->elems < kernel_cache->max_elems);
}
  
/************************** Compute estimates ******************************/

void compute_xa_estimates(MODEL *model, long int *label, 
			  long int *unlabeled, long int totdoc, 
			  DOC **docs, double *lin, double *a, 
			  KERNEL_PARM *kernel_parm, 
			  LEARN_PARM *learn_parm, double *error, 
			  double *recall, double *precision) 
     /* Computes xa-estimate of error rate, recall, and precision. See
        T. Joachims, Estimating the Generalization Performance of an SVM
        Efficiently, IMCL, 2000. */
{
  long i,looerror,looposerror,loonegerror;
  long totex,totposex;
  double xi,r_delta,r_delta_sq,sim=0;
  long *sv2dnum=NULL,*sv=NULL,svnum;

  r_delta=estimate_r_delta(docs,totdoc,kernel_parm); 
  r_delta_sq=r_delta*r_delta;

  looerror=0;
  looposerror=0;
  loonegerror=0;
  totex=0;
  totposex=0;
  svnum=0;

  if(learn_parm->xa_depth > 0) {
    sv = (long *)my_malloc(sizeof(long)*(totdoc+11));
    for(i=0;i<totdoc;i++) 
      sv[i]=0;
    for(i=1;i<model->sv_num;i++) 
      if(a[model->supvec[i]->docnum] 
	 < (learn_parm->svm_cost[model->supvec[i]->docnum]
	    -learn_parm->epsilon_a)) {
	sv[model->supvec[i]->docnum]=1;
	svnum++;
      }
    sv2dnum = (long *)my_malloc(sizeof(long)*(totdoc+11));
    clear_index(sv2dnum);
    compute_index(sv,totdoc,sv2dnum);
  }

  for(i=0;i<totdoc;i++) {
    if(unlabeled[i]) {
      /* ignore it */
    }
    else {
      xi=1.0-((lin[i]-model->b)*(double)label[i]);
      if(xi<0) xi=0;
      if(label[i]>0) {
	totposex++;
      }
      if((learn_parm->rho*a[i]*r_delta_sq+xi) >= 1.0) {
	if(learn_parm->xa_depth > 0) {  /* makes assumptions */
	  sim=distribute_alpha_t_greedily(sv2dnum,svnum,docs,a,i,label,
					  kernel_parm,learn_parm,
		            (double)((1.0-xi-a[i]*r_delta_sq)/(2.0*a[i])));
	}
	if((learn_parm->xa_depth == 0) || 
	   ((a[i]*kernel(kernel_parm,docs[i],docs[i])+a[i]*2.0*sim+xi) >= 1.0)) { 
	  looerror++;
	  if(label[i]>0) {
	    looposerror++;
	  }
	  else {
	    loonegerror++;
	  }
	}
      }
      totex++;
    }
  }

  (*error)=((double)looerror/(double)totex)*100.0;
  (*recall)=(1.0-(double)looposerror/(double)totposex)*100.0;
  (*precision)=(((double)totposex-(double)looposerror)
    /((double)totposex-(double)looposerror+(double)loonegerror))*100.0;

  free(sv);
  free(sv2dnum);
}


double distribute_alpha_t_greedily(long int *sv2dnum, long int svnum, 
				   DOC **docs, double *a, 
				   long int docnum, 
				   long int *label, 
				   KERNEL_PARM *kernel_parm, 
				   LEARN_PARM *learn_parm, double thresh)
     /* Experimental Code improving plain XiAlpha Estimates by
	computing a better bound using a greedy optimzation strategy. */
{
  long best_depth=0;
  long i,j,k,d,skip,allskip;
  double best,best_val[101],val,init_val_sq,init_val_lin;
  long best_ex[101];
  CFLOAT *cache,*trow;

  cache=(CFLOAT *)my_malloc(sizeof(CFLOAT)*learn_parm->xa_depth*svnum);
  trow = (CFLOAT *)my_malloc(sizeof(CFLOAT)*svnum);

  for(k=0;k<svnum;k++) {
    trow[k]=kernel(kernel_parm,docs[docnum],docs[sv2dnum[k]]);
  }

  init_val_sq=0;
  init_val_lin=0;
  best=0;

  for(d=0;d<learn_parm->xa_depth;d++) {
    allskip=1;
    if(d>=1) {
      init_val_sq+=cache[best_ex[d-1]+svnum*(d-1)]; 
      for(k=0;k<d-1;k++) {
        init_val_sq+=2.0*cache[best_ex[k]+svnum*(d-1)]; 
      }
      init_val_lin+=trow[best_ex[d-1]]; 
    }
    for(i=0;i<svnum;i++) {
      skip=0;
      if(sv2dnum[i] == docnum) skip=1;
      for(j=0;j<d;j++) {
	if(i == best_ex[j]) skip=1;
      }

      if(!skip) {
	val=init_val_sq;
	if(kernel_parm->kernel_type == LINEAR) 
	  val+=docs[sv2dnum[i]]->fvec->twonorm_sq;
	else
	  val+=kernel(kernel_parm,docs[sv2dnum[i]],docs[sv2dnum[i]]);
	for(j=0;j<d;j++) {
	  val+=2.0*cache[i+j*svnum];
	}
	val*=(1.0/(2.0*(d+1.0)*(d+1.0)));
	val-=((init_val_lin+trow[i])/(d+1.0));

	if(allskip || (val < best_val[d])) {
	  best_val[d]=val;
	  best_ex[d]=i;
	}
	allskip=0;
	if(val < thresh) {
	  i=svnum;
	  /*	  printf("EARLY"); */
	}
      }
    }
    if(!allskip) {
      for(k=0;k<svnum;k++) {
	  cache[d*svnum+k]=kernel(kernel_parm,
				  docs[sv2dnum[best_ex[d]]],
				  docs[sv2dnum[k]]);
      }
    }
    if((!allskip) && ((best_val[d] < best) || (d == 0))) {
      best=best_val[d];
      best_depth=d;
    }
    if(allskip || (best < thresh)) {
      d=learn_parm->xa_depth;
    }
  }    

  free(cache);
  free(trow);

  /*  printf("Distribute[%ld](%ld)=%f, ",docnum,best_depth,best); */
  return(best);
}


void estimate_transduction_quality(MODEL *model, long int *label, 
				   long int *unlabeled, 
				   long int totdoc, DOC **docs, double *lin) 
     /* Loo-bound based on observation that loo-errors must have an
	equal distribution in both training and test examples, given
	that the test examples are classified correctly. Compare
	chapter "Constraints on the Transductive Hyperplane" in my
	Dissertation. */
{
  long i,j,l=0,ulab=0,lab=0,labpos=0,labneg=0,ulabpos=0,ulabneg=0,totulab=0;
  double totlab=0,totlabpos=0,totlabneg=0,labsum=0,ulabsum=0;
  double r_delta,r_delta_sq,xi,xisum=0,asum=0;

  r_delta=estimate_r_delta(docs,totdoc,&(model->kernel_parm)); 
  r_delta_sq=r_delta*r_delta;

  for(j=0;j<totdoc;j++) {
    if(unlabeled[j]) {
      totulab++;
    }
    else {
      totlab++;
      if(label[j] > 0) 
	totlabpos++;
      else 
	totlabneg++;
    }
  }
  for(j=1;j<model->sv_num;j++) {
    i=model->supvec[j]->docnum;
    xi=1.0-((lin[i]-model->b)*(double)label[i]);
    if(xi<0) xi=0;

    xisum+=xi;
    asum+=fabs(model->alpha[j]);
    if(unlabeled[i]) {
      ulabsum+=(fabs(model->alpha[j])*r_delta_sq+xi);
    }
    else {
      labsum+=(fabs(model->alpha[j])*r_delta_sq+xi);
    }
    if((fabs(model->alpha[j])*r_delta_sq+xi) >= 1) { 
      l++;
      if(unlabeled[model->supvec[j]->docnum]) {
	ulab++;
	if(model->alpha[j] > 0) 
	  ulabpos++;
	else 
	  ulabneg++;
      }
      else {
	lab++;
	if(model->alpha[j] > 0) 
	  labpos++;
	else 
	  labneg++;
      }
    }
  }
  printf("xacrit>=1: labeledpos=%.5f labeledneg=%.5f default=%.5f\n",(double)labpos/(double)totlab*100.0,(double)labneg/(double)totlab*100.0,(double)totlabpos/(double)(totlab)*100.0);
  printf("xacrit>=1: unlabelpos=%.5f unlabelneg=%.5f\n",(double)ulabpos/(double)totulab*100.0,(double)ulabneg/(double)totulab*100.0);
  printf("xacrit>=1: labeled=%.5f unlabled=%.5f all=%.5f\n",(double)lab/(double)totlab*100.0,(double)ulab/(double)totulab*100.0,(double)l/(double)(totdoc)*100.0);
  printf("xacritsum: labeled=%.5f unlabled=%.5f all=%.5f\n",(double)labsum/(double)totlab*100.0,(double)ulabsum/(double)totulab*100.0,(double)(labsum+ulabsum)/(double)(totdoc)*100.0);
  printf("r_delta_sq=%.5f xisum=%.5f asum=%.5f\n",r_delta_sq,xisum,asum);
}

double estimate_margin_vcdim(MODEL *model, double w, double R, 
			     KERNEL_PARM *kernel_parm) 
     /* optional: length of model vector in feature space */
     /* optional: radius of ball containing the data */
{
  double h;

  /* follows chapter 5.6.4 in [Vapnik/95] */

  if(w<0) {
    w=model_length_s(model,kernel_parm);
  }
  if(R<0) {
    R=estimate_sphere(model,kernel_parm); 
  }
  h = w*w * R*R +1; 
  return(h);
}

double estimate_sphere(MODEL *model, KERNEL_PARM *kernel_parm) 
                          /* Approximates the radius of the ball containing */
                          /* the support vectors by bounding it with the */
{                         /* length of the longest support vector. This is */
  register long j;        /* pretty good for text categorization, since all */
  double xlen,maxxlen=0;  /* documents have feature vectors of length 1. It */
  DOC *nulldoc;           /* assumes that the center of the ball is at the */
  WORD nullword;          /* origin of the space. */

  nullword.wnum=0;
  nulldoc=create_example(-2,0,0,0.0,create_svector(&nullword,"",1.0)); 

  for(j=1;j<model->sv_num;j++) {
    xlen=sqrt(kernel(kernel_parm,model->supvec[j],model->supvec[j])
	      -2*kernel(kernel_parm,model->supvec[j],nulldoc)
	      +kernel(kernel_parm,nulldoc,nulldoc));
    if(xlen>maxxlen) {
      maxxlen=xlen;
    }
  }

  free_example(nulldoc,1);
  return(maxxlen);
}

double estimate_r_delta(DOC **docs, long int totdoc, KERNEL_PARM *kernel_parm)
{
  long i;
  double maxxlen,xlen;
  DOC *nulldoc;           /* assumes that the center of the ball is at the */
  WORD nullword;          /* origin of the space. */

  nullword.wnum=0;
  nulldoc=create_example(-2,0,0,0.0,create_svector(&nullword,"",1.0)); 

  maxxlen=0;
  for(i=0;i<totdoc;i++) {
    xlen=sqrt(kernel(kernel_parm,docs[i],docs[i])
	      -2*kernel(kernel_parm,docs[i],nulldoc)
	      +kernel(kernel_parm,nulldoc,nulldoc));
    if(xlen>maxxlen) {
      maxxlen=xlen;
    }
  }

  free_example(nulldoc,1);
  return(maxxlen);
}

double estimate_r_delta_average(DOC **docs, long int totdoc, 
				KERNEL_PARM *kernel_parm)
{
  long i;
  double avgxlen;
  DOC *nulldoc;           /* assumes that the center of the ball is at the */
  WORD nullword;          /* origin of the space. */

  nullword.wnum=0;
  nulldoc=create_example(-2,0,0,0.0,create_svector(&nullword,"",1.0)); 

  avgxlen=0;
  for(i=0;i<totdoc;i++) {
    avgxlen+=sqrt(kernel(kernel_parm,docs[i],docs[i])
		  -2*kernel(kernel_parm,docs[i],nulldoc)
		  +kernel(kernel_parm,nulldoc,nulldoc));
  }

  free_example(nulldoc,1);
  return(avgxlen/totdoc);
}

double length_of_longest_document_vector(DOC **docs, long int totdoc, 
					 KERNEL_PARM *kernel_parm)
{
  long i;
  double maxxlen,xlen;

  maxxlen=0;
  for(i=0;i<totdoc;i++) {
    xlen=sqrt(kernel(kernel_parm,docs[i],docs[i]));
    if(xlen>maxxlen) {
      maxxlen=xlen;
    }
  }

  return(maxxlen);
}

/****************************** IO-handling **********************************/

void write_prediction(char *predfile, MODEL *model, double *lin, 
		      double *a, long int *unlabeled, 
		      long int *label, long int totdoc, 
		      LEARN_PARM *learn_parm)
{
  FILE *predfl;
  long i;
  double dist,a_max;

  if(verbosity>=1) {
    printf("Writing prediction file..."); fflush(stdout);
  }
  if ((predfl = fopen (predfile, "w")) == NULL)
  { perror (predfile); exit (1); }
  a_max=learn_parm->epsilon_a;
  for(i=0;i<totdoc;i++) {
    if((unlabeled[i]) && (a[i]>a_max)) {
      a_max=a[i];
    }
  }
  for(i=0;i<totdoc;i++) {
    if(unlabeled[i]) {
      if((a[i]>(learn_parm->epsilon_a))) {
	dist=(double)label[i]*(1.0-learn_parm->epsilon_crit-a[i]/(a_max*2.0));
      }
      else {
	dist=(lin[i]-model->b);
      }
      if(dist>0) {
	fprintf(predfl,"%.8g:+1 %.8g:-1\n",dist,-dist);
      }
      else {
	fprintf(predfl,"%.8g:-1 %.8g:+1\n",-dist,dist);
      }
    }
  }
  fclose(predfl);
  if(verbosity>=1) {
    printf("done\n");
  }
}

void write_alphas(char *alphafile, double *a, 
		  long int *label, long int totdoc)
{
  FILE *alphafl;
  long i;

  if(verbosity>=1) {
    printf("Writing alpha file..."); fflush(stdout);
  }
  if ((alphafl = fopen (alphafile, "w")) == NULL)
  { perror (alphafile); exit (1); }
  for(i=0;i<totdoc;i++) {
    fprintf(alphafl,"%.18g\n",a[i]*(double)label[i]);
  }
  fclose(alphafl);
  if(verbosity>=1) {
    printf("done\n");
  }
}