Line data    Source code

       1             : #include "LandweberIteration.h"
       2             : #include "DataContainer.h"
       3             : #include "Functional.h"
       4             : #include "LinearOperator.h"
       5             : #include "Logger.h"
       6             : #include "TypeCasts.hpp"
       7             : #include "PowerIterations.h"
       8             : #include <iostream>
       9             : 
      10             : namespace elsa
      11             : {
      12             :     template <typename data_t>
      13             :     LandweberIteration<data_t>::LandweberIteration(const LinearOperator<data_t>& A,
      14             :                                                    const DataContainer<data_t>& b, data_t stepSize)
      15             :         : Solver<data_t>(), A_(A.clone()), b_(b), residual_(emptylike(b)), stepSize_{stepSize}
      16           8 :     {
      17           8 :     }
      18             : 
      19             :     template <typename data_t>
      20             :     LandweberIteration<data_t>::LandweberIteration(const LinearOperator<data_t>& A,
      21             :                                                    const DataContainer<data_t>& b)
      22             :         : Solver<data_t>(), A_(A.clone()), b_(b), residual_(emptylike(b))
      23           0 :     {
      24           0 :     }
      25             : 
      26             :     template <typename data_t>
      27             :     DataContainer<data_t> LandweberIteration<data_t>::setup(std::optional<DataContainer<data_t>> x0)
      28           4 :     {
      29           4 :         auto x = extract_or(x0, A_->getDomainDescriptor());
      30             : 
      31           4 :         if (!tam_) {
      32           4 :             tam_ = setupOperators(*A_);
      33           4 :         }
      34             : 
      35           4 :         if (!stepSize_.isInitialized()) {
      36             :             // Choose step length to be just below \f$\frac{2}{\sigma^2}\f$, where \f$\sigma\f$
      37             :             // is the largest eigenvalue of \f$T * A^T * M * A\f$. This is computed using the power
      38             :             // iterations.
      39           0 :             auto Anorm = powerIterations(*tam_ * *A_);
      40           0 :             stepSize_ = 0.9 * (2. / Anorm);
      41           0 :         }
      42             : 
      43           4 :         Logger::get("LandweberIterations")->info("Using Steplength: {}", *stepSize_);
      44             : 
      45             :         // setup done!
      46           4 :         this->configured_ = true;
      47             : 
      48           4 :         return x;
      49           4 :     }
      50             : 
      51             :     template <typename data_t>
      52             :     DataContainer<data_t> LandweberIteration<data_t>::step(DataContainer<data_t> x)
      53           4 :     {
      54             :         // Compute Ax - b memory efficient
      55           4 :         A_->apply(x, residual_);
      56           4 :         residual_ -= b_;
      57             : 
      58           4 :         x -= (*stepSize_) * tam_->apply(residual_);
      59             : 
      60           4 :         projection_(x);
      61             : 
      62           4 :         return x;
      63           4 :     }
      64             : 
      65             :     template <typename data_t>
      66             :     std::string LandweberIteration<data_t>::formatHeader() const
      67           4 :     {
      68           4 :         return fmt::format("{:^12} | {:^12} |", "x norm", "Residual");
      69           4 :     }
      70             : 
      71             :     template <typename data_t>
      72             :     std::string LandweberIteration<data_t>::formatStep(const DataContainer<data_t>& x) const
      73           4 :     {
      74           4 :         return fmt::format("{:>12.5f} | {:>12.5f} |", x.l2Norm(), residual_.l2Norm());
      75           4 :     }
      76             : 
      77             :     template <typename data_t>
      78             :     bool LandweberIteration<data_t>::isEqual(const Solver<data_t>& other) const
      79           4 :     {
      80           4 :         auto landweber = downcast_safe<LandweberIteration<data_t>>(&other);
      81           4 :         return landweber && *A_ == *landweber->A_ && b_ == landweber->b_
      82           4 :                && stepSize_ == landweber->stepSize_;
      83           4 :     }
      84             : 
      85             :     template <class data_t>
      86             :     void LandweberIteration<data_t>::setProjection(
      87             :         const std::function<void(DataContainer<data_t>&)> projection)
      88           0 :     {
      89           0 :         projection_ = projection;
      90           0 :     }
      91             : 
      92             :     // ------------------------------------------
      93             :     // explicit template instantiation
      94             :     template class LandweberIteration<float>;
      95             :     template class LandweberIteration<double>;
      96             : } // namespace elsa