Line data    Source code

       1             : /**
       2             : 
       3             :    Bojan Nikolic <bojan@bnikolic.co.uk>, <b.nikolic@mrao.cam.ac.uk>
       4             :    Initial version 2009
       5             : 
       6             :    \file twoerrline.cxx
       7             :    Renamed to twoerrline.cc 2023
       8             : 
       9             :    Fit lines to data which have errors in both coordinates
      10             : */
      11             : 
      12             : #include <memory>
      13             : #include "twoerrline.h"
      14             : #include "twoerrline_ml.h"
      15             : #include "gradientminim.h"
      16             : #include "lmmin.h"
      17             : #include "lmminutils.h"
      18             : 
      19             : namespace Minim {
      20             : 
      21           0 :   void LFit_MaxL(const std::vector<double> &xvals,
      22             :                  const std::vector<double> &yvals,
      23             :                  double sigma_x,
      24             :                  double sigma_y,
      25             :                  LineFit &res)
      26             :   {
      27             :     // We dont get the covariance matrix using the BFGS2 algorithm
      28           0 :     res.cv[0]=res.cv[3]=-1;
      29             : 
      30             :     Minim::LineTwoErrML lml(xvals, 
      31             :                             yvals,
      32             :                             sigma_x, 
      33           0 :                             sigma_y);
      34           0 :     lml.a=res.a;
      35           0 :     lml.b=res.b;
      36             : 
      37           0 :     Minim::BFGS2Minim m(lml);
      38           0 :     m.solve();
      39             :     
      40           0 :     res.a=lml.a;
      41           0 :     res.b=lml.b;
      42             :     
      43           0 :   }
      44             : 
      45           0 :   void LFit_LM(const std::vector<double> &xvals,
      46             :                const std::vector<double> &yvals,
      47             :                double sigma_x,
      48             :                double sigma_y,
      49             :                LineFit &res)
      50             :   {
      51             :     Minim::LineTwoErr_LavMarq lml(xvals, 
      52             :                                   yvals,
      53             :                                   sigma_x,
      54           0 :                                   sigma_y);
      55           0 :     Minim::LMMin m(lml);
      56             : 
      57           0 :     m.getbyname("a")->setp(res.a);
      58           0 :     m.getbyname("b")->setp(res.b);
      59             : 
      60           0 :     m.solve();
      61             :     
      62           0 :     res.a=m.getbyname("a")->getp();
      63           0 :     res.b=m.getbyname("b")->getp();
      64             :     
      65           0 :     std::unique_ptr<std::vector<double> > cv(CVMatrix(m,1e-5));
      66             : 
      67           0 :     std::copy(cv->begin(),
      68             :               cv->end(),
      69             :               res.cv.begin());
      70             :     
      71             : 
      72           0 :   }
      73             : 
      74             : }