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