Line data    Source code

       1             : // This file is part of the Mamico project. For conditions of distribution
       2             : // and use, please see the copyright notice in Mamico's main folder, or at
       3             : // www5.in.tum.de/mamico
       4             : 
       5             : /** Implementation of Gauss.h
       6             :  *  @author Felix Maurer
       7             :  */
       8             : 
       9             : // Member functions of Gauss.h:
      10             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
      11             : // TODO: This is very much hardcoded for kernel radius = 1
      12           0 : void coupling::filtering::Gauss<dim, scope...>::operator()() {
      13             : #ifdef DEBUG_GAUSS
      14             :   std::cout << "             GAUSS(Dim: " << _dim << "): Applying filter..." << std::endl;
      15             : #endif
      16             : 
      17             :   using namespace coupling::indexing;
      18             : 
      19             :   ScalarIndex index;
      20             :   ScalarIndex indexAbove;
      21             :   ScalarIndex indexBelow;
      22             : 
      23           0 :   for (unsigned int raw_index = 0; raw_index < _inputCells.size(); raw_index++) {
      24             : 
      25             :     // construct CellIndex using raw index
      26           0 :     index = ScalarIndex{raw_index};
      27             : 
      28             :     // make use of implicit conversion between linear <-> vector indices
      29           0 :     indexAbove = getIndexAbove(index, _dim);
      30           0 :     indexBelow = getIndexBelow(index, _dim);
      31             : 
      32             :     /*
      33             :     std::cout << "INDEX: " << index << std::endl;
      34             :     std::cout << "ABOVE: " << indexAbove << std::endl;
      35             :     std::cout << "BELOW: " << indexBelow << std::endl;
      36             :     */
      37             : 
      38             :     //[0] = below, [1] = at index, [2] = above
      39           0 :     double weights[3] = {_kernel[0], _kernel[1], _kernel[2]};
      40             : 
      41             :     // only one of these two cases can occur at once
      42           0 :     if (indexBelow == index) {
      43             :       // set weight of out of bounds index to zero
      44           0 :       weights[0] = 0;
      45             : 
      46           0 :       switch (_extrapolationStrategy) {
      47             :       case NONE:
      48             :         break;
      49           0 :       case MIRROR:
      50           0 :         weights[2] *= 2;
      51           0 :         break;
      52           0 :       case REFLECT:
      53           0 :         weights[1] *= 2;
      54           0 :         break;
      55             :       }
      56             : 
      57             :       // normalize the other two weights
      58           0 :       double weightsSum = weights[1] + weights[2];
      59           0 :       weights[1] /= weightsSum;
      60           0 :       weights[2] /= weightsSum;
      61             :     }
      62             : 
      63             :     // Since the domain has to be larger than one unit in all directions, not
      64             :     // both the index above and below can be illegal at the same time.
      65           0 :     else if (indexAbove == index) {
      66             :       // set weight of out of bounds index to zero
      67           0 :       weights[2] = 0;
      68             : 
      69           0 :       switch (_extrapolationStrategy) {
      70             :       case NONE:
      71             :         break;
      72           0 :       case MIRROR:
      73           0 :         weights[0] *= 2;
      74           0 :         break;
      75           0 :       case REFLECT:
      76           0 :         weights[1] *= 2;
      77           0 :         break;
      78             :       }
      79             : 
      80             :       // normalize the other two weights
      81           0 :       double weightsSum = weights[0] + weights[1];
      82           0 :       weights[0] /= weightsSum;
      83           0 :       weights[1] /= weightsSum;
      84             :     }
      85             : 
      86             :     // apply to scalars
      87           0 :     for (const auto scalarProperty : _scalarAccessFunctionPairs) {
      88           0 :       (_outputCells[raw_index]->*scalarProperty.set)((_inputCells[indexBelow.get()]->*scalarProperty.get)() * weights[0] +
      89           0 :                                                      (_inputCells[raw_index]->*scalarProperty.get)() * weights[1] +
      90           0 :                                                      (_inputCells[indexAbove.get()]->*scalarProperty.get)() * weights[2]);
      91             :     }
      92             :     // apply to vectors
      93           0 :     for (const auto vectorProperty : _vectorAccessFunctionPairs) {
      94           0 :       (_outputCells[raw_index]->*vectorProperty.set)((_inputCells[indexBelow.get()]->*vectorProperty.get)() * weights[0] +
      95           0 :                                                      (_inputCells[raw_index]->*vectorProperty.get)() * weights[1] +
      96           0 :                                                      (_inputCells[indexAbove.get()]->*vectorProperty.get)() * weights[2]);
      97             :     }
      98             :   } // index
      99             : 
     100             : #ifdef DEBUG_GAUSS
     101             :   std::cout << "done." << std::endl;
     102             : #endif
     103           0 : }
     104             : 
     105             : // Private functions of Gauss.h:
     106             : 
     107             : /**
     108             :  * @brief Computes the one-dimensional Gaussian kernel.
     109             :  * Kernel size is 1 + 2 * GAUSS_KERNEL_RADIUS, because there is one central cell a
     110             :  *  and GAUSS_KERNEL_RADIUS neighboring cells before and after that.
     111             :  *
     112             :  * @returns an array with kernel weights, normalized to sum(weights)=1
     113             :  */
     114             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
     115           0 : std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> coupling::filtering::Gauss<dim, scope...>::generateKernel() {
     116             :   // std::cout << "gauss: kernel radius: " << GAUSS_KERNEL_RADIUS << std::endl;
     117             : 
     118             :   std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> kernel;
     119           0 :   unsigned int i = 0; // index in kernel
     120           0 :   double sum = 0;     // used for normalization
     121             : 
     122             :   // fill kernel using Gaussian error function
     123           0 :   for (int x = -GAUSS_KERNEL_RADIUS; x <= GAUSS_KERNEL_RADIUS; x++) {
     124           0 :     kernel[i] = gaussianDensityFunction(x);
     125             : 
     126           0 :     sum += kernel[i];
     127           0 :     i++;
     128             :   }
     129             : 
     130             :   // normalize kernel. by definition, sum cannot be 0.
     131           0 :   for (auto& k : kernel) {
     132           0 :     k *= 1 / sum;
     133             : 
     134             : #ifdef DEBUG_GAUSS
     135             :     std::cout << "           GAUSS(Dim: " << _dim << "): Add kernel entry: " << k << std::endl;
     136             : #endif
     137             :   }
     138             : 
     139           0 :   return kernel;
     140             : }
     141             : 
     142             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
     143           0 : constexpr double coupling::filtering::Gauss<dim, scope...>::gaussianDensityFunction(int x) {
     144           0 :   double ans = (1 / (_sigma * sqrt(2 * M_PI))) * exp(-0.5 * (x / _sigma) * (x / _sigma));
     145           0 :   return ans;
     146             : }
     147             : 
     148             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
     149             : typename coupling::filtering::Gauss<dim, scope...>::VectorIndex
     150           0 : coupling::filtering::Gauss<dim, scope...>::getIndexAbove(const coupling::filtering::Gauss<dim, scope...>::VectorIndex index, unsigned int d) {
     151             :   // check for border indices
     152           0 :   const coupling::indexing::BaseIndex<dim> boundary = coupling::filtering::Gauss<dim, scope...>::VectorIndex::upperBoundary;
     153           0 :   for (unsigned int d = 0; d < dim; d++) {
     154             :     // convert index to BaseIndex and compare with boundary.
     155           0 :     if ((coupling::indexing::BaseIndex<dim>{index}).get()[d] >= boundary.get()[d])
     156             :       // if index is out of bounds, return it unmodified. this gets detected in
     157             :       // operator()'s boundary handling.
     158           0 :       return index;
     159             :   }
     160             : 
     161             :   // get raw vector index
     162           0 :   tarch::la::Vector<dim, int> ans = index.get();
     163             :   // get upper index in dimension d
     164           0 :   ans[d] += 1;
     165             : 
     166           0 :   return VectorIndex{ans};
     167             : }
     168             : 
     169             : template <unsigned int dim, coupling::indexing::IndexTrait... scope>
     170             : typename coupling::filtering::Gauss<dim, scope...>::VectorIndex
     171           0 : coupling::filtering::Gauss<dim, scope...>::getIndexBelow(const coupling::filtering::Gauss<dim, scope...>::VectorIndex index, unsigned int d) {
     172             :   // check for border indices
     173           0 :   const coupling::indexing::BaseIndex<dim> boundary = coupling::filtering::Gauss<dim, scope...>::VectorIndex::lowerBoundary;
     174           0 :   for (unsigned int d = 0; d < dim; d++) {
     175             :     // convert index to BaseIndex and compare with boundary.
     176           0 :     if ((coupling::indexing::BaseIndex<dim>{index}).get()[d] <= boundary.get()[d])
     177             :       // if index is out of bounds, return it unmodified. this gets detected in
     178             :       // operator()'s boundary handling.
     179           0 :       return index;
     180             :   }
     181             : 
     182             :   // copy input index
     183           0 :   tarch::la::Vector<dim, int> ans = index.get();
     184             :   // get lower index in dimension d
     185           0 :   ans[d] -= 1;
     186             : 
     187           0 :   return VectorIndex{ans};
     188             : }