LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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 <algorithm> /* for sort */
26#include <climits>
27#include <cmath>
28#include <cstring> /* for memset*/
29#include <iostream>
30
31#include "lsst/log/Log.h"
33
34namespace {
35LOG_LOGGER _log = LOG_GET("lsst.jointcal.Histo4d");
36}
37
38namespace lsst {
39namespace jointcal {
40
41SparseHisto4d::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 : _data(nEntries, 0) {
45 double indexMax = n1 * n2 * n3 * n4;
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 _dataSize = nEntries;
63 _ndata = 0;
64 _sorted = false;
65}
66
67int SparseHisto4d::code_value(const double x[4]) const {
68 int index = 0;
69 for (int idim = 0; idim < 4; ++idim) {
70 int i = (int)std::floor((x[idim] - _minVal[idim]) * _scale[idim]);
71 if (i < 0 || i >= _n[idim]) return -1;
72 index = index * _n[idim] + i;
73 }
74 return index;
75}
76
77void SparseHisto4d::inverse_code(int code, double x[4]) const {
78 for (int i = 3; i >= 0; --i) {
79 int bin = code % _n[i];
80 code /= _n[i];
81 x[i] = _minVal[i] + ((double)bin + 0.5) / _scale[i];
82 }
83}
84
86 if (!_sorted) {
87 std::sort(_data.begin(), _data.end());
88 _sorted = true;
89 }
90}
91
92void SparseHisto4d::fill(const double x[4])
93
94{
95 int code = code_value(x);
96 if (code < 0) return;
97 if (_ndata == _dataSize) {
98 _dataSize *= 2;
99 _data.resize(_dataSize);
100 }
101 _data[_ndata++] = code;
102 _sorted = false;
103}
104
105void SparseHisto4d::fill(const double x1, const double x2, const double x3, const double x4) {
106 static double x[4];
107 x[0] = x1;
108 x[1] = x2;
109 x[2] = x3;
110 x[3] = x4;
111 fill(x);
112}
113
114int SparseHisto4d::maxBin(double x[4]) {
115 sort();
116 if (_ndata == 0) return 0;
117 int maxval = _data[0];
118 int maxCount = 1;
119 int oldval = _data[0];
120 int currentCount = 1;
121 for (int i = 1; i < _ndata; ++i) {
122 if (_data[i] == oldval)
123 currentCount++;
124 else
125 currentCount = 1;
126 if (currentCount > maxCount) {
127 maxCount = currentCount;
128 maxval = _data[i];
129 }
130 oldval = _data[i];
131 }
132 inverse_code(maxval, x);
133 return maxCount;
134}
135
136void SparseHisto4d::zeroBin(double x[4]) {
137 sort();
138 int code = code_value(x);
139 // inefficient locator...
140 int start = 0;
141 while ((_data[start] < code) && start < _ndata) start++;
142 // find how many identical entries we have
143 int end = std::min(start + 1, _ndata);
144 while (end < _ndata && _data[start] == _data[end]) end++;
145 int shift = end - start;
146 int lastShift = _ndata - (end - start);
147 for (int i = start; i < lastShift; ++i) _data[i] = _data[i + shift];
148 _ndata -= shift;
149}
150
151void SparseHisto4d::binLimits(const double x[4], const int iDim, double &xMin, double &xMax) const {
152 int code = code_value(x);
153 double xCenter[4];
154 inverse_code(code, xCenter);
155 xMin = xCenter[iDim] - 0.5 / _scale[iDim];
156 xMax = xCenter[iDim] + 0.5 / _scale[iDim];
157}
158
160 for (int i = 0; i < _ndata; ++i) // DEBUG
161 std::cout << _data[i] << ' ';
163}
164} // namespace jointcal
165} // namespace lsst
int end
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:659
#define LOG_GET(logger)
Returns a Log object associated with logger.
Definition Log.h:75
#define LOG_LOGGER
Definition Log.h:714
T begin(T... args)
void zeroBin(double x[4])
Definition Histo4d.cc:136
int maxBin(double x[4])
Definition Histo4d.cc:114
void inverse_code(int code, double x[4]) const
Definition Histo4d.cc:77
void fill(const double x[4])
Definition Histo4d.cc:92
void binLimits(const double x[4], int idim, double &xMin, double &xMax) const
return the bin limits of dimension idim (0<=idim<4), around point X.
Definition Histo4d.cc:151
int code_value(const double x[4]) const
Definition Histo4d.cc:67
T end(T... args)
T endl(T... args)
T floor(T... args)
T min(T... args)
T resize(T... args)
T sort(T... args)