Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <bojan@bnikolic.co.uk>, <b.nikolic@mrao.cam.ac.uk>
       3             :    Initial version 2008
       4             : 
       5             :    \file robustline.cxx
       6             :    Renamed to robustline.cc 2023
       7             : 
       8             : */
       9             : 
      10             : #include "robustline.h"
      11             : 
      12             : #include <cmath>
      13             : #include <iostream>
      14             : 
      15             : namespace Minim {
      16             : 
      17           0 :   RobustLineObsMod::RobustLineObsMod(const std::vector<double> &xobs,
      18           0 :                                      const std::vector<double> &yobs):
      19           0 :     nobs(xobs.size()),
      20           0 :     xobs(xobs),
      21           0 :     yobs(yobs),
      22           0 :     ysigma(xobs.size(),1.0)
      23             :   {
      24           0 :     a=1;
      25           0 :     b=0;
      26           0 :   }
      27             : 
      28           0 :   void RobustLineObsMod::residuals(std::vector<double> &res) const
      29             :   {
      30           0 :     res.resize(nobs);
      31           0 :     for(size_t i =0; i<nobs; ++i)
      32             :     {
      33           0 :       res[i]=(yobs[i]-xobs[i]*a-b)/ysigma[i];
      34             :     }
      35           0 :   }
      36             : 
      37           0 :   void RobustLineObsMod::dres_da(std::vector<double> &res) const
      38             :   {
      39           0 :     res.resize(nobs);
      40           0 :     for(size_t i =0; i<nobs; ++i)
      41             :     {
      42           0 :       res[i]= (-1.0 * xobs[i]/ysigma[i]);
      43             :     }
      44           0 :   }
      45             : 
      46           0 :   void RobustLineObsMod::dres_db(std::vector<double> &res) const
      47             :   {
      48           0 :     res.resize(nobs);
      49           0 :     for(size_t i =0; i<nobs; ++i)
      50             :     {
      51           0 :       res[i]= (-1.0/ysigma[i]) ;
      52             :     }
      53           0 :   }
      54             : 
      55           0 :   double RobustLineObsMod::lLikely(void) const
      56             :   {
      57           0 :     std::vector<double> res;
      58           0 :     residuals(res);
      59             : 
      60           0 :     double tot=0;
      61           0 :     for(size_t i=0; i<nobs; ++i)
      62             :     {
      63           0 :       if (fabs(res[i])  < ysigma[i])
      64             :       {
      65           0 :         tot += pow(res[i],2);
      66             :       }
      67             :       else
      68             :       {
      69           0 :         tot+=fabs(res[i]);
      70             :       }
      71             :     }
      72           0 :     return tot;
      73           0 :   }
      74             : 
      75           0 :   void RobustLineObsMod::lGrd(std::vector< double > &out) const
      76             :   {
      77           0 :     out.resize(2);
      78           0 :     out[0]=out[1]=0;
      79             : 
      80           0 :     std::vector<double> res;
      81           0 :     residuals(res);
      82             : 
      83           0 :     std::vector<double> da,db;
      84           0 :     dres_da(da);
      85           0 :     dres_db(db);    
      86             : 
      87           0 :     for(size_t i=0; i<nobs; ++i)
      88             :     {
      89           0 :       if(fabs(res[i]) < ysigma[i])
      90             :       {
      91           0 :         out[0]+= -2 * res[i] * xobs[i]/ pow(ysigma[i],2);
      92           0 :         out[1]+= -2 * res[i] / pow(ysigma[i],2);
      93             :       }
      94           0 :       else if(res[i]>0)
      95             :       {
      96           0 :         out[0] += da[i];
      97           0 :         out[1] += db[i];
      98             :       }
      99             :       else
     100             :       {
     101           0 :         out[0] -= da[i];
     102           0 :         out[1] -= db[i];
     103             :       }
     104             :     }
     105           0 :   }
     106             : 
     107             : }