Line data    Source code

       1             : #include "EmissionLogLikelihood.h"
       2             : #include "Scaling.h"
       3             : #include "TypeCasts.hpp"
       4             : 
       5             : #include <cmath>
       6             : #include <stdexcept>
       7             : 
       8             : namespace elsa
       9             : {
      10             :     template <typename data_t>
      11           0 :     EmissionLogLikelihood<data_t>::EmissionLogLikelihood(const DataDescriptor& domainDescriptor,
      12             :                                                          const DataContainer<data_t>& y,
      13             :                                                          const DataContainer<data_t>& r)
      14             :         : Functional<data_t>(domainDescriptor),
      15           0 :           _y{std::make_unique<DataContainer<data_t>>(y)},
      16           0 :           _r{std::make_unique<DataContainer<data_t>>(r)}
      17             :     {
      18             :         // sanity check
      19           0 :         if (domainDescriptor != y.getDataDescriptor() || domainDescriptor != r.getDataDescriptor())
      20           0 :             throw InvalidArgumentError(
      21             :                 "EmissionLogLikelihood: descriptor and y/r not matching in size.");
      22           0 :     }
      23             : 
      24             :     template <typename data_t>
      25           0 :     EmissionLogLikelihood<data_t>::EmissionLogLikelihood(const DataDescriptor& domainDescriptor,
      26             :                                                          const DataContainer<data_t>& y)
      27           0 :         : Functional<data_t>(domainDescriptor), _y{std::make_unique<DataContainer<data_t>>(y)}
      28             :     {
      29             :         // sanity check
      30           0 :         if (domainDescriptor != y.getDataDescriptor())
      31           0 :             throw InvalidArgumentError(
      32             :                 "EmissionLogLikelihood: descriptor and y not matching in size.");
      33           0 :     }
      34             : 
      35             :     template <typename data_t>
      36           0 :     EmissionLogLikelihood<data_t>::EmissionLogLikelihood(const Residual<data_t>& residual,
      37             :                                                          const DataContainer<data_t>& y,
      38             :                                                          const DataContainer<data_t>& r)
      39             :         : Functional<data_t>(residual),
      40           0 :           _y{std::make_unique<DataContainer<data_t>>(y)},
      41           0 :           _r{std::make_unique<DataContainer<data_t>>(r)}
      42             :     {
      43             :         // sanity check
      44           0 :         if (residual.getRangeDescriptor() != y.getDataDescriptor()
      45           0 :             || residual.getRangeDescriptor() != r.getDataDescriptor())
      46           0 :             throw InvalidArgumentError(
      47             :                 "EmissionLogLikelihood: residual and y/r not matching in size.");
      48           0 :     }
      49             : 
      50             :     template <typename data_t>
      51           0 :     EmissionLogLikelihood<data_t>::EmissionLogLikelihood(const Residual<data_t>& residual,
      52             :                                                          const DataContainer<data_t>& y)
      53           0 :         : Functional<data_t>(residual), _y{std::make_unique<DataContainer<data_t>>(y)}
      54             :     {
      55             :         // sanity check
      56           0 :         if (residual.getRangeDescriptor() != y.getDataDescriptor())
      57           0 :             throw InvalidArgumentError(
      58             :                 "EmissionLogLikelihood: residual and y not matching in size.");
      59           0 :     }
      60             : 
      61             :     template <typename data_t>
      62           0 :     data_t EmissionLogLikelihood<data_t>::evaluateImpl(const DataContainer<data_t>& Rx)
      63             :     {
      64           0 :         auto result = static_cast<data_t>(0.0);
      65             : 
      66           0 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      67           0 :             data_t temp = Rx[i];
      68           0 :             if (_r)
      69           0 :                 temp += (*_r)[i];
      70             : 
      71           0 :             result += temp - (*_y)[i] * std::log(temp);
      72             :         }
      73             : 
      74           0 :         return result;
      75             :     }
      76             : 
      77             :     template <typename data_t>
      78           0 :     void EmissionLogLikelihood<data_t>::getGradientInPlaceImpl(DataContainer<data_t>& Rx)
      79             :     {
      80           0 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      81           0 :             data_t temp = Rx[i];
      82           0 :             if (_r)
      83           0 :                 temp += (*_r)[i];
      84             : 
      85           0 :             Rx[i] = 1 - (*_y)[i] / temp;
      86             :         }
      87           0 :     }
      88             : 
      89             :     template <typename data_t>
      90             :     LinearOperator<data_t>
      91           0 :         EmissionLogLikelihood<data_t>::getHessianImpl(const DataContainer<data_t>& Rx)
      92             :     {
      93           0 :         DataContainer<data_t> scaleFactors(Rx.getDataDescriptor());
      94           0 :         for (index_t i = 0; i < Rx.getSize(); ++i) {
      95           0 :             data_t temp = Rx[i];
      96           0 :             if (_r)
      97           0 :                 temp += (*_r)[i];
      98             : 
      99           0 :             scaleFactors[i] = (*_y)[i] / (temp * temp);
     100             :         }
     101             : 
     102           0 :         return leaf(Scaling<data_t>(Rx.getDataDescriptor(), scaleFactors));
     103           0 :     }
     104             : 
     105             :     template <typename data_t>
     106           0 :     EmissionLogLikelihood<data_t>* EmissionLogLikelihood<data_t>::cloneImpl() const
     107             :     {
     108           0 :         if (_r)
     109           0 :             return new EmissionLogLikelihood<data_t>(this->getResidual(), *_y, *_r);
     110             :         else
     111           0 :             return new EmissionLogLikelihood<data_t>(this->getResidual(), *_y);
     112             :     }
     113             : 
     114             :     template <typename data_t>
     115           0 :     bool EmissionLogLikelihood<data_t>::isEqual(const Functional<data_t>& other) const
     116             :     {
     117           0 :         if (!Functional<data_t>::isEqual(other))
     118           0 :             return false;
     119             : 
     120           0 :         auto otherELL = downcast_safe<EmissionLogLikelihood>(&other);
     121           0 :         if (!otherELL)
     122           0 :             return false;
     123             : 
     124           0 :         if (*_y != *otherELL->_y)
     125           0 :             return false;
     126             : 
     127           0 :         if (_r && *_r != *otherELL->_r)
     128           0 :             return false;
     129             : 
     130           0 :         return true;
     131             :     }
     132             : 
     133             :     // ------------------------------------------
     134             :     // explicit template instantiation
     135             :     template class EmissionLogLikelihood<float>;
     136             :     template class EmissionLogLikelihood<double>;
     137             :     // no complex instantiations, they make no sense
     138             : 
     139             : } // namespace elsa