Line data    Source code

       1             : /**
       2             :    Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
       3             : 
       4             :    This file is part of LibAIR and is licensed under GNU Public
       5             :    License Version 2
       6             :    
       7             :    \file numalgo.hpp
       8             :    Renamed numalgo.h
       9             :    
      10             :    Numerical algorithms used in LibAIR
      11             : */
      12             : #ifndef _LIBAIR_NUMALGO_HPP__
      13             : #define _LIBAIR_NUMALGO_HPP__
      14             : 
      15             : #include <cmath>
      16             : #include <vector>
      17             : #include <array>
      18             : #include <functional>
      19             : #include <casacore/scimath/Mathematics/Combinatorics.h> // was <boost/math/special_functions/binomial.hpp>
      20             : 
      21             : namespace LibAIR2 {
      22             : 
      23             :   inline double 
      24    77889705 :   GaussError( std::vector<double>::const_iterator obs,
      25             :               const std::vector<double> & mod,
      26             :               const std::vector<double> & noise)
      27             :   {
      28    77889705 :     double res= 0 ;
      29    77889705 :     double pref=0;
      30   389448525 :     for (size_t i =0 ; i < noise.size() ; ++i)
      31             :     {
      32   311558820 :       res +=   std::pow(   ( *(obs+i) - mod[i])/noise[i], 2)*0.5;
      33   311558820 :       pref +=  0.5*std::log(2*M_PI*pow(noise[i],2));
      34             :     }
      35             : 
      36    77889705 :     return pref+res;
      37             : 
      38             :   }
      39             : 
      40             :   /** Compute finite differences up to (N-1)-th order of function f
      41             : 
      42             :       This should be accellerated by storing the previous evaluations
      43             :       
      44             :       \param c The central point around which to compute the difference
      45             :    */
      46             :   template<size_t N>
      47             :   std::array<double, N>
      48             :   CenFD(const std::function<double (double)> &f,
      49             :         double c,
      50             :         double delt)
      51             :   {
      52             :     std::array<double, N> res;
      53             :     for (size_t n=0; n<N; ++n)
      54             :     {
      55             :       res[n]=0;
      56             :       if (n==0)
      57             :       {
      58             :         res[n]=f(c);
      59             :       }
      60             :       else
      61             :       {
      62             :         for (size_t i=0; i<=n; ++i)
      63             :           res[n]+=std::pow(-1.0,(int)i)*casacore::Combinatorics::choose(n,i)*f(c+(0.5*n-i)*delt);
      64             :         res[n]/=std::pow(delt, (int)n);
      65             :       }
      66             :     }
      67             :     return res;
      68             :   }
      69             : 
      70             :   /** Compute the central difference derivative approximation of an
      71             :       array-valued function
      72             :       
      73             :       \param N Derivative order up to N-1 is calculated
      74             : 
      75             :       \param K Size of the array produced by the function
      76             : 
      77             :       \param f Function to differentiate
      78             : 
      79             :       \param c Central point around which to calculate the derivative
      80             :       
      81             :       \param delt Displacement to use for the differnetiation
      82             :    */
      83             :   template<size_t N, size_t K>
      84             :   std::array< std::array<double, K>, N>
      85           0 :   CenFDV(const std::function< std::array<double, K> (double)> &f,
      86             :          double c,
      87             :          double delt)
      88             :   {
      89             :     std::array< std::array<double, K>, N> res;
      90           0 :     for (size_t n=0; n<N; ++n)
      91             :     {
      92           0 :       for (size_t k=0; k<K; ++k)
      93             :       {
      94           0 :         res[n][k]=0;
      95             :       }
      96             : 
      97           0 :       if (n==0)
      98             :       {
      99           0 :         res[n]=f(c);
     100             :       }
     101             :       else
     102             :       {
     103           0 :         for (size_t i=0; i<=n; ++i)
     104             :         {
     105           0 :           double coeff=std::pow(-1.0,(int)i)*casacore::Combinatorics::choose(n,i);
     106           0 :           std::array<double, K> ir =f(c+(0.5*n-i)*delt);
     107           0 :           for(size_t k=0; k<K; ++k)
     108             :           {
     109           0 :             res[n][k]+=coeff*ir[k];
     110             :           }
     111             :         }         
     112           0 :         for(size_t k=0; k<K; ++k)
     113             :         {
     114           0 :           res[n][k]/=std::pow(delt, (int)n);
     115             :         }
     116             :       }
     117             :     }
     118           0 :     return res;
     119             :   }
     120             : 
     121             : 
     122             : }
     123             : #endif
     124             :