{

dst = new T[n];

memcpy((void *)dst,(void *)src,sizeof(T)*n);

{

double tmp = base, ret = 1.0;

for(int t=times; t>0; t/=2)

{

if(t%2==1) ret*=tmp;

tmp = tmp * tmp;

}

return ret;

{

/*

Rprintf(s);

{

char buf[BUFSIZ];

va_list ap;

va_start(ap,fmt);

vsprintf(buf,fmt,ap);

va_end(ap);

(*svm_print_string)(buf);

{

Cache(int l,long int size);

~Cache();

// request data [0,len)

// return some position p where [p,len) need to be filled

// (p >= len if nothing needs to be filled)

int get_data(const int index, Qfloat **data, int len);

void swap_index(int i, int j);

int l;

long int size;

struct head_t

{

head_t *prev, *next; // a circular list

Qfloat *data;

int len; // data[0,len) is cached in this entry

};

head_t *head;

head_t lru_head;

void lru_delete(head_t *h);

void lru_insert(head_t *h);

{

head = (head_t *)calloc(l,sizeof(head_t)); // initialized to 0

size /= sizeof(Qfloat);

size -= l * sizeof(head_t) / sizeof(Qfloat);

size = max(size, 2 * (long int) l); // cache must be large enough for two columns

lru_head.next = lru_head.prev = &lru_head;

{

// delete from current location

h->prev->next = h->next;

h->next->prev = h->prev;

{

// insert to last position

h->next = &lru_head;

h->prev = lru_head.prev;

h->prev->next = h;

h->next->prev = h;

{

head_t *h = &head[index];

if(h->len) lru_delete(h);

int more = len - h->len;

if(more > 0)

{

// free old space

while(size < more)

{

head_t *old = lru_head.next;

lru_delete(old);

free(old->data);

size += old->len;

old->data = 0;

old->len = 0;

}

// allocate new space

h->data = (Qfloat *)realloc(h->data,sizeof(Qfloat)*len);

size -= more;

swap(h->len,len);

}

lru_insert(h);

*data = h->data;

return len;

{

if(i==j) return;

if(head[i].len) lru_delete(&head[i]);

if(head[j].len) lru_delete(&head[j]);

swap(head[i].data,head[j].data);

swap(head[i].len,head[j].len);

if(head[i].len) lru_insert(&head[i]);

if(head[j].len) lru_insert(&head[j]);

if(i>j) swap(i,j);

for(head_t *h = lru_head.next; h!=&lru_head; h=h->next)

{

if(h->len > i)

{

if(h->len > j)

swap(h->data[i],h->data[j]);

else

{

// give up

lru_delete(h);

free(h->data);

size += h->len;

h->data = 0;

h->len = 0;

}

}

}

virtual Qfloat *get_Q(int column, int len) const = 0;

virtual double *get_QD() const = 0;

virtual void swap_index(int i, int j) const = 0;

virtual ~QMatrix() {}

#ifdef _DENSE_REP

Kernel(int l, svm_node * x, const svm_parameter& param);

#else

Kernel(int l, svm_node * const * x, const svm_parameter& param);

#endif

virtual ~Kernel();

static double k_function(const svm_node *x, const svm_node *y,

const svm_parameter& param);

virtual Qfloat *get_Q(int column, int len) const = 0;

virtual double *get_QD() const = 0;

virtual void swap_index(int i, int j) const // no so const...

{

swap(x[i],x[j]);

if(x_square) swap(x_square[i],x_square[j]);

}

double (Kernel::*kernel_function)(int i, int j) const;

#ifdef _DENSE_REP

svm_node *x;

#else

const svm_node **x;

#endif

double *x_square;

// svm_parameter

const int kernel_type;

const int degree;

const double gamma;

const double coef0;

static double dot(const svm_node *px, const svm_node *py);

#ifdef _DENSE_REP

static double dot(const svm_node &px, const svm_node &py);

#endif

double kernel_linear(int i, int j) const

{

return dot(x[i],x[j]);

}

double kernel_poly(int i, int j) const

{

return powi(gamma*dot(x[i],x[j])+coef0,degree);

}

double kernel_rbf(int i, int j) const

{

return exp(-gamma*(x_square[i]+x_square[j]-2*dot(x[i],x[j])));

}

double kernel_sigmoid(int i, int j) const

{

return tanh(gamma*dot(x[i],x[j])+coef0);

}

double kernel_precomputed(int i, int j) const

{

#ifdef _DENSE_REP

return (x+i)->values[(int)((x+j)->values[0])];

#else

return x[i][(int)(x[j][0].value)].value;

#endif

}

#endif

:kernel_type(param.kernel_type), degree(param.degree),

gamma(param.gamma), coef0(param.coef0)

{

switch(kernel_type)

{

case LINEAR:

kernel_function = &Kernel::kernel_linear;

break;

case POLY:

kernel_function = &Kernel::kernel_poly;

break;

case RBF:

kernel_function = &Kernel::kernel_rbf;

break;

case SIGMOID:

kernel_function = &Kernel::kernel_sigmoid;

break;

case PRECOMPUTED:

kernel_function = &Kernel::kernel_precomputed;

break;

}

clone(x,x_,l);

if(kernel_type == RBF)

{

x_square = new double[l];

for(int i=0;i<l;i++)

x_square[i] = dot(x[i],x[i]);

}

else

x_square = 0;

{

double sum = 0;

int dim = min(px->dim, py->dim);

for (int i = 0; i < dim; i++)

sum += (px->values)[i] * (py->values)[i];

return sum;

{

double sum = 0;

int dim = min(px.dim, py.dim);

for (int i = 0; i < dim; i++)

sum += px.values[i] * py.values[i];

return sum;

{

double sum = 0;

while(px->index != -1 && py->index != -1)

{

if(px->index == py->index)

{

sum += px->value * py->value;

++px;

++py;

}

else

{

if(px->index > py->index)

++py;

else

++px;

}

}

return sum;

const svm_parameter& param)

{

switch(param.kernel_type)

{

case LINEAR:

return dot(x,y);

case POLY:

return powi(param.gamma*dot(x,y)+param.coef0,param.degree);

case RBF:

{

double sum = 0;

#ifdef _DENSE_REP

int dim = min(x->dim, y->dim), i;

for (i = 0; i < dim; i++)

{

double d = x->values[i] - y->values[i];

sum += d*d;

}

for (; i < x->dim; i++)

sum += x->values[i] * x->values[i];

for (; i < y->dim; i++)

sum += y->values[i] * y->values[i];

#else

while(x->index != -1 && y->index !=-1)

{

if(x->index == y->index)

{

double d = x->value - y->value;

sum += d*d;

++x;

++y;

}

else

{

if(x->index > y->index)

{

sum += y->value * y->value;

++y;

}

else

{

sum += x->value * x->value;

++x;

}

}

}

while(x->index != -1)

{

sum += x->value * x->value;

++x;

}

while(y->index != -1)

{

sum += y->value * y->value;

++y;

}

#endif

return exp(-param.gamma*sum);

}

case SIGMOID:

return tanh(param.gamma*dot(x,y)+param.coef0);

case PRECOMPUTED: //x: test (validation), y: SV

#ifdef _DENSE_REP

return x->values[(int)(y->values[0])];

#else

return x[(int)(y->value)].value;

#endif

default:

return 0; // Unreachable

}

Solver() {};

virtual ~Solver() {};

struct SolutionInfo {

double obj;

double rho;

double upper_bound_p;

double upper_bound_n;

double r; // for Solver_NU

};

void Solve(int l, const QMatrix& Q, const double *p_, const schar *y_,

double *alpha_, double Cp, double Cn, double eps,

SolutionInfo* si, int shrinking);

void Solve_b(int l, const QMatrix& Q, const double *p_, const schar *y_,

double *alpha_, double Cp, double Cn, double *upbound_, double eps,

SolutionInfo* si, int shrinking);

int active_size;

schar *y;

double *G; // gradient of objective function

enum { LOWER_BOUND, UPPER_BOUND, FREE };

char *alpha_status; // LOWER_BOUND, UPPER_BOUND, FREE

double *alpha;

const QMatrix *Q;

const double *QD;

double eps;

double Cp,Cn;

double *p;

int *active_set;

double *G_bar; // gradient, if we treat free variables as 0

int l;

bool unshrink; // XXX

bool ubound;

double *upbound;

double get_C(int i)

{

if (ubound) {

return upbound[i];

} else {

return (y[i] > 0)? Cp : Cn;

}

}

void update_alpha_status(int i)

{

if(alpha[i] >= get_C(i))

alpha_status[i] = UPPER_BOUND;

else if(alpha[i] <= 0)

alpha_status[i] = LOWER_BOUND;

else alpha_status[i] = FREE;

}

bool is_upper_bound(int i) { return alpha_status[i] == UPPER_BOUND; }

bool is_lower_bound(int i) { return alpha_status[i] == LOWER_BOUND; }

bool is_free(int i) { return alpha_status[i] == FREE; }

void swap_index(int i, int j);

void reconstruct_gradient();

virtual int select_working_set(int &i, int &j);

virtual double calculate_rho();

virtual void do_shrinking();

bool be_shrunk(int i, double Gmax1, double Gmax2);

{

Q->swap_index(i,j);

swap(y[i],y[j]);

swap(G[i],G[j]);

swap(alpha_status[i],alpha_status[j]);

swap(alpha[i],alpha[j]);

swap(p[i],p[j]);

swap(active_set[i],active_set[j]);

swap(G_bar[i],G_bar[j]);

if (ubound) {

swap(upbound[i],upbound[j]);

}

{

// reconstruct inactive elements of G from G_bar and free variables

if(active_size == l) return;

int i,j;

int nr_free = 0;

for(j=active_size;j<l;j++)

G[j] = G_bar[j] + p[j];

for(j=0;j<active_size;j++)

if(is_free(j))

nr_free++;

if(2*nr_free < active_size)

info("\nWARNING: using -h 0 may be faster\n");

if (nr_free*l > 2*active_size*(l-active_size))

{

for(i=active_size;i<l;i++)

{

const Qfloat *Q_i = Q->get_Q(i,active_size);

for(j=0;j<active_size;j++)

if(is_free(j))

G[i] += alpha[j] * Q_i[j];

}

}

else

{

for(i=0;i<active_size;i++)

if(is_free(i))

{

const Qfloat *Q_i = Q->get_Q(i,l);

double alpha_i = alpha[i];

for(j=active_size;j<l;j++)

G[j] += alpha_i * Q_i[j];

}

}

double *alpha_, double Cp, double Cn, double eps,

SolutionInfo* si, int shrinking)

{

this->l = l;

this->Q = &Q;

QD=Q.get_QD();

clone(p, p_,l);

clone(y, y_,l);

clone(alpha,alpha_,l);

this->Cp = Cp;

this->Cn = Cn;

this->eps = eps;

unshrink = false;

this->ubound = 0;

// initialize alpha_status

{

alpha_status = new char[l];

for(int i=0;i<l;i++)

update_alpha_status(i);

}

// initialize active set (for shrinking)

{

active_set = new int[l];

for(int i=0;i<l;i++)

active_set[i] = i;

active_size = l;

}

// initialize gradient

{

G = new double[l];

G_bar = new double[l];

int i;

for(i=0;i<l;i++)

{

G[i] = p[i];

G_bar[i] = 0;

}

for(i=0;i<l;i++)

if(!is_lower_bound(i))

{

const Qfloat *Q_i = Q.get_Q(i,l);

double alpha_i = alpha[i];

int j;

for(j=0;j<l;j++)

G[j] += alpha_i*Q_i[j];

if(is_upper_bound(i))

for(j=0;j<l;j++)

G_bar[j] += get_C(i) * Q_i[j];

}

}

// optimization step

int iter = 0;

int max_iter = max(10000000, l>INT_MAX/100 ? INT_MAX : 100*l);

int counter = min(l,1000)+1;

while(iter < max_iter)

{

// show progress and do shrinking

if(--counter == 0)

{

counter = min(l,1000);

if(shrinking) do_shrinking();

info(".");

}

int i,j;

if(select_working_set(i,j)!=0)

{

// reconstruct the whole gradient

reconstruct_gradient();

// reset active set size and check

active_size = l;

info("*");

if(select_working_set(i,j)!=0)

break;

else

counter = 1; // do shrinking next iteration

}

++iter;

// update alpha[i] and alpha[j], handle bounds carefully

const Qfloat *Q_i = Q.get_Q(i,active_size);

const Qfloat *Q_j = Q.get_Q(j,active_size);

double C_i = get_C(i);

double C_j = get_C(j);

double old_alpha_i = alpha[i];

double old_alpha_j = alpha[j];

if(y[i]!=y[j])

{

double quad_coef = QD[i]+QD[j]+2*Q_i[j];

if (quad_coef <= 0)

quad_coef = TAU;

double delta = (-G[i]-G[j])/quad_coef;

double diff = alpha[i] - alpha[j];

alpha[i] += delta;

alpha[j] += delta;

if(diff > 0)

{

if(alpha[j] < 0)

{

alpha[j] = 0;

alpha[i] = diff;

}

}

else

{

if(alpha[i] < 0)

{

alpha[i] = 0;

alpha[j] = -diff;

}

}

if(diff > C_i - C_j)

{

if(alpha[i] > C_i)

{

alpha[i] = C_i;

alpha[j] = C_i - diff;

}

}

else

{

if(alpha[j] > C_j)

{

alpha[j] = C_j;

alpha[i] = C_j + diff;

}

}

}

else

{

double quad_coef = QD[i]+QD[j]-2*Q_i[j];

if (quad_coef <= 0)

quad_coef = TAU;

double delta = (G[i]-G[j])/quad_coef;

double sum = alpha[i] + alpha[j];

alpha[i] -= delta;

alpha[j] += delta;

if(sum > C_i)

{

if(alpha[i] > C_i)

{

alpha[i] = C_i;

alpha[j] = sum - C_i;

}

}

else

{

if(alpha[j] < 0)

{

alpha[j] = 0;

alpha[i] = sum;

}

}

if(sum > C_j)

{

if(alpha[j] > C_j)

{

alpha[j] = C_j;

alpha[i] = sum - C_j;

}

}

else

{

if(alpha[i] < 0)

{

alpha[i] = 0;

alpha[j] = sum;

}

}

}

// update G

double delta_alpha_i = alpha[i] - old_alpha_i;

double delta_alpha_j = alpha[j] - old_alpha_j;

for(int k=0;k<active_size;k++)

{

G[k] += Q_i[k]*delta_alpha_i + Q_j[k]*delta_alpha_j;

}

// update alpha_status and G_bar

{

bool ui = is_upper_bound(i);

bool uj = is_upper_bound(j);

update_alpha_status(i);

update_alpha_status(j);

int k;

if(ui != is_upper_bound(i))

{

Q_i = Q.get_Q(i,l);

if(ui)

for(k=0;k<l;k++)

G_bar[k] -= C_i * Q_i[k];

else

for(k=0;k<l;k++)

G_bar[k] += C_i * Q_i[k];

}

if(uj != is_upper_bound(j))

{

Q_j = Q.get_Q(j,l);

if(uj)

for(k=0;k<l;k++)

G_bar[k] -= C_j * Q_j[k];

else

for(k=0;k<l;k++)

G_bar[k] += C_j * Q_j[k];

}

}

}

if(iter >= max_iter)

{

if(active_size < l)

{

// reconstruct the whole gradient to calculate objective value

reconstruct_gradient();

active_size = l;

info("*");

}

REprintf("\nWARNING: reaching max number of iterations\n");

}

// calculate rho

si->rho = calculate_rho();

// calculate objective value

{

double v = 0;

int i;

for(i=0;i<l;i++)

v += alpha[i] * (G[i] + p[i]);

si->obj = v/2;

}

// put back the solution

{

for(int i=0;i<l;i++)

alpha_[active_set[i]] = alpha[i];

}

// juggle everything back

/*{

si->upper_bound_p = Cp;

si->upper_bound_n = Cn;

info("\noptimization finished, #iter = %d\n",iter);

delete[] p;

delete[] y;

delete[] alpha;

delete[] alpha_status;

delete[] active_set;

delete[] G;

delete[] G_bar;

double *alpha_, double Cp, double Cn, double *upbound_, double eps,

SolutionInfo* si, int shrinking)

{

this->l = l;

this->Q = &Q;

QD=Q.get_QD();

clone(p, p_,l);

clone(y, y_,l);

clone(alpha,alpha_,l);

clone(upbound,upbound_,l);

this->Cp = Cp;

this->Cn = Cn;

this->eps = eps;

unshrink = false;

this->ubound = 1;

// initialize alpha_status

{

alpha_status = new char[l];

for(int i=0;i<l;i++)

update_alpha_status(i);

}

// initialize active set (for shrinking)

{

active_set = new int[l];

for(int i=0;i<l;i++)

active_set[i] = i;

active_size = l;

}

// initialize gradient

{

G = new double[l];

G_bar = new double[l];

int i;

for(i=0;i<l;i++)

{

G[i] = p[i];

G_bar[i] = 0;

}

for(i=0;i<l;i++)

if(!is_lower_bound(i))

{

const Qfloat *Q_i = Q.get_Q(i,l);

double alpha_i = alpha[i];

int j;

for(j=0;j<l;j++)

G[j] += alpha_i*Q_i[j];

if(is_upper_bound(i))

for(j=0;j<l;j++)

G_bar[j] += get_C(i) * Q_i[j];

}

}

// optimization step

int iter = 0;

int counter = min(l,1000)+1;

while(1)

{

// show progress and do shrinking

if(--counter == 0)

{

counter = min(l,1000);

if(shrinking) do_shrinking();

info(".");

}

int i,j;

if(select_working_set(i,j)!=0)

{

// reconstruct the whole gradient

reconstruct_gradient();

// reset active set size and check

active_size = l;

info("*");

if(select_working_set(i,j)!=0)

break;

else

counter = 1; // do shrinking next iteration

}

++iter;

// update alpha[i] and alpha[j], handle bounds carefully

const Qfloat *Q_i = Q.get_Q(i,active_size);

const Qfloat *Q_j = Q.get_Q(j,active_size);

double C_i = get_C(i);

double C_j = get_C(j);

double old_alpha_i = alpha[i];

double old_alpha_j = alpha[j];

if(y[i]!=y[j])

{

double quad_coef = Q_i[i]+Q_j[j]+2*Q_i[j];

if (quad_coef <= 0)

quad_coef = TAU;

double delta = (-G[i]-G[j])/quad_coef;

double diff = alpha[i] - alpha[j];

alpha[i] += delta;

alpha[j] += delta;

if(diff > 0)

{

if(alpha[j] < 0)

{

alpha[j] = 0;

alpha[i] = diff;

}

}

else

{

if(alpha[i] < 0)

{

alpha[i] = 0;

alpha[j] = -diff;

}

}

if(diff > C_i - C_j)

{

if(alpha[i] > C_i)

{

alpha[i] = C_i;

alpha[j] = C_i - diff;

}

}

else

{