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