Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <bojan@bnikolic.co.uk>, <b.nikolic@mrao.cam.ac.uk>
       3             : 
       4             :    \file priors.cxx
       5             :    Renamed to priors.cc 2023
       6             : 
       7             : */
       8             : 
       9             : #include "priors.h"
      10             : 
      11             : #include "paramalgo.h"
      12             : #include <cmath>
      13             : 
      14             : namespace Minim {
      15             : 
      16             :   const double PriorNLikelihood::lkl_h=1e9;
      17           0 :   PriorNLikelihood::PriorNLikelihood(MLikelihood * mod):
      18           0 :     _mod(mod)
      19             :   {
      20           0 :   }
      21             : 
      22           0 :   PriorNLikelihood::~PriorNLikelihood()
      23             :   {
      24           0 :   }
      25             : 
      26           0 :   void  PriorNLikelihood::AddParams (std::vector< Minim::DParamCtr > &pars)
      27             :   {
      28           0 :     _mod->AddParams(pars);
      29           0 :   }
      30             : 
      31           0 :   IndependentPriors::IndependentPriors(MLikelihood * mod):
      32           0 :     PriorNLikelihood(mod)    
      33             :   {
      34           0 :     AddParams(_mpars);
      35           0 :   }
      36             : 
      37           0 :   IndependentPriors::~IndependentPriors(void)
      38             :   {
      39           0 :   }
      40             : 
      41           0 :   void IndependentPriors::AddPrior( const std::string & pname,
      42             :                                     double low,
      43             :                                     double high)
      44             :   {
      45             :     fprior_t pr;
      46             :     
      47           0 :     DParamCtr * par=findName(_mpars,pname);
      48           0 :     if (par == NULL)
      49           0 :       throw ParamError(pname+" not found");
      50           0 :     pr.p=par->p;
      51           0 :     pr.pmin=low;
      52           0 :     pr.pmax=high;
      53           0 :     priorlist.push_back(pr);
      54           0 :   }
      55             : 
      56           0 :   double IndependentFlatPriors::pprob(void) const
      57             :   {
      58           0 :     double lpriorprop = 0 ;
      59           0 :     for (priorlist_t::const_iterator i = pbegin();
      60           0 :          i != pend();
      61           0 :          ++i)
      62             :     {
      63           0 :       const double p = * (i->p);
      64           0 :       if (p < i->pmin or p > i->pmax)
      65           0 :         lpriorprop += lkl_h;
      66             :     }
      67           0 :     return lpriorprop;
      68             :   }
      69             : 
      70           0 :   IndependentFlatPriors::IndependentFlatPriors(MLikelihood * mod):
      71           0 :     IndependentPriors(mod)
      72             :   {
      73           0 :   }
      74             :   
      75           0 :   IndependentFlatPriors::~IndependentFlatPriors(void)
      76             :   {
      77           0 :   }
      78             : 
      79           0 :   LogFlatPriors::LogFlatPriors(MLikelihood * mod):
      80           0 :     IndependentPriors(mod)
      81             :   {
      82           0 :   }
      83             :   
      84           0 :   LogFlatPriors::~LogFlatPriors(void)
      85             :   {
      86           0 :   }
      87             : 
      88           0 :   double LogFlatPriors::pprob(void) const
      89             :   {
      90           0 :     double lpriorprop = 0 ;
      91           0 :     for (priorlist_t::const_iterator i = pbegin();
      92           0 :          i != pend();
      93           0 :          ++i)
      94             :     {
      95           0 :       const double p = * (i->p);
      96           0 :       if (p < i->pmin or p > i->pmax)
      97             :       {
      98           0 :         lpriorprop += lkl_h;
      99             :       }
     100             :       else
     101             :       {
     102           0 :         lpriorprop += std::log(p);
     103             :       }
     104             :     }
     105           0 :     return lpriorprop;
     106             :   }
     107             : 
     108             : 
     109             : }
     110             :