Line data    Source code

       1             : /** 
       2             :     \file twoerrline_imp.cxx
       3             : 
       4             :     Bojan Nikolic <bojan@bnikolic.co.uk>, <b.nikolic@mrao.cam.ac.uk>
       5             :     Initial version 2009    
       6             : 
       7             :     Renamed to twoerrline_ml.cc 2023
       8             : 
       9             : */
      10             : 
      11             : #include "twoerrline_ml.h"
      12             : 
      13             : #include <casacore/casa/Arrays/MatrixMath.h>
      14             : 
      15             : namespace Minim {
      16             : 
      17           0 :   LineTwoErrML::LineTwoErrML(const std::vector<double> &xvals,
      18             :                              const std::vector<double> &yvals,
      19             :                              double sigmax,
      20           0 :                              double sigmay):
      21           0 :     xobs(xvals.size()),
      22           0 :     yobs(yvals.size()),
      23           0 :     sigmax(sigmax),
      24           0 :     sigmay(sigmay),
      25           0 :     nobs(xvals.size())
      26             :   {
      27           0 :     for (size_t i=0; i<xobs.size(); ++i)
      28             :     {
      29           0 :       xobs[i]=xvals[i];
      30           0 :       yobs[i]=yvals[i];
      31             :     }
      32           0 :   }
      33             : 
      34           0 :   void LineTwoErrML::residuals(casacore::Vector<double> &res) const
      35             :   {
      36           0 :     casacore::Vector<double> ub(xobs.size(),b); // initialize with all vector elements == b
      37           0 :     res=yobs-xobs*a-ub;
      38           0 :   }
      39             : 
      40           0 :   double LineTwoErrML::lLikely(void) const
      41             :   {
      42           0 :     casacore::Vector<double> ub(xobs.size(),b);
      43           0 :     const double ressq=std::pow(casacore::norm(yobs-xobs*a-ub),2);
      44           0 :     const double r=0.5*ressq/(pow(sigmay,2)+ pow(sigmax*a,2));
      45           0 :     return r;
      46           0 :   }
      47             : 
      48           0 :   void LineTwoErrML::lGrd(std::vector< double > &res) const
      49             :   {
      50           0 :     res.resize(2);
      51           0 :     casacore::Vector<double> ub(xobs.size(),b);
      52           0 :     const casacore::Vector<double> rr(yobs-xobs*a-ub);
      53             :     
      54           0 :     const double st=(pow(sigmay,2)+ pow(sigmax*a,2));
      55             : 
      56           0 :     res[0]= -1.0* casacore::innerProduct(xobs,rr) / st - casacore::innerProduct(rr,rr)/ pow(st,2)* a* pow(sigmax,2);
      57           0 :     res[1]= -1.0* casacore::innerProduct(rr,rr) / st;
      58           0 :   }
      59             : 
      60           0 :   LineTwoErr_LavMarq::LineTwoErr_LavMarq(const std::vector<double> &xvals,
      61             :                                          const std::vector<double> &yvals,
      62             :                                          double sigmax,
      63           0 :                                          double sigmay):
      64           0 :     m(xvals,
      65             :       yvals,
      66             :       sigmax,
      67           0 :       sigmay)
      68             :   {
      69           0 :     m.a=1.0;
      70           0 :     m.b=1.0;
      71           0 :   }
      72             : 
      73           0 :   unsigned LineTwoErr_LavMarq::nres (void) const
      74             :   {
      75           0 :     return m.nobs;
      76             :   }
      77             : 
      78             : 
      79           0 :   void  LineTwoErr_LavMarq::residuals (std::vector<double> &res) const
      80             :   {
      81           0 :     res.resize(m.nobs);
      82           0 :     casacore::Vector<double> rr(m.nobs);
      83           0 :     m.residuals(rr);
      84           0 :     const double st= pow((pow(m.sigmay,2)+ pow(m.sigmax*m.a,2)), 0.5);
      85           0 :     rr=rr/st;
      86           0 :     std::copy(rr.begin(),
      87           0 :               rr.end(),
      88             :               res.begin());
      89           0 :   }
      90             : 
      91           0 :   void LineTwoErr_LavMarq::AddParams(std::vector< Minim::DParamCtr > &pars)
      92             :   {
      93           0 :     m.AddParams(pars);
      94           0 :   }
      95             : 
      96             : 
      97             : 
      98             : }
      99             : 
     100             : 
     101             : 
     102             :