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             :    \file almaabs_i.cpp
       7             :    Renamed almaabs_i.cc 2023
       8             :    
       9             : */
      10             : #include <iostream>
      11             : 
      12             : #include "almaabs_i.h"
      13             : 
      14             : #include "../model_make.h"
      15             : #include "../dtdltools.h"
      16             : #include "bnmin1/src/nestedinitial.h"
      17             : 
      18             : 
      19             : namespace LibAIR2 {
      20             : 
      21             :   const double iALMAAbsRetLL::thermNoise=1.0;
      22             :   const size_t iALMAAbsRet::n_ss=200;
      23             : 
      24          81 :   iALMAAbsRetLL::iALMAAbsRetLL(const std::vector<double> &TObs,
      25             :                                double el,
      26          81 :                                const ALMAWVRCharacter &WVRChar):
      27          81 :     cm(new CouplingModel(mkSingleLayerWater(WVRChar,
      28             :                                             PartTable,
      29          81 :                                             AirCont))),
      30          81 :     m(cm),
      31          81 :     ll(new AbsNormMeasure(m))
      32             :   {
      33          81 :     cm->setSpill(0.98, 275);
      34             :     
      35          81 :     m.setZA(0.5 * M_PI -el);
      36          81 :     ll->obs=TObs;
      37         162 :     ll->thermNoise=std::vector<double>(TObs.size(),
      38          81 :                                        thermNoise);
      39             : 
      40          81 :   }
      41             : 
      42          81 :   iALMAAbsRet::iALMAAbsRet(const std::vector<double> &TObs,
      43             :                            double el,
      44             :                            const ALMAWVRCharacter &WVRChar,
      45          81 :                            unsigned rseed):
      46          81 :     ls(TObs, 
      47             :        el, 
      48             :        WVRChar),
      49          81 :     pll(ls.ll),
      50          81 :     rseed(rseed),
      51         162 :     evidence()
      52             :   {
      53          81 :     pll.AddPrior("n", 0, 10);
      54          81 :     pll.AddPrior("T", 250, 295);
      55          81 :     pll.AddPrior("P", 300, 550);
      56          81 :   }
      57             : 
      58          81 :   bool iALMAAbsRet::sample(void)
      59             :   {
      60             :     // Create starting set
      61          81 :     std::list<Minim::MCPoint> ss;
      62          81 :     startSetDirect(pll,
      63             :                    n_ss,
      64             :                    ss,
      65             :                    rseed);
      66             : 
      67             :     // Create the nested sampler
      68          81 :     ns.reset(new Minim::NestedS(pll));
      69          81 :     (*ns)["coupling"]->dofit=false;
      70          81 :     ns->reset(ss);
      71             : 
      72             :     // So far not obvious it is necessary to enable this
      73             :     //ns->InitalS(new Minim::InitialRandom(n_ss));
      74             : 
      75          81 :     evidence=ns->sample(10000);
      76          81 :     post=ns->g_post();
      77             : 
      78          81 :     if (post.size() < 10000 )
      79             :     {
      80           4 :       std::cout<<"Terminated after "<<post.size()<<std::endl;
      81             :     }
      82          81 :     if(evidence == 0.){
      83           4 :       return false;
      84             :     }
      85          77 :     return true;
      86             : 
      87          81 :   }
      88             : 
      89          77 :   void iALMAAbsRet::g_Pars(ALMAResBase &res)
      90             :   {
      91          77 :     res.ev=evidence;
      92             : 
      93         154 :     std::vector<double> m1(3), m2(3);
      94          77 :     moment1(post,
      95             :             evidence,
      96             :             m1);
      97             : 
      98          77 :     moment2(post,
      99             :             m1,
     100             :             evidence,
     101             :             m2);
     102             : 
     103          77 :     res.c=m1[0];
     104          77 :     res.c_err=std::pow(m2[0], 0.5);
     105          77 :   }
     106             : 
     107          77 :   void iALMAAbsRet::g_Coeffs(ALMAResBase &r)
     108             :   {
     109          77 :     r.ev=evidence;
     110          77 :     dTdLMom1(post,
     111          77 :              *ns,
     112             :              ls.m,
     113             :              evidence,
     114             :              1e-10,
     115          77 :              r.dTdL);
     116             : 
     117          77 :     dTdLMom2(post,
     118          77 :              *ns,
     119             :              ls.m,
     120          77 :              r.dTdL,
     121             :              evidence,
     122             :              1e-10,
     123          77 :              r.dTdL_err);
     124          77 :   }
     125             : 
     126             : 
     127             : }
     128             : 
     129             : 
     130             :