Line data    Source code

       1             : #include "LinearResidual.h"
       2             : #include "Identity.h"
       3             : #include "TypeCasts.hpp"
       4             : 
       5             : #include <stdexcept>
       6             : 
       7             : namespace elsa
       8             : {
       9             :     template <typename data_t>
      10           0 :     LinearResidual<data_t>::LinearResidual(const DataDescriptor& descriptor)
      11           0 :         : Residual<data_t>(descriptor, descriptor)
      12             :     {
      13           0 :     }
      14             : 
      15             :     template <typename data_t>
      16           0 :     LinearResidual<data_t>::LinearResidual(const DataContainer<data_t>& b)
      17           0 :         : Residual<data_t>(b.getDataDescriptor(), b.getDataDescriptor()), _dataVector{b}
      18             :     {
      19           0 :     }
      20             : 
      21             :     template <typename data_t>
      22           0 :     LinearResidual<data_t>::LinearResidual(const LinearOperator<data_t>& A)
      23           0 :         : Residual<data_t>(A.getDomainDescriptor(), A.getRangeDescriptor()), _operator{A.clone()}
      24             :     {
      25           0 :     }
      26             : 
      27             :     template <typename data_t>
      28           0 :     LinearResidual<data_t>::LinearResidual(const LinearOperator<data_t>& A,
      29             :                                            const DataContainer<data_t>& b)
      30             :         : Residual<data_t>(A.getDomainDescriptor(), A.getRangeDescriptor()),
      31           0 :           _operator{A.clone()},
      32           0 :           _dataVector{b}
      33             :     {
      34           0 :         if (A.getRangeDescriptor().getNumberOfCoefficients() != b.getSize())
      35           0 :             throw InvalidArgumentError("LinearResidual: A and b do not match");
      36           0 :     }
      37             : 
      38             :     template <typename data_t>
      39           0 :     bool LinearResidual<data_t>::hasOperator() const
      40             :     {
      41           0 :         return static_cast<bool>(_operator);
      42             :     }
      43             : 
      44             :     template <typename data_t>
      45           0 :     bool LinearResidual<data_t>::hasDataVector() const
      46             :     {
      47           0 :         return _dataVector.has_value();
      48             :     }
      49             : 
      50             :     template <typename data_t>
      51           0 :     const LinearOperator<data_t>& LinearResidual<data_t>::getOperator() const
      52             :     {
      53           0 :         if (!_operator)
      54           0 :             throw Error("LinearResidual::getOperator: operator not present");
      55             : 
      56           0 :         return *_operator;
      57             :     }
      58             : 
      59             :     template <typename data_t>
      60           0 :     const DataContainer<data_t>& LinearResidual<data_t>::getDataVector() const
      61             :     {
      62           0 :         if (!_dataVector)
      63           0 :             throw Error("LinearResidual::getDataVector: data vector not present");
      64             : 
      65           0 :         return *_dataVector;
      66             :     }
      67             : 
      68             :     template <typename data_t>
      69           0 :     LinearResidual<data_t>* LinearResidual<data_t>::cloneImpl() const
      70             :     {
      71           0 :         if (hasOperator() && hasDataVector())
      72           0 :             return new LinearResidual<data_t>(getOperator(), getDataVector());
      73             : 
      74           0 :         if (hasOperator())
      75           0 :             return new LinearResidual<data_t>(getOperator());
      76             : 
      77           0 :         if (hasDataVector())
      78           0 :             return new LinearResidual<data_t>(getDataVector());
      79             : 
      80           0 :         return new LinearResidual<data_t>(this->getDomainDescriptor());
      81             :     }
      82             : 
      83             :     template <typename data_t>
      84           0 :     bool LinearResidual<data_t>::isEqual(const Residual<data_t>& other) const
      85             :     {
      86           0 :         if (!Residual<data_t>::isEqual(other))
      87           0 :             return false;
      88             : 
      89           0 :         auto otherLinearResidual = downcast_safe<LinearResidual>(&other);
      90           0 :         if (!otherLinearResidual)
      91           0 :             return false;
      92             : 
      93           0 :         if (hasOperator() != otherLinearResidual->hasOperator()
      94           0 :             || hasDataVector() != otherLinearResidual->hasDataVector())
      95           0 :             return false;
      96             : 
      97           0 :         if ((_operator && !otherLinearResidual->_operator)
      98           0 :             || (!_operator && otherLinearResidual->_operator)
      99           0 :             || (_dataVector && !otherLinearResidual->_dataVector)
     100           0 :             || (!_dataVector && otherLinearResidual->_dataVector))
     101           0 :             return false;
     102             : 
     103           0 :         if (_operator && otherLinearResidual->_operator
     104           0 :             && *_operator != *otherLinearResidual->_operator)
     105           0 :             return false;
     106             : 
     107           0 :         if (_dataVector && otherLinearResidual->_dataVector
     108           0 :             && *_dataVector != *otherLinearResidual->_dataVector)
     109           0 :             return false;
     110             : 
     111           0 :         return true;
     112             :     }
     113             : 
     114             :     template <typename data_t>
     115           0 :     void LinearResidual<data_t>::evaluateImpl(const DataContainer<data_t>& x,
     116             :                                               DataContainer<data_t>& result) const
     117             :     {
     118           0 :         if (hasOperator())
     119           0 :             _operator->apply(x, result);
     120             :         else
     121           0 :             result = x;
     122             : 
     123           0 :         if (hasDataVector())
     124           0 :             result -= *_dataVector;
     125           0 :     }
     126             : 
     127             :     template <typename data_t>
     128             :     LinearOperator<data_t>
     129           0 :         LinearResidual<data_t>::getJacobianImpl([[maybe_unused]] const DataContainer<data_t>& x)
     130             :     {
     131           0 :         if (hasOperator())
     132           0 :             return leaf(*_operator);
     133             :         else
     134           0 :             return leaf(Identity<data_t>(this->getRangeDescriptor()));
     135             :     }
     136             : 
     137             :     // ------------------------------------------
     138             :     // explicit template instantiation
     139             :     template class LinearResidual<float>;
     140             :     template class LinearResidual<double>;
     141             :     template class LinearResidual<complex<float>>;
     142             :     template class LinearResidual<complex<double>>;
     143             : 
     144             : } // namespace elsa