LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
Histo4d.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 /*
3  * This file is part of jointcal.
4  *
5  * Developed for the LSST Data Management System.
6  * This product includes software developed by the LSST Project
7  * (https://www.lsst.org).
8  * See the COPYRIGHT file at the top-level directory of this distribution
9  * for details of code ownership.
10  *
11  * This program is free software: you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation, either version 3 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with this program. If not, see <https://www.gnu.org/licenses/>.
23  */
24 
25 #include <iostream>
26 #include <cmath>
27 #include <string.h> /* for memset*/
28 #include <algorithm> /* for sort */
29 #include <limits.h>
30 
31 #include "lsst/log/Log.h"
32 #include "lsst/jointcal/Histo4d.h"
33 
34 namespace {
35 LOG_LOGGER _log = LOG_GET("jointcal.Histo4d");
36 }
37 
38 namespace lsst {
39 namespace jointcal {
40 
41 SparseHisto4d::SparseHisto4d(const int n1, double min1, double max1, const int n2, double min2, double max2,
42  const int n3, double min3, double max3, const int n4, double min4, double max4,
43  const int nEntries) {
44  double indexMax = n1 * n2 * n3 * n4;
45  _data.reset();
46  if (indexMax > double(INT_MAX))
47  LOGLS_WARN(_log, "Cannot hold a 4D histo with more than " << INT_MAX << " values.");
48  _n[0] = n1;
49  _n[1] = n2;
50  _n[2] = n3;
51  _n[3] = n4;
52  _minVal[0] = min1;
53  _minVal[1] = min2;
54  _minVal[2] = min3;
55  _minVal[3] = min4;
56  _maxVal[0] = max1;
57  _maxVal[1] = max2;
58  _maxVal[2] = max3;
59  _maxVal[3] = max4;
60 
61  for (int i = 0; i < 4; ++i) _scale[i] = _n[i] / (_maxVal[i] - _minVal[i]);
62  _data.reset(new int[nEntries]);
63  _dataSize = nEntries;
64  _ndata = 0;
65  _sorted = false;
66 }
67 
68 int SparseHisto4d::code_value(const double x[4]) const {
69  int index = 0;
70  for (int idim = 0; idim < 4; ++idim) {
71  int i = (int)std::floor((x[idim] - _minVal[idim]) * _scale[idim]);
72  if (i < 0 || i >= _n[idim]) return -1;
73  index = index * _n[idim] + i;
74  }
75  return index;
76 }
77 
78 void SparseHisto4d::inverse_code(int code, double x[4]) const {
79  for (int i = 3; i >= 0; --i) {
80  int bin = code % _n[i];
81  code /= _n[i];
82  x[i] = _minVal[i] + ((double)bin + 0.5) / _scale[i];
83  }
84 }
85 
87  if (!_sorted) {
88  std::sort(_data.get(), _data.get() + _ndata);
89  _sorted = true;
90  }
91 }
92 
93 void SparseHisto4d::fill(const double x[4])
94 
95 {
96  int code = code_value(x);
97  if (code < 0) return;
98  if (_ndata == _dataSize) {
99  std::unique_ptr<int[]> newData(new int[_dataSize * 2]);
100  memcpy(newData.get(), _data.get(), _dataSize * sizeof(_data[0]));
101  _data.swap(newData);
102  _dataSize *= 2;
103  }
104  _data[_ndata++] = code;
105  _sorted = false;
106 }
107 
108 void SparseHisto4d::fill(const double x1, const double x2, const double x3, const double x4) {
109  static double x[4];
110  x[0] = x1;
111  x[1] = x2;
112  x[2] = x3;
113  x[3] = x4;
114  fill(x);
115 }
116 
117 int SparseHisto4d::maxBin(double x[4]) {
118  sort();
119  if (_ndata == 0) return 0;
120  int maxval = _data[0];
121  int maxCount = 1;
122  int oldval = _data[0];
123  int currentCount = 1;
124  for (int i = 1; i < _ndata; ++i) {
125  if (_data[i] == oldval)
126  currentCount++;
127  else
128  currentCount = 1;
129  if (currentCount > maxCount) {
130  maxCount = currentCount;
131  maxval = _data[i];
132  }
133  oldval = _data[i];
134  }
135  inverse_code(maxval, x);
136  return maxCount;
137 }
138 
139 void SparseHisto4d::zeroBin(double x[4]) {
140  sort();
141  int code = code_value(x);
142  // inefficient locator...
143  int start = 0;
144  while ((_data[start] < code) && start < _ndata) start++;
145  // find how many identical entries we have
146  int end = std::min(start + 1, _ndata);
147  while (end < _ndata && _data[start] == _data[end]) end++;
148  int shift = end - start;
149  int lastShift = _ndata - (end - start);
150  for (int i = start; i < lastShift; ++i) _data[i] = _data[i + shift];
151  _ndata -= shift;
152 }
153 
154 void SparseHisto4d::binLimits(const double x[4], const int iDim, double &xMin, double &xMax) const {
155  int code = code_value(x);
156  double xCenter[4];
157  inverse_code(code, xCenter);
158  xMin = xCenter[iDim] - 0.5 / _scale[iDim];
159  xMax = xCenter[iDim] + 0.5 / _scale[iDim];
160 }
161 
162 void SparseHisto4d::print() const {
163  for (int i = 0; i < _ndata; ++i) // DEBUG
164  std::cout << _data[i] << ' ';
165  std::cout << std::endl;
166 }
167 } // namespace jointcal
168 } // namespace lsst
int end
double x
LSST DM logging module built on log4cxx.
#define LOGLS_WARN(logger, message)
Log a warn-level message using an iostream-based interface.
Definition: Log.h:648
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition: Log.h:75
#define LOG_LOGGER
Definition: Log.h:703
void zeroBin(double x[4])
Definition: Histo4d.cc:139
void inverse_code(const int code, double x[4]) const
Definition: Histo4d.cc:78
int maxBin(double x[4])
Definition: Histo4d.cc:117
void binLimits(const double x[4], const int idim, double &xMin, double &xMax) const
return the bin limits of dimension idim (0<=idim<4), around point X.
Definition: Histo4d.cc:154
void fill(const double x[4])
Definition: Histo4d.cc:93
int code_value(const double x[4]) const
Definition: Histo4d.cc:68
T endl(T... args)
T floor(T... args)
T get(T... args)
T memcpy(T... args)
T min(T... args)
A base class for image defects.
T reset(T... args)
T sort(T... args)
T swap(T... args)