Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "LinearOperator.h"
       4             : #include "BlockDescriptor.h"
       5             : 
       6             : #include <vector>
       7             : 
       8             : namespace elsa
       9             : {
      10             :     /**
      11             :      * @brief Class representing a block operator matrix
      12             :      *
      13             :      * @author Matthias Wieczorek - initial code
      14             :      * @author David Frank - rewrite
      15             :      * @author Nikola Dinev - automatic descriptor generation, rewrite
      16             :      *
      17             :      * @tparam data_t  data type for the domain and range of the operator, defaulting to real_t
      18             :      *
      19             :      * A block linear operator represents a block operator matrix \f$ B \f$ consisting of multiple
      20             :      * matrices \f$ A_i, i=1,\ldots,n \f$ stacked:
      21             :      * - rowwise \f[ B = \begin{bmatrix}
      22             :      *                      A_{1}\\
      23             :      *                      A_{2}\\
      24             :      *                      \vdots\\
      25             :      *                      A_{n}
      26             :      *                   \end{bmatrix} \f]
      27             :      *
      28             :      * - columnwise  \f[ B = \begin{bmatrix} A_{1} & A_{2} & \hdots & A_{n} \end{bmatrix} \f]
      29             :      */
      30             :     template <typename data_t = real_t>
      31             :     class BlockLinearOperator : public LinearOperator<data_t>
      32             :     {
      33             :     public:
      34             :         /// possible arrangements of the blocks
      35             :         enum BlockType {
      36             :             ROW,
      37             :             COL,
      38             :         };
      39             : 
      40             :         /// convenience typedef for a vector of pointers to LinearOperator
      41             :         using OperatorList = typename std::vector<std::unique_ptr<LinearOperator<data_t>>>;
      42             : 
      43             :         /**
      44             :          * @brief Construct a BlockLinearOperator of the given BlockType from the list of operators
      45             :          *
      46             :          * @param[in] ops the list of operators
      47             :          * @param[in] blockType the fashion in which the operators are to be stacked
      48             :          *
      49             :          * @throw InvalidArgumentError if ops is empty
      50             :          *
      51             :          * The domain and range descriptors of the BlockLinearOperator are generated automatically
      52             :          * based on the descriptors of the operators in the list. For the block descriptor, a
      53             :          * PartitionDescriptor is preferentially generated, if not possible a RandomBlocksDescriptor
      54             :          * is generated instead. For the non-block descriptor the best common descriptor is chosen
      55             :          * (see DataDescriptor::bestCommon()).
      56             :          */
      57             :         BlockLinearOperator(const OperatorList& ops, BlockType blockType);
      58             : 
      59             :         /**
      60             :          * @brief Construct a BlockLinearOperator of the given BlockType from the list of operators,
      61             :          * and additionally manually set the domain and range descriptors of the operator
      62             :          *
      63             :          * @param[in] domainDescriptor descriptor of the domain of the operator
      64             :          * @param[in] rangeDescriptor descriptor of the range of the operator
      65             :          * @param[in] ops the list of operators
      66             :          * @param[in] blockType the fashion in which the operators are to be stacked
      67             :          *
      68             :          * @throw InvalidArgumentError if the passed in descriptors are not suitable for the
      69             :          * BlockLinearOperator
      70             :          */
      71             :         BlockLinearOperator(const DataDescriptor& domainDescriptor,
      72             :                             const DataDescriptor& rangeDescriptor, const OperatorList& ops,
      73             :                             BlockType blockType);
      74             : 
      75             :         /// default destructor
      76           0 :         ~BlockLinearOperator() override = default;
      77             : 
      78             :         BlockLinearOperator& operator=(BlockLinearOperator&) = delete;
      79             : 
      80             :         /// return the operator corresponding to the i-th block of the matrix
      81             :         const LinearOperator<data_t>& getIthOperator(index_t i) const;
      82             : 
      83             :         /// return the total number of blocks
      84             :         index_t numberOfOps() const;
      85             : 
      86             :     protected:
      87             :         /// protected copy constructor; used for cloning
      88             :         BlockLinearOperator(const BlockLinearOperator& other);
      89             : 
      90             :         /// apply the block linear operator
      91             :         void applyImpl(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) const override;
      92             : 
      93             :         /// apply the adjoint of the block linear operator
      94             :         void applyAdjointImpl(const DataContainer<data_t>& y,
      95             :                               DataContainer<data_t>& Aty) const override;
      96             : 
      97             :         /// implement the polymorphic clone operation
      98             :         BlockLinearOperator<data_t>* cloneImpl() const override;
      99             : 
     100             :         /// implement the polymorphic comparison operation
     101             :         bool isEqual(const LinearOperator<data_t>& other) const override;
     102             : 
     103             :     private:
     104             :         /// list specifying the individual operators corresponding to each block
     105             :         OperatorList _operatorList;
     106             : 
     107             :         /// determines in which fashion the operators are concatenated - rowwise or columnwise
     108             :         BlockType _blockType;
     109             : 
     110             :         /// lift from base class
     111             :         using LinearOperator<data_t>::_domainDescriptor;
     112             :         using LinearOperator<data_t>::_rangeDescriptor;
     113             : 
     114             :         /// returns the best fitting domain descriptor based on the operator list and block type
     115             :         static std::unique_ptr<DataDescriptor>
     116             :             determineDomainDescriptor(const OperatorList& operatorList, BlockType blockType);
     117             : 
     118             :         /// returns the best fitting range descriptor based on the operator list and block type
     119             :         static std::unique_ptr<DataDescriptor>
     120             :             determineRangeDescriptor(const OperatorList& operatorList, BlockType blockType);
     121             : 
     122             :         /// finds the best fitting block descriptor, such that each block is described by the
     123             :         /// corresponding descriptor in the list
     124             :         static std::unique_ptr<BlockDescriptor>
     125             :             bestBlockDescriptor(const std::vector<const DataDescriptor*>&);
     126             :     };
     127             : } // namespace elsa