SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Cplex.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 1999-2009 Soeren Sonnenburg
8  * Copyright (C) 1999-2009 Fraunhofer Institute FIRST and Max-Planck-Society
9  */
10 
11 #include <shogun/lib/config.h>
12 
13 #ifdef USE_CPLEX
14 #include <unistd.h>
15 
17 #include <shogun/io/SGIO.h>
18 #include <shogun/lib/Signal.h>
20 
21 using namespace shogun;
22 
24 : CSGObject(), env(NULL), lp(NULL), lp_initialized(false)
25 {
26 }
27 
29 {
30  cleanup();
31 }
32 
33 bool CCplex::init(E_PROB_TYPE typ, int32_t timeout)
34 {
35  problem_type=typ;
36 
37  while (env==NULL)
38  {
39  int status = 1; /* for calling external lib */
40  env = CPXopenCPLEX (&status);
41 
42  if ( env == NULL )
43  {
44  char errmsg[1024];
45  SG_WARNING("Could not open CPLEX environment.\n")
46  CPXgeterrorstring (env, status, errmsg);
47  SG_WARNING("%s", errmsg)
48  SG_WARNING("retrying in %d seconds\n", timeout)
49  sleep(timeout);
50  }
51  else
52  {
53  /* Turn on output to the screen */
54 
55  status = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_OFF);
56  if (status)
57  SG_ERROR("Failure to turn off screen indicator, error %d.\n", status)
58 
59  {
60  status = CPXsetintparam (env, CPX_PARAM_DATACHECK, CPX_ON);
61  if (status)
62  SG_ERROR("Failure to turn on data checking, error %d.\n", status)
63  else
64  {
65  lp = CPXcreateprob (env, &status, "shogun");
66 
67  if ( lp == NULL )
68  SG_ERROR("Failed to create optimization problem.\n")
69  else
70  CPXchgobjsen (env, lp, CPX_MIN); /* Problem is minimization */
71 
72  if (problem_type==E_QP)
73  status = CPXsetintparam (env, CPX_PARAM_QPMETHOD, 0);
74  else if (problem_type == E_LINEAR)
75  status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 0);
76  if (status)
77  SG_ERROR("Failure to select dual lp/qp optimization, error %d.\n", status)
78 
79  }
80  }
81  }
82  }
83 
84  return (lp != NULL) && (env != NULL);
85 }
86 
89  int32_t* idx_bound, int32_t num_bound, int32_t* w_zero, int32_t num_zero,
90  float64_t* vee, int32_t num_dim, bool use_bias)
91 {
92  const int cmatsize=10*1024*1024; //no more than 10mil. elements
93  bool result=false;
94  int32_t num_variables=num_dim + num_bound+num_zero; // xi, beta
95 
96  ASSERT(num_dim>0)
97  ASSERT(num_dim>0)
98 
99  // setup LP part
100  float64_t* lb=SG_MALLOC(float64_t, num_variables);
101  float64_t* ub=SG_MALLOC(float64_t, num_variables);
102  float64_t* obj=SG_MALLOC(float64_t, num_variables);
103  char* sense = SG_MALLOC(char, num_dim);
104  int* cmatbeg=SG_MALLOC(int, num_variables);
105  int* cmatcnt=SG_MALLOC(int, num_variables);
106  int* cmatind=SG_MALLOC(int, cmatsize);
107  double* cmatval=SG_MALLOC(double, cmatsize);
108 
109  for (int32_t i=0; i<num_variables; i++)
110  {
111  obj[i]=0;
112 
113  if (i<num_dim) //xi part
114  {
115  lb[i]=-CPX_INFBOUND;
116  ub[i]=+CPX_INFBOUND;
117  }
118  else //beta part
119  {
120  lb[i]=0.0;
121  ub[i]=1.0;
122  }
123  }
124 
125  int32_t offs=0;
126  for (int32_t i=0; i<num_variables; i++)
127  {
128  if (i<num_dim) //sum_xi
129  {
130  sense[i]='E';
131  cmatbeg[i]=offs;
132  cmatcnt[i]=1;
133  cmatind[offs]=offs;
134  cmatval[offs]=1.0;
135  offs++;
136  ASSERT(offs<cmatsize)
137  }
138  else if (i<num_dim+num_zero) // Z_w*beta_w
139  {
140  cmatbeg[i]=offs;
141  cmatcnt[i]=1;
142  cmatind[offs]=w_zero[i-num_dim];
143  cmatval[offs]=-1.0;
144  offs++;
145  ASSERT(offs<cmatsize)
146  }
147  else // Z_x*beta_x
148  {
149  int32_t idx=idx_bound[i-num_dim-num_zero];
150  int32_t vlen=0;
151  bool vfree=false;
152  //SG_PRINT("idx=%d\n", idx)
154  //SG_PRINT("vlen=%d\n", vlen)
155 
156  cmatbeg[i]=offs;
157  cmatcnt[i]=vlen;
158 
159  float64_t val= -C*labels->get_confidence(idx);
160 
161  if (vlen>0)
162  {
163  for (int32_t j=0; j<vlen; j++)
164  {
165  cmatind[offs]=vec.features[j].feat_index;
166  cmatval[offs]=-val*vec.features[j].entry;
167  offs++;
168  ASSERT(offs<cmatsize)
169  //SG_PRINT("vec[%d]=%10.10f\n", j, vec.features[j].entry)
170  }
171 
172  if (use_bias)
173  {
174  cmatcnt[i]++;
175  cmatind[offs]=num_dim-1;
176  cmatval[offs]=-val;
177  offs++;
178  ASSERT(offs<cmatsize)
179  }
180  }
181  else
182  {
183  if (use_bias)
184  {
185  cmatcnt[i]++;
186  cmatind[offs]=num_dim-1;
187  cmatval[offs]=-val;
188  }
189  else
190  cmatval[offs]=0.0;
191  offs++;
192  ASSERT(offs<cmatsize)
193  }
194 
195  features->free_feature_vector(vec, idx);
196  }
197  }
198 
199  result = CPXcopylp(env, lp, num_variables, num_dim, CPX_MIN,
200  obj, vee, sense, cmatbeg, cmatcnt, cmatind, cmatval, lb, ub, NULL) == 0;
201 
202  if (!result)
203  SG_ERROR("CPXcopylp failed.\n")
204 
205  //write_problem("problem.lp");
206 
207  SG_FREE(sense);
208  SG_FREE(lb);
209  SG_FREE(ub);
210  SG_FREE(obj);
211  SG_FREE(cmatbeg);
212  SG_FREE(cmatcnt);
213  SG_FREE(cmatind);
214  SG_FREE(cmatval);
215 
217  int* qmatbeg=SG_MALLOC(int, num_variables);
218  int* qmatcnt=SG_MALLOC(int, num_variables);
219  int* qmatind=SG_MALLOC(int, num_variables);
220  double* qmatval=SG_MALLOC(double, num_variables);
221 
222  float64_t diag=2.0;
223 
224  for (int32_t i=0; i<num_variables; i++)
225  {
226  if (i<num_dim) //|| (!use_bias && i<num_dim)) //xi
227  //if ((i<num_dim-1) || (!use_bias && i<num_dim)) //xi
228  {
229  qmatbeg[i]=i;
230  qmatcnt[i]=1;
231  qmatind[i]=i;
232  qmatval[i]=diag;
233  }
234  else
235  {
236  //qmatbeg[i]= (use_bias) ? (num_dim-1) : (num_dim);
237  qmatbeg[i]= num_dim;
238  qmatcnt[i]=0;
239  qmatind[i]=0;
240  qmatval[i]=0;
241  }
242  }
243 
244  if (result)
245  result = CPXcopyquad(env, lp, qmatbeg, qmatcnt, qmatind, qmatval) == 0;
246 
247  SG_FREE(qmatbeg);
248  SG_FREE(qmatcnt);
249  SG_FREE(qmatind);
250  SG_FREE(qmatval);
251 
252  if (!result)
253  SG_ERROR("CPXcopyquad failed.\n")
254 
255  //write_problem("problem.lp");
256  //write_Q("problem.qp");
257 
258  return result;
259 }
260 
261 bool CCplex::setup_lpboost(float64_t C, int32_t num_cols)
262 {
263  init(E_LINEAR);
264  int32_t status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1); //primal simplex
265  if (status)
266  SG_ERROR("Failure to select dual lp optimization, error %d.\n", status)
267 
268  double* obj=SG_MALLOC(double, num_cols);
269  double* lb=SG_MALLOC(double, num_cols);
270  double* ub=SG_MALLOC(double, num_cols);
271 
272  for (int32_t i=0; i<num_cols; i++)
273  {
274  obj[i]=-1;
275  lb[i]=0;
276  ub[i]=C;
277  }
278 
279  status = CPXnewcols(env, lp, num_cols, obj, lb, ub, NULL, NULL);
280  if ( status )
281  {
282  char errmsg[1024];
283  CPXgeterrorstring (env, status, errmsg);
284  SG_ERROR("%s", errmsg)
285  }
286  SG_FREE(obj);
287  SG_FREE(lb);
288  SG_FREE(ub);
289  return status==0;
290 }
291 
293  float64_t factor, SGSparseVectorEntry<float64_t>* h, int32_t len, int32_t ulen,
294  CBinaryLabels* label)
295 {
296  int amatbeg[1]; /* for calling external lib */
297  int amatind[len+1]; /* for calling external lib */
298  double amatval[len+1]; /* for calling external lib */
299  double rhs[1]; /* for calling external lib */
300  char sense[1];
301 
302  amatbeg[0]=0;
303  rhs[0]=1;
304  sense[0]='L';
305 
306  for (int32_t i=0; i<len; i++)
307  {
308  int32_t idx=h[i].feat_index;
309  float64_t val=factor*h[i].entry;
310  amatind[i]=idx;
311  amatval[i]=label->get_confidence(idx)*val;
312  }
313 
314  int32_t status = CPXaddrows (env, lp, 0, 1, len, rhs, sense, amatbeg, amatind, amatval, NULL, NULL);
315 
316  if ( status )
317  SG_ERROR("Failed to add the new row.\n")
318 
319  return status == 0;
320 }
321 
323  float64_t C, CSparseFeatures<float64_t>* x, CBinaryLabels* y, bool use_bias)
324 {
325  ASSERT(x)
326  ASSERT(y)
327  int32_t num_vec=y->get_num_labels();
328  int32_t num_feat=x->get_num_features();
329  int64_t nnz=x->get_num_nonzero_entries();
330 
331  //number of variables: b,w+,w-,xi concatenated
332  int32_t num_dims=1+2*num_feat+num_vec;
333  int32_t num_constraints=num_vec;
334 
335  float64_t* lb=SG_MALLOC(float64_t, num_dims);
336  float64_t* ub=SG_MALLOC(float64_t, num_dims);
337  float64_t* f=SG_MALLOC(float64_t, num_dims);
338  float64_t* b=SG_MALLOC(float64_t, num_dims);
339 
340  //number of non zero entries in A (b,w+,w-,xi)
341  int64_t amatsize=((int64_t) num_vec)+nnz+nnz+num_vec;
342 
343  int* amatbeg=SG_MALLOC(int, num_dims); /* for calling external lib */
344  int* amatcnt=SG_MALLOC(int, num_dims); /* for calling external lib */
345  int* amatind=SG_MALLOC(int, amatsize); /* for calling external lib */
346  double* amatval= SG_MALLOC(double, amatsize); /* for calling external lib */
347 
348  for (int32_t i=0; i<num_dims; i++)
349  {
350  if (i==0) //b
351  {
352  if (use_bias)
353  {
354  lb[i]=-CPX_INFBOUND;
355  ub[i]=+CPX_INFBOUND;
356  }
357  else
358  {
359  lb[i]=0;
360  ub[i]=0;
361  }
362  f[i]=0;
363  }
364  else if (i<2*num_feat+1) //w+,w-
365  {
366  lb[i]=0;
367  ub[i]=CPX_INFBOUND;
368  f[i]=1;
369  }
370  else //xi
371  {
372  lb[i]=0;
373  ub[i]=CPX_INFBOUND;
374  f[i]=C;
375  }
376  }
377 
378  for (int32_t i=0; i<num_constraints; i++)
379  b[i]=-1;
380 
381  char* sense=SG_MALLOC(char, num_constraints);
382  memset(sense,'L',sizeof(char)*num_constraints);
383 
384  //construct A
385  int64_t offs=0;
386 
387  //b part of A
388  amatbeg[0]=offs;
389  amatcnt[0]=num_vec;
390 
391  for (int32_t i=0; i<num_vec; i++)
392  {
393  amatind[offs]=i;
394  amatval[offs]=-y->get_confidence(i);
395  offs++;
396  }
397 
398  //w+ and w- part of A
399  int32_t num_sfeat=0;
400  int32_t num_svec=0;
401  SGSparseVector<float64_t>* sfeat= x->get_transposed(num_sfeat, num_svec);
402  ASSERT(sfeat)
403 
404  for (int32_t i=0; i<num_svec; i++)
405  {
406  amatbeg[i+1]=offs;
407  amatcnt[i+1]=sfeat[i].num_feat_entries;
408 
409  for (int32_t j=0; j<sfeat[i].num_feat_entries; j++)
410  {
411  int32_t row=sfeat[i].features[j].feat_index;
412  float64_t val=sfeat[i].features[j].entry;
413 
414  amatind[offs]=row;
415  amatval[offs]=-y->get_confidence(row)*val;
416  offs++;
417  }
418  }
419 
420  for (int32_t i=0; i<num_svec; i++)
421  {
422  amatbeg[num_svec+i+1]=offs;
423  amatcnt[num_svec+i+1]=sfeat[i].num_feat_entries;
424 
425  for (int32_t j=0; j<sfeat[i].num_feat_entries; j++)
426  {
427  int32_t row=sfeat[i].features[j].feat_index;
428  float64_t val=sfeat[i].features[j].entry;
429 
430  amatind[offs]=row;
431  amatval[offs]=y->get_confidence(row)*val;
432  offs++;
433  }
434  }
435 
436  x->clean_tsparse(sfeat, num_svec);
437 
438  //xi part of A
439  for (int32_t k=0; k<num_vec; k++)
440  {
441  amatbeg[1+2*num_feat+k]=offs;
442  amatcnt[1+2*num_feat+k]=1;
443  amatind[offs]=k;
444  amatval[offs]=-1;
445  offs++;
446  }
447 
448  int32_t status = CPXsetintparam (env, CPX_PARAM_LPMETHOD, 1); //barrier
449  if (status)
450  SG_ERROR("Failure to select barrier optimization, error %d.\n", status)
451  CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON);
452 
453  bool result = CPXcopylp(env, lp, num_dims, num_constraints, CPX_MIN,
454  f, b, sense, amatbeg, amatcnt, amatind, amatval, lb, ub, NULL) == 0;
455 
456 
457  SG_FREE(amatval);
458  SG_FREE(amatcnt);
459  SG_FREE(amatind);
460  SG_FREE(amatbeg);
461  SG_FREE(b);
462  SG_FREE(f);
463  SG_FREE(ub);
464  SG_FREE(lb);
465 
466  return result;
467 }
468 
470 {
471  bool result=false;
472 
473  if (lp)
474  {
475  int32_t status = CPXfreeprob(env, &lp);
476  lp = NULL;
477  lp_initialized = false;
478 
479  if (status)
480  SG_WARNING("CPXfreeprob failed, error code %d.\n", status)
481  else
482  result = true;
483  }
484 
485  if (env)
486  {
487  int32_t status = CPXcloseCPLEX (&env);
488  env=NULL;
489 
490  if (status)
491  {
492  char errmsg[1024];
493  SG_WARNING("Could not close CPLEX environment.\n")
494  CPXgeterrorstring (env, status, errmsg);
495  SG_WARNING("%s", errmsg)
496  }
497  else
498  result = true;
499  }
500  return result;
501 }
502 
504  float64_t* H, int32_t rows, int32_t cols, int* &qmatbeg, int* &qmatcnt,
505  int* &qmatind, double* &qmatval)
506 {
507  qmatbeg=SG_MALLOC(int, cols);
508  qmatcnt=SG_MALLOC(int, cols);
509  qmatind=SG_MALLOC(int, cols*rows);
510  qmatval = H;
511 
512  if (!(qmatbeg && qmatcnt && qmatind))
513  {
514  SG_FREE(qmatbeg);
515  SG_FREE(qmatcnt);
516  SG_FREE(qmatind);
517  return false;
518  }
519 
520  for (int32_t i=0; i<cols; i++)
521  {
522  qmatcnt[i]=rows;
523  qmatbeg[i]=i*rows;
524  for (int32_t j=0; j<rows; j++)
525  qmatind[i*rows+j]=j;
526  }
527 
528  return true;
529 }
530 
532  float64_t* objective, float64_t* constraints_mat, int32_t rows,
533  int32_t cols, float64_t* rhs, float64_t* lb, float64_t* ub)
534 {
535  bool result=false;
536 
537  int* qmatbeg=NULL;
538  int* qmatcnt=NULL;
539  int* qmatind=NULL;
540  double* qmatval=NULL;
541 
542  char* sense = NULL;
543 
544  if (constraints_mat==0)
545  {
546  ASSERT(rows==0)
547  rows=1;
548  float64_t dummy=0;
549  rhs=&dummy;
550  sense=SG_MALLOC(char, rows);
551  memset(sense,'L',sizeof(char)*rows);
552  constraints_mat=SG_MALLOC(float64_t, cols);
553  memset(constraints_mat, 0, sizeof(float64_t)*cols);
554 
555  result=dense_to_cplex_sparse(constraints_mat, 0, cols, qmatbeg, qmatcnt, qmatind, qmatval);
556  ASSERT(result)
557  result = CPXcopylp(env, lp, cols, rows, CPX_MIN,
558  objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
559  SG_FREE(constraints_mat);
560  }
561  else
562  {
563  sense=SG_MALLOC(char, rows);
564  memset(sense,'L',sizeof(char)*rows);
565  result=dense_to_cplex_sparse(constraints_mat, rows, cols, qmatbeg, qmatcnt, qmatind, qmatval);
566  result = CPXcopylp(env, lp, cols, rows, CPX_MIN,
567  objective, rhs, sense, qmatbeg, qmatcnt, qmatind, qmatval, lb, ub, NULL) == 0;
568  }
569 
570  SG_FREE(sense);
571  SG_FREE(qmatbeg);
572  SG_FREE(qmatcnt);
573  SG_FREE(qmatind);
574 
575  if (!result)
576  SG_WARNING("CPXcopylp failed.\n")
577 
578  return result;
579 }
580 
581 bool CCplex::setup_qp(float64_t* H, int32_t dim)
582 {
583  int* qmatbeg=NULL;
584  int* qmatcnt=NULL;
585  int* qmatind=NULL;
586  double* qmatval=NULL;
587  bool result=dense_to_cplex_sparse(H, dim, dim, qmatbeg, qmatcnt, qmatind, qmatval);
588  if (result)
589  result = CPXcopyquad(env, lp, qmatbeg, qmatcnt, qmatind, qmatval) == 0;
590 
591  SG_FREE(qmatbeg);
592  SG_FREE(qmatcnt);
593  SG_FREE(qmatind);
594 
595  if (!result)
596  SG_WARNING("CPXcopyquad failed.\n")
597 
598  return result;
599 }
600 
602 {
603  int solnstat; /* for calling external lib */
604  double objval;
605  int status=1; /* for calling external lib */
606 
607  if (problem_type==E_QP)
608  status = CPXqpopt (env, lp);
609  else if (problem_type == E_LINEAR)
610  status = CPXlpopt (env, lp);
611 
612  if (status)
613  SG_WARNING("Failed to optimize QP.\n")
614 
615  status = CPXsolution (env, lp, &solnstat, &objval, sol, lambda, NULL, NULL);
616 
617  //SG_PRINT("obj:%f\n", objval)
618 
619  return (status==0);
620 }
621 #endif

SHOGUN Machine Learning Toolbox - Documentation