CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
VoronoiAlgorithm.cc
Go to the documentation of this file.
1 #include "RecoHI/HiJetAlgos/interface/VoronoiAlgorithm.h"
2 
3 extern "C" {
4 
5  extern struct ascal_s {
6  double objnor;
7  double rhsnor;
8  double scdiff;
9  int scpass;
10  int scalmet;
11  } ascal_;
12 
13  extern struct compl_s {
14  double climit;
15  double ccorr;
16  } compl_;
17 
18  extern struct dims_s {
19  int n;
20  int n1;
21  int m;
22  int mn;
23  int nz;
24  int cfree;
25  int pivotn;
26  int denwin;
27  int rfree;
28  } dims_;
29 
30  extern struct drop_s {
31  double tfixvar;
32  double tfixslack;
33  double slklim;
34  } drop_;
35 
36  extern struct factor_s {
37  double tpiv1;
38  double tpiv2;
39  double tabs;
40  double trabs;
41  double lam;
42  double tfind;
43  double order;
44  double supdens;
45  } factor_;
46 
47  extern struct initv_s {
48  double prmin;
49  double upmax;
50  double dumin;
51  int stamet;
52  int safmet;
53  int premet;
54  int regul;
55  } initv_;
56 
57  extern struct itref_s {
58  double tresx;
59  double tresy;
60  int maxref;
61  } itref_;
62 
63  extern struct mscal_s {
64  double varadd;
65  double slkadd;
66  double scfree;
67  } mscal_;
68 
69  extern struct numer_s {
70  double tplus;
71  double tzer;
72  } numer_;
73 
74  extern struct param_s {
75  double palpha;
76  double dalpha;
77  } param_;
78 
79  extern struct predc_s {
80  double target;
81  double tsmall;
82  double tlarge;
83  double center;
84  double corstp;
85  int mincc;
86  int maxcc;
87  } predc_;
88 
89  extern struct predp_s {
90  double ccstop;
91  double barset;
92  double bargrw;
93  double barmin;
94  int mincor;
95  int maxcor;
96  int inibar;
97  } predp_;
98 
99  extern struct setden_s {
100  double maxdense;
101  double densgap;
102  int setlam;
103  int denslen;
104  } setden_;
105 
106  extern struct sprnod_s {
107  int psupn;
108  int ssupn;
109  int maxsnz;
110  } sprnod_;
111 
112  extern struct toler_s {
113  double tsdir;
114  double topt1;
115  double topt2;
116  double tfeas1;
117  double tfeas2;
118  double feas1;
119  double feas2;
120  double pinfs;
121  double dinfs;
122  double inftol;
123  int maxiter;
124  } toler_;
125 
126  extern int solver_(
127  double *obj, double *rhs, double *lbound, double *ubound,
128  /* Size dims_.mn working arrays */
129  double *diag, double *odiag,
130  /* Size dims_.mn primal values */
131  double *xs,
132  /* Size dims_.mn working arrays */
133  double *dxs, double *dxsn, double *up,
134  /* Size dims_.mn dual residuals */
135  double *dspr,
136  /* Size dims_.mn working arrays */
137  double *ddspr, double *ddsprn, double *dsup, double *ddsup,
138  double *ddsupn,
139  /* Size dims_.m dual values */
140  double *dv,
141  /* Size dims_.m working arrays */
142  double *ddv, double *ddvn, double *prinf,
143  /* Size dims_.mn working arrays */
144  double *upinf, double *duinf, double *scale,
145  /* Size dims_.cfree nonzero values */
146  double *nonzeros,
147  /* Size dims_.n working array */
148  int *vartyp,
149  /* Size dims_.m working array */
150  int *slktyp,
151  /* Size dims_.n1 column pointer */
152  int *colpnt,
153  /* Size dims_.mn working arrays */
154  int *ecolpnt, int *count, int *vcstat, int *pivots,
155  int *invprm, int *snhead, int *nodtyp, int *inta1,
156  int *prehis,
157  /* Size dims_.cfree row index */
158  int *rowidx,
159  /* Size dims_.rfree working array */
160  int *rindex,
161  /* Scalar termination code */
162  int *code,
163  /* Scalar optimum value */
164  double *opt,
165  int *iter, int *corect, int *fixn, int *dropn, int *fnzmax,
166  int *fnzmin, double *addobj, double *bigbou,
167  /* Practical +inf */
168  double *big,
169  int *ft);
170 }
171 
172 namespace {
173 
174  class problem_t {
175  public:
176  static constexpr double infinity = INFINITY;
177  // These are equivalent to the constants used by IBM ILOG
178  // CPLEX
179  static const int minimize = 1;
180  static const int maximize = -1;
181  static const char less_equal = 'L';
182  static const char equal = 'E';
183  static const char greater_equal = 'G';
184  static const char range = 'R';
185  protected:
186  bool _optimized;
187  bool _variable_named;
188 
189  std::vector<double> _rhs;
190  std::vector<char> _constraint_sense;
191  std::vector<double> _range_value;
192  std::vector<char *> _row_name;
193 
194  std::vector<double> _objective;
195  std::vector<double> _lower_bound;
196  std::vector<double> _upper_bound;
197  std::vector<char *> _column_name;
198 
199  std::vector<int> _coefficient_row;
200  std::vector<int> _coefficient_column;
201  std::vector<double> _coefficient_value;
202 
203  void clear_row(void)
204  {
205  _rhs.clear();
206  _constraint_sense.clear();
207  for (std::vector<char *>::const_iterator iterator =
208  _row_name.begin();
209  iterator != _row_name.end(); iterator++) {
210  delete [] *iterator;
211  }
212  _row_name.clear();
213  }
214  void clear_column(void)
215  {
216  _objective.clear();
217  _lower_bound.clear();
218  _upper_bound.clear();
219  for (std::vector<char *>::const_iterator iterator =
220  _column_name.begin();
221  iterator != _column_name.end(); iterator++) {
222  free(*iterator);
223  }
224  _column_name.clear();
225  }
226  void clear_coefficient(void)
227  {
228  _coefficient_row.clear();
229  _coefficient_column.clear();
230  _coefficient_value.clear();
231  }
232  void to_csr(std::vector<int> &column_pointer,
233  std::vector<int> &row_index,
234  std::vector<double> &nonzero_value,
235  const int index_start = 1)
236  {
237  // Convert from coordinate (COO) storage into compressed
238  // sparse row (CSR)
239 
240  std::vector<std::vector<std::pair<int, double> > >
241  column_index(_objective.size(),
242  std::vector<std::pair<int, double> >());
243 
244  std::vector<int>::const_iterator iterator_row =
245  _coefficient_row.begin();
246  std::vector<int>::const_iterator iterator_column =
247  _coefficient_column.begin();
248  std::vector<double>::const_iterator iterator_value =
249  _coefficient_value.begin();
250 
251  for (; iterator_value != _coefficient_value.end();
252  iterator_row++, iterator_column++, iterator_value++) {
253  column_index[*iterator_column].push_back(
254  std::pair<int, double>(
255  // Conversion into Fortran indexing
256  *iterator_row + index_start, *iterator_value));
257  }
258 
259  for (std::vector<std::vector<std::pair<int, double> > >::
260  iterator iterator_outer = column_index.begin();
261  iterator_outer != column_index.end(); iterator_outer++) {
262  // Conversion into Fortran indexing
263  column_pointer.push_back(row_index.size() + index_start);
264  std::sort(iterator_outer->begin(), iterator_outer->end());
265  for (std::vector<std::pair<int, double> >::const_iterator
266  iterator_inner = iterator_outer->begin();
267  iterator_inner != iterator_outer->end();
268  iterator_inner++) {
269  row_index.push_back(iterator_inner->first);
270  nonzero_value.push_back(iterator_inner->second);
271  }
272  }
273  // Conversion into Fortran indexing
274  column_pointer.push_back(row_index.size() + index_start);
275  }
276  public:
277  problem_t(const bool variable_named)
278  : _optimized(false), _variable_named(variable_named)
279  {
280  }
281  ~problem_t(void)
282  {
283  clear();
284  }
285  void clear(void)
286  {
287  clear_row();
288  clear_column();
289  clear_coefficient();
290  }
291  virtual int populate(void) = 0;
292  size_t nrow(void) const
293  {
294  return _rhs.size();
295  }
296  size_t ncolumn(void) const
297  {
298  return _objective.size();
299  }
300  void push_back_row(const char constraint_sense,
301  const double rhs)
302  {
303  _rhs.push_back(rhs);
304  _constraint_sense.push_back(constraint_sense);
305 
306  if (_variable_named) {
307  static const size_t name_length = 24;
308  char *name = new char[name_length];
309 
310  snprintf(name, name_length, "c%llu",
311  static_cast<unsigned long long>(_rhs.size()));
312  _row_name.push_back(name);
313  }
314  }
315  void push_back_row(const std::string &constraint_sense,
316  const double rhs)
317  {
318  char cplex_sense;
319 
320  if (constraint_sense == "<=") {
321  cplex_sense = 'L';
322  push_back_row(rhs, cplex_sense);
323  }
324  else if (constraint_sense == "==") {
325  cplex_sense = 'E';
326  push_back_row(rhs, cplex_sense);
327  }
328  else if (constraint_sense == ">=") {
329  cplex_sense = 'G';
330  push_back_row(rhs, cplex_sense);
331  }
332  else {
333  fprintf(stderr, "%s:%d: illegal sense (`%s')\n",
334  __FILE__, __LINE__, constraint_sense.c_str());
335  }
336  }
337  void push_back_column(const double objective,
338  const double lower_bound,
339  const double upper_bound)
340  {
341  _objective.push_back(objective);
342  _lower_bound.push_back(lower_bound);
343  _upper_bound.push_back(upper_bound);
344 
345  if (_variable_named) {
346  static const size_t name_length = 24;
347  char *name = new char[name_length];
348 
349  snprintf(name, name_length, "x%llu",
350  static_cast<unsigned long long>(
351  _objective.size()));
352  _column_name.push_back(name);
353  }
354  }
355  void push_back_coefficient(
356  const int row, const int column, const double value)
357  {
358  _coefficient_row.push_back(row);
359  _coefficient_column.push_back(column);
360  _coefficient_value.push_back(value);
361  }
362  virtual int optimize(void) = 0;
363  int optimize_primal_simplex(void)
364  {
365  return optimize();
366  }
367  int optimize_dual_simplex(void)
368  {
369  return optimize();
370  }
371  int optimize_barrier(void)
372  {
373  return optimize();
374  }
375  int optimize_network(void)
376  {
377  return optimize();
378  }
379  int optimize_sifting(void)
380  {
381  return optimize();
382  }
383  int optimize_concurrent(void)
384  {
385  return optimize();
386  }
387  int optimize_deterministic_concurrent(void)
388  {
389  return optimize();
390  }
391  //
392  virtual int solve(
393  int &solution_status, double &objective_value,
394  std::vector<double> &variable_primal,
395  std::vector<double> &variable_dual,
396  std::vector<double> &variable_slack_surplus,
397  std::vector<double> &reduced_cost) = 0;
398  int solve(
399  double &objective_value,
400  std::vector<double> &variable_primal,
401  std::vector<double> &variable_dual,
402  std::vector<double> &variable_slack_surplus,
403  std::vector<double> &reduced_cost)
404  {
405  int solution_status;
406 
407  return solve(solution_status, objective_value,
408  variable_primal, variable_dual,
409  variable_slack_surplus, reduced_cost);
410  }
411  int solve(
412  std::vector<double> &variable_primal,
413  std::vector<double> &variable_dual,
414  std::vector<double> &variable_slack_surplus,
415  std::vector<double> &reduced_cost)
416  {
417  int solution_status;
418  double objective_value;
419 
420  return solve(solution_status, objective_value,
421  variable_primal, variable_dual,
422  variable_slack_surplus, reduced_cost);
423  }
424  virtual int solve(
425  int &solution_status, double &objective_value,
426  std::vector<double> &variable_primal,
427  std::vector<double> &variable_dual)
428  {
429  std::vector<double> variable_slack_surplus;
430  std::vector<double> reduced_cost;
431 
432  return solve(solution_status, objective_value,
433  variable_primal, variable_dual,
434  variable_slack_surplus, reduced_cost);
435  }
436  int solve(
437  std::vector<double> &variable_primal,
438  std::vector<double> &variable_dual)
439  {
440  int solution_status;
441  double objective_value;
442 
443  return solve(solution_status, objective_value, variable_primal,
444  variable_dual);
445  }
446  int solve(
447  std::vector<double> &variable_primal)
448  {
449  std::vector<double> variable_dual;
450 
451  return solve(variable_primal, variable_dual);
452  }
453  };
454 
455  class environment_t {
456  protected:
457  std::vector<problem_t *> _problem;
458  };
459 
460 #ifndef BPMPD_INFINITY_BOUND
461 #define BPMPD_INFINITY_BOUND 1e+30
462 #endif // BPMPD_INFINITY_BOUND
463 
464 #ifndef BPMPD_VERSION_STRING
465 #define BPMPD_VERSION_STRING "Version 2.11 (1996 December)"
466 #endif // BPMPD_VERSION_STRING
467 
468  class bpmpd_environment_t;
469 
470  class bpmpd_problem_t : public problem_t {
471  public:
472  static constexpr double infinity = BPMPD_INFINITY_BOUND;
473  protected:
474  double _objective_sense;
475  double _objective_value;
476  std::vector<double> _variable_primal;
477  std::vector<double> _variable_dual;
478  int _solution_status;
479  private:
480  int ampl_solution_status(const int termination_code)
481  {
482  switch (termination_code) {
483  // General memory limit (no solution)
484  case -2: return 520;
485  // Memory limit during iterations
486  case -1: return 510;
487  // No optimum
488  case 1: return 506;
489  // Optimal solution
490  case 2: return 0;
491  // Primal infeasible
492  case 3: return 503;
493  // Dual infeasible
494  case 4: return 503;
495  default:
496  fprintf(stderr, "%s:%d: error: unknown termination code "
497  "%d\n", __FILE__, __LINE__, termination_code);
498  return 500;
499  }
500  fprintf(stderr, "%s:%d: %d\n", __FILE__, __LINE__,
501  termination_code);
502  }
503  void set_bpmpd_parameter(void)
504  {
505  // Set the parameter as in subroutine BPMAIN
506 
507  setden_.maxdense = 0.15;
508  setden_.densgap = 3.00;
509  setden_.setlam = 0;
510  factor_.supdens = 250;
511  setden_.denslen = 10;
512 
513  mscal_.varadd = 1e-12;
514  mscal_.slkadd = 1e+16;
515  mscal_.scfree = 1.1e-6;
516 
517  factor_.tpiv1 = 1e-3;
518  factor_.tpiv2 = 1e-8;
519  factor_.tabs = 1e-12;
520  factor_.trabs = 1e-15;
521  factor_.lam = 1e-5;
522  factor_.tfind = 25;
523  factor_.order = 2;
524 
525  sprnod_.psupn = 3;
526  sprnod_.ssupn = 4;
527  sprnod_.maxsnz = 9999999;
528 
529  toler_.tsdir = 1e-16;
530  toler_.topt1 = 1e-8;
531  toler_.topt2 = 1e-16;
532  toler_.tfeas1 = 1e-7;
533  toler_.tfeas2 = 1e-7;
534  toler_.feas1 = 1e-2;
535  toler_.feas2 = 1e-2;
536  toler_.pinfs = 1e-6;
537  toler_.dinfs = 1e-6;
538  toler_.inftol = 1e+4;
539  toler_.maxiter = 99;
540 
541  numer_.tplus = 1e-11;
542  numer_.tzer = 1e-35;
543 
544  itref_.tresx = 1e-9;
545  itref_.tresy = 1e-9;
546  itref_.maxref = 5;
547 
548  ascal_.objnor = 1e+2;
549  ascal_.rhsnor = 0;
550  ascal_.scdiff = 1;
551  ascal_.scalmet = 2;
552  ascal_.scpass = 5;
553 
554  predp_.ccstop = 1.01;
555  predp_.barset = 2.50e-1;
556  predp_.bargrw = 1.00e+2;
557  predp_.barmin = 1.00e-10;
558  predp_.mincor = 1;
559  predp_.maxcor = 1;
560  predp_.inibar = 0;
561 
562  predc_.target = 9e-2;
563  predc_.tsmall = 2e-1;
564  predc_.tlarge = 2e+1;
565  predc_.center = 5;
566  predc_.corstp = 1.01;
567  predc_.mincc = 0;
568  predc_.maxcc = 9;
569 
570  param_.palpha = 0.999;
571  param_.dalpha = 0.999;
572 
573  drop_.tfixvar = 1e-16;
574  drop_.tfixslack = 1e-16;
575  drop_.slklim = 1e-16;
576 
577  initv_.prmin = 100;
578  initv_.upmax = 50000;
579  initv_.dumin = 100;
580  initv_.stamet = 2;
581  initv_.safmet = -3;
582  initv_.premet = 511;
583  initv_.regul = 0;
584 
585  compl_.climit = 1;
586  compl_.ccorr = 1e-5;
587  }
588  void append_constraint_sense_bound(
589  std::vector<double> &lower_bound_bpmpd,
590  std::vector<double> &upper_bound_bpmpd)
591  {
592  for (std::vector<char>::const_iterator iterator =
593  _constraint_sense.begin();
594  iterator != _constraint_sense.end(); iterator++) {
595  switch (*iterator) {
596  case less_equal:
597  lower_bound_bpmpd.push_back(-BPMPD_INFINITY_BOUND);
598  upper_bound_bpmpd.push_back(0);
599  break;
600  case equal:
601  lower_bound_bpmpd.push_back(0);
602  upper_bound_bpmpd.push_back(0);
603  break;
604  case greater_equal:
605  lower_bound_bpmpd.push_back(0);
606  upper_bound_bpmpd.push_back(BPMPD_INFINITY_BOUND);
607  break;
608  case range:
609  {
610  const size_t index =
611  iterator - _constraint_sense.begin();
612 
613  lower_bound_bpmpd.push_back(0);
614  upper_bound_bpmpd.push_back(_range_value[index] -
615  _rhs[index]);
616  }
617  break;
618  }
619  }
620  lower_bound_bpmpd.reserve(dims_.mn);
621  lower_bound_bpmpd.insert(lower_bound_bpmpd.end(),
622  _lower_bound.begin(),
623  _lower_bound.end());
624  upper_bound_bpmpd.reserve(dims_.mn);
625  upper_bound_bpmpd.insert(upper_bound_bpmpd.end(),
626  _upper_bound.begin(),
627  _upper_bound.end());
628  }
629  protected:
630  int populate(void)
631  {
632  return 0;
633  }
634  public:
635  bpmpd_problem_t(void)
636  : problem_t(false), _objective_sense(1.0),
637  _objective_value(NAN)
638  {
639  }
640  inline size_t nrow(void) const
641  {
642  return dims_.m;
643  }
644  inline size_t ncolumn(void) const
645  {
646  return dims_.n;
647  }
648  inline void set_objective_sense(int sense)
649  {
650  _objective_sense = sense;
651 
652  // This will be multiplied to the objective coefficients
653  // (i.e. sign flip when sense = -1 for maximization)
654  }
655  inline void save(std::string filename)
656  {
657  // MPS save?
658  }
659  void push_back_column(const double objective,
660  const double lower_bound = 0,
661  const double upper_bound = BPMPD_INFINITY_BOUND)
662  {
663  problem_t::push_back_column(objective, lower_bound,
664  upper_bound);
665  }
666  void set_size(const size_t nonzero_value_size,
667  const double mf = 6.0, const size_t ma = 0)
668  {
669  dims_.n = _objective.size();
670  dims_.m = _rhs.size();
671  dims_.mn = dims_.m + dims_.n;
672  dims_.n1 = dims_.n + 1;
673 
674  dims_.nz = nonzero_value_size;
675 
676  // Adapted from the AMPL interface
677  // http://www.netlib.org/ampl/solvers/bpmpd/
678 
679  const size_t rp_23 = 4 * dims_.m + 18 * dims_.mn;
680 
681  dims_.cfree = static_cast<int>(
682  std::max(2.0, mf) * (rp_23 + dims_.nz)) +
683  std::max(static_cast<size_t>(0), ma);
684  dims_.rfree = ((dims_.cfree + rp_23) >> 1) + 11 * dims_.m +
685  8 * dims_.n;
686  }
687  int optimize(void)
688  {
689  // Usually referred to as the variable "IA" for the CSR format
690  std::vector<int> column_pointer;
691  // Usually referred to as the variable "JA" for the CSR format
692  std::vector<int> row_index;
693  // Usually referred to as the variable "A" for the CSR format
694  std::vector<double> nonzero_value;
695 
696  to_csr(column_pointer, row_index, nonzero_value);
697  std::transform(_objective.begin(), _objective.end(),
698  _objective.begin(),
699  std::bind1st(std::multiplies<double>(),
700  _objective_sense));
701 
702  // Try 1x, 2x, and 4x the default memory allocation
703  for (size_t i = 0; i < 3; i++) {
704  set_size(nonzero_value.size(), 6.0 * (1 << i));
705 
706  nonzero_value.resize(dims_.cfree);
707  row_index.resize(dims_.cfree);
708 
709  set_bpmpd_parameter();
710 
711  std::vector<double> diag(dims_.mn, 0);
712  std::vector<double> odiag(dims_.mn, 0);
713  std::vector<double> dxs(dims_.mn, 0);
714  std::vector<double> dxsn(dims_.mn, 0);
715  std::vector<double> up(dims_.mn, 0);
716  std::vector<double> dual_residual(dims_.mn, 0);
717  std::vector<double> ddspr(dims_.mn, 0);
718  std::vector<double> ddsprn(dims_.mn, 0);
719  std::vector<double> dsup(dims_.mn, 0);
720  std::vector<double> ddsup(dims_.mn, 0);
721  std::vector<double> ddsupn(dims_.mn, 0);
722  std::vector<double> ddv(dims_.m, 0);
723  std::vector<double> ddvn(dims_.m, 0);
724  std::vector<double> prinf(dims_.m, 0);
725  std::vector<double> upinf(dims_.mn, 0);
726  std::vector<double> duinf(dims_.mn, 0);
727  std::vector<double> scale(dims_.mn, 0);
728  std::vector<int> vartyp(dims_.n, 0);
729  std::vector<int> slktyp(dims_.m, 0);
730  std::vector<int> colpnt(dims_.n1, 0);
731  std::vector<int> ecolpnt(dims_.mn, 0);
732  std::vector<int> count(dims_.mn, 0);
733  std::vector<int> vcstat(dims_.mn, 0);
734  std::vector<int> pivots(dims_.mn, 0);
735  std::vector<int> invprm(dims_.mn, 0);
736  std::vector<int> snhead(dims_.mn, 0);
737  std::vector<int> nodtyp(dims_.mn, 0);
738  std::vector<int> inta1(dims_.mn, 0);
739  std::vector<int> prehis(dims_.mn, 0);
740  std::vector<int> rindex(dims_.rfree, 0);
741  int termination_code;
742  int iter;
743  int correct;
744  int fixn;
745  int dropn;
746  int fnzmax = 0;
747  int fnzmin = -1;
748  double addobj = 0;
749  // Numerical limit of bounds to be ignored
750  double bigbou = 1e+15;
751  double infinity_copy = BPMPD_INFINITY_BOUND;
752  int ft;
753 
754  _variable_primal.resize(dims_.mn);
755  _variable_dual.resize(dims_.m);
756 
757  std::vector<double> lower_bound_bpmpd = _lower_bound;
758  std::vector<double> upper_bound_bpmpd = _upper_bound;
759 
760  append_constraint_sense_bound(lower_bound_bpmpd,
761  upper_bound_bpmpd);
762 
763  solver_(&_objective[0], &_rhs[0], &lower_bound_bpmpd[0],
764  &upper_bound_bpmpd[0], &diag[0], &odiag[0],
765  &_variable_primal[0], &dxs[0], &dxsn[0], &up[0],
766  &dual_residual[0], &ddspr[0], &ddsprn[0],
767  &dsup[0], &ddsup[0], &ddsupn[0],
768  &_variable_dual[0], &ddv[0], &ddvn[0], &prinf[0],
769  &upinf[0], &duinf[0], &scale[0],
770  &nonzero_value[0], &vartyp[0], &slktyp[0],
771  &column_pointer[0], &ecolpnt[0], &count[0],
772  &vcstat[0], &pivots[0], &invprm[0], &snhead[0],
773  &nodtyp[0], &inta1[0], &prehis[0], &row_index[0],
774  &rindex[0], &termination_code, &_objective_value,
775  &iter, &correct, &fixn, &dropn, &fnzmax, &fnzmin,
776  &addobj, &bigbou, &infinity_copy, &ft);
777 
778  _objective_value *= _objective_sense;
779  _solution_status = ampl_solution_status(termination_code);
780  if (termination_code != -2) {
781  // No out-of-memory errors
782  break;
783  }
784  }
785 
786  return _solution_status == 0 ? 0 : 1;
787  }
788  int solve(
789  int &solution_status, double &objective_value,
790  std::vector<double> &variable_primal,
791  std::vector<double> &variable_dual,
792  std::vector<double> &variable_slack_surplus,
793  std::vector<double> &reduced_cost)
794  {
795  // This set of solution is not implemented yet (or readily
796  // available from BPMPD)
797 
798  return 1;
799  }
800  int solve(
801  int &solution_status, double &objective_value,
802  std::vector<double> &variable_primal,
803  std::vector<double> &variable_dual)
804  {
805  variable_primal = std::vector<double>(
806  _variable_primal.begin(),
807  _variable_primal.begin() + _objective.size());
808  variable_dual = _variable_dual;
809 
810  return 0;
811  }
812  friend class bpmpd_environment_t;
813  };
814 
815  class bpmpd_environment_t : public environment_t {
816  public:
817  bpmpd_environment_t(void)
818  {
819  }
820  ~bpmpd_environment_t(void)
821  {
822  }
823  int set_verbose(void);
824  int set_data_checking(void);
825  inline std::string version_str(void) const
826  {
827  return BPMPD_VERSION_STRING;
828  }
829  bpmpd_problem_t problem(std::string name = "")
830  {
831  return bpmpd_problem_t();
832  }
833  };
834 
835 }
836 
837 namespace {
838 
839  double angular_range_reduce(const double x)
840  {
841  if (!std::isfinite(x)) {
842  return NAN;
843  }
844 
845  static const double cody_waite_x_max = 1608.4954386379741381;
846  static const double two_pi_0 = 6.2831853071795649157;
847  static const double two_pi_1 = 2.1561211432631314669e-14;
848  static const double two_pi_2 = 1.1615423895917441336e-27;
849  double ret;
850 
851  if (x >= -cody_waite_x_max && x <= cody_waite_x_max) {
852  static const double inverse_two_pi =
853  0.15915494309189534197;
854  const double k = rint(x * inverse_two_pi);
855  ret = ((x - (k * two_pi_0)) - k * two_pi_1) -
856  k * two_pi_2;
857  }
858  else {
859  long double sin_x;
860  long double cos_x;
861 
862  sincosl(x, &sin_x, &cos_x);
863  ret = (double)atan2l(sin_x, cos_x);
864  }
865  if (ret == -M_PI) {
866  ret = M_PI;
867  }
868 
869  return ret;
870  }
871 
872 
873  size_t pf_id_reduce(const int pf_id)
874  {
875  // Particle::pdgId_ PFCandidate::particleId_
876  // PFCandidate::ParticleType Particle
877  // 0 0 X unknown, or dummy
878  // +211, -211 1 h charged hadron
879  // +11, -11 2 e electron
880  // +13, -13 3 mu muon
881  // 22 4 gamma photon
882  // 130 5 h0 neutral hadron
883  // 130 6 h_HF hadronic energy in an HF tower
884  // 22 7 egamma_HF electromagnetic energy in an HF tower
885 
886  if (pf_id == 4) {
887  return 1;
888  }
889  else if (pf_id >= 5 && pf_id <= 7) {
890  return 2;
891  }
892 
893  return 0;
894  }
895 
896 }
897 
898  void VoronoiAlgorithm::initialize_geometry(void)
899  {
900  static const size_t ncms_hcal_edge_pseudorapidity = 82 + 1;
901  static const double cms_hcal_edge_pseudorapidity[
902  ncms_hcal_edge_pseudorapidity] = {
903  -5.191, -4.889, -4.716, -4.538, -4.363, -4.191, -4.013,
904  -3.839, -3.664, -3.489, -3.314, -3.139, -2.964, -2.853,
905  -2.650, -2.500, -2.322, -2.172, -2.043, -1.930, -1.830,
906  -1.740, -1.653, -1.566, -1.479, -1.392, -1.305, -1.218,
907  -1.131, -1.044, -0.957, -0.879, -0.783, -0.696, -0.609,
908  -0.522, -0.435, -0.348, -0.261, -0.174, -0.087,
909  0.000,
910  0.087, 0.174, 0.261, 0.348, 0.435, 0.522, 0.609,
911  0.696, 0.783, 0.879, 0.957, 1.044, 1.131, 1.218,
912  1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830,
913  1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853,
914  2.964, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013,
915  4.191, 4.363, 4.538, 4.716, 4.889, 5.191
916  };
917 
918  _cms_hcal_edge_pseudorapidity = std::vector<double>(
919  cms_hcal_edge_pseudorapidity,
920  cms_hcal_edge_pseudorapidity +
921  ncms_hcal_edge_pseudorapidity);
922 
923  static const size_t ncms_ecal_edge_pseudorapidity = 344 + 1;
924 
925  for (size_t i = 0; i < ncms_ecal_edge_pseudorapidity; i++) {
926  _cms_ecal_edge_pseudorapidity.push_back(
927  i * (2 * 2.9928 /
928  (ncms_ecal_edge_pseudorapidity - 1)) -
929  2.9928);
930  }
931  }
932  void VoronoiAlgorithm::allocate(void)
933  {
934  _perp_fourier = new boost::multi_array<double, 4>(
935  boost::extents[_edge_pseudorapidity.size() - 1]
936  [nreduced_particle_flow_id][nfourier][2]);
937  }
938  void VoronoiAlgorithm::deallocate(void)
939  {
940  delete _perp_fourier;
941  }
942  void VoronoiAlgorithm::event_fourier(void)
943  {
944  std::fill(_perp_fourier->data(),
945  _perp_fourier->data() +
946  _perp_fourier->num_elements(),
947  0);
948 
949  for (std::vector<particle_t>::const_iterator iterator =
950  _event.begin();
951  iterator != _event.end(); iterator++) {
952  const unsigned int reduced_id =
953  iterator->reduced_particle_flow_id;
954 
955  for (size_t k = 1; k < _edge_pseudorapidity.size();
956  k++) {
957  if (iterator->momentum.Eta() >=
958  _edge_pseudorapidity[k - 1] &&
959  iterator->momentum.Eta() <
960  _edge_pseudorapidity[k]) {
961  const double azimuth =
962  iterator->momentum.Phi();
963 
964  for (size_t l = 0; l < nfourier; l++) {
965  (*_perp_fourier)[k - 1][reduced_id]
966  [l][0] +=
967  iterator->momentum.Pt() *
968  cos(l * azimuth);
969  (*_perp_fourier)[k - 1][reduced_id]
970  [l][1] +=
971  iterator->momentum.Pt() *
972  sin(l * azimuth);
973  }
974  }
975  }
976  }
977  }
978  void VoronoiAlgorithm::feature_extract(void)
979  {
980  const size_t nfeature = 2 * nfourier - 1;
981 
982  _feature.resize(nfeature);
983 
984  // Scale factor to get 95% of the coefficient below 1.0
985  // (where however one order of magnitude tolerance is
986  // acceptable). This is valid for nfourier < 18 (where
987  // interference behavior with the HF geometry starts to
988  // appear)
989 
990  std::vector<double> scale(nfourier, 1.0 / 200.0);
991 
992  if (nfourier >= 1) {
993  scale[0] = 1.0 / 5400.0;
994  }
995  if (nfourier >= 2) {
996  scale[1] = 1.0 / 130.0;
997  }
998  if (nfourier >= 3) {
999  scale[2] = 1.0 / 220.0;
1000  }
1001 
1002  const size_t index_edge_end =
1003  _edge_pseudorapidity.size() - 2;
1004 
1005  _feature[0] = 0;
1006  for (size_t j = 0; j < nreduced_particle_flow_id; j++) {
1007  _feature[0] += scale[0] *
1008  ((*_perp_fourier)[0 ][j][0][0] +
1009  (*_perp_fourier)[index_edge_end][j][0][0]);
1010  }
1011 
1012  for (size_t k = 1; k < nfourier; k++) {
1013  _feature[2 * k - 1] = 0;
1014  for (size_t j = 0; j < nreduced_particle_flow_id; j++) {
1015  _feature[2 * k - 1] += scale[k] *
1016  ((*_perp_fourier)[0 ][j][k][0] +
1017  (*_perp_fourier)[index_edge_end][j][k][0]);
1018  }
1019  _feature[2 * k] = 0;
1020  for (size_t j = 0; j < nreduced_particle_flow_id; j++) {
1021  _feature[2 * k] += scale[k] *
1022  ((*_perp_fourier)[0 ][j][k][1] +
1023  (*_perp_fourier)[index_edge_end][j][k][1]);
1024  }
1025  }
1026 
1027 
1028 #if 0
1029  const double event_plane = atan2(_feature[4], _feature[3]);
1030  const double v2 =
1031  sqrt(_feature[3] * _feature[3] +
1032  _feature[4] * _feature[4]) / _feature[0];
1033 #endif
1034  }
1035  void VoronoiAlgorithm::voronoi_area_incident(void)
1036  {
1037  // Make the Voronoi diagram
1038 
1039  voronoi_diagram_t diagram;
1040 
1041  // Reverse Voronoi face lookup
1042 #ifdef HAVE_SPARSEHASH
1043  // The "empty" or default value of the hash table
1044  const voronoi_diagram_t::Face face_empty;
1045  google::dense_hash_map<voronoi_diagram_t::Face_handle,
1046  size_t, hash<voronoi_diagram_t::Face_handle> >
1047  face_index;
1048 
1049  face_index.set_empty_key(face_empty);
1050 #else // HAVE_SPARSEHASH
1051  std::map<voronoi_diagram_t::Face_handle, size_t>
1052  face_index;
1053 #endif // HAVE_SPARSEHASH
1054 
1055  for (std::vector<particle_t>::const_iterator iterator =
1056  _event.begin();
1057  iterator != _event.end(); iterator++) {
1058  // Make two additional replicas with azimuth +/- 2 pi
1059  // (and use only the middle) to mimick the azimuthal
1060  // cyclicity
1061  for (int k = -1; k <= 1; k++) {
1062  const point_2d_t p(
1063  iterator->momentum.Eta(),
1064  iterator->momentum.Phi() +
1065  k * (2 * M_PI));
1066  const voronoi_diagram_t::Face_handle handle =
1067  diagram.insert(p);
1068 
1069  face_index[handle] = iterator - _event.begin();
1070  }
1071  }
1072 
1073  // Extract the Voronoi cells as polygon and calculate the
1074  // area associated with individual particles
1075 
1076  for (std::vector<particle_t>::iterator iterator =
1077  _event.begin();
1078  iterator != _event.end(); iterator++) {
1079  const voronoi_diagram_t::Locate_result result =
1080  diagram.locate(*iterator);
1081  const voronoi_diagram_t::Face_handle *face =
1082  boost::get<voronoi_diagram_t::Face_handle>(
1083  &result);
1084  double polygon_area;
1085 
1086  if (face != NULL) {
1087  voronoi_diagram_t::Ccb_halfedge_circulator
1088  circulator_start = (*face)->outer_ccb();
1089  bool unbounded = false;
1090  polygon_t polygon;
1091 
1092  voronoi_diagram_t::Ccb_halfedge_circulator
1093  circulator = circulator_start;
1094 
1095  // Circle around the edges and extract the polygon
1096  // vertices
1097  do {
1098  if (circulator->has_target()) {
1099  polygon.push_back(
1100  circulator->target()->point());
1101  _event[face_index[*face]].incident.
1102  insert(
1103  _event.begin() +
1104  face_index[circulator->twin()->
1105  face()]);
1106  }
1107  else {
1108  unbounded = true;
1109  break;
1110  }
1111  }
1112  while (++circulator != circulator_start);
1113  if (unbounded) {
1114  polygon_area = INFINITY;
1115  }
1116  else {
1117  polygon_area = polygon.area();
1118  }
1119  }
1120  else {
1121  polygon_area = NAN;
1122  }
1123  iterator->area = fabs(polygon_area);
1124  }
1125  }
1126  void VoronoiAlgorithm::subtract_momentum(void)
1127  {
1128  for (std::vector<particle_t>::iterator iterator =
1129  _event.begin();
1130  iterator != _event.end(); iterator++) {
1131  int predictor_index = -1;
1132  int interpolation_index = -1;
1133  double density = 0;
1134  double pred_0 = 0;
1135 
1136  for (size_t l = 1; l < _edge_pseudorapidity.size(); l++) {
1137  if (iterator->momentum.Eta() >=
1138  _edge_pseudorapidity[l - 1] &&
1139  iterator->momentum.Eta() <
1140  _edge_pseudorapidity[l]) {
1141  predictor_index = l - 1;
1142  }
1143  }
1144 
1145  for (size_t j = 0; j < nreduced_particle_flow_id; j++) {
1146  if (j == 2) {
1147  // HCAL
1148  for (size_t l = 1;
1149  l < _cms_hcal_edge_pseudorapidity.size(); l++) {
1150  if (iterator->momentum.Eta() >=
1151  _cms_hcal_edge_pseudorapidity[l - 1] &&
1152  iterator->momentum.Eta() <
1153  _cms_hcal_edge_pseudorapidity[l]) {
1154  interpolation_index = l - 1;
1155  }
1156  }
1157  }
1158  else {
1159  // Tracks or ECAL clusters
1160  for (size_t l = 1;
1161  l < _cms_ecal_edge_pseudorapidity.size(); l++) {
1162  if (iterator->momentum.Eta() >=
1163  _cms_ecal_edge_pseudorapidity[l - 1] &&
1164  iterator->momentum.Eta() <
1165  _cms_ecal_edge_pseudorapidity[l]) {
1166  interpolation_index = l - 1;
1167  }
1168  }
1169  }
1170 
1171  if (predictor_index >= 0 && interpolation_index >= 0) {
1172  // Calculate the aggregated prediction and
1173  // interpolation for the pseudorapidity segment
1174 
1175  const double azimuth = iterator->momentum.Phi();
1176  const float (*p)[2][82] =
1177 #ifdef STANDALONE
1178  ue_predictor_pf[j][predictor_index]
1179 #else // STANDALONE
1180  ue->ue_predictor_pf[j][predictor_index]
1181 #endif // STANDALONE
1182  ;
1183  double pred = 0;
1184 
1185  for (size_t l = 0; l < nfourier; l++) {
1186  for (size_t m = 0; m < 2; m++) {
1187  float u = p[l][m][0];
1188 
1189  for (size_t n = 0; n < 2 * nfourier - 1; n++) {
1190  u += (((((((((p[l][m][9 * n + 9]) *
1191  _feature[n] +
1192  p[l][m][9 * n + 8]) *
1193  _feature[n] +
1194  p[l][m][9 * n + 7]) *
1195  _feature[n] +
1196  p[l][m][9 * n + 6]) *
1197  _feature[n] +
1198  p[l][m][9 * n + 5]) *
1199  _feature[n] +
1200  p[l][m][9 * n + 4]) *
1201  _feature[n] +
1202  p[l][m][9 * n + 3]) *
1203  _feature[n] +
1204  p[l][m][9 * n + 2]) *
1205  _feature[n] +
1206  p[l][m][9 * n + 1]) *
1207  _feature[n];
1208  }
1209 
1210  pred += u * (l == 0 ? 1.0 : 2.0) *
1211  (m == 0 ? cos(l * azimuth) :
1212  sin(l * azimuth));
1213  if (l == 0 && m == 0) {
1214  pred_0 += u /
1215  (2.0 * M_PI *
1216  (_edge_pseudorapidity[predictor_index + 1] -
1217  _edge_pseudorapidity[predictor_index]));
1218  }
1219  }
1220  }
1221 
1222  double interp;
1223 
1224 #ifdef STANDALONE
1225  if (j == 0) {
1226  interp =
1227  ue_interpolation_pf0[predictor_index][
1228  interpolation_index];
1229  }
1230  else if (j == 1) {
1231  interp =
1232  ue_interpolation_pf1[predictor_index][
1233  interpolation_index];
1234  }
1235  else if (j == 2) {
1236  interp =
1237  ue_interpolation_pf2[predictor_index][
1238  interpolation_index];
1239  }
1240 #else // STANDALONE
1241  if (j == 0) {
1242  interp =
1243  ue->ue_interpolation_pf0[predictor_index][
1244  interpolation_index];
1245  }
1246  else if (j == 1) {
1247  interp =
1248  ue->ue_interpolation_pf1[predictor_index][
1249  interpolation_index];
1250  }
1251  else if (j == 2) {
1252  interp =
1253  ue->ue_interpolation_pf2[predictor_index][
1254  interpolation_index];
1255  }
1256 #endif // STANDALONE
1257  // Interpolate down to the finely binned
1258  // pseudorapidity
1259 
1260  density += pred /
1261  (2.0 * M_PI *
1262  (_edge_pseudorapidity[predictor_index + 1] -
1263  _edge_pseudorapidity[predictor_index])) *
1264  interp;
1265  }
1266  }
1267 
1268  if (std::isfinite(iterator->area) && density >= 0) {
1269  // Subtract the PF candidate by density times
1270  // Voronoi cell area
1271  iterator->momentum_perp_subtracted =
1272  iterator->momentum.Pt() -
1273  density * iterator->area;
1274  }
1275  else {
1276  iterator->momentum_perp_subtracted =
1277  iterator->momentum.Pt();
1278  }
1279  iterator->momentum_perp_subtracted_unequalized =
1280  iterator->momentum_perp_subtracted;
1281  }
1282  }
1283  void VoronoiAlgorithm::recombine_link(void)
1284  {
1285  boost::multi_array<double, 2> radial_distance_square(
1286  boost::extents[_event.size()][_event.size()]);
1287 
1288  for (std::vector<particle_t>::const_iterator
1289  iterator_outer = _event.begin();
1290  iterator_outer != _event.end(); iterator_outer++) {
1291  radial_distance_square
1292  [iterator_outer - _event.begin()]
1293  [iterator_outer - _event.begin()] = 0;
1294 
1295  for (std::vector<particle_t>::const_iterator
1296  iterator_inner = _event.begin();
1297  iterator_inner != iterator_outer;
1298  iterator_inner++) {
1299  const double deta = iterator_outer->momentum.Eta() -
1300  iterator_inner->momentum.Eta();
1301  const double dphi = angular_range_reduce(
1302  iterator_outer->momentum.Phi() -
1303  iterator_inner->momentum.Phi());
1304 
1305  radial_distance_square
1306  [iterator_outer - _event.begin()]
1307  [iterator_inner - _event.begin()] =
1308  deta * deta + dphi * dphi;
1309  radial_distance_square
1310  [iterator_inner - _event.begin()]
1311  [iterator_outer - _event.begin()] =
1312  radial_distance_square
1313  [iterator_outer - _event.begin()]
1314  [iterator_inner - _event.begin()];
1315  }
1316  }
1317 
1318  _active.clear();
1319 
1320  for (std::vector<particle_t>::const_iterator
1321  iterator_outer = _event.begin();
1322  iterator_outer != _event.end(); iterator_outer++) {
1323  double incident_area_sum = iterator_outer->area;
1324 
1325  for (std::set<std::vector<particle_t>::iterator>::
1326  const_iterator iterator_inner =
1327  iterator_outer->incident.begin();
1328  iterator_inner !=
1329  iterator_outer->incident.end();
1330  iterator_inner++) {
1331  incident_area_sum += (*iterator_inner)->area;
1332  }
1333  _active.push_back(incident_area_sum < 2.0);
1334  }
1335 
1336  _recombine.clear();
1337  _recombine_index = std::vector<std::vector<size_t> >(
1338  _event.size(), std::vector<size_t>());
1339  _recombine_unsigned = std::vector<std::vector<size_t> >(
1340  _event.size(), std::vector<size_t>());
1341  _recombine_tie.clear();
1342 
1343  // 36 cells corresponds to ~ 3 layers, note that for
1344  // hexagonal tiling, cell in proximity = 3 * layer *
1345  // (layer + 1)
1346  static const size_t npair_max = 36;
1347 
1348  for (size_t i = 0; i < _event.size(); i++) {
1349  for (size_t j = 0; j < _event.size(); j++) {
1350  const bool active_i_j = _active[i] && _active[j];
1351  const size_t incident_count =
1352  _event[i].incident.count(_event.begin() + j) +
1353  _event[j].incident.count(_event.begin() + i);
1354 
1355  if (active_i_j &&
1356  (radial_distance_square[i][j] <
1357  _radial_distance_square_max ||
1358  incident_count > 0)) {
1359  _recombine_unsigned[i].push_back(j);
1360  }
1361  }
1362 
1363  if (_event[i].momentum_perp_subtracted < 0) {
1364  std::vector<double> radial_distance_square_list;
1365 
1366  for (std::vector<size_t>::const_iterator iterator =
1367  _recombine_unsigned[i].begin();
1368  iterator != _recombine_unsigned[i].end();
1369  iterator++) {
1370  const size_t j = *iterator;
1371 
1372  if (_event[j].momentum_perp_subtracted > 0) {
1373  radial_distance_square_list.push_back(
1374  radial_distance_square[i][j]);
1375  }
1376  }
1377 
1378  double radial_distance_square_max_equalization_cut =
1379  _radial_distance_square_max;
1380 
1381  if (radial_distance_square_list.size() >= npair_max) {
1382  std::sort(radial_distance_square_list.begin(),
1383  radial_distance_square_list.end());
1384  radial_distance_square_max_equalization_cut =
1385  radial_distance_square_list[npair_max - 1];
1386  }
1387 
1388  for (std::vector<size_t>::const_iterator iterator =
1389  _recombine_unsigned[i].begin();
1390  iterator != _recombine_unsigned[i].end();
1391  iterator++) {
1392  const size_t j = *iterator;
1393 
1394  if (_event[j].momentum_perp_subtracted > 0 &&
1395  radial_distance_square[i][j] <
1396  radial_distance_square_max_equalization_cut) {
1397  _recombine_index[j].push_back(
1398  _recombine.size());
1399  _recombine_index[i].push_back(
1400  _recombine.size());
1401  _recombine.push_back(
1402  std::pair<size_t, size_t>(i, j));
1403  _recombine_tie.push_back(
1404  radial_distance_square[i][j] /
1405  _radial_distance_square_max);
1406  }
1407  }
1408  }
1409  }
1410  }
1411  void VoronoiAlgorithm::lp_populate(void *lp_problem)
1412  {
1413  bpmpd_problem_t *p = reinterpret_cast<bpmpd_problem_t *>(lp_problem);
1414 
1415  // The minimax problem is transformed into the LP notation
1416  // using the cost variable trick:
1417  //
1418  // Minimize c
1419  // Subject to:
1420  // c + sum_l t_kl + n_k >= 0 for negative cells n_k
1421  // c - sum_k t_kl + p_l >= 0 for positive cells p_l
1422 
1423  // Common LP mistakes during code development and their
1424  // CPLEX errors when running CPLEX in data checking mode:
1425  //
1426  // Error 1201 (column index ... out of range): Bad column
1427  // indexing, usually index_column out of bound for the
1428  // cost variables.
1429  //
1430  // Error 1222 (duplicate entry): Forgetting to increment
1431  // index_row, or index_column out of bound for the cost
1432  // variables.
1433 
1434  p->set_objective_sense(bpmpd_problem_t::minimize);
1435 
1436  // Rows (RHS of the constraints) of the LP problem
1437 
1438  static const size_t nsector_azimuth = 12;
1439 
1440  // Approximatively 2 pi / nsector_azimuth segmentation of
1441  // the CMS HCAL granularity
1442 
1443  static const size_t ncms_hcal_edge_pseudorapidity = 19 + 1;
1444  static const double cms_hcal_edge_pseudorapidity[
1445  ncms_hcal_edge_pseudorapidity] = {
1446  -5.191, -4.538, -4.013,
1447  -3.489, -2.853, -2.322, -1.830, -1.305, -0.783, -0.261,
1448  0.261, 0.783, 1.305, 1.830, 2.322, 2.853, 3.489,
1449  4.013, 4.538, 5.191
1450  };
1451 
1452  size_t nedge_pseudorapidity;
1453  const double *edge_pseudorapidity;
1454 
1455  nedge_pseudorapidity = ncms_hcal_edge_pseudorapidity;
1456  edge_pseudorapidity = cms_hcal_edge_pseudorapidity;
1457 
1458  const size_t nsuperblock = (nedge_pseudorapidity - 2) *
1459  nsector_azimuth;
1460 
1461  size_t index_row = 0;
1462  for (size_t index_pseudorapidity = 0;
1463  index_pseudorapidity < nedge_pseudorapidity - 2;
1464  index_pseudorapidity++) {
1465  for (size_t index_azimuth = 0;
1466  index_azimuth < nsector_azimuth - 1;
1467  index_azimuth++) {
1468  const size_t index_column =
1469  index_pseudorapidity * nsector_azimuth +
1470  index_azimuth;
1471  p->push_back_row(
1472  bpmpd_problem_t::greater_equal, 0);
1473  p->push_back_coefficient(
1474  index_row, index_column, 1);
1475  p->push_back_coefficient(
1476  index_row, nsuperblock + index_column, -1);
1477  index_row++;
1478  p->push_back_row(
1479  bpmpd_problem_t::greater_equal, 0);
1480  p->push_back_coefficient(
1481  index_row, index_column, 1);
1482  p->push_back_coefficient(
1483  index_row, nsuperblock + index_column + 1, -1);
1484  index_row++;
1485  p->push_back_row(
1486  bpmpd_problem_t::greater_equal, 0);
1487  p->push_back_coefficient(
1488  index_row, index_column, 1);
1489  p->push_back_coefficient(
1490  index_row,
1491  nsuperblock + index_column + nsector_azimuth, -1);
1492  index_row++;
1493  p->push_back_row(
1494  bpmpd_problem_t::greater_equal, 0);
1495  p->push_back_coefficient(
1496  index_row, index_column, 1);
1497  p->push_back_coefficient(
1498  index_row,
1499  nsuperblock + index_column + nsector_azimuth + 1,
1500  -1);
1501  index_row++;
1502  }
1503  const size_t index_column =
1504  index_pseudorapidity * nsector_azimuth +
1505  nsector_azimuth - 1;
1506  p->push_back_row(
1507  bpmpd_problem_t::greater_equal, 0);
1508  p->push_back_coefficient(
1509  index_row, index_column, 1);
1510  p->push_back_coefficient(
1511  index_row, nsuperblock + index_column, -1);
1512  index_row++;
1513  p->push_back_row(
1514  bpmpd_problem_t::greater_equal, 0);
1515  p->push_back_coefficient(
1516  index_row, index_column, 1);
1517  p->push_back_coefficient(
1518  index_row,
1519  nsuperblock + index_column - (nsector_azimuth - 1),
1520  -1);
1521  index_row++;
1522  p->push_back_row(
1523  bpmpd_problem_t::greater_equal, 0);
1524  p->push_back_coefficient(
1525  index_row, index_column, 1);
1526  p->push_back_coefficient(
1527  index_row,
1528  nsuperblock + index_column + nsector_azimuth, -1);
1529  index_row++;
1530  p->push_back_row(
1531  bpmpd_problem_t::greater_equal, 0);
1532  p->push_back_coefficient(
1533  index_row, index_column, 1);
1534  p->push_back_coefficient(
1535  index_row,
1536  nsuperblock + index_column + nsector_azimuth -
1537  (nsector_azimuth - 1),
1538  -1);
1539  index_row++;
1540  }
1541 
1542  const size_t nstaggered_block =
1543  (nedge_pseudorapidity - 1) * nsector_azimuth;
1544  const size_t nblock = nsuperblock + 2 * nstaggered_block;
1545 
1546  _nblock_subtract = std::vector<size_t>(_event.size(), 0);
1547 
1548  std::vector<size_t>
1549  positive_index(_event.size(), _event.size());
1550  size_t positive_count = 0;
1551 
1552  for (std::vector<particle_t>::const_iterator iterator =
1553  _event.begin();
1554  iterator != _event.end(); iterator++) {
1555  if (iterator->momentum_perp_subtracted >= 0) {
1556  positive_index[iterator - _event.begin()] =
1557  positive_count;
1558  positive_count++;
1559  }
1560  }
1561 
1562  _ncost = nblock + positive_count;
1563 
1564  const double sum_unequalized_0 = _equalization_threshold.first;
1565  const double sum_unequalized_1 = (2.0 / 3.0) * _equalization_threshold.first + (1.0 / 3.0) * _equalization_threshold.second;
1566  const double sum_unequalized_2 = (1.0 / 3.0) * _equalization_threshold.first + (2.0 / 3.0) * _equalization_threshold.second;
1567  const double sum_unequalized_3 = _equalization_threshold.second;
1568 
1569  std::vector<particle_t>::const_iterator
1570  iterator_particle = _event.begin();
1571  std::vector<bool>::const_iterator iterator_active =
1572  _active.begin();
1573  std::vector<std::vector<size_t> >::const_iterator
1574  iterator_recombine_index_outer =
1575  _recombine_index.begin();
1576  std::vector<std::vector<size_t> >::const_iterator
1577  iterator_recombine_unsigned_outer =
1578  _recombine_unsigned.begin();
1579  size_t index_column_max = _ncost - 1;
1580  for (; iterator_particle != _event.end();
1581  iterator_particle++, iterator_active++,
1582  iterator_recombine_index_outer++,
1583  iterator_recombine_unsigned_outer++) {
1584  if (*iterator_active) {
1585  int index_pseudorapidity = -1;
1586 
1588  for (size_t i = 1; i < nedge_pseudorapidity; i++) {
1589  if (iterator_particle->momentum.Eta() >= edge_pseudorapidity[i - 1] &&
1590  iterator_particle->momentum.Eta() < edge_pseudorapidity[i]) {
1591  index_pseudorapidity = i - 1;
1592  }
1593  }
1594 
1595  const int index_azimuth = floor(
1596  (iterator_particle->momentum.Phi() + M_PI) *
1597  ((nsector_azimuth >> 1) / M_PI));
1598 
1599  if (index_pseudorapidity != -1) {
1600  // p_i - sum t - u = c_i
1601  // or: c_i + u + sum_t = p_i
1602  // n_i + sum t - u <= 0
1603  // or: u - sum_t >= n_i
1604 
1605  // Inequality RHS
1606  p->push_back_row(
1607  iterator_particle->momentum_perp_subtracted >= 0 ?
1608  bpmpd_problem_t::equal :
1609  bpmpd_problem_t::greater_equal,
1610  iterator_particle->momentum_perp_subtracted);
1611 
1612  // Energy transfer coefficients t_kl
1613  const double sign = iterator_particle->momentum_perp_subtracted >= 0 ? 1 : -1;
1614  const size_t index_column_block_subtract =
1615  nsuperblock +
1616  (nedge_pseudorapidity - 1) * nsector_azimuth +
1617  index_pseudorapidity * nsector_azimuth +
1618  index_azimuth;
1619 
1620  _nblock_subtract[iterator_particle - _event.begin()] =
1621  index_column_block_subtract;
1622 
1623  if (iterator_particle->momentum_perp_subtracted >= 0) {
1624  const size_t index_column_cost =
1625  nblock + positive_index[iterator_particle - _event.begin()];
1626 
1627  p->push_back_coefficient(
1628  index_row, index_column_cost, 1);
1629  index_column_max =
1630  std::max(index_column_max, index_column_cost);
1631  }
1632  p->push_back_coefficient(
1633  index_row, index_column_block_subtract, 1);
1634  index_column_max =
1635  std::max(index_column_max, index_column_block_subtract);
1636 
1637  for (std::vector<size_t>::const_iterator
1638  iterator_recombine_index_inner =
1639  iterator_recombine_index_outer->begin();
1640  iterator_recombine_index_inner !=
1641  iterator_recombine_index_outer->end();
1642  iterator_recombine_index_inner++) {
1643  const size_t index_column =
1644  *iterator_recombine_index_inner +
1645  _ncost;
1646 
1647  p->push_back_coefficient(
1648  index_row, index_column, sign);
1649  index_column_max =
1650  std::max(index_column_max, index_column);
1651  }
1652  index_row++;
1653 
1654  const size_t index_column_block =
1655  nsuperblock +
1656  index_pseudorapidity * nsector_azimuth +
1657  index_azimuth;
1658 
1659  // sum_R c_i - o_i >= -d
1660  // or: d + sum_R c_i >= o_i
1661  // sum_R c_i - o_i <= d
1662  // or: d - sum_R c_i >= -o_i
1663 
1664  double sum_unequalized;
1665 
1666  sum_unequalized = 0;
1667  for (std::vector<size_t>::const_iterator
1668  iterator_recombine_unsigned_inner =
1669  iterator_recombine_unsigned_outer->begin();
1670  iterator_recombine_unsigned_inner !=
1671  iterator_recombine_unsigned_outer->end();
1672  iterator_recombine_unsigned_inner++) {
1673  sum_unequalized +=
1674  _event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted;
1675  }
1676  sum_unequalized = std::max(0.0, sum_unequalized);
1677 
1678  if (sum_unequalized >= sum_unequalized_3 ||
1679  (sum_unequalized >= sum_unequalized_2 &&
1680  (iterator_particle - _event.begin()) % 2 == 0) ||
1681  (sum_unequalized >= sum_unequalized_1 &&
1682  (iterator_particle - _event.begin()) % 4 == 0) ||
1683  (sum_unequalized >= sum_unequalized_0 &&
1684  (iterator_particle - _event.begin()) % 8 == 0)) {
1685 
1686  const double weight = sum_unequalized *
1687  std::min(1.0, std::max(1e-3,
1688  iterator_particle->area));
1689 
1690  if (weight > 0) {
1691  p->push_back_row(
1692  bpmpd_problem_t::greater_equal,
1693  sum_unequalized);
1694 
1695  p->push_back_coefficient(
1696  index_row, index_column_block, 1.0 / weight);
1697 
1698  for (std::vector<size_t>::const_iterator
1699  iterator_recombine_unsigned_inner =
1700  iterator_recombine_unsigned_outer->begin();
1701  iterator_recombine_unsigned_inner !=
1702  iterator_recombine_unsigned_outer->end();
1703  iterator_recombine_unsigned_inner++) {
1704  if (_event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted >= 0) {
1705  const size_t index_column_cost =
1706  nblock +
1707  positive_index[*iterator_recombine_unsigned_inner];
1708 
1709  p->push_back_coefficient(
1710  index_row, index_column_cost, 1);
1711  index_column_max =
1712  std::max(index_column_max, index_column_cost);
1713  }
1714  }
1715  index_row++;
1716 
1717  p->push_back_row(
1718  bpmpd_problem_t::greater_equal,
1719  -sum_unequalized);
1720 
1721  p->push_back_coefficient(
1722  index_row, index_column_block, _positive_bound_scale / weight);
1723 
1724  for (std::vector<size_t>::const_iterator iterator_recombine_unsigned_inner = iterator_recombine_unsigned_outer->begin();
1725  iterator_recombine_unsigned_inner != iterator_recombine_unsigned_outer->end();
1726  iterator_recombine_unsigned_inner++) {
1727  if (_event[*iterator_recombine_unsigned_inner].momentum_perp_subtracted >= 0) {
1728  const size_t index_column_cost =
1729  nblock +
1730  positive_index[*iterator_recombine_unsigned_inner];
1731 
1732  p->push_back_coefficient(
1733  index_row, index_column_cost, -1);
1734  index_column_max =
1735  std::max(index_column_max, index_column_cost);
1736  }
1737  }
1738  index_row++;
1739  }
1740 
1741  }
1742 
1743  }
1744  }
1745  }
1746 
1747  // Epsilon that breaks the degeneracy, in the same units
1748  // as the pT of the event (i.e. GeV)
1749  static const double epsilon_degeneracy = 1e-2;
1750 
1751  // Columns (variables and the objective coefficients) of
1752  // the LP problem
1753  //
1754  // Cost variables (objective coefficient 1)
1755  for (size_t i = 0; i < nsuperblock; i++) {
1756  p->push_back_column(
1757  1, 0, bpmpd_problem_t::infinity);
1758  }
1759  for (size_t i = nsuperblock; i < nsuperblock + nstaggered_block; i++) {
1760  p->push_back_column(
1761  0, 0, bpmpd_problem_t::infinity);
1762  }
1763  for (size_t i = nsuperblock + nstaggered_block; i < nsuperblock + 2 * nstaggered_block; i++) {
1764  p->push_back_column(
1765  0, 0, bpmpd_problem_t::infinity);
1766  }
1767  for (size_t i = nsuperblock + 2 * nstaggered_block; i < _ncost; i++) {
1768  p->push_back_column(
1769  0, 0, bpmpd_problem_t::infinity);
1770  }
1771  //fprintf(stderr, "%s:%d: %lu %lu\n", __FILE__, __LINE__, index_column_max, recombine_tie.size());
1772  // Energy transfer coefficients t_kl (objective
1773  // coefficient 0 + epsilon)
1774  for (size_t i = _ncost; i <= index_column_max; i++) {
1775  p->push_back_column(
1776  epsilon_degeneracy * _recombine_tie[i - _ncost],
1777  0, bpmpd_problem_t::infinity);
1778  }
1779  }
1780  void VoronoiAlgorithm::equalize(void)
1781  {
1782  bpmpd_problem_t lp_problem = reinterpret_cast<bpmpd_environment_t *>(_lp_environment)->problem();
1783 
1784  recombine_link();
1785  lp_populate(&lp_problem);
1786  lp_problem.optimize();
1787 
1788  int solution_status;
1789  double objective_value;
1790  std::vector<double> x;
1791  std::vector<double> pi;
1792 
1793  lp_problem.solve(solution_status, objective_value,
1794  x, pi);
1795 
1796  for (size_t k = _ncost; k < x.size(); k++) {
1797  if (_event[_recombine[k - _ncost].first].
1798  momentum_perp_subtracted < 0 &&
1799  _event[_recombine[k - _ncost].second].
1800  momentum_perp_subtracted >= 0 && x[k] >= 0) {
1801  _event[_recombine[k - _ncost].first].
1802  momentum_perp_subtracted += x[k];
1803  _event[_recombine[k - _ncost].second].
1804  momentum_perp_subtracted -= x[k];
1805  }
1806  }
1807  for (size_t k = 0; k < _event.size(); k++) {
1808  if (_nblock_subtract[k] != 0 &&
1809  x[_nblock_subtract[k]] >= 0) {
1810  _event[k].momentum_perp_subtracted -=
1811  x[_nblock_subtract[k]];
1812  }
1813  }
1814  }
1815  void VoronoiAlgorithm::remove_nonpositive(void)
1816  {
1817  for (std::vector<particle_t>::iterator iterator =
1818  _event.begin();
1819  iterator != _event.end(); iterator++) {
1820  iterator->momentum_perp_subtracted = std::max(
1821  0.0, iterator->momentum_perp_subtracted);
1822  }
1823  }
1824  void VoronoiAlgorithm::subtract_if_necessary(void)
1825  {
1826  if (!_subtracted) {
1827  event_fourier();
1828  feature_extract();
1829  voronoi_area_incident();
1830  subtract_momentum();
1831  if (_remove_nonpositive) {
1832  equalize();
1833  remove_nonpositive();
1834  }
1835  _subtracted = true;
1836  }
1837  }
1838 
1839  VoronoiAlgorithm::VoronoiAlgorithm(
1840  const double dr_max,
1841  const bool isRealData,
1842  const bool isCalo,
1843  const std::pair<double, double> equalization_threshold,
1844  const bool remove_nonpositive)
1845  : _remove_nonpositive(remove_nonpositive),
1846  _equalization_threshold(equalization_threshold),
1847  _radial_distance_square_max(dr_max * dr_max),
1848  _positive_bound_scale(0.2),
1849  _subtracted(false),
1850  ue(NULL)
1851  {
1852  initialize_geometry();
1853  ue = new UECalibration(isRealData, isCalo);
1854  static const size_t nedge_pseudorapidity = 15 + 1;
1855  static const double edge_pseudorapidity[nedge_pseudorapidity] = {
1856  -5.191, -2.650, -2.043, -1.740, -1.479, -1.131, -0.783, -0.522, 0.522, 0.783, 1.131, 1.479, 1.740, 2.043, 2.650, 5.191
1857  };
1858 
1859  _edge_pseudorapidity = std::vector<double>(
1860  edge_pseudorapidity,
1861  edge_pseudorapidity + nedge_pseudorapidity);
1862  allocate();
1863  }
1864 
1865  VoronoiAlgorithm::~VoronoiAlgorithm(void)
1866  {
1867  deallocate();
1868  }
1869 
1879  void VoronoiAlgorithm::push_back_particle(
1880  const double perp, const double pseudorapidity,
1881  const double azimuth,
1882  const unsigned int reduced_particle_flow_id)
1883  {
1884  math::PtEtaPhiELorentzVector p(perp, pseudorapidity, azimuth, NAN);
1885 
1886  p.SetE(p.P());
1887  _event.push_back(particle_t(p, reduced_particle_flow_id));
1888  }
1892  void VoronoiAlgorithm::clear(void)
1893  {
1894  _event.clear();
1895  _subtracted = false;
1896  }
1902  std::vector<double> VoronoiAlgorithm::subtracted_equalized_perp(void)
1903  {
1904  subtract_if_necessary();
1905 
1906  std::vector<double> ret;
1907 
1908  for (std::vector<particle_t>::const_iterator iterator =
1909  _event.begin();
1910  iterator != _event.end(); iterator++) {
1911  ret.push_back(iterator->momentum_perp_subtracted);
1912  }
1913 
1914  return ret;
1915  }
1916  std::vector<double> VoronoiAlgorithm::subtracted_unequalized_perp(void)
1917  {
1918  subtract_if_necessary();
1919 
1920  std::vector<double> ret;
1921 
1922  for (std::vector<particle_t>::const_iterator iterator =
1923  _event.begin();
1924  iterator != _event.end(); iterator++) {
1925  ret.push_back(iterator->momentum_perp_subtracted_unequalized);
1926  }
1927 
1928  return ret;
1929  }
1936  std::vector<double> VoronoiAlgorithm::particle_area(void)
1937  {
1938  subtract_if_necessary();
1939 
1940  std::vector<double> ret;
1941 
1942  for (std::vector<particle_t>::const_iterator iterator =
1943  _event.begin();
1944  iterator != _event.end(); iterator++) {
1945  ret.push_back(iterator->area);
1946  }
1947 
1948  return ret;
1949  }
1958  std::vector<std::set<size_t> > VoronoiAlgorithm::particle_incident(void)
1959  {
1960  subtract_if_necessary();
1961 
1962  std::vector<std::set<size_t> > ret;
1963 
1964  for (std::vector<particle_t>::const_iterator
1965  iterator_outer = _event.begin();
1966  iterator_outer != _event.end(); iterator_outer++) {
1967  std::set<size_t> e;
1968 
1969  for (std::set<std::vector<particle_t>::iterator>::
1970  const_iterator iterator_inner =
1971  iterator_outer->incident.begin();
1972  iterator_inner != iterator_outer->incident.begin();
1973  iterator_inner++) {
1974  e.insert(*iterator_inner - _event.begin());
1975  }
1976  ret.push_back(e);
1977  }
1978 
1979  return ret;
1980  }
1981  std::vector<double> VoronoiAlgorithm::perp_fourier(void)
1982  {
1983  subtract_if_necessary();
1984 
1985  return std::vector<double>(
1986  _perp_fourier->data(),
1987  _perp_fourier->data() +
1988  _perp_fourier->num_elements());
1989  }
1990  size_t VoronoiAlgorithm::nedge_pseudorapidity(void) const
1991  {
1992  return _edge_pseudorapidity.size();
1993  }
predc_
struct predc_s predc_
param_
struct @392 param_
predp_s::mincor
int mincor
Definition: VoronoiAlgorithm.cc:94
factor_s
Definition: VoronoiAlgorithm.cc:36
VoronoiAlgorithm::_active
std::vector< bool > _active
Definition: VoronoiAlgorithm.h:76
param_s
Definition: VoronoiAlgorithm.cc:74
VoronoiAlgorithm::_recombine_index
std::vector< std::vector< size_t > > _recombine_index
Definition: VoronoiAlgorithm.h:78
toler_
struct toler_s toler_
initv_s::stamet
int stamet
Definition: VoronoiAlgorithm.cc:51
i
int i
Definition: DBlmapReader.cc:9
sprnod_s::ssupn
int ssupn
Definition: VoronoiAlgorithm.cc:108
initv_s::premet
int premet
Definition: VoronoiAlgorithm.cc:53
drop_s::tfixvar
double tfixvar
Definition: VoronoiAlgorithm.cc:31
lumiContext.fill
string fill
Definition: lumiContext.py:319
toler_s::feas2
double feas2
Definition: VoronoiAlgorithm.cc:119
setden_
struct setden_s setden_
param_s::dalpha
double dalpha
Definition: VoronoiAlgorithm.cc:76
predc_s::center
double center
Definition: VoronoiAlgorithm.cc:83
toler_s::tfeas2
double tfeas2
Definition: VoronoiAlgorithm.cc:117
initv_s::safmet
int safmet
Definition: VoronoiAlgorithm.cc:52
toler_s::inftol
double inftol
Definition: VoronoiAlgorithm.cc:122
VoronoiAlgorithm::perp_fourier
std::vector< double > perp_fourier(void)
Definition: VoronoiAlgorithm.cc:1981
VoronoiAlgorithm::subtract_momentum
void subtract_momentum(void)
Definition: VoronoiAlgorithm.cc:1126
compl_s::ccorr
double ccorr
Definition: VoronoiAlgorithm.cc:15
sprnod_s::psupn
int psupn
Definition: VoronoiAlgorithm.cc:107
VoronoiAlgorithm::event_fourier
void event_fourier(void)
Definition: VoronoiAlgorithm.cc:942
predc_s::target
double target
Definition: VoronoiAlgorithm.cc:80
run_regression.ret
int ret
Definition: run_regression.py:388
factor_s::tabs
double tabs
Definition: VoronoiAlgorithm.cc:39
factor_s::tpiv1
double tpiv1
Definition: VoronoiAlgorithm.cc:37
VoronoiAlgorithm::particle_t
Definition: VoronoiAlgorithm.h:26
mscal_s::scfree
double scfree
Definition: VoronoiAlgorithm.cc:66
predc_s::maxcc
int maxcc
Definition: VoronoiAlgorithm.cc:86
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
sprnod_
struct sprnod_s sprnod_
predc_s::mincc
int mincc
Definition: VoronoiAlgorithm.cc:85
toler_s::dinfs
double dinfs
Definition: VoronoiAlgorithm.cc:121
predc_s::corstp
double corstp
Definition: VoronoiAlgorithm.cc:84
itref_s
Definition: VoronoiAlgorithm.cc:57
toler_s
Definition: VoronoiAlgorithm.cc:112
VoronoiAlgorithm::voronoi_diagram_t
CGAL::Voronoi_diagram_2< delaunay_triangulation_t, CGAL::Delaunay_triangulation_adaptation_traits_2< delaunay_triangulation_t >, CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2< delaunay_triangulation_t > > voronoi_diagram_t
Definition: VoronoiAlgorithm.h:57
VoronoiAlgorithm::_positive_bound_scale
double _positive_bound_scale
Definition: VoronoiAlgorithm.h:71
dims_s::cfree
int cfree
Definition: VoronoiAlgorithm.cc:24
dims_s::rfree
int rfree
Definition: VoronoiAlgorithm.cc:27
VoronoiAlgorithm::remove_nonpositive
void remove_nonpositive(void)
Definition: VoronoiAlgorithm.cc:1815
setden_s::densgap
double densgap
Definition: VoronoiAlgorithm.cc:101
NULL
#define NULL
Definition: scimark2.h:8
sprnod_s::maxsnz
int maxsnz
Definition: VoronoiAlgorithm.cc:109
dims_s::nz
int nz
Definition: VoronoiAlgorithm.cc:23
min
#define min(a, b)
Definition: mlp_lapack.h:161
VoronoiAlgorithm::_edge_pseudorapidity
std::vector< double > _edge_pseudorapidity
Definition: VoronoiAlgorithm.h:65
itref_s::tresy
double tresy
Definition: VoronoiAlgorithm.cc:59
predc_s
Definition: VoronoiAlgorithm.cc:79
VoronoiAlgorithm::initialize_geometry
void initialize_geometry(void)
Definition: VoronoiAlgorithm.cc:898
numer_s::tplus
double tplus
Definition: VoronoiAlgorithm.cc:70
factor_s::supdens
double supdens
Definition: VoronoiAlgorithm.cc:44
numer_
struct numer_s numer_
predp_s::bargrw
double bargrw
Definition: VoronoiAlgorithm.cc:92
factor_s::tfind
double tfind
Definition: VoronoiAlgorithm.cc:42
setden_s::setlam
int setlam
Definition: VoronoiAlgorithm.cc:102
dims_s
Definition: VoronoiAlgorithm.cc:18
VoronoiAlgorithm::_recombine_tie
std::vector< double > _recombine_tie
Definition: VoronoiAlgorithm.h:80
predp_s
Definition: VoronoiAlgorithm.cc:89
initv_s::upmax
double upmax
Definition: VoronoiAlgorithm.cc:49
ascal_s::rhsnor
double rhsnor
Definition: VoronoiAlgorithm.cc:7
pi
const Double_t pi
Definition: trackSplitPlot.h:35
solver_
int solver_(double *obj, double *rhs, double *lbound, double *ubound, double *diag, double *odiag, double *xs, double *dxs, double *dxsn, double *up, double *dspr, double *ddspr, double *ddsprn, double *dsup, double *ddsup, double *ddsupn, double *dv, double *ddv, double *ddvn, double *prinf, double *upinf, double *duinf, double *scale, double *nonzeros, int *vartyp, int *slktyp, int *colpnt, int *ecolpnt, int *count, int *vcstat, int *pivots, int *invprm, int *snhead, int *nodtyp, int *inta1, int *prehis, int *rowidx, int *rindex, int *code, double *opt, int *iter, int *corect, int *fixn, int *dropn, int *fnzmax, int *fnzmin, double *addobj, double *bigbou, double *big, int *ft)
ascal_s::scalmet
int scalmet
Definition: VoronoiAlgorithm.cc:10
compl_s
Definition: VoronoiAlgorithm.cc:13
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:234
setden_s
Definition: VoronoiAlgorithm.cc:99
VoronoiAlgorithm::point_2d_t
delaunay_triangulation_t::Point point_2d_t
Definition: VoronoiAlgorithm.h:25
VoronoiAlgorithm::_cms_ecal_edge_pseudorapidity
std::vector< double > _cms_ecal_edge_pseudorapidity
Definition: VoronoiAlgorithm.h:67
VoronoiAlgorithm::particle_area
std::vector< double > particle_area(void)
Definition: VoronoiAlgorithm.cc:1936
toler_s::maxiter
int maxiter
Definition: VoronoiAlgorithm.cc:123
toler_s::tsdir
double tsdir
Definition: VoronoiAlgorithm.cc:113
VoronoiAlgorithm::lp_populate
void lp_populate(void *lp_problem)
Definition: VoronoiAlgorithm.cc:1411
VoronoiAlgorithm::subtract_if_necessary
void subtract_if_necessary(void)
Definition: VoronoiAlgorithm.cc:1824
initv_s::prmin
double prmin
Definition: VoronoiAlgorithm.cc:48
VoronoiAlgorithm::push_back_particle
void push_back_particle(const double perp, const double pseudorapidity, const double azimuth, const unsigned int reduced_particle_flow_id)
Definition: VoronoiAlgorithm.cc:1879
itref_s::maxref
int maxref
Definition: VoronoiAlgorithm.cc:60
compl_
struct compl_s compl_
VoronoiAlgorithm::_recombine
std::vector< std::pair< size_t, size_t > > _recombine
Definition: VoronoiAlgorithm.h:77
numer_s
Definition: VoronoiAlgorithm.cc:69
initv_
struct initv_s initv_
itref_
struct itref_s itref_
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
pileupReCalc_HLTpaths.scale
int scale
Definition: pileupReCalc_HLTpaths.py:111
BPMPD_VERSION_STRING
#define BPMPD_VERSION_STRING
Definition: VoronoiAlgorithm.cc:465
factor_s::order
double order
Definition: VoronoiAlgorithm.cc:43
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
hitfit::clear
void clear(CLHEP::HepGenMatrix &m)
Helper function: Reset all elements of a matrix to 0.
Definition: matutil.cc:168
VoronoiAlgorithm::nreduced_particle_flow_id
static const size_t nreduced_particle_flow_id
Definition: VoronoiAlgorithm.h:62
VoronoiAlgorithm::particle_incident
std::vector< std::set< size_t > > particle_incident(void)
Definition: VoronoiAlgorithm.cc:1958
query.result
tuple result
Definition: query.py:137
predc_s::tsmall
double tsmall
Definition: VoronoiAlgorithm.cc:81
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
VoronoiAlgorithm::nfourier
static const size_t nfourier
Definition: VoronoiAlgorithm.h:63
infinity
const double infinity
Definition: CSCChamberFitter.cc:10
patZpeak.handle
tuple handle
Definition: patZpeak.py:22
j
int j
Definition: DBlmapReader.cc:9
dims_s::n1
int n1
Definition: VoronoiAlgorithm.cc:20
compl_s::climit
double climit
Definition: VoronoiAlgorithm.cc:14
VoronoiAlgorithm::_nblock_subtract
std::vector< size_t > _nblock_subtract
Definition: VoronoiAlgorithm.h:82
dims_s::denwin
int denwin
Definition: VoronoiAlgorithm.cc:26
factor_s::lam
double lam
Definition: VoronoiAlgorithm.cc:41
VoronoiAlgorithm::_feature
std::vector< double > _feature
Definition: VoronoiAlgorithm.h:75
dims_s::n
int n
Definition: VoronoiAlgorithm.cc:19
initv_s::regul
int regul
Definition: VoronoiAlgorithm.cc:54
drop_s::tfixslack
double tfixslack
Definition: VoronoiAlgorithm.cc:32
first
bool first
Definition: L1TdeRCT.cc:94
ascal_s
Definition: VoronoiAlgorithm.cc:5
initv_s
Definition: VoronoiAlgorithm.cc:47
mscal_s::varadd
double varadd
Definition: VoronoiAlgorithm.cc:64
param_s::palpha
double palpha
Definition: VoronoiAlgorithm.cc:75
drop_s
Definition: VoronoiAlgorithm.cc:30
dims_s::m
int m
Definition: VoronoiAlgorithm.cc:21
ascal_s::objnor
double objnor
Definition: VoronoiAlgorithm.cc:6
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
predc_s::tlarge
double tlarge
Definition: VoronoiAlgorithm.cc:82
toler_s::feas1
double feas1
Definition: VoronoiAlgorithm.cc:118
ascal_s::scpass
int scpass
Definition: VoronoiAlgorithm.cc:9
toler_s::topt2
double topt2
Definition: VoronoiAlgorithm.cc:115
M_PI
#define M_PI
Definition: BFit3D.cc:3
VoronoiAlgorithm::_radial_distance_square_max
double _radial_distance_square_max
Definition: VoronoiAlgorithm.h:70
setden_s::denslen
int denslen
Definition: VoronoiAlgorithm.cc:103
VoronoiAlgorithm::_perp_fourier
boost::multi_array< double, 4 > * _perp_fourier
Definition: VoronoiAlgorithm.h:74
dims_s::pivotn
int pivotn
Definition: VoronoiAlgorithm.cc:25
math::PtEtaPhiELorentzVector
PtEtaPhiELorentzVectorD PtEtaPhiELorentzVector
Lorentz vector with cartesian internal representation.
Definition: LorentzVector.h:28
dims_s::mn
int mn
Definition: VoronoiAlgorithm.cc:22
mscal_s::slkadd
double slkadd
Definition: VoronoiAlgorithm.cc:65
VoronoiAlgorithm::_recombine_unsigned
std::vector< std::vector< size_t > > _recombine_unsigned
Definition: VoronoiAlgorithm.h:79
VoronoiAlgorithm::ue
UECalibration * ue
Definition: VoronoiAlgorithm.h:86
UECalibration::ue_interpolation_pf0
float ue_interpolation_pf0[15][344]
Definition: UECalibration.h:90
predp_s::ccstop
double ccstop
Definition: VoronoiAlgorithm.cc:90
initv_s::dumin
double dumin
Definition: VoronoiAlgorithm.cc:50
VoronoiAlgorithm::feature_extract
void feature_extract(void)
Definition: VoronoiAlgorithm.cc:978
factor_s::tpiv2
double tpiv2
Definition: VoronoiAlgorithm.cc:38
mscal_
struct mscal_s mscal_
mscal_s
Definition: VoronoiAlgorithm.cc:63
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
predp_s::inibar
int inibar
Definition: VoronoiAlgorithm.cc:96
numer_s::tzer
double tzer
Definition: VoronoiAlgorithm.cc:71
perp
T perp() const
Magnitude of transverse component.
Definition: Basic3DVectorLD.h:148
VoronoiAlgorithm::voronoi_area_incident
void voronoi_area_incident(void)
Definition: VoronoiAlgorithm.cc:1035
begin
#define begin
Definition: vmac.h:31
create_public_pileup_plots.transform
transform
Definition: create_public_pileup_plots.py:254
VoronoiAlgorithm::_equalization_threshold
std::pair< double, double > _equalization_threshold
Definition: VoronoiAlgorithm.h:69
big
Definition: big.h:6
VoronoiAlgorithm::nedge_pseudorapidity
size_t nedge_pseudorapidity(void) const
Definition: VoronoiAlgorithm.cc:1990
drop_
struct drop_s drop_
VoronoiAlgorithm::VoronoiAlgorithm
VoronoiAlgorithm(const double dr_max, const bool isRealData=true, const bool isCalo=false, const std::pair< double, double > equalization_threshold=std::pair< double, double >(5.0, 35.0), const bool remove_nonpositive=true)
Definition: VoronoiAlgorithm.cc:1839
lut2db_cfg.filename
tuple filename
Definition: lut2db_cfg.py:20
itref_s::tresx
double tresx
Definition: VoronoiAlgorithm.cc:58
UECalibration::ue_predictor_pf
float ue_predictor_pf[3][15][5][2][82]
Definition: UECalibration.h:90
factor_s::trabs
double trabs
Definition: VoronoiAlgorithm.cc:40
predp_s::barset
double barset
Definition: VoronoiAlgorithm.cc:91
ascal_s::scdiff
double scdiff
Definition: VoronoiAlgorithm.cc:8
dims_
struct dims_s dims_
VoronoiAlgorithm::_cms_hcal_edge_pseudorapidity
std::vector< double > _cms_hcal_edge_pseudorapidity
Definition: VoronoiAlgorithm.h:66
VoronoiAlgorithm::polygon_t
CGAL::Polygon_2< CGAL::Exact_predicates_inexact_constructions_kernel > polygon_t
Definition: VoronoiAlgorithm.h:60
VoronoiAlgorithm::subtracted_unequalized_perp
std::vector< double > subtracted_unequalized_perp(void)
Definition: VoronoiAlgorithm.cc:1916
vdt::x
x
Definition: VDTMath.h:216
UECalibration::ue_interpolation_pf1
float ue_interpolation_pf1[15][344]
Definition: UECalibration.h:90
predp_
struct predp_s predp_
ascal_
struct ascal_s ascal_
toler_s::topt1
double topt1
Definition: VoronoiAlgorithm.cc:114
UECalibration::ue_interpolation_pf2
float ue_interpolation_pf2[15][82]
Definition: UECalibration.h:90
edm::eventsetup::heterocontainer::insert
bool insert(Storage &, ItemType *, const IdTag &)
toler_s::pinfs
double pinfs
Definition: VoronoiAlgorithm.cc:120
toler_s::tfeas1
double tfeas1
Definition: VoronoiAlgorithm.cc:116
predp_s::maxcor
int maxcor
Definition: VoronoiAlgorithm.cc:95
drop_s::slklim
double slklim
Definition: VoronoiAlgorithm.cc:33
factor_
struct factor_s factor_
setden_s::maxdense
double maxdense
Definition: VoronoiAlgorithm.cc:100
VoronoiAlgorithm::subtracted_equalized_perp
std::vector< double > subtracted_equalized_perp(void)
Definition: VoronoiAlgorithm.cc:1902
predp_s::barmin
double barmin
Definition: VoronoiAlgorithm.cc:93
constexpr
#define constexpr
Definition: GCC11Compatibility.h:29
sprnod_s
Definition: VoronoiAlgorithm.cc:106
BPMPD_INFINITY_BOUND
#define BPMPD_INFINITY_BOUND
Definition: VoronoiAlgorithm.cc:461
runtimedef::set
void set(const std::string &name, int value)
set the flag, with a run-time name
Definition: ProfilingTools.cc:67
Generated for CMSSW Reference Manual by doxygen 1.8.5