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