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