Line data Source code
1 : /**
2 : Bojan Nikolic <bojan@bnikolic.co.uk>
3 : Initial version 2009
4 :
5 : This file is part of BNMin1 and is licensed under GNU General
6 : Public License version 2
7 :
8 : \file mcpoint.cxx
9 : Renamed to mcpoint.cc 2023
10 : */
11 :
12 : #include <cmath>
13 : #include <algorithm>
14 :
15 : #include <gsl/gsl_math.h>
16 : #include <gsl/gsl_eigen.h>
17 :
18 : #include "mcpoint.h"
19 : #include "bnmin_main.h"
20 :
21 : namespace Minim {
22 :
23 771024 : MCPoint::MCPoint(void):
24 771024 : p(0),
25 771024 : ll(-9999),
26 771024 : fval(0)
27 : {
28 771024 : }
29 :
30 16200 : MCPoint::MCPoint(const std::vector<double> &p):
31 16200 : p(p),
32 16200 : ll(-9999),
33 16200 : fval(0)
34 : {
35 16200 : }
36 :
37 81 : MCPoint::MCPoint(size_t np):
38 81 : p(np),
39 81 : ll(-9999),
40 81 : fval(0)
41 : {
42 81 : }
43 :
44 3116480 : MCPoint::MCPoint(const MCPoint &other):
45 3116480 : p(other.p),
46 3116480 : ll(other.ll),
47 3116480 : fval(other.fval)
48 : {
49 3116480 : }
50 :
51 0 : MCPoint &MCPoint::operator=(const MCPoint &other)
52 : {
53 0 : p=other.p;
54 0 : ll=other.ll;
55 0 : fval=other.fval;
56 0 : return *this;
57 : }
58 :
59 77 : void moment1(const std::list<WPPoint> &l,
60 : std::vector<double> &res)
61 : {
62 77 : const size_t n=l.begin()->p.size();
63 77 : res=std::vector<double>(n, 0.0);
64 77 : for(std::list<WPPoint>::const_iterator i=l.begin();
65 770077 : i!= l.end();
66 770000 : ++i)
67 : {
68 3080000 : for (size_t j=0; j<n; ++j)
69 : {
70 2310000 : res[j]+= (i->p[j] * i->w * exp(- i->ll));
71 : }
72 : }
73 77 : }
74 :
75 77 : void moment1(const std::list<WPPoint> &l,
76 : double Z,
77 : std::vector<double> &res)
78 : {
79 77 : moment1(l,res);
80 308 : for(size_t j=0; j<res.size(); ++j)
81 231 : res[j] /= Z;
82 77 : }
83 :
84 77 : void moment2(const std::list<WPPoint> &l,
85 : const std::vector<double> &m1,
86 : std::vector<double> &res)
87 : {
88 77 : const size_t n=m1.size();
89 77 : res=std::vector<double>(n, 0.0);
90 77 : for(std::list<WPPoint>::const_iterator i=l.begin();
91 770077 : i!= l.end();
92 770000 : ++i)
93 : {
94 3080000 : for (size_t j=0; j<n; ++j)
95 : {
96 2310000 : res[j]+= ( pow(i->p[j]-m1[j],2.0) * i->w * exp(- i->ll));
97 : }
98 : }
99 77 : }
100 :
101 77 : void moment2(const std::list<WPPoint> &l,
102 : const std::vector<double> &m1,
103 : double Z,
104 : std::vector<double> &res)
105 : {
106 77 : moment2(l, m1, res);
107 308 : for(size_t j=0; j<res.size(); ++j)
108 231 : res[j] /= Z;
109 77 : }
110 :
111 771024 : void moment1(const std::set<MCPoint> &s,
112 : std::vector<double> &res)
113 : {
114 771024 : const size_t n=s.begin()->p.size();
115 771024 : res=std::vector<double>(n, 0.0);
116 :
117 771024 : size_t N=0;
118 771024 : for(std::set<MCPoint>::const_iterator i=s.begin();
119 154975824 : i!= s.end();
120 154204800 : ++i)
121 : {
122 154204800 : if(i->p.size() != n)
123 : {
124 0 : std::string errmsg = "In function moment1 expected "+std::to_string(n)+" but received "+std::to_string(i->p.size())+" pars ";
125 0 : throw BaseErr(errmsg); // NParsErr
126 0 : }
127 616819200 : for (size_t j=0; j<n; ++j)
128 : {
129 462614400 : res[j]+= (i->p[j]);
130 : }
131 154204800 : ++N;
132 : }
133 :
134 3084096 : for(size_t j=0; j<res.size(); ++j)
135 : {
136 2313072 : res[j]/=N;
137 : }
138 771024 : }
139 :
140 0 : void moment2(const std::set<MCPoint> &s,
141 : const std::vector<double> &m1,
142 : std::vector<double> &res)
143 : {
144 0 : const size_t n=m1.size();
145 0 : res=std::vector<double>(n, 0.0);
146 :
147 0 : size_t N=0;
148 0 : for(std::set<MCPoint>::const_iterator i=s.begin();
149 0 : i!= s.end();
150 0 : ++i)
151 : {
152 0 : for (size_t j=0; j<n; ++j)
153 : {
154 0 : res[j]+= pow(i->p[j]-m1[j], 2);
155 : }
156 0 : ++N;
157 : }
158 :
159 0 : for(size_t j=0; j<res.size(); ++j)
160 : {
161 0 : res[j]/=N;
162 : }
163 0 : }
164 :
165 771024 : void omoment2(const std::set<MCPoint> &s,
166 : const std::vector<double> &m1,
167 : std::vector<double> &res)
168 : {
169 771024 : const size_t n=m1.size();
170 771024 : res=std::vector<double>(n*n, 0.0);
171 :
172 771024 : size_t N=0;
173 771024 : for(std::set<MCPoint>::const_iterator i=s.begin();
174 154975824 : i!= s.end();
175 154204800 : ++i)
176 : {
177 616819200 : for (size_t j=0; j<n; ++j)
178 : {
179 1850457600 : for(size_t k=0; k<n; ++k)
180 : {
181 1387843200 : res[j*n+k] += (i->p[j]-m1[j])*(i->p[k]-m1[k]);
182 : }
183 : }
184 154204800 : ++N;
185 : }
186 :
187 7710240 : for(size_t j=0; j<res.size(); ++j)
188 : {
189 6939216 : res[j]/=N;
190 : }
191 :
192 771024 : }
193 771024 : void omoment2(const std::set<MCPoint> &s,
194 : std::vector<double> &res)
195 : {
196 771024 : std::vector<double> m1;
197 771024 : moment1(s, m1);
198 771024 : omoment2(s, m1, res);
199 771024 : }
200 :
201 0 : void StdDev(const std::set<MCPoint> &s,
202 : std::vector<double> &res)
203 : {
204 0 : std::vector<double> m1, m2;
205 0 : moment1(s, m1);
206 0 : moment2(s, m1, m2);
207 0 : res.resize(m2.size());
208 0 : for(size_t j=0; j<res.size(); ++j)
209 : {
210 0 : res[j]=pow(m2[j],0.5);
211 : }
212 0 : }
213 :
214 771024 : void principalCV(const std::vector<double> &cv,
215 : std::vector<double> &eigvals,
216 : std::vector<double> &eigvects)
217 : {
218 771024 : const size_t n=sqrt(cv.size());
219 : gsl_matrix_view m
220 771024 : = gsl_matrix_view_array (const_cast<double*>(&cv[0]), n, n);
221 :
222 771024 : gsl_vector *eval = gsl_vector_alloc (n);
223 771024 : gsl_matrix *evec = gsl_matrix_alloc (n, n);
224 :
225 : gsl_eigen_symmv_workspace * w =
226 771024 : gsl_eigen_symmv_alloc (n);
227 :
228 771024 : gsl_eigen_symmv (&m.matrix,
229 : eval,
230 : evec,
231 : w);
232 :
233 771024 : gsl_eigen_symmv_free (w);
234 :
235 771024 : gsl_eigen_symmv_sort (eval,
236 : evec,
237 : GSL_EIGEN_SORT_ABS_ASC);
238 :
239 771024 : eigvals.resize(n);
240 771024 : eigvects.resize(n*n);
241 3084096 : for(size_t j=0; j<n; ++j)
242 : {
243 2313072 : eigvals[j]=gsl_vector_get (eval, j);
244 9252288 : for(size_t i=0; i<n; ++i)
245 : {
246 6939216 : eigvects[j*n+i]= gsl_matrix_get(evec, i,j);
247 : }
248 : }
249 :
250 771024 : gsl_vector_free (eval);
251 771024 : gsl_matrix_free (evec);
252 771024 : }
253 :
254 0 : void postHist(const std::list<WPPoint> &l,
255 : double /*Z*/,
256 : const std::vector<double> &low,
257 : const std::vector<double> &high,
258 : size_t nbins,
259 : std::vector<double> &res)
260 : {
261 0 : const size_t ndim=low.size();
262 :
263 0 : res.resize(pow(nbins,ndim));
264 0 : std::fill(res.begin(), res.end(), 0.0);
265 :
266 :
267 0 : std::vector<double> deltas(ndim);
268 0 : for(size_t i=0; i<ndim; ++i)
269 : {
270 0 : deltas[i]=(high[i]-low[i])/nbins;
271 : }
272 :
273 0 : for(std::list<WPPoint>::const_iterator i=l.begin();
274 0 : i!= l.end();
275 0 : ++i)
276 : {
277 0 : bool inside=true;
278 0 : size_t k=0;
279 0 : for (size_t j=0; j<ndim; ++j)
280 : {
281 0 : int dimi = int((i->p[j]-low[j])/deltas[j]);
282 0 : if (dimi >= 0 and dimi < (int)nbins)
283 : {
284 0 : k+= dimi * pow(nbins, ndim-j-1);
285 : }
286 : else
287 : {
288 0 : inside=false;
289 : }
290 : }
291 0 : if (inside)
292 : {
293 0 : res[k]+= i->w * exp(- i->ll);
294 : }
295 : }
296 0 : }
297 :
298 0 : void marginHist(const std::list<WPPoint> &l,
299 : size_t pi,
300 : double Z,
301 : double low,
302 : double high,
303 : size_t nbins,
304 : std::vector<double> &res)
305 : {
306 0 : res.resize(nbins);
307 0 : std::fill(res.begin(), res.end(),
308 0 : 0.0);
309 :
310 0 : const double d=(high-low)/nbins;
311 0 : for(std::list<WPPoint>::const_iterator i=l.begin();
312 0 : i!= l.end();
313 0 : ++i)
314 : {
315 0 : int k=int((i->p[pi]-low)/d);
316 0 : if (k > 0 and k < (int)nbins)
317 : {
318 0 : res[k]+= i->w * exp(- i->ll);
319 : }
320 : }
321 :
322 0 : for(size_t i=0; i<res.size(); ++i)
323 : {
324 0 : res[i]/=Z;
325 : }
326 0 : }
327 :
328 0 : void marginHist2D(const std::list<WPPoint> &l,
329 : double /*Z*/,
330 : size_t i,
331 : double ilow,
332 : double ihigh,
333 : size_t j,
334 : double jlow,
335 : double jhigh,
336 : size_t nbins,
337 : std::vector<double> &res)
338 : {
339 : // Two dimensions only
340 0 : res.resize(pow(nbins,2));
341 0 : std::fill(res.begin(), res.end(),
342 0 : 0.0);
343 0 : const double idelta=(ihigh-ilow)/nbins;
344 0 : const double jdelta=(jhigh-jlow)/nbins;
345 :
346 0 : for(std::list<WPPoint>::const_iterator p=l.begin();
347 0 : p!= l.end();
348 0 : ++p)
349 : {
350 :
351 0 : int dimi = int((p->p[i]-ilow)/idelta);
352 0 : int dimj = int((p->p[j]-jlow)/jdelta);
353 :
354 0 : if (dimi >= 0 and dimi<((int)nbins) and dimj >= 0 and dimj < ((int)nbins))
355 : {
356 0 : const size_t k= dimi*nbins + dimj;
357 0 : res[k]+= p->w * exp(- p->ll);
358 : }
359 :
360 : }
361 0 : }
362 : }
363 :
364 :
|
