Line data Source code
1 : /**
2 : Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
3 : Initial version January 2010.
4 : Maintained by ESO since 2013.
5 :
6 : This file is part of LibAIR and is licensed under GNU Public
7 : License Version 2
8 :
9 : \file arraygains.cpp
10 : Renamed arraygains.cc 2023
11 :
12 : Structure to hold gains derived from WVR data
13 : */
14 :
15 : #include "arraygains.h"
16 : #include "arraydata.h"
17 : #include "dtdlcoeffs.h"
18 : #include <cassert>
19 : #include <iostream>
20 :
21 : namespace LibAIR2 {
22 :
23 :
24 0 : ArrayGains::ArrayGains(const std::vector<double> &time,
25 : const std::vector<double> &el,
26 : const std::vector<size_t> &state,
27 : const std::vector<size_t> &field,
28 : const std::vector<size_t> &source,
29 0 : size_t nAnt):
30 0 : time(time),
31 0 : el(el),
32 0 : state(state), field(field), source(source),
33 0 : path(time.size(),nAnt,0.),
34 0 : nAnt(nAnt)
35 : {
36 0 : }
37 :
38 0 : void ArrayGains::calc(const InterpArrayData &wvrdata,
39 : const std::vector<double> &coeffs)
40 : {
41 0 : std::vector<double> c;
42 0 : reweight_thermal(coeffs, c);
43 :
44 0 : const size_t ntimes=wvrdata.nTimes();
45 0 : for (size_t i=0; i<ntimes; ++i)
46 : {
47 0 : for(size_t j=0; j<wvrdata.nAnts; ++j)
48 : {
49 0 : double cpath=0;
50 0 : for(size_t k=0; k<4; ++k)
51 : {
52 0 : double T = wvrdata.g_wvrdata()(i,j,k);
53 :
54 0 : if(T>0. && coeffs[k]>0)
55 : {
56 0 : cpath+=T*c[k];
57 : }
58 : }
59 0 : path(i,j)=cpath;
60 : }
61 : }
62 0 : }
63 :
64 0 : void ArrayGains::calc(const InterpArrayData &wvrdata,
65 : const std::vector<double> &coeffs,
66 : const std::vector<double> &coeffs2,
67 : const std::vector<double> &TRef)
68 : {
69 0 : std::vector<double> c, c2;
70 0 : reweight_thermal(coeffs, coeffs2, c, c2);
71 :
72 0 : const size_t ntimes=wvrdata.nTimes();
73 0 : for (size_t i=0; i<ntimes; ++i)
74 : {
75 0 : for(size_t j=0; j<wvrdata.nAnts; ++j)
76 : {
77 0 : double cpath=0;
78 0 : for(size_t k=0; k<4; ++k)
79 : {
80 0 : const double T = wvrdata.g_wvrdata()(i,j,k);
81 0 : if(T>0. && coeffs[k]>0)
82 : {
83 0 : cpath+=(T-TRef[k])*c[k]+ 0.5*std::pow(T-TRef[k], 2)*c2[k];
84 : }
85 : }
86 0 : path(i,j)=cpath;
87 : }
88 : }
89 0 : }
90 :
91 0 : void ArrayGains::calc(const InterpArrayData &wvrdata,
92 : const std::vector<double> &coeffs,
93 : const std::vector<double> &weights)
94 : {
95 0 : std::vector<double> c(4);
96 0 : for(size_t i=0; i<4; ++i)
97 0 : c[i]=weights[i]/coeffs[i];
98 :
99 0 : const size_t ntimes=wvrdata.nTimes();
100 0 : for (size_t i=0; i<ntimes; ++i)
101 : {
102 0 : for(size_t j=0; j<wvrdata.nAnts; ++j)
103 : {
104 0 : double cpath=0;
105 0 : for(size_t k=0; k<4; ++k)
106 : {
107 0 : const double T = wvrdata.g_wvrdata()(i,j,k);
108 0 : if(T>0. && coeffs[k]>0)
109 : {
110 0 : cpath+=T*c[k];
111 : }
112 : }
113 0 : path(i,j)=cpath;
114 : }
115 : }
116 0 : }
117 :
118 0 : void ArrayGains::calc(const InterpArrayData &wvrdata,
119 : const dTdLCoeffsBase &coeffs)
120 : {
121 0 : const size_t ntimes=wvrdata.nTimes();
122 :
123 : // c are the original coefficients, cw are the reweighted
124 : // coefficients
125 0 : std::vector<double> c, c2, cw;
126 :
127 0 : for (size_t i=0; i<ntimes; ++i)
128 : {
129 0 : for(size_t j=0; j<wvrdata.nAnts; ++j)
130 : {
131 0 : double cpath=0;
132 0 : coeffs.get(j,
133 0 : wvrdata.g_time()[i],
134 : M_PI/2.,
135 : c,
136 : c2);
137 0 : reweight_thermal(c, cw);
138 0 : for(size_t k=0; k<4; ++k)
139 : {
140 0 : const double T = wvrdata.g_wvrdata()(i,j,k);
141 0 : if(T>0. && c[k]!=0)
142 : {
143 : //std::cout << "i j k cw[k] T " << i << " " << j << " " << k << " " << cw[k] << " " << T << std::endl;
144 0 : cpath+=T*cw[k];
145 : }
146 : }
147 0 : path(i,j)=cpath;
148 : }
149 : }
150 0 : }
151 :
152 0 : void ArrayGains::calcLinI(const ArrayGains &o)
153 : {
154 0 : const std::vector<double> &otime=o.g_time();
155 0 : const path_t &opath=o.g_path();
156 0 : const size_t ntimes=time.size();
157 0 : const size_t notimes=otime.size();
158 :
159 0 : size_t oi=0;
160 :
161 0 : for (size_t i=0; i<ntimes; ++i)
162 : {
163 0 : while(oi < notimes and otime[oi]<time[i])
164 : {
165 0 : ++oi;
166 : }
167 :
168 0 : if (oi==0)
169 : {
170 0 : for(size_t j=0; j<nAnt; ++j)
171 0 : path[i][j]=opath(0,j);
172 : }
173 0 : else if (oi == notimes)
174 : {
175 0 : for(size_t j=0; j<nAnt; ++j)
176 0 : path(i,j)=opath(notimes-1,j);
177 : }
178 : else
179 : {
180 0 : const double tdelta=otime[oi]-otime[oi-1];
181 0 : const double w=(time[i]-otime[oi-1])/tdelta;
182 0 : for(size_t j=0; j<nAnt; ++j)
183 : {
184 0 : path(i,j)=opath(oi-1,j)*(1-w) + opath(oi,j)*w;
185 : }
186 : }
187 : }
188 0 : }
189 :
190 0 : void ArrayGains::scale(double s)
191 : {
192 : //const size_t ntimes=time.size();
193 : //for (size_t i=0; i<ntimes; ++i)
194 : //{
195 : // for(size_t j=0; j<nAnt; ++j)
196 : // {
197 : // path[i][j]=path[i][j]*s;
198 : // }
199 : //}
200 0 : path *= s;
201 0 : }
202 :
203 0 : double ArrayGains::deltaPath(size_t timei,
204 : size_t i,
205 : size_t j) const
206 : {
207 0 : return path(timei,i)-path(timei,j);
208 : }
209 :
210 0 : double ArrayGains::absPath(size_t timei,
211 : size_t i) const
212 : {
213 0 : return path(timei,i);
214 : }
215 :
216 0 : double ArrayGains::greatestRMSBl(const std::vector<std::pair<double, double> > &tmask) const
217 : {
218 0 : const size_t ntimes=time.size();
219 0 : double maxrms=0;
220 0 : for(size_t i=0; i<nAnt; ++i)
221 0 : for(size_t j=i+1; j<nAnt;++j)
222 : {
223 0 : double sx=0, sx2=0;
224 0 : size_t csegment=0;
225 0 : size_t npoints=0;
226 0 : if(tmask.size()>0)
227 : {
228 0 : for (size_t k=0; k<ntimes; ++k)
229 : {
230 0 : if (time[k]<tmask[csegment].first)
231 0 : continue;
232 0 : if (time[k]>=tmask[csegment].first && time[k]<=tmask[csegment].second)
233 : {
234 0 : if(absPath(k,i)>0. && absPath(k,j)>0.)
235 : {
236 0 : const double d=deltaPath(k, i, j);
237 0 : sx+=d;
238 0 : sx2+= (d*d);
239 0 : ++npoints;
240 : }
241 : }
242 0 : if (time[k]>=tmask[csegment].second && csegment<tmask.size()-1)
243 0 : ++csegment;
244 : }
245 : }
246 0 : if (npoints>0)
247 : {
248 0 : const double rms= pow(sx2/((double)npoints)-pow(sx/((double)npoints), 2), 0.5);
249 0 : if (rms>maxrms)
250 : {
251 0 : maxrms=rms;
252 : }
253 : }
254 : }
255 0 : return maxrms;
256 : }
257 :
258 0 : void ArrayGains::pathRMSAnt(std::vector<double> &res) const
259 : {
260 0 : std::vector<std::pair<double, double> > tmask;
261 0 : tmask.push_back(std::pair<double, double>(g_time()[0], g_time()[g_time().size()-1]));
262 0 : pathRMSAnt(tmask, res);
263 0 : }
264 :
265 0 : void ArrayGains::pathRMSAnt(const std::vector<std::pair<double, double> > &tmask,
266 : std::vector<double> &res) const
267 : {
268 :
269 0 : res.resize(nAnt);
270 :
271 0 : for(size_t i=0; i<nAnt; ++i)
272 : {
273 0 : double sx=0, sx2=0;
274 0 : size_t csegment=0;
275 0 : size_t npoints=0;
276 :
277 0 : if(tmask.size()>0)
278 : {
279 0 : for (size_t k=0; k<time.size(); ++k)
280 : {
281 0 : if(time[k]<tmask[csegment].first)
282 0 : continue;
283 :
284 0 : if (time[k]>=tmask[csegment].first && time[k]<=tmask[csegment].second)
285 : {
286 0 : if(absPath(k, i)>0.)
287 : {
288 0 : const double d=absPath(k, i) * std::sin(el[k]);
289 : //std::cout << "i k d " << i << " " << k << " " << d << std::endl;
290 0 : sx+=d;
291 0 : sx2+= (d*d);
292 0 : ++npoints;
293 : }
294 : }
295 0 : if (time[k]>=tmask[csegment].second && csegment<tmask.size()-1)
296 0 : ++csegment;
297 : }
298 : }
299 0 : if(npoints>0)
300 : {
301 0 : double dnpoints = npoints;
302 0 : const double rms= pow(sx2/dnpoints - pow(sx/dnpoints, 2), 0.5);
303 : //std::cout << "sx2 dnpoints sx rms " << sx2 << " " << dnpoints << " " << sx << " " << rms << std::endl;
304 0 : res[i]=rms;
305 : }
306 : else
307 : {
308 0 : res[i]=0.;
309 : }
310 : }
311 0 : }
312 :
313 0 : void ArrayGains::pathDiscAnt(const ArrayGains &other,
314 : std::vector<double> &res) const
315 : {
316 0 : res.resize(nAnt);
317 0 : size_t ntimes=time.size();
318 0 : for(size_t i=0; i<nAnt; ++i)
319 : {
320 0 : double sx=0, sx2=0;
321 0 : double np = 0;
322 0 : for (size_t k=0; k<ntimes; ++k)
323 : {
324 0 : if(absPath(k, i)>0. && other.absPath(k,i)>0.)
325 : {
326 0 : const double d=(absPath(k, i)-other.absPath(k,i));
327 0 : sx+=d;
328 0 : sx2+= (d*d);
329 0 : np+=1.;
330 : }
331 : }
332 0 : if(np>0)
333 : {
334 0 : const double rms= pow(sx2/np-pow(sx/np, 2), 0.5);
335 0 : res[i]=rms;
336 : }
337 : else
338 : {
339 0 : res[i] = 0.;
340 : }
341 : }
342 0 : }
343 :
344 0 : void ArrayGains::pathDiscAnt(const ArrayGains &other,
345 : const std::vector<std::pair<double, double> > &tmask,
346 : std::vector<double> &res) const
347 : {
348 0 : res.resize(nAnt);
349 0 : const size_t ntimes=time.size();
350 : //std::cerr << ntimes << std::endl;
351 0 : for(size_t i=0; i<nAnt; ++i)
352 : {
353 0 : double sx=0, sx2=0;
354 0 : size_t np=0;
355 0 : size_t ns=0;
356 0 : if (tmask.size()>0)
357 : {
358 0 : for (size_t k=0; k<ntimes; ++k)
359 : {
360 0 : if (time[k]<tmask[ns].first)
361 0 : continue;
362 0 : if (time[k]>=tmask[ns].first && time[k]<=tmask[ns].second)
363 : {
364 : //std::cerr << "k i absPath(k, i) other.absPath(k,i) " << k << " " << i << " " << absPath(k, i) << " " << other.absPath(k,i) << std:: endl;
365 0 : if(absPath(k, i)>0. && other.absPath(k,i)>0.)
366 : {
367 0 : const double d=(absPath(k, i)-other.absPath(k,i));
368 0 : sx+=d;
369 0 : sx2+= (d*d);
370 0 : ++np;
371 : }
372 : }
373 0 : if (time[k]>=tmask[ns].second && ns<tmask.size()-1)
374 0 : ++ns;
375 : }
376 : }
377 0 : res[i]=0.;
378 0 : if(np>0){
379 0 : double dnp = np;
380 0 : double rms2 = sx2/dnp - pow(sx/dnp, 2);
381 0 : if(rms2>0.){
382 0 : double rms= pow(rms2, 0.5); //the standard deviation
383 0 : res[i]=rms;
384 : }
385 : }
386 : }
387 0 : }
388 :
389 0 : void ArrayGains::blankSources(std::set<size_t> &flagsrc)
390 : {
391 0 : const size_t ntimes=time.size();
392 0 : for (size_t i=0; i<ntimes; ++i)
393 : {
394 0 : if (flagsrc.count(source[i]))
395 : {
396 0 : for(size_t j=0; j<nAnt; ++j)
397 : {
398 0 : path(i,j)=0;
399 : }
400 : }
401 : }
402 0 : }
403 :
404 0 : void ArrayGains::blankTimes(std::vector<double> &blanktimes)
405 : {
406 0 : const size_t ntimes=time.size();
407 0 : size_t k = 0;
408 0 : for (size_t i=0; i<ntimes; ++i){
409 0 : if (time[i]==blanktimes[k]){
410 0 : for(size_t j=0; j<nAnt; ++j){
411 0 : path(i,j)=0;
412 : }
413 0 : k++;
414 0 : if(k>=blanktimes.size()){
415 0 : break;
416 : }
417 : }
418 : }
419 0 : }
420 :
421 :
422 0 : void reweight(const std::vector<double> &coeffs,
423 : std::vector<double> &res )
424 : {
425 0 : const size_t n=coeffs.size();
426 0 : res.resize(n);
427 :
428 0 : double sum=0;
429 0 : for(size_t i=0; i<n; ++i)
430 : {
431 0 : if (coeffs[i] != 0)
432 : {
433 0 : res[i]=1.0/coeffs[i];
434 0 : sum+=1.0;
435 : }
436 : else
437 : {
438 0 : res[i]=0;
439 : }
440 : }
441 :
442 0 : for(size_t i=0; i<n; ++i)
443 : {
444 0 : res[i]/=sum;
445 : }
446 0 : }
447 :
448 0 : void reweight_thermal(const std::vector<double> &coeffs,
449 : std::vector<double> &res)
450 : {
451 0 : std::vector<double> c2(4,0.0);
452 0 : std::vector<double> res2(4,0.0);
453 0 : reweight_thermal(coeffs, c2, res, res2);
454 0 : }
455 :
456 0 : void reweight_thermal(const std::vector<double> &coeffs,
457 : const std::vector<double> &c2,
458 : std::vector<double> &res,
459 : std::vector<double> &res2)
460 : {
461 0 : std::array<double, 4> noise = { 0.1, 0.08, 0.08, 0.09 };
462 0 : const size_t n=coeffs.size(); // should be <=4!
463 0 : assert(n<=4);
464 :
465 0 : res.resize(n);
466 0 : res2.resize(n);
467 :
468 0 : double sum=0;
469 0 : for(size_t i=0; i<n; ++i)
470 : {
471 0 : if (coeffs[i] != 0)
472 : {
473 0 : const double dLdT = 1.0/coeffs[i];
474 0 : const double w=std::pow(noise[i] * dLdT, -2);
475 0 : const double d2LdT2 = -1.0 * c2[i] * std::pow(coeffs[i], -3);
476 0 : res[i]=w*dLdT;
477 0 : res2[i]=w*d2LdT2;
478 0 : sum+=w;
479 : }
480 : else
481 : {
482 0 : res[i]=0;
483 0 : res2[i]=0;
484 : }
485 : }
486 :
487 0 : for(size_t i=0; i<n; ++i)
488 : {
489 0 : if (sum>0.)
490 : {
491 0 : res[i]/=sum;
492 0 : res2[i]/=sum;
493 : }
494 0 : else if (res[i]>0. || res2[i]>0.)
495 : {
496 0 : std::cerr << "Cannot perform thermal reweighting: sum of all weights is zero." << std::endl;
497 : }
498 : }
499 0 : }
500 :
501 0 : double thermal_error(const std::vector<double> &coeffs)
502 : {
503 0 : std::vector<double> rwc;
504 0 : reweight_thermal(coeffs, rwc);
505 0 : std::array<double, 4> noise = {0.1, 0.08, 0.08, 0.09};
506 0 : double sum=0;
507 0 : for(size_t i=0; i<4; ++i)
508 : {
509 0 : if (rwc[i] != 0)
510 : {
511 0 : sum+=std::pow( noise[i]*rwc[i], 2);
512 : }
513 : }
514 0 : return pow(sum, 0.5);
515 0 : }
516 :
517 :
518 : }
519 :
|
