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          81 :   PriorNLikelihood::PriorNLikelihood(MLikelihood * mod):
      18          81 :     _mod(mod)
      19             :   {
      20          81 :   }
      21             : 
      22          81 :   PriorNLikelihood::~PriorNLikelihood()
      23             :   {
      24          81 :   }
      25             : 
      26         162 :   void  PriorNLikelihood::AddParams (std::vector< Minim::DParamCtr > &pars)
      27             :   {
      28         162 :     _mod->AddParams(pars);
      29         162 :   }
      30             : 
      31          81 :   IndependentPriors::IndependentPriors(MLikelihood * mod):
      32          81 :     PriorNLikelihood(mod)    
      33             :   {
      34          81 :     AddParams(_mpars);
      35          81 :   }
      36             : 
      37          81 :   IndependentPriors::~IndependentPriors(void)
      38             :   {
      39          81 :   }
      40             : 
      41         243 :   void IndependentPriors::AddPrior( const std::string & pname,
      42             :                                     double low,
      43             :                                     double high)
      44             :   {
      45             :     fprior_t pr;
      46             :     
      47         243 :     DParamCtr * par=findName(_mpars,pname);
      48         243 :     if (par == NULL)
      49           0 :       throw ParamError(pname+" not found");
      50         243 :     pr.p=par->p;
      51         243 :     pr.pmin=low;
      52         243 :     pr.pmax=high;
      53         243 :     priorlist.push_back(pr);
      54         243 :   }
      55             : 
      56    77873505 :   double IndependentFlatPriors::pprob(void) const
      57             :   {
      58    77873505 :     double lpriorprop = 0 ;
      59    77873505 :     for (priorlist_t::const_iterator i = pbegin();
      60   311494020 :          i != pend();
      61   233620515 :          ++i)
      62             :     {
      63   233620515 :       const double p = * (i->p);
      64   233620515 :       if (p < i->pmin or p > i->pmax)
      65     2856434 :         lpriorprop += lkl_h;
      66             :     }
      67    77873505 :     return lpriorprop;
      68             :   }
      69             : 
      70          81 :   IndependentFlatPriors::IndependentFlatPriors(MLikelihood * mod):
      71          81 :     IndependentPriors(mod)
      72             :   {
      73          81 :   }
      74             :   
      75          81 :   IndependentFlatPriors::~IndependentFlatPriors(void)
      76             :   {
      77          81 :   }
      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             :