Line data    Source code

       1             : #include "DataHandlerMapCPU.h"
       2             : #include "DataHandlerCPU.h"
       3             : #include "TypeCasts.hpp"
       4             : 
       5             : namespace elsa
       6             : {
       7             :     template <typename data_t>
       8           0 :     DataHandlerMapCPU<data_t>::DataHandlerMapCPU(DataHandlerCPU<data_t>* dataOwner, data_t* data,
       9             :                                                  index_t n)
      10           0 :         : _map(data, n), _dataOwner{dataOwner}
      11             :     {
      12             :         // sanity checks performed in getBlock()
      13           0 : #pragma omp critical
      14             :         {
      15             :             // add self to list of Maps referring to the _dataOwner
      16           0 :             _dataOwner->_associatedMaps.push_front(this);
      17           0 :             _handle = _dataOwner->_associatedMaps.begin();
      18             :         }
      19           0 :     }
      20             : 
      21             :     template <typename data_t>
      22           0 :     DataHandlerMapCPU<data_t>::DataHandlerMapCPU(const DataHandlerMapCPU<data_t>& other)
      23           0 :         : _map{other._map}, _dataOwner{other._dataOwner}
      24             :     {
      25           0 : #pragma omp critical
      26             :         {
      27             :             // add self to list of Maps referring to the _dataOwner
      28           0 :             _dataOwner->_associatedMaps.push_front(this);
      29           0 :             _handle = _dataOwner->_associatedMaps.begin();
      30             :         }
      31           0 :     }
      32             : 
      33             :     template <typename data_t>
      34           0 :     DataHandlerMapCPU<data_t>::~DataHandlerMapCPU()
      35             :     {
      36             :         // remove self from list of Maps referring to the _dataOwner
      37           0 :         if (_dataOwner) {
      38           0 : #pragma omp critical
      39           0 :             _dataOwner->_associatedMaps.erase(_handle);
      40             :         }
      41           0 :     }
      42             : 
      43             :     template <typename data_t>
      44           0 :     index_t DataHandlerMapCPU<data_t>::getSize() const
      45             :     {
      46           0 :         return static_cast<index_t>(_map.size());
      47             :     }
      48             : 
      49             :     template <typename data_t>
      50           0 :     data_t& DataHandlerMapCPU<data_t>::operator[](index_t index)
      51             :     {
      52           0 :         _dataOwner->detach();
      53           0 :         return _map[index];
      54             :     }
      55             : 
      56             :     template <typename data_t>
      57           0 :     const data_t& DataHandlerMapCPU<data_t>::operator[](index_t index) const
      58             :     {
      59           0 :         return _map[index];
      60             :     }
      61             : 
      62             :     template <typename data_t>
      63           0 :     data_t DataHandlerMapCPU<data_t>::dot(const DataHandler<data_t>& v) const
      64             :     {
      65           0 :         if (v.getSize() != getSize())
      66           0 :             throw InvalidArgumentError("DataHandlerMapCPU: dot product argument has wrong size");
      67             : 
      68             :         // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
      69           0 :         if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
      70           0 :             return _map.dot(*otherHandler->_data);
      71           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&v)) {
      72           0 :             return _map.dot(otherHandler->_map);
      73             :         } else {
      74           0 :             return this->slowDotProduct(v);
      75             :         }
      76             :     }
      77             : 
      78             :     template <typename data_t>
      79           0 :     GetFloatingPointType_t<data_t> DataHandlerMapCPU<data_t>::squaredL2Norm() const
      80             :     {
      81           0 :         return _map.squaredNorm();
      82             :     }
      83             : 
      84             :     template <typename data_t>
      85           0 :     GetFloatingPointType_t<data_t> DataHandlerMapCPU<data_t>::l2Norm() const
      86             :     {
      87           0 :         return _map.norm();
      88             :     }
      89             : 
      90             :     template <typename data_t>
      91           0 :     index_t DataHandlerMapCPU<data_t>::l0PseudoNorm() const
      92             :     {
      93             :         using FloatType = GetFloatingPointType_t<data_t>;
      94           0 :         return (_map.array().cwiseAbs() >= std::numeric_limits<FloatType>::epsilon()).count();
      95             :     }
      96             : 
      97             :     template <typename data_t>
      98           0 :     GetFloatingPointType_t<data_t> DataHandlerMapCPU<data_t>::l1Norm() const
      99             :     {
     100           0 :         return _map.array().abs().sum();
     101             :     }
     102             : 
     103             :     template <typename data_t>
     104           0 :     GetFloatingPointType_t<data_t> DataHandlerMapCPU<data_t>::lInfNorm() const
     105             :     {
     106           0 :         return _map.array().abs().maxCoeff();
     107             :     }
     108             : 
     109             :     template <typename data_t>
     110           0 :     data_t DataHandlerMapCPU<data_t>::sum() const
     111             :     {
     112           0 :         return _map.sum();
     113             :     }
     114             : 
     115             :     template <typename data_t>
     116           0 :     data_t DataHandlerMapCPU<data_t>::minElement() const
     117             :     {
     118             :         if constexpr (isComplex<data_t>) {
     119           0 :             throw LogicError("DataHandlerCPU: minElement of complex type not supported");
     120             :         } else {
     121           0 :             return _map.minCoeff();
     122             :         }
     123             :     }
     124             : 
     125             :     template <typename data_t>
     126           0 :     data_t DataHandlerMapCPU<data_t>::maxElement() const
     127             :     {
     128             :         if constexpr (isComplex<data_t>) {
     129           0 :             throw LogicError("DataHandlerCPU: maxElement of complex type not supported");
     130             :         } else {
     131           0 :             return _map.maxCoeff();
     132             :         }
     133             :     }
     134             : 
     135             :     template <typename data_t>
     136           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::fft(const DataDescriptor& source_desc,
     137             :                                                         FFTNorm norm)
     138             :     {
     139             :         // detaches internally
     140           0 :         this->_dataOwner->fft(source_desc, norm);
     141           0 :         return *this;
     142             :     }
     143             : 
     144             :     template <typename data_t>
     145           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::ifft(const DataDescriptor& source_desc,
     146             :                                                          FFTNorm norm)
     147             :     {
     148             :         // detaches internally
     149           0 :         this->_dataOwner->ifft(source_desc, norm);
     150           0 :         return *this;
     151             :     }
     152             : 
     153             :     template <typename data_t>
     154           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator+=(const DataHandler<data_t>& v)
     155             :     {
     156           0 :         if (v.getSize() != getSize())
     157           0 :             throw InvalidArgumentError("DataHandler: addition argument has wrong size");
     158             : 
     159           0 :         _dataOwner->detach();
     160             : 
     161             :         // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
     162           0 :         if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
     163           0 :             _map += *otherHandler->_data;
     164           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&v)) {
     165           0 :             _map += otherHandler->_map;
     166             :         } else {
     167           0 :             this->slowAddition(v);
     168             :         }
     169             : 
     170           0 :         return *this;
     171             :     }
     172             : 
     173             :     template <typename data_t>
     174           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator-=(const DataHandler<data_t>& v)
     175             :     {
     176           0 :         if (v.getSize() != getSize())
     177           0 :             throw InvalidArgumentError("DataHandler: subtraction argument has wrong size");
     178             : 
     179           0 :         _dataOwner->detach();
     180             : 
     181             :         // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
     182           0 :         if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
     183           0 :             _map -= *otherHandler->_data;
     184           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&v)) {
     185           0 :             _map -= otherHandler->_map;
     186             :         } else {
     187           0 :             this->slowSubtraction(v);
     188             :         }
     189             : 
     190           0 :         return *this;
     191             :     }
     192             : 
     193             :     template <typename data_t>
     194           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator*=(const DataHandler<data_t>& v)
     195             :     {
     196           0 :         if (v.getSize() != getSize())
     197           0 :             throw InvalidArgumentError("DataHandler: multiplication argument has wrong size");
     198             : 
     199           0 :         _dataOwner->detach();
     200             : 
     201             :         // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
     202           0 :         if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
     203           0 :             _map.array() *= otherHandler->_data->array();
     204           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&v)) {
     205           0 :             _map.array() *= otherHandler->_map.array();
     206             :         } else {
     207           0 :             this->slowMultiplication(v);
     208             :         }
     209             : 
     210           0 :         return *this;
     211             :     }
     212             : 
     213             :     template <typename data_t>
     214           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator/=(const DataHandler<data_t>& v)
     215             :     {
     216           0 :         if (v.getSize() != getSize())
     217           0 :             throw InvalidArgumentError("DataHandler: division argument has wrong size");
     218             : 
     219           0 :         _dataOwner->detach();
     220             : 
     221             :         // use Eigen if the other handler is CPU or Map, otherwise use the slow fallback version
     222           0 :         if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&v)) {
     223           0 :             _map.array() /= otherHandler->_data->array();
     224           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&v)) {
     225           0 :             _map.array() /= otherHandler->_map.array();
     226             :         } else {
     227           0 :             this->slowDivision(v);
     228             :         }
     229             : 
     230           0 :         return *this;
     231             :     }
     232             : 
     233             :     template <typename data_t>
     234             :     DataHandlerMapCPU<data_t>&
     235           0 :         DataHandlerMapCPU<data_t>::operator=(const DataHandlerMapCPU<data_t>& v)
     236             :     {
     237           0 :         if (v.getSize() != getSize())
     238           0 :             throw InvalidArgumentError("DataHandler: assignment argument has wrong size");
     239             : 
     240           0 :         if (getSize() == _dataOwner->getSize() && v.getSize() == v._dataOwner->getSize()) {
     241           0 :             _dataOwner->attach(v._dataOwner->_data);
     242             :         } else {
     243           0 :             _dataOwner->detachWithUninitializedBlock(_map.data() - _dataOwner->_data->data(),
     244           0 :                                                      getSize());
     245           0 :             _map = v._map;
     246             :         }
     247             : 
     248           0 :         return *this;
     249             :     }
     250             : 
     251             :     template <typename data_t>
     252           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator+=(data_t scalar)
     253             :     {
     254           0 :         _dataOwner->detach();
     255             : 
     256           0 :         _map.array() += scalar;
     257           0 :         return *this;
     258             :     }
     259             : 
     260             :     template <typename data_t>
     261           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator-=(data_t scalar)
     262             :     {
     263           0 :         _dataOwner->detach();
     264             : 
     265           0 :         _map.array() -= scalar;
     266           0 :         return *this;
     267             :     }
     268             : 
     269             :     template <typename data_t>
     270           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator*=(data_t scalar)
     271             :     {
     272           0 :         _dataOwner->detach();
     273             : 
     274           0 :         _map.array() *= scalar;
     275           0 :         return *this;
     276             :     }
     277             : 
     278             :     template <typename data_t>
     279           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator/=(data_t scalar)
     280             :     {
     281           0 :         _dataOwner->detach();
     282             : 
     283           0 :         _map.array() /= scalar;
     284           0 :         return *this;
     285             :     }
     286             : 
     287             :     template <typename data_t>
     288           0 :     DataHandler<data_t>& DataHandlerMapCPU<data_t>::operator=(data_t scalar)
     289             :     {
     290           0 :         _dataOwner->detach();
     291             : 
     292           0 :         _map.setConstant(scalar);
     293           0 :         return *this;
     294             :     }
     295             : 
     296             :     template <typename data_t>
     297             :     std::unique_ptr<DataHandler<data_t>>
     298           0 :         DataHandlerMapCPU<data_t>::getBlock(index_t startIndex, index_t numberOfElements)
     299             :     {
     300           0 :         if (startIndex >= getSize() || numberOfElements > getSize() - startIndex)
     301           0 :             throw InvalidArgumentError("DataHandler: requested block out of bounds");
     302             : 
     303           0 :         return std::unique_ptr<DataHandlerMapCPU<data_t>>(
     304           0 :             new DataHandlerMapCPU{_dataOwner, _map.data() + startIndex, numberOfElements});
     305             :     }
     306             : 
     307             :     template <typename data_t>
     308             :     std::unique_ptr<const DataHandler<data_t>>
     309           0 :         DataHandlerMapCPU<data_t>::getBlock(index_t startIndex, index_t numberOfElements) const
     310             :     {
     311           0 :         if (startIndex >= getSize() || numberOfElements > getSize() - startIndex)
     312           0 :             throw InvalidArgumentError("DataHandler: requested block out of bounds");
     313             : 
     314             :         // using a const_cast here is fine as long as the DataHandlers never expose the internal
     315             :         // Eigen objects
     316           0 :         auto mutableData = const_cast<data_t*>(_map.data() + startIndex);
     317           0 :         return std::unique_ptr<const DataHandlerMapCPU<data_t>>(
     318           0 :             new DataHandlerMapCPU{_dataOwner, mutableData, numberOfElements});
     319             :     }
     320             : 
     321             :     template <typename data_t>
     322           0 :     DataHandlerCPU<data_t>* DataHandlerMapCPU<data_t>::cloneImpl() const
     323             :     {
     324           0 :         if (getSize() == _dataOwner->getSize()) {
     325           0 :             return new DataHandlerCPU<data_t>{*_dataOwner};
     326             :         } else {
     327           0 :             return new DataHandlerCPU<data_t>{_map};
     328             :         }
     329             :     }
     330             : 
     331             :     template <typename data_t>
     332           0 :     bool DataHandlerMapCPU<data_t>::isEqual(const DataHandler<data_t>& other) const
     333             :     {
     334           0 :         if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&other)) {
     335             : 
     336           0 :             if (_map.size() != otherHandler->_map.size())
     337           0 :                 return false;
     338             : 
     339           0 :             if (_map.data() != otherHandler->_map.data() && _map != otherHandler->_map) {
     340           0 :                 return false;
     341             :             }
     342             : 
     343           0 :             return true;
     344           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&other)) {
     345             : 
     346           0 :             if (_map.size() != otherHandler->_data->size())
     347           0 :                 return false;
     348             : 
     349           0 :             if (_map.data() != otherHandler->_data->data() && _map != *otherHandler->_data)
     350           0 :                 return false;
     351             : 
     352           0 :             return true;
     353             :         } else
     354           0 :             return false;
     355             :     }
     356             : 
     357             :     template <typename data_t>
     358           0 :     void DataHandlerMapCPU<data_t>::assign(const DataHandler<data_t>& other)
     359             :     {
     360             : 
     361           0 :         if (auto otherHandler = downcast_safe<DataHandlerMapCPU>(&other)) {
     362           0 :             if (getSize() == _dataOwner->getSize()
     363           0 :                 && otherHandler->getSize() == otherHandler->_dataOwner->getSize()) {
     364           0 :                 _dataOwner->attach(otherHandler->_dataOwner->_data);
     365             :             } else {
     366           0 :                 _dataOwner->detachWithUninitializedBlock(_map.data() - _dataOwner->_data->data(),
     367           0 :                                                          getSize());
     368           0 :                 _map = otherHandler->_map;
     369             :             }
     370           0 :         } else if (auto otherHandler = downcast_safe<DataHandlerCPU<data_t>>(&other)) {
     371           0 :             if (getSize() == _dataOwner->getSize()) {
     372           0 :                 _dataOwner->attach(otherHandler->_data);
     373             :             } else {
     374           0 :                 _dataOwner->detachWithUninitializedBlock(_map.data() - _dataOwner->_data->data(),
     375           0 :                                                          getSize());
     376           0 :                 _map = *otherHandler->_data;
     377             :             }
     378             :         } else
     379           0 :             this->slowAssign(other);
     380           0 :     }
     381             : 
     382             :     template <typename data_t>
     383           0 :     void DataHandlerMapCPU<data_t>::assign(DataHandler<data_t>&& other)
     384             :     {
     385           0 :         assign(other);
     386           0 :     }
     387             : 
     388             :     template <typename data_t>
     389           0 :     typename DataHandlerMapCPU<data_t>::DataMap_t DataHandlerMapCPU<data_t>::accessData()
     390             :     {
     391           0 :         _dataOwner->detach();
     392           0 :         return _map;
     393             :     }
     394             : 
     395             :     template <typename data_t>
     396           0 :     typename DataHandlerMapCPU<data_t>::DataMap_t DataHandlerMapCPU<data_t>::accessData() const
     397             :     {
     398           0 :         return _map;
     399             :     }
     400             : 
     401             :     // ------------------------------------------
     402             :     // explicit template instantiation
     403             :     template class DataHandlerMapCPU<float>;
     404             :     template class DataHandlerMapCPU<complex<float>>;
     405             :     template class DataHandlerMapCPU<double>;
     406             :     template class DataHandlerMapCPU<complex<double>>;
     407             :     template class DataHandlerMapCPU<index_t>;
     408             : 
     409             : } // namespace elsa