SHOGUN  v2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GNPPLib.cpp
Go to the documentation of this file.
1 /*-----------------------------------------------------------------------
2  *
3  * This program is free software; you can redistribute it and/or modify
4  * it under the terms of the GNU General Public License as published by
5  * the Free Software Foundation; either version 3 of the License, or
6  * (at your option) any later version.
7  *
8  * Library of solvers for Generalized Nearest Point Problem (GNPP).
9  *
10  * Written (W) 1999-2008 Vojtech Franc, xfrancv@cmp.felk.cvut.cz
11  * Copyright (C) 1999-2008 Center for Machine Perception, CTU FEL Prague
12  *
13 -------------------------------------------------------------------- */
14 
15 #include <math.h>
16 #include <limits.h>
17 #include <shogun/lib/common.h>
18 #include <shogun/io/SGIO.h>
19 #include <shogun/mathematics/Math.h>
20 
21 #include <shogun/classifier/svm/GNPPLib.h>
22 #include <shogun/kernel/Kernel.h>
23 
24 using namespace shogun;
25 
26 #define HISTORY_BUF 1000000
27 
28 #define MINUS_INF INT_MIN
29 #define PLUS_INF INT_MAX
30 
31 #define INDEX(ROW,COL,DIM) ((COL*DIM)+ROW)
32 
33 CGNPPLib::CGNPPLib()
34 {
35  SG_UNSTABLE("CGNPPLib::CGNPPLib()", "\n");
36 
37  kernel_columns = NULL;
38  cache_index = NULL;
39  first_kernel_inx = 0;
40  Cache_Size = 0;
41  m_num_data = 0;
42  m_reg_const = 0.0;
43  m_vector_y = NULL;
44  m_kernel = NULL;
45 }
46 
47 CGNPPLib::CGNPPLib(
48  float64_t* vector_y, CKernel* kernel, int32_t num_data, float64_t reg_const)
49 : CSGObject()
50 {
51  m_reg_const = reg_const;
52  m_num_data = num_data;
53  m_vector_y = vector_y;
54  m_kernel = kernel;
55 
56  Cache_Size = ((int64_t) kernel->get_cache_size())*1024*1024/(sizeof(float64_t)*num_data);
57  Cache_Size = CMath::min(Cache_Size, (int64_t) num_data);
58 
59  SG_INFO("using %d kernel cache lines\n", Cache_Size);
60  ASSERT(Cache_Size>=2);
61 
62  /* allocates memory for kernel cache */
63  kernel_columns = SG_MALLOC(float64_t*, Cache_Size);
64  cache_index = SG_MALLOC(float64_t, Cache_Size);
65 
66  for(int32_t i = 0; i < Cache_Size; i++ )
67  {
68  kernel_columns[i] = SG_MALLOC(float64_t, num_data);
69  cache_index[i] = -2;
70  }
71  first_kernel_inx = 0;
72 
73 }
74 
75 CGNPPLib::~CGNPPLib()
76 {
77  for(int32_t i = 0; i < Cache_Size; i++ )
78  SG_FREE(kernel_columns[i]);
79 
80  SG_FREE(cache_index);
81  SG_FREE(kernel_columns);
82 }
83 
84 /* --------------------------------------------------------------
85  QP solver based on Mitchell-Demyanov-Malozemov algorithm.
86 
87  Usage: exitflag = gnpp_mdm( diag_H, vector_c, vector_y,
88  dim, tmax, tolabs, tolrel, th, &alpha, &t, &aHa11, &aHa22, &History );
89 -------------------------------------------------------------- */
90 int8_t CGNPPLib::gnpp_mdm(float64_t *diag_H,
91  float64_t *vector_c,
92  float64_t *vector_y,
93  int32_t dim,
94  int32_t tmax,
95  float64_t tolabs,
96  float64_t tolrel,
97  float64_t th,
98  float64_t *alpha,
99  int32_t *ptr_t,
100  float64_t *ptr_aHa11,
101  float64_t *ptr_aHa22,
102  float64_t **ptr_History,
103  int32_t verb)
104 {
105  float64_t LB;
106  float64_t UB;
107  float64_t aHa11, aHa12, aHa22, ac1, ac2;
108  float64_t tmp;
109  float64_t Huu, Huv, Hvv;
110  float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
111  float64_t lambda;
112  float64_t delta1, delta2;
113  float64_t *History;
114  float64_t *Ha1;
115  float64_t *Ha2;
116  float64_t *tmp_ptr;
117  float64_t *col_u, *col_v;
118  float64_t *col_v1, *col_v2;
119  int64_t u1=0, u2=0;
120  int64_t v1, v2;
121  int64_t i;
122  int64_t t;
123  int64_t History_size;
124  int8_t exitflag;
125 
126  /* ------------------------------------------------------------ */
127  /* Initialization */
128  /* ------------------------------------------------------------ */
129 
130  Ha1 = SG_MALLOC(float64_t, dim);
131  if( Ha1 == NULL ) SG_ERROR("Not enough memory.\n");
132  Ha2 = SG_MALLOC(float64_t, dim);
133  if( Ha2 == NULL ) SG_ERROR("Not enough memory.\n");
134 
135  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
136  History = SG_MALLOC(float64_t, History_size*2);
137  if( History == NULL ) SG_ERROR("Not enough memory.\n");
138 
139  /* inx1 = firts of find( y ==1 ), inx2 = firts of find( y ==2 ) */
140  v1 = -1; v2 = -1; i = 0;
141  while( (v1 == -1 || v2 == -1) && i < dim ) {
142  if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
143  if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; }
144  i++;
145  }
146 
147  col_v1 = (float64_t*)get_col(v1,-1);
148  col_v2 = (float64_t*)get_col(v2,v1);
149 
150  aHa12 = col_v1[v2];
151  aHa11 = diag_H[v1];
152  aHa22 = diag_H[v2];
153  ac1 = vector_c[v1];
154  ac2 = vector_c[v2];
155 
156  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
157  for( i = 0; i < dim; i++ )
158  {
159  alpha[i] = 0;
160  Ha1[i] = col_v1[i];
161  Ha2[i] = col_v2[i];
162 
163  beta = Ha1[i] + Ha2[i] + vector_c[i];
164 
165  if( vector_y[i] == 1 && min_beta1 > beta ) {
166  u1 = i;
167  min_beta1 = beta;
168  }
169 
170  if( vector_y[i] == 2 && min_beta2 > beta ) {
171  u2 = i;
172  min_beta2 = beta;
173  }
174  }
175 
176  alpha[v1] = 1;
177  alpha[v2] = 1;
178 
179  UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
180  LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
181 
182  delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
183  delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
184 
185  t = 0;
186  History[INDEX(0,0,2)] = LB;
187  History[INDEX(1,0,2)] = UB;
188 
189  if( verb ) {
190  SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
191  UB, LB, UB-LB,(UB-LB)/UB);
192  }
193 
194  /* Stopping conditions */
195  if( UB-LB <= tolabs ) exitflag = 1;
196  else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
197  else if(LB > th) exitflag = 3;
198  else exitflag = -1;
199 
200  /* ------------------------------------------------------------ */
201  /* Main optimization loop */
202  /* ------------------------------------------------------------ */
203 
204  while( exitflag == -1 )
205  {
206  t++;
207 
208  if( delta1 > delta2 )
209  {
210  col_u = (float64_t*)get_col(u1,-1);
211  col_v = (float64_t*)get_col(v1,u1);
212 
213  Huu = diag_H[u1];
214  Hvv = diag_H[v1];
215  Huv = col_u[v1];
216 
217  lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
218  lambda = CMath::min(1.0,lambda);
219 
220  tmp = lambda*alpha[v1];
221 
222  aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
223  aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
224  ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
225 
226  alpha[u1] = alpha[u1] + tmp;
227  alpha[v1] = alpha[v1] - tmp;
228 
229  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
230  max_beta1 = MINUS_INF; max_beta2 = MINUS_INF;
231  for( i = 0; i < dim; i ++ )
232  {
233  Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
234 
235  beta = Ha1[i] + Ha2[i] + vector_c[i];
236  if( vector_y[i] == 1 )
237  {
238  if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
239  if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
240  }
241  else
242  {
243  if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
244  if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
245  }
246  }
247  }
248  else
249  {
250  col_u = (float64_t*)get_col(u2,-1);
251  col_v = (float64_t*)get_col(v2,u2);
252 
253  Huu = diag_H[u2];
254  Hvv = diag_H[v2];
255  Huv = col_u[v2];
256 
257  lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
258  lambda = CMath::min(1.0,lambda);
259 
260  tmp = lambda*alpha[v2];
261  aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
262  aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
263  ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
264 
265  alpha[u2] = alpha[u2] + tmp;
266  alpha[v2] = alpha[v2] - tmp;
267 
268  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
269  max_beta1 = MINUS_INF; max_beta2 = MINUS_INF;
270  for(i = 0; i < dim; i++ )
271  {
272  Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
273 
274  beta = Ha1[i] + Ha2[i] + vector_c[i];
275 
276  if( vector_y[i] == 1 )
277  {
278  if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
279  if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
280  }
281  else
282  {
283  if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
284  if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
285  }
286  }
287  }
288 
289  UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
290  LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
291 
292  delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
293  delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
294 
295  /* Stopping conditions */
296  if( UB-LB <= tolabs ) exitflag = 1;
297  else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
298  else if(LB > th) exitflag = 3;
299  else if(t >= tmax) exitflag = 0;
300 
301  if( verb && (t % verb) == 0) {
302  SG_PRINT("%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
303  t, UB, LB, UB-LB,(UB-LB)/UB);
304  }
305 
306  /* Store selected values */
307  if( t < History_size ) {
308  History[INDEX(0,t,2)] = LB;
309  History[INDEX(1,t,2)] = UB;
310  }
311  else {
312  tmp_ptr = SG_MALLOC(float64_t, (History_size+HISTORY_BUF)*2);
313  if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.\n");
314  for( i = 0; i < History_size; i++ ) {
315  tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
316  tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
317  }
318  tmp_ptr[INDEX(0,t,2)] = LB;
319  tmp_ptr[INDEX(1,t,2)] = UB;
320 
321  History_size += HISTORY_BUF;
322  SG_FREE(History);
323  History = tmp_ptr;
324  }
325  }
326 
327  /* print info about last iteration*/
328  if(verb && (t % verb) ) {
329  SG_PRINT("Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
330  UB, LB, UB-LB,(UB-LB)/UB);
331  }
332 
333  /*------------------------------------------------------- */
334  /* Set outputs */
335  /*------------------------------------------------------- */
336  (*ptr_t) = t;
337  (*ptr_aHa11) = aHa11;
338  (*ptr_aHa22) = aHa22;
339  (*ptr_History) = History;
340 
341  /* Free memory */
342  SG_FREE(Ha1);
343  SG_FREE(Ha2);
344 
345  return( exitflag );
346 }
347 
348 
349 /* --------------------------------------------------------------
350  QP solver based on Improved MDM algorithm (u fixed v optimized)
351 
352  Usage: exitflag = gnpp_imdm( diag_H, vector_c, vector_y,
353  dim, tmax, tolabs, tolrel, th, &alpha, &t, &aHa11, &aHa22, &History );
354 -------------------------------------------------------------- */
355 int8_t CGNPPLib::gnpp_imdm(float64_t *diag_H,
356  float64_t *vector_c,
357  float64_t *vector_y,
358  int32_t dim,
359  int32_t tmax,
360  float64_t tolabs,
361  float64_t tolrel,
362  float64_t th,
363  float64_t *alpha,
364  int32_t *ptr_t,
365  float64_t *ptr_aHa11,
366  float64_t *ptr_aHa22,
367  float64_t **ptr_History,
368  int32_t verb)
369 {
370  float64_t LB;
371  float64_t UB;
372  float64_t aHa11, aHa12, aHa22, ac1, ac2;
373  float64_t tmp;
374  float64_t Huu, Huv, Hvv;
375  float64_t min_beta1, max_beta1, min_beta2, max_beta2, beta;
376  float64_t lambda;
377  float64_t delta1, delta2;
378  float64_t improv, max_improv;
379  float64_t *History;
380  float64_t *Ha1;
381  float64_t *Ha2;
382  float64_t *tmp_ptr;
383  float64_t *col_u, *col_v;
384  float64_t *col_v1, *col_v2;
385  int64_t u1=0, u2=0;
386  int64_t v1, v2;
387  int64_t i;
388  int64_t t;
389  int64_t History_size;
390  int8_t exitflag;
391  int8_t which_case;
392 
393  /* ------------------------------------------------------------ */
394  /* Initialization */
395  /* ------------------------------------------------------------ */
396 
397  Ha1 = SG_MALLOC(float64_t, dim);
398  if( Ha1 == NULL ) SG_ERROR("Not enough memory.\n");
399  Ha2 = SG_MALLOC(float64_t, dim);
400  if( Ha2 == NULL ) SG_ERROR("Not enough memory.\n");
401 
402  History_size = (tmax < HISTORY_BUF ) ? tmax+1 : HISTORY_BUF;
403  History = SG_MALLOC(float64_t, History_size*2);
404  if( History == NULL ) SG_ERROR("Not enough memory.\n");
405 
406  /* inx1 = firts of find( y ==1 ), inx2 = firts of find( y ==2 ) */
407  v1 = -1; v2 = -1; i = 0;
408  while( (v1 == -1 || v2 == -1) && i < dim ) {
409  if( v1 == -1 && vector_y[i] == 1 ) { v1 = i; }
410  if( v2 == -1 && vector_y[i] == 2 ) { v2 = i; }
411  i++;
412  }
413 
414  col_v1 = (float64_t*)get_col(v1,-1);
415  col_v2 = (float64_t*)get_col(v2,v1);
416 
417  aHa12 = col_v1[v2];
418  aHa11 = diag_H[v1];
419  aHa22 = diag_H[v2];
420  ac1 = vector_c[v1];
421  ac2 = vector_c[v2];
422 
423  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
424  for( i = 0; i < dim; i++ )
425  {
426  alpha[i] = 0;
427  Ha1[i] = col_v1[i];
428  Ha2[i] = col_v2[i];
429 
430  beta = Ha1[i] + Ha2[i] + vector_c[i];
431 
432  if( vector_y[i] == 1 && min_beta1 > beta ) {
433  u1 = i;
434  min_beta1 = beta;
435  }
436 
437  if( vector_y[i] == 2 && min_beta2 > beta ) {
438  u2 = i;
439  min_beta2 = beta;
440  }
441  }
442 
443  alpha[v1] = 1;
444  alpha[v2] = 1;
445 
446  UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
447  LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
448 
449  delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
450  delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
451 
452  t = 0;
453  History[INDEX(0,0,2)] = LB;
454  History[INDEX(1,0,2)] = UB;
455 
456  if( verb ) {
457  SG_PRINT("Init: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
458  UB, LB, UB-LB,(UB-LB)/UB);
459  }
460 
461  if( delta1 > delta2 )
462  {
463  which_case = 1;
464  col_u = (float64_t*)get_col(u1,v1);
465  col_v = col_v1;
466  }
467  else
468  {
469  which_case = 2;
470  col_u = (float64_t*)get_col(u2,v2);
471  col_v = col_v2;
472  }
473 
474  /* Stopping conditions */
475  if( UB-LB <= tolabs ) exitflag = 1;
476  else if(UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
477  else if(LB > th) exitflag = 3;
478  else exitflag = -1;
479 
480  /* ------------------------------------------------------------ */
481  /* Main optimization loop */
482  /* ------------------------------------------------------------ */
483 
484  while( exitflag == -1 )
485  {
486  t++;
487 
488  if( which_case == 1 )
489  {
490  Huu = diag_H[u1];
491  Hvv = diag_H[v1];
492  Huv = col_u[v1];
493 
494  lambda = delta1/(alpha[v1]*(Huu - 2*Huv + Hvv ));
495  lambda = CMath::min(1.0,lambda);
496 
497  tmp = lambda*alpha[v1];
498 
499  aHa11 = aHa11 + 2*tmp*(Ha1[u1]-Ha1[v1])+tmp*tmp*( Huu - 2*Huv + Hvv );
500  aHa12 = aHa12 + tmp*(Ha2[u1]-Ha2[v1]);
501  ac1 = ac1 + tmp*(vector_c[u1]-vector_c[v1]);
502 
503  alpha[u1] = alpha[u1] + tmp;
504  alpha[v1] = alpha[v1] - tmp;
505 
506  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
507  max_beta1 = MINUS_INF; max_beta2 = MINUS_INF;
508  for( i = 0; i < dim; i ++ )
509  {
510  Ha1[i] = Ha1[i] + tmp*(col_u[i] - col_v[i]);
511 
512  beta = Ha1[i] + Ha2[i] + vector_c[i];
513  if( vector_y[i] == 1 )
514  {
515  if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
516  if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
517  }
518  else
519  {
520  if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
521  if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
522  }
523  }
524  }
525  else
526  {
527  Huu = diag_H[u2];
528  Hvv = diag_H[v2];
529  Huv = col_u[v2];
530 
531  lambda = delta2/(alpha[v2]*( Huu - 2*Huv + Hvv ));
532  lambda = CMath::min(1.0,lambda);
533 
534  tmp = lambda*alpha[v2];
535  aHa22 = aHa22 + 2*tmp*( Ha2[u2]-Ha2[v2]) + tmp*tmp*( Huu - 2*Huv + Hvv);
536  aHa12 = aHa12 + tmp*(Ha1[u2]-Ha1[v2]);
537  ac2 = ac2 + tmp*( vector_c[u2]-vector_c[v2] );
538 
539  alpha[u2] = alpha[u2] + tmp;
540  alpha[v2] = alpha[v2] - tmp;
541 
542  min_beta1 = PLUS_INF; min_beta2 = PLUS_INF;
543  max_beta1 = MINUS_INF; max_beta2 = MINUS_INF;
544  for(i = 0; i < dim; i++ )
545  {
546  Ha2[i] = Ha2[i] + tmp*( col_u[i] - col_v[i] );
547 
548  beta = Ha1[i] + Ha2[i] + vector_c[i];
549 
550  if( vector_y[i] == 1 )
551  {
552  if( min_beta1 > beta ) { u1 = i; min_beta1 = beta; }
553  if( max_beta1 < beta && alpha[i] > 0 ) { v1 = i; max_beta1 = beta; }
554  }
555  else
556  {
557  if( min_beta2 > beta ) { u2 = i; min_beta2 = beta; }
558  if( max_beta2 < beta && alpha[i] > 0) { v2 = i; max_beta2 = beta; }
559  }
560  }
561  }
562 
563  UB = 0.5*(aHa11 + 2*aHa12 + aHa22) + ac1 + ac2;
564  LB = min_beta1 + min_beta2 - 0.5*(aHa11 + 2*aHa12 + aHa22);
565 
566  delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
567  delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
568 
569  if( delta1 > delta2 )
570  {
571  col_u = (float64_t*)get_col(u1,-1);
572 
573  /* search for optimal v while u is fixed */
574  for( max_improv = MINUS_INF, i = 0; i < dim; i++ ) {
575 
576  if( vector_y[i] == 1 && alpha[i] != 0 ) {
577 
578  beta = Ha1[i] + Ha2[i] + vector_c[i];
579 
580  if( beta >= min_beta1 ) {
581 
582  tmp = diag_H[u1] - 2*col_u[i] + diag_H[i];
583  if( tmp != 0 ) {
584  improv = (0.5*(beta-min_beta1)*(beta-min_beta1))/tmp;
585 
586  if( improv > max_improv ) {
587  max_improv = improv;
588  v1 = i;
589  }
590  }
591  }
592  }
593  }
594  col_v = (float64_t*)get_col(v1,u1);
595  delta1 = Ha1[v1] + Ha2[v1] + vector_c[v1] - min_beta1;
596  which_case = 1;
597 
598  }
599  else
600  {
601  col_u = (float64_t*)get_col(u2,-1);
602 
603  /* search for optimal v while u is fixed */
604  for( max_improv = MINUS_INF, i = 0; i < dim; i++ ) {
605 
606  if( vector_y[i] == 2 && alpha[i] != 0 ) {
607 
608  beta = Ha1[i] + Ha2[i] + vector_c[i];
609 
610  if( beta >= min_beta2 ) {
611 
612  tmp = diag_H[u2] - 2*col_u[i] + diag_H[i];
613  if( tmp != 0 ) {
614  improv = (0.5*(beta-min_beta2)*(beta-min_beta2))/tmp;
615 
616  if( improv > max_improv ) {
617  max_improv = improv;
618  v2 = i;
619  }
620  }
621  }
622  }
623  }
624 
625  col_v = (float64_t*)get_col(v2,u2);
626  delta2 = Ha1[v2] + Ha2[v2] + vector_c[v2] - min_beta2;
627  which_case = 2;
628  }
629 
630 
631  /* Stopping conditions */
632  if( UB-LB <= tolabs ) exitflag = 1;
633  else if( UB-LB <= CMath::abs(UB)*tolrel ) exitflag = 2;
634  else if(LB > th) exitflag = 3;
635  else if(t >= tmax) exitflag = 0;
636 
637  if( verb && (t % verb) == 0) {
638  SG_PRINT("%d: UB=%f,LB=%f,UB-LB=%f,(UB-LB)/|UB|=%f\n",
639  t, UB, LB, UB-LB,(UB-LB)/UB);
640  }
641 
642  /* Store selected values */
643  if( t < History_size ) {
644  History[INDEX(0,t,2)] = LB;
645  History[INDEX(1,t,2)] = UB;
646  }
647  else {
648  tmp_ptr = SG_MALLOC(float64_t, (History_size+HISTORY_BUF)*2);
649  if( tmp_ptr == NULL ) SG_ERROR("Not enough memory.\n");
650  for( i = 0; i < History_size; i++ ) {
651  tmp_ptr[INDEX(0,i,2)] = History[INDEX(0,i,2)];
652  tmp_ptr[INDEX(1,i,2)] = History[INDEX(1,i,2)];
653  }
654  tmp_ptr[INDEX(0,t,2)] = LB;
655  tmp_ptr[INDEX(1,t,2)] = UB;
656 
657  History_size += HISTORY_BUF;
658  SG_FREE(History);
659  History = tmp_ptr;
660  }
661  }
662 
663  /* print info about last iteration*/
664  if(verb && (t % verb) ) {
665  SG_PRINT("Exit: UB=%f, LB=%f, UB-LB=%f, (UB-LB)/|UB|=%f \n",
666  UB, LB, UB-LB,(UB-LB)/UB);
667  }
668 
669  /*------------------------------------------------------- */
670  /* Set outputs */
671  /*------------------------------------------------------- */
672  (*ptr_t) = t;
673  (*ptr_aHa11) = aHa11;
674  (*ptr_aHa22) = aHa22;
675  (*ptr_History) = History;
676 
677  /* Free memory */
678  SG_FREE(Ha1);
679  SG_FREE(Ha2);
680 
681  return( exitflag );
682 }
683 
684 
685 float64_t* CGNPPLib::get_col(int64_t a, int64_t b)
686 {
687  float64_t *col_ptr;
688  float64_t y;
689  int64_t i;
690  int64_t inx;
691 
692  inx = -1;
693  for( i=0; i < Cache_Size; i++ ) {
694  if( cache_index[i] == a ) { inx = i; break; }
695  }
696 
697  if( inx != -1 ) {
698  col_ptr = kernel_columns[inx];
699  return( col_ptr );
700  }
701 
702  col_ptr = kernel_columns[first_kernel_inx];
703  cache_index[first_kernel_inx] = a;
704 
705  first_kernel_inx++;
706  if( first_kernel_inx >= Cache_Size ) first_kernel_inx = 0;
707 
708  y = m_vector_y[a];
709  for( i=0; i < m_num_data; i++ ) {
710  if( m_vector_y[i] == y )
711  {
712  col_ptr[i] = 2*m_kernel->kernel(i,a);
713  }
714  else
715  {
716  col_ptr[i] = -2*m_kernel->kernel(i,a);
717  }
718  }
719 
720  col_ptr[a] = col_ptr[a] + m_reg_const;
721 
722  return( col_ptr );
723 }
SHOGUN Machine Learning Toolbox - Documentation