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

SHOGUN Machine Learning Toolbox - Documentation