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           0 :   iALMAAbsRetLL::iALMAAbsRetLL(const std::vector<double> &TObs,
      25             :                                double el,
      26           0 :                                const ALMAWVRCharacter &WVRChar):
      27           0 :     cm(new CouplingModel(mkSingleLayerWater(WVRChar,
      28             :                                             PartTable,
      29           0 :                                             AirCont))),
      30           0 :     m(cm),
      31           0 :     ll(new AbsNormMeasure(m))
      32             :   {
      33           0 :     cm->setSpill(0.98, 275);
      34             :     
      35           0 :     m.setZA(0.5 * M_PI -el);
      36           0 :     ll->obs=TObs;
      37           0 :     ll->thermNoise=std::vector<double>(TObs.size(),
      38           0 :                                        thermNoise);
      39             : 
      40           0 :   }
      41             : 
      42           0 :   iALMAAbsRet::iALMAAbsRet(const std::vector<double> &TObs,
      43             :                            double el,
      44             :                            const ALMAWVRCharacter &WVRChar,
      45           0 :                            unsigned rseed):
      46           0 :     ls(TObs, 
      47             :        el, 
      48             :        WVRChar),
      49           0 :     pll(ls.ll),
      50           0 :     rseed(rseed),
      51           0 :     evidence()
      52             :   {
      53           0 :     pll.AddPrior("n", 0, 10);
      54           0 :     pll.AddPrior("T", 250, 295);
      55           0 :     pll.AddPrior("P", 300, 550);
      56           0 :   }
      57             : 
      58           0 :   bool iALMAAbsRet::sample(void)
      59             :   {
      60             :     // Create starting set
      61           0 :     std::list<Minim::MCPoint> ss;
      62           0 :     startSetDirect(pll,
      63             :                    n_ss,
      64             :                    ss,
      65             :                    rseed);
      66             : 
      67             :     // Create the nested sampler
      68           0 :     ns.reset(new Minim::NestedS(pll));
      69           0 :     (*ns)["coupling"]->dofit=false;
      70           0 :     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           0 :     evidence=ns->sample(10000);
      76           0 :     post=ns->g_post();
      77             : 
      78           0 :     if (post.size() < 10000 )
      79             :     {
      80           0 :       std::cout<<"Terminated after "<<post.size()<<std::endl;
      81             :     }
      82           0 :     if(evidence == 0.){
      83           0 :       return false;
      84             :     }
      85           0 :     return true;
      86             : 
      87           0 :   }
      88             : 
      89           0 :   void iALMAAbsRet::g_Pars(ALMAResBase &res)
      90             :   {
      91           0 :     res.ev=evidence;
      92             : 
      93           0 :     std::vector<double> m1(3), m2(3);
      94           0 :     moment1(post,
      95             :             evidence,
      96             :             m1);
      97             : 
      98           0 :     moment2(post,
      99             :             m1,
     100             :             evidence,
     101             :             m2);
     102             : 
     103           0 :     res.c=m1[0];
     104           0 :     res.c_err=std::pow(m2[0], 0.5);
     105           0 :   }
     106             : 
     107           0 :   void iALMAAbsRet::g_Coeffs(ALMAResBase &r)
     108             :   {
     109           0 :     r.ev=evidence;
     110           0 :     dTdLMom1(post,
     111           0 :              *ns,
     112             :              ls.m,
     113             :              evidence,
     114             :              1e-10,
     115           0 :              r.dTdL);
     116             : 
     117           0 :     dTdLMom2(post,
     118           0 :              *ns,
     119             :              ls.m,
     120           0 :              r.dTdL,
     121             :              evidence,
     122             :              1e-10,
     123           0 :              r.dTdL_err);
     124           0 :   }
     125             : 
     126             : 
     127             : }
     128             : 
     129             : 
     130             :