loo.c

留一模型选择法leave-one-out model selection

#include <math.h>
#include <ctype.h>
#include <float.h>
#include <string.h>
#include <assert.h>
#include <limits.h>
#include <signal.h>
#include "bsvm.h"
#include "bsvm_util.h"
#include "loo.h"

svm_param_t svm_param;
kernel_param_t kernel_param;
static int N;
static int Nsv=0, Nfree=0;
static sample_t *sample;
static int *workingset, *chosenflag;
static int iter=0, old_iter=0;
static double *x, *Qx, *sortbase;
static double *loo_init_x, *loo_init_Qx;
static int *SV;
static BQP qp;

static double *positive_max=0, *negative_max=0;
static int *positive_set=0, *negative_set=0;

/*		min 0.5*x'*Qx - e'*x
subject to 	0 <= x <= svm_param.C
		Qij = sample[i].a*sample[j].a*kernel_function(sample+i, sample+j)
		Qx = Q*x     */

static double **cache, **QB;
static int *Q2cache, *cache2Q, usingcacherows;
static double *cache2G;

/* kernel_eval may be greater than 2^32 */
static double kernel_eval;

/*	QB: Q[workingset]
	Q2cache: row Q[i] is cached in cache[Q2cache[i]] (-1 means not cached)
	cache2Q: cache[i] save row Q[cache2Q[i]] */

double (*kernel_function)(sample_t*,sample_t*);

static int compare_by_Qx(const void *first, const void *second)
{
	double key1, key2;
	key1 = Qx[*((int*)first)];
	key2 = Qx[*((int*)second)];
	if(key1>key2)
	{
		return 1;
	}
	else if(key1<key2)
	{
		return -1;
	}
	else
	{
		return 0;
	}
}

static void sort_SV_by_Qx()
{
	qsort(SV, Nsv, sizeof(int), compare_by_Qx);
}

static void store_Qx_and_x()
{
	int i;
	
	for(i=0; i<N; i++)
	{
		loo_init_x[i] = x[i];
		loo_init_Qx[i] = Qx[i];
	}
}

static void update_Qx_by_remove_sv(int testv)
{
	// The following does this: Qx[i] = (Qx - x[testv]*Q[testv])[i]
	int i;
	
	if(Q2cache[testv] != -1)
	{
		double g = cache2G[Q2cache[testv]];
		if(g == kernel_param.gamma)
		{
			for(i=0; i<N; i++)
			{
				Qx[i] = loo_init_Qx[i] - loo_init_x[testv]*cache[Q2cache[testv]][i];
			}
		}
		else
		{
			double new_cache, old_cache;
			double r = kernel_param.gamma/g;
			cache2G[Q2cache[testv]] = kernel_param.gamma;
			
			for(i=0; i<N; i++)
			{
				old_cache = cache[Q2cache[testv]][i];
				if(sample[i].a!=sample[testv].a)
					old_cache = -old_cache-1;
				else
					old_cache--;
				
				assert(old_cache>=0);
				new_cache = pow(old_cache, r)+1;
				
				if(sample[i].a!=sample[testv].a)
					new_cache = -new_cache;

				cache[Q2cache[testv]][i] = new_cache;
				Qx[i] = loo_init_Qx[i] - loo_init_x[testv]*new_cache;
			}
		}
	}
	else
	{
		double v;
		if(svm_param.verbosity >= 3) printf("miss!\n");
		for(i=0; i<N; i++)
		{
			//if(i==testv) continue; // Qx[testv] won't be used
			v = kernel_function(&sample[testv], &sample[i]);
			if(sample[testv].a != sample[i].a)
				v = -v;
			Qx[i] = loo_init_Qx[i] - loo_init_x[testv]*v;
		}

		kernel_eval += N;
	}
}

double decision_function(sample_t *testsample)
{
	int i;
	double f = 0;

	for (i=0;i<N;i++)
		if(x[i]>ZERO)
			f += sample[i].a*x[i]*kernel_function(testsample, &sample[i]);
	
	return f;
}

/* least-recently-used strategy */
static int LRU()
{
	int i, j;
	
	if (usingcacherows < svm_param.cachesize)
		i = usingcacherows++;
	else
	{
		for (i=0,j=1;j<usingcacherows;j++)
			if (chosenflag[cache2Q[i]] > chosenflag[cache2Q[j]])
				i = j;
		Q2cache[cache2Q[i]] = -1;
	}
	return i;
}

static void get_QB(int q)
{
	int i, j;
	
	for (i=0;i<q;i++)
	{
		int wi = workingset[i];
		if (Q2cache[wi] != -1)
		{
			double g = cache2G[Q2cache[wi]];
			if(g != kernel_param.gamma)
			{
				double new_cache, old_cache;
				double r = kernel_param.gamma/g;
				cache2G[Q2cache[wi]] = kernel_param.gamma;
				
				for(j=0; j<N; j++)
				{
					old_cache = cache[Q2cache[wi]][j];
					if(sample[j].a!=sample[wi].a)
						old_cache = -old_cache-1;
					else
						old_cache--;
				
					assert(old_cache>=0);
					new_cache = pow(old_cache, r)+1;
				
					if(sample[j].a!=sample[wi].a)
						new_cache = -new_cache;

					cache[Q2cache[wi]][j] = new_cache;
				}
			}

			QB[i] = cache[Q2cache[wi]];
		}
		else
		{
			/* miss */
			int cache_index = LRU();
			for (j=0;j<N;j++)
			{
				double v = kernel_function(&sample[wi], &sample[j]);
				if(sample[wi].a != sample[j].a)
					v = -v;
				cache[cache_index][j] = v;
			}
			kernel_eval += N;
			Q2cache[wi] = cache_index;
			cache2Q[cache_index] = wi;
			cache2G[cache_index] = kernel_param.gamma;
			QB[i] = cache[cache_index];
		}
	}
}

void initialize_cache()
{
	int i;

	svm_param.cachesize = svm_param.cachesize*(1L << 20)/sizeof(double)/N;
	if (svm_param.cachesize > N)
		svm_param.cachesize = N;
	if (svm_param.cachesize < svm_param.qpsize)
		myerror("Cache size is too small\n");
	QB = (double **) xmalloc(sizeof(double *)*svm_param.qpsize);
	cache = (double **) xmalloc(sizeof(double *)*svm_param.cachesize);
	for (i=0;i<svm_param.cachesize;i++)
		cache[i] = (double *) xmalloc(sizeof(double)*N);
	Q2cache = (int *) xmalloc(sizeof(int)*N);
	cache2G = (double *) xmalloc(sizeof(double)*svm_param.cachesize);
	for (i=0;i<N;i++)
		Q2cache[i] = -1;
	cache2Q = (int *) xmalloc(sizeof(int)*svm_param.cachesize);
	usingcacherows = 0;
	kernel_eval = 0;
}

void free_cache()
{
	int i;
	
	free(QB);
	for (i=0;i<svm_param.cachesize;i++)
		free(cache[i]);
	free(cache);
	free(Q2cache);
	free(cache2Q);
}

double infnorm_of_grad()
{
	int i;
	double maxdiff = 0, diff;

	for (i=0;i<N;i++)
	{
//		if(i==testv) continue;
		diff = 1 - Qx[i];
		if ((diff > 0 && x[i] >= COST) || (diff < 0 && x[i] <= ZERO))
			diff = 0;
		
		diff = fabs(diff);
		sortbase[i] = diff;
		if (diff > maxdiff)
			maxdiff = diff; 	 
	}
	return maxdiff;
}

int select_workingset(int it)
{
	int i, j, sz = 0, sz0;

	if(Nfree == 0)	/* choose q/2 + and q/2 - */
	{
		for (i=0;i<svm_param.qpsize/2;i++)
			positive_max[i] = negative_max[i] = -DBL_MAX;

		for(i=0;i<N;i++)
		{
//			if(i==testv) continue;
			if(sample[i].a == 1)
			{
				double v = sortbase[i];
				if(v>positive_max[0])
				{
					for(j=1;j<svm_param.qpsize/2;j++)
					{
						if(v<=positive_max[j]) break;
						positive_max[j-1]=positive_max[j];
						positive_set[j-1]=positive_set[j];
					}
					positive_max[j-1]=v;
					positive_set[j-1]=i;
				}
			}
			else
			{
				double v = sortbase[i];
				if(v>negative_max[0])
				{
					for(j=1;j<svm_param.qpsize/2;j++)
					{
						if(v<=negative_max[j]) break;
						negative_max[j-1]=negative_max[j];
						negative_set[j-1]=negative_set[j];
					}
					negative_max[j-1]=v;
					negative_set[j-1]=i;
				}
			}
		}
		
		for(i=0;i<svm_param.qpsize/2;i++)
		{
			//if(positive_max[i]!=-DBL_MAX)
			if(positive_max[i]>0)
			{
				workingset[sz++] = positive_set[i];
				chosenflag[positive_set[i]]=it;
			}
			//if(negative_max[i]!=-DBL_MAX)
			if(negative_max[i]>0)
			{
				workingset[sz++] = negative_set[i];
				chosenflag[negative_set[i]]=it;
			}
		}
	}
	else
	{
		/* choose other free variables */
		for(i=0;i<svm_param.qpsize/2;i++)
			positive_max[i] = DBL_MAX;
		for(i=0;i<N;i++)
		{
			double v = sortbase[i];
//			if(i==testv) continue;
			if(x[i]>ZERO && x[i]<COST &&
			v<positive_max[0] )
			{
				for(j=1;j<svm_param.qpsize/2;j++)
				{
					if(v>=positive_max[j]) break;
					positive_max[j-1]=positive_max[j];
					positive_set[j-1]=positive_set[j];
				}
				positive_max[j-1]=v;
				positive_set[j-1]=i;
			}
		}
		
		for(i=0;i<svm_param.qpsize/2;i++)
		{
			//if(positive_max[i]!=DBL_MAX)
			if(positive_max[i]>ZERO)
			{
				workingset[sz++] = positive_set[i];
				chosenflag[positive_set[i]]=it;
			}
		}

		/* choose q/2 according to gradient */
		sz0 = sz;
		for (i=0;i<svm_param.qpsize - sz0;i++)
			positive_max[i] = -DBL_MAX;
		for(i=0;i<N;i++)
		{
			double v = sortbase[i];
//			if(i==testv) continue;
			if(v>positive_max[0] && chosenflag[i]<it)
			{
				for(j=1;j<svm_param.qpsize - sz0;j++)
				{
					if(v<=positive_max[j]) break;
					positive_max[j-1]=positive_max[j];
					positive_set[j-1]=positive_set[j];
				}
				positive_max[j-1]=v;
				positive_set[j-1]=i;
			}
		}
		
		for(i=0;i<svm_param.qpsize - sz0;i++)
		{
			//if(positive_max[i]!=-DBL_MAX)
			if(positive_max[i]>0)
			{
				workingset[sz++] = positive_set[i];
				chosenflag[positive_set[i]]=it;
			}
		}
	}

	/* return workingset size */
	return sz;
}
void construct_subQP(int qpn)
{
	int i, j, Bi, Bj;

	qp.C = svm_param.C;
	qp.n = qpn;
	for (i=0;i<qp.n;i++)
		qp.p[i] = Qx[workingset[i]] - 1;
	for (i=0;i<qp.n;i++)
	{
		Bi = workingset[i];
		qp.x[i] = x[Bi];
		qp.Q[i*qp.n+i] = QB[i][Bi];
		qp.p[i] -= qp.Q[i*qp.n+i]*x[Bi];
		for (j=i+1;j<qp.n;j++)