Line data Source code
1 : /**
2 : Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
3 : Initial version November 2009
4 : Maintained by ESO since 2013.
5 :
6 : Renamed almaabs.cc 2023
7 :
8 : */
9 :
10 : #include "almaabs.h"
11 :
12 : #include <iostream>
13 : #include <stdio.h> /*** for sprintf(...) ***/
14 : #include <cassert>
15 :
16 : #include "almaabs_i.h"
17 : #include "almaopts.h"
18 : #include "arraydata.h"
19 : #include "dtdlcoeffs.h"
20 : #include "segmentation.h"
21 : #include "../radiometermeasure.h"
22 : #include "../dtdltools.h"
23 :
24 : #include "../model_make.h"
25 : #include "bnmin1/src/nestedsampler.h"
26 : #include "bnmin1/src/nestedinitial.h"
27 :
28 :
29 : #ifdef ELLIPSESAMPLE
30 : #include "bnmin1/src/prior_sampler.h"
31 : #include "bnmin1/src/mcmonitor.h"
32 : #endif
33 :
34 : namespace LibAIR2 {
35 :
36 0 : ALMAAbsRet::ALMAAbsRet(const std::vector<double> &TObs,
37 : double el,
38 : const ALMAWVRCharacter &WVRChar,
39 0 : unsigned rseed):
40 0 : i(new iALMAAbsRet(TObs,
41 : el,
42 : WVRChar,
43 0 : rseed)),
44 0 : valid(true)
45 : {
46 0 : if( ! i->sample()){
47 0 : valid = false;
48 : }
49 0 : }
50 :
51 0 : ALMAAbsRet::~ALMAAbsRet()
52 : {
53 0 : }
54 :
55 0 : bool ALMAAbsRet::g_Res(ALMAResBase &res)
56 : {
57 0 : if(valid){
58 0 : i->g_Pars(res);
59 0 : i->g_Coeffs(res);
60 0 : return true;
61 : }
62 : else{
63 0 : return false;
64 : }
65 : }
66 :
67 0 : std::ostream & operator<<(std::ostream &os,
68 : const ALMAAbsInput &i)
69 : {
70 : char tstr[120];
71 0 : sprintf(tstr, "[%5.2f, %5.2f, %5.2f, %5.2f]", i.TObs[0], i.TObs[1], i.TObs[2], i.TObs[3]);
72 :
73 0 : os<<i.time<<"\t"
74 0 : <<i.antno<<"\t"
75 : <<tstr
76 0 : <<"\t"
77 0 : <<i.el
78 0 : <<"\t"
79 0 : <<i.state;
80 0 : return os;
81 : }
82 :
83 0 : std::ostream & operator<<(std::ostream &os,
84 : const ALMAAbsInpL &i)
85 : {
86 0 : os<<" Input WVR data information: "<<std::endl
87 0 : <<"-----------------------------------"<<std::endl;
88 : os<<"Timestamp"<<"\t"
89 : <<"Ant. #"<<"\t"
90 : <<"["<<"Chn 0"<<", "
91 : <<"Chn 1"<<", "
92 : <<"Chn 2"<<", "
93 : <<"Chn 3"<<"]"
94 : <<"\t"
95 : <<"Elevation"
96 : <<"\t"
97 0 : <<"STATE_ID"<<std::endl;
98 0 : for(const ALMAAbsInput &x: i)
99 : {
100 0 : os<<x<<std::endl;
101 : }
102 0 : return os;
103 : }
104 :
105 : /** \brief Check that the supplied temperatures are reasonable
106 : */
107 0 : static void checkTObs(const std::vector<double> &TObs)
108 : {
109 0 : int above = 0;
110 0 : int below = 0;
111 0 : for(size_t i=0; i<TObs.size(); ++i)
112 : {
113 0 : if (TObs[i]<2.7){
114 0 : below++;
115 : }
116 0 : else if (TObs[i]>350.){
117 0 : above++;
118 : }
119 : }
120 0 : if (above>0 || below>0){
121 : char tstr[120];
122 0 : sprintf(tstr, "Values out of range: Out of %d TObs values, %d were below 2.7 K and %d were above 350 K.",
123 0 : (int)TObs.size(), below, above);
124 0 : throw std::runtime_error(tstr);
125 : }
126 0 : }
127 :
128 0 : static void getMidPointData(const InterpArrayData &d,
129 : int refant,
130 : std::vector<double> &TObs,
131 : double &el,
132 : double &time,
133 : size_t &state)
134 : {
135 0 : TObs.resize(4);
136 0 : const size_t midpoint=static_cast<size_t>(d.g_time().size() *0.5);
137 0 : for(size_t k=0; k<4; ++k)
138 : {
139 0 : TObs[k]=d.g_wvrdata()(midpoint, refant, k);
140 : }
141 0 : checkTObs(TObs);
142 0 : el=d.g_el()[midpoint];
143 0 : time=d.g_time()[midpoint];
144 0 : state=d.g_state()[midpoint];
145 0 : }
146 :
147 : // Convert to units of K/meter
148 0 : static void convertKm(std::vector<double> &dTdL)
149 : {
150 0 : for(size_t i=0; i<4; ++i)
151 : {
152 0 : dTdL[i]*=1e3;
153 : }
154 0 : }
155 0 : ALMAAbsInpL SimpleSingleI(const InterpArrayData &d,
156 : int refant)
157 : {
158 0 : ALMAAbsInpL res;
159 : ALMAAbsInput a;
160 :
161 0 : std::vector<double> TObs(4);
162 : double el, time;
163 : size_t state;
164 0 : getMidPointData(d, refant,
165 : TObs,
166 : el, time,
167 : state);
168 :
169 0 : a.antno=refant;
170 0 : for(size_t i=0; i<4; ++i)
171 0 : a.TObs[i]=TObs[i];
172 0 : a.el=el;
173 0 : a.time=time;
174 0 : a.state=state;
175 0 : res.push_back(a);
176 0 : return res;
177 0 : }
178 :
179 0 : ALMAAbsInpL MultipleUniformI(const InterpArrayData &d,
180 : size_t n,
181 : const std::set<size_t> &states,
182 : int refant)
183 : {
184 0 : ALMAAbsInpL res;
185 0 : const size_t nrows=d.g_time().size();
186 0 : const size_t ndelta=nrows / (n+1);
187 :
188 0 : for(size_t i=0; i<n; ++i)
189 : {
190 : ALMAAbsInput a;
191 0 : size_t row=ndelta*(i+1);
192 0 : while (row<(nrows-1) and states.count(d.g_state()[row])==0)
193 : {
194 0 : ++row;
195 : }
196 :
197 0 : if (states.count(d.g_state()[row]) == 0)
198 : {
199 0 : throw std::runtime_error("Could not find a row with a sky state");
200 : }
201 0 : a.antno=refant;
202 0 : for(size_t k=0; k<4; ++k)
203 0 : a.TObs[k]=d.g_wvrdata()(row, refant, k);
204 0 : a.el=d.g_el()[row];
205 0 : a.time=d.g_time()[row];
206 0 : a.state=d.g_state()[row];
207 :
208 0 : res.push_back(a);
209 : }
210 0 : return res;
211 0 : }
212 :
213 0 : ALMAAbsInpL FieldMidPointI(const InterpArrayData &d,
214 : const std::vector<std::pair<double, double> > &fb,
215 : const std::set<size_t> &states,
216 : int refant)
217 : {
218 0 : ALMAAbsInpL res;
219 :
220 0 : const size_t nrows=d.g_time().size();
221 0 : size_t row=0;
222 0 : for(size_t i=0; i<fb.size(); ++i)
223 : {
224 : ALMAAbsInput a;
225 0 : const double ftime=fb[i].first;
226 0 : const double ltime=fb[i].second;
227 0 : const double mtime=(ftime+ltime)/2;
228 :
229 0 : while (row<(nrows-1) and d.g_time()[row] < mtime)
230 : {
231 0 : ++row;
232 : }
233 :
234 : // Now find the next row with a good state
235 0 : while (row<(nrows-1) and states.count(d.g_state()[row])==0)
236 : {
237 0 : ++row;
238 : }
239 0 : if (states.count(d.g_state()[row]) == 0)
240 : {
241 0 : throw std::runtime_error("Could not find a row with a sky state");
242 : }
243 0 : a.antno=refant;
244 0 : for(size_t k=0; k<4; ++k)
245 : {
246 0 : a.TObs[k]=d.g_wvrdata()(row, refant, k);
247 : }
248 0 : a.el=d.g_el()[row];
249 0 : a.time=d.g_time()[row];
250 0 : a.state=d.g_state()[row];
251 0 : a.source=d.g_source()[row];
252 0 : res.push_back(a);
253 : }
254 0 : return res;
255 0 : }
256 :
257 : dTdLCoeffsBase *
258 0 : SimpleSingle(const InterpArrayData &d, int refant)
259 : {
260 0 : std::vector<double> TObs(4);
261 : double el, time;
262 : size_t state;
263 0 : getMidPointData(d, refant, TObs, el, time, state);
264 :
265 0 : ALMAWVRCharacter wvrchar;
266 : ALMAAbsRet ar(TObs,
267 : M_PI/2.0,
268 0 : wvrchar);
269 0 : ALMAResBase res;
270 0 : ar.g_Res(res);
271 :
272 : // Convert to units of K/meter
273 0 : std::vector<double> dTdL(4);
274 0 : std::vector<double> dTdL_err(4);
275 0 : for(size_t i=0; i<4; ++i)
276 : {
277 0 : dTdL[i]=res.dTdL[i];
278 0 : dTdL_err[i]=res.dTdL_err[i];
279 : }
280 :
281 0 : convertKm(dTdL);
282 0 : convertKm(dTdL_err);
283 :
284 : return new dTdLCoeffsSingle(dTdL,
285 : dTdL_err
286 0 : );
287 :
288 0 : }
289 :
290 : // There is too much repetition here, need to merge with other
291 : // functions here
292 : dTdLCoeffsBase *
293 0 : SimpleSingleCont(const InterpArrayData &d, int refant)
294 : {
295 :
296 0 : std::vector<double> TObs(4);
297 : double el, time;
298 : size_t state;
299 0 : getMidPointData(d, refant, TObs, el, time, state);
300 :
301 0 : LibAIR2::ALMARetOpts opts;
302 0 : LibAIR2::ALMAContRes res;
303 0 : ALMAWVRCharacter wvrchar;
304 :
305 0 : ALMAAbsContRetrieve(TObs,
306 : el,
307 : wvrchar,
308 : res,
309 : opts);
310 :
311 : // Convert to units of K/meter
312 0 : std::vector<double> dTdL(4);
313 0 : std::vector<double> dTdL_err(4);
314 0 : for(size_t i=0; i<4; ++i)
315 : {
316 0 : dTdL[i]=res.dTdL[i];
317 0 : dTdL_err[i]=res.dTdL_err[i];
318 : }
319 :
320 0 : convertKm(dTdL);
321 0 : convertKm(dTdL_err);
322 :
323 : return new dTdLCoeffsSingle(dTdL,
324 0 : dTdL_err);
325 :
326 0 : }
327 :
328 0 : ALMAResBaseList doALMAAbsRet(ALMAAbsInpL &il,
329 : std::vector<std::pair<double, double> > &fb,
330 : AntSet& problemAnts,
331 : unsigned rseed)
332 : {
333 :
334 0 : problemAnts.clear();
335 :
336 0 : AntSet resZeroAnts;
337 0 : AntSet resNonZeroAnts;
338 :
339 0 : ALMAResBaseList res;
340 :
341 0 : ALMAAbsInpL newil;
342 0 : std::vector<std::pair<double, double> > newfb;
343 0 : bool fbFilled = (fb.size()>0);
344 :
345 0 : if (rseed==0){ // if rseed is set to zero, use a default value
346 0 : rseed=43;
347 : }
348 :
349 0 : size_t count=0;
350 :
351 0 : for(const ALMAAbsInput &x: il)
352 : {
353 0 : bool problematic = false;
354 0 : std::vector<double> TObs(4);
355 0 : for(size_t i=0; i<4; ++i){
356 0 : TObs[i]=x.TObs[i];
357 : }
358 : try {
359 0 : checkTObs(TObs);
360 : }
361 0 : catch(const std::runtime_error rE){
362 0 : std::cout << std::endl << "WARNING: problem with Tobs of antenna " << x.antno
363 0 : << std::endl << " LibAIR2::checkTObs: " << rE.what() << std::endl;
364 0 : std::cerr << std::endl << "WARNING: problem with Tobs of antenna " << x.antno
365 0 : << std::endl << " LibAIR2::checkTObs: " << rE.what() << std::endl;
366 0 : problematic = true;
367 0 : }
368 0 : ALMAWVRCharacter wvrchar;
369 :
370 0 : std::cout << "Using random seed " << rseed << std::endl;
371 :
372 : ALMAAbsRet ar(TObs,
373 0 : x.el,
374 : wvrchar,
375 0 : rseed);
376 :
377 0 : ALMAResBase *ares=new ALMAResBase;
378 0 : if(!ar.g_Res(*ares)){
379 0 : std::cout << "WARNING: Bayesian evidence was zero for antenna " << x.antno << std::endl
380 0 : << " TObs was " << TObs[0] << " " << TObs[1] << " " << TObs[2] << " " <<TObs[3]
381 0 : << " K, elevation " << x.el/M_PI*180. << " deg" << std::endl;
382 0 : std::cerr << "WARNING: Bayesian evidence was zero for antenna " << x.antno << std::endl
383 0 : << " TObs was " << TObs[0] << " " << TObs[1] << " " << TObs[2] << " " <<TObs[3]
384 0 : << " K, elevation " << x.el/M_PI*180. << " deg" << std::endl;
385 : }
386 : else{
387 0 : res.ptr_list.push_back(ares);
388 0 : newil.push_back(x);
389 0 : if(fbFilled){
390 0 : newfb.push_back(fb[count++]);
391 : }
392 : }
393 :
394 0 : if(problematic){
395 0 : problemAnts.insert(x.antno);
396 : }
397 :
398 0 : if(res.ptr_list.size()==0){
399 0 : resZeroAnts.insert(x.antno);
400 : }
401 : else{
402 0 : resNonZeroAnts.insert(x.antno);
403 : }
404 :
405 0 : }
406 :
407 0 : for(AntSet::const_iterator it=resZeroAnts.begin(); it!=resZeroAnts.end(); it++){
408 0 : if(resNonZeroAnts.count(*it)==0){ // in none of the retrievals was there a non-zero evidence
409 0 : problemAnts.insert(*it);
410 : }
411 : }
412 :
413 0 : il = newil;
414 0 : fb = newfb;
415 :
416 0 : return res;
417 0 : }
418 :
419 0 : static void ALMAAbsRetP(ALMAResBase *i,
420 : std::array<double, 4> &dTdL,
421 : std::array<double, 4> &dTdL_err)
422 : {
423 0 : for(size_t k=0; k<4; ++k)
424 : {
425 0 : dTdL[k]= i->dTdL[k]*1e3;
426 0 : dTdL_err[k]= i->dTdL_err[k]*1e3;
427 : }
428 0 : }
429 :
430 : dTdLCoeffsBase *
431 0 : ALMAAbsProcessor(const ALMAAbsInpL &inp,
432 : ALMAResBaseList &r)
433 : {
434 0 : assert(inp.size() == r.ptr_list.size());
435 :
436 0 : std::unique_ptr<dTdLCoeffsBase> res;
437 0 : if (inp.size()>0){
438 0 : dTdLCoeffsSingleInterpolated *rr=new dTdLCoeffsSingleInterpolated();
439 0 : res=std::unique_ptr<dTdLCoeffsBase>(rr);
440 :
441 0 : for(ALMAAbsInpL::const_iterator i=inp.begin();
442 0 : i!=inp.end();
443 0 : ++i)
444 : {
445 0 : ALMAResBase *rp = *(r.ptr_list.begin());
446 :
447 : std::array<double, 4> dTdL;
448 : std::array<double, 4> dTdL_err;
449 0 : ALMAAbsRetP(rp, dTdL, dTdL_err);
450 :
451 0 : rr->insert(i->time,
452 : dTdL,
453 : dTdL_err);
454 :
455 0 : r.ptr_list.pop_front(); // remove this almaresbase pointer from the list
456 0 : delete rp; // and free it
457 : }
458 :
459 : }
460 0 : return res.release();
461 0 : }
462 :
463 :
464 : dTdLCoeffsSingleInterpolated *
465 0 : SimpleMultiple(const std::vector<std::pair<double, double> > &fb,
466 : ALMAResBaseList &r)
467 : {
468 :
469 0 : assert(fb.size()==r.ptr_list.size());
470 :
471 0 : std::unique_ptr<dTdLCoeffsSingleInterpolated> res(new dTdLCoeffsSingleInterpolated());
472 :
473 0 : std::list<ALMAResBase*>::iterator ir=r.ptr_list.begin();
474 0 : for(size_t i=0; i<fb.size(); ++i)
475 : {
476 :
477 : // Convert to units of K/meter
478 : std::array<double, 4> dTdL;
479 : std::array<double, 4> dTdL_err;
480 0 : ALMAAbsRetP(*ir, dTdL, dTdL_err);
481 0 : res->insert(fb[i].first,
482 : dTdL,
483 : dTdL_err);
484 0 : res->insert(fb[i].second,
485 : dTdL,
486 : dTdL_err);
487 0 : ++ir;
488 : }
489 :
490 0 : return res.release();
491 :
492 0 : }
493 :
494 0 : void ALMAAbsContRetrieve(const std::vector<double> &TObs,
495 : double el,
496 : const ALMAWVRCharacter &WVRChar,
497 : ALMAContRes &res,
498 : const ALMARetOpts &opts)
499 : {
500 0 : CouplingModel *cm= new CouplingModel(mkCloudy(WVRChar,
501 : PartTable,
502 0 : AirCont));
503 0 : PPDipModel m(cm);
504 0 : AbsNormMeasure *ll=new AbsNormMeasure(m);
505 0 : Minim::IndependentFlatPriors pll(ll);
506 :
507 0 : cm->setSpill(0.98, 275);
508 0 : m.setZA(0.5 * M_PI -el);
509 0 : ll->obs=TObs;
510 0 : ll->thermNoise=std::vector<double>(TObs.size(),
511 0 : 1.0);
512 :
513 0 : if (opts.OSFPriors)
514 : {
515 0 : pll.AddPrior("n", 0, 15);
516 0 : pll.AddPrior("T", 250, 295);
517 0 : pll.AddPrior("P", 550, 750);
518 0 : pll.AddPrior("tau183", 0, 0.2);
519 : }
520 : else
521 : {
522 0 : pll.AddPrior("n", 0, 10);
523 0 : pll.AddPrior("T", 250, 295);
524 0 : pll.AddPrior("P", 300, 550);
525 0 : pll.AddPrior("tau183", 0, 0.2);
526 : }
527 :
528 :
529 0 : std::list<Minim::MCPoint> ss;
530 0 : startSetDirect(pll,
531 : 200,
532 : ss);
533 :
534 : // Create the nested sampler
535 0 : std::unique_ptr<Minim::NestedS> ns;
536 0 : ns.reset(new Minim::NestedS(pll));
537 0 : (*ns)["coupling"]->dofit=false;
538 :
539 : #ifdef ELLIPSESAMPLE
540 : Minim::EllipsoidCPSampler *cps=new Minim::EllipsoidCPSampler(pll,
541 : *ns);
542 : cps->reshape_maxp=50;
543 :
544 : ns->reset(ss,
545 : cps);
546 : cps->reshape();
547 : #else
548 0 : ns->reset(ss);
549 : #endif
550 :
551 : #if 0
552 : Minim::SOutMCMon *pp=new Minim::SOutMCMon();
553 : ns->mon=pp;
554 : cps->mon=pp;
555 : #endif
556 : //ns->InitalS(new Minim::InitialRandom(200));
557 :
558 0 : double evidence=ns->sample(3000);
559 : /// The posterior
560 0 : std::list<Minim::WPPoint> post;
561 0 : post=ns->g_post();
562 :
563 0 : if (post.size() < 10000 )
564 : {
565 0 : std::cout<<"Terminated after "<<post.size()
566 0 : <<std::endl;
567 : }
568 :
569 0 : res.ev=evidence;
570 :
571 0 : std::vector<double> m1(4), m2(4);
572 0 : moment1(post,
573 : evidence,
574 : m1);
575 :
576 0 : moment2(post,
577 : m1,
578 : evidence,
579 : m2);
580 :
581 0 : res.c=m1[0];
582 0 : res.c_err=std::pow(m2[0], 0.5);
583 :
584 0 : res.tau183=m1[3];
585 0 : res.tau183_err=std::pow(m2[3], 0.5);
586 :
587 0 : dTdLMom1(post,
588 0 : *ns,
589 : m,
590 : evidence,
591 : 1e-10,
592 0 : res.dTdL);
593 :
594 0 : dTdLMom2(post,
595 0 : *ns,
596 : m,
597 0 : res.dTdL,
598 : evidence,
599 : 1e-10,
600 0 : res.dTdL_err);
601 :
602 0 : }
603 :
604 : }
605 :
606 :
|
