C++ GridStorage类代码示例

  /**
   * Descends on all dimensions beside dim_sweep. Class functor for dim_sweep.
   * Boundaries are regarded
   *
   * @param source coefficients of the sparse grid
   * @param result coefficients of the function computed by sweep
   * @param index current grid position
   * @param dim_list list of dimensions, that should be handled
   * @param dim_rem number of remaining dims
   * @param dim_sweep static dimension, in this dimension the functor is executed
   */
  void sweep_Boundary_rec(DataMatrix& source, DataMatrix& result,
                          grid_iterator& index,
                          std::vector<size_t>& dim_list, size_t dim_rem,
                          size_t dim_sweep) {
    if (dim_rem == 0) {
      functor(source, result, index, dim_sweep);
    } else {
      typedef GridStorage::index_type::level_type level_type;
      typedef GridStorage::index_type::index_type index_type;

      level_type current_level;
      index_type current_index;

      index.get(dim_list[dim_rem - 1], current_level, current_index);

      // handle level greater zero
      if (current_level > 0) {
        // given current point to next dim
        sweep_Boundary_rec(source, result, index, dim_list, dim_rem - 1,
                           dim_sweep);

        if (!index.hint()) {
          index.leftChild(dim_list[dim_rem - 1]);

          if (!storage->end(index.seq())) {
            sweep_Boundary_rec(source, result, index, dim_list, dim_rem,
                               dim_sweep);
          }

          index.stepRight(dim_list[dim_rem - 1]);

          if (!storage->end(index.seq())) {
            sweep_Boundary_rec(source, result, index, dim_list, dim_rem,
                               dim_sweep);
          }

          index.up(dim_list[dim_rem - 1]);
        }
      } else {  // handle level zero
        sweep_Boundary_rec(source, result, index, dim_list, dim_rem - 1,
                           dim_sweep);

        index.resetToRightLevelZero(dim_list[dim_rem - 1]);
        sweep_Boundary_rec(source, result, index, dim_list, dim_rem - 1,
                           dim_sweep);

        if (!index.hint()) {
          index.resetToLevelOne(dim_list[dim_rem - 1]);

          if (!storage->end(index.seq())) {
            sweep_Boundary_rec(source, result, index, dim_list, dim_rem,
                               dim_sweep);
          }
        }

        index.resetToLeftLevelZero(dim_list[dim_rem - 1]);
      }
    }
  }

  /**
   * Recursive traversal of the "tree" of basis functions for evaluation, used in operator().
   * For a given evaluation point \f$x\f$, it stores tuples (std::pair) of
   * \f$(i,\phi_i(x))\f$ in the result vector for all basis functions that are non-zero.
   *
   * @param basis a sparse grid basis
   * @param point evaluation point within the domain
   * @param current_dim the dimension currently looked at (recursion parameter)
   * @param value the value of the evaluation of the current basis function up to (excluding) dimension current_dim (product of the evaluations of the one-dimensional ones)
   * @param working iterator working on the GridStorage of the basis
   * @param source array of indices for each dimension (identifying the indices of the current grid point)
   * @param alpha the spars grid's ansatzfunctions coefficients
   * @param result reference to a float_t into which the result should be stored
   */
  void rec(BASIS& basis, const DataVector& point, size_t current_dim,
           float_t value, GridStorage::grid_iterator& working,
           GridStorage::index_type::index_type* source, const DataVector& alpha,
           float_t& result) {
    typedef GridStorage::index_type::level_type level_type;
    typedef GridStorage::index_type::index_type index_type;

    const unsigned int BITS_IN_BYTE = 8;
    // maximum possible level for the index type
    const level_type max_level = static_cast<level_type> (
                                   sizeof(index_type) * BITS_IN_BYTE - 1);
    index_type src_index = source[current_dim];

    level_type work_level = 1;

    while (true) {
      size_t seq = working.seq();

      if (storage->end(seq)) {
        break;
      } else {
        index_type work_index;
        level_type temp;

        working.get(current_dim, temp, work_index);

        float_t new_value = basis.eval(work_level, work_index,
                                       point[current_dim]);
        new_value *= value;

        if (current_dim == storage->dim() - 1) {
          result += (alpha[seq] * new_value);
        } else {
          rec(basis, point, current_dim + 1, new_value,
              working, source, alpha, result);
        }
      }

      if (working.hint()) {
        break;
      }

      // this decides in which direction we should descend by evaluating
      // the corresponding bit
      // the bits are coded from left to right starting with level 1
      // being in position max_level
      bool right = (src_index & (1 << (max_level - work_level))) > 0;
      ++work_level;

      if (right) {
        working.rightChild(current_dim);
      } else {
        working.leftChild(current_dim);
      }
    }

    working.resetToLevelOne(current_dim);
  }

/**
 * If during the refinement one or more points of the shadow register are added
 * to the actual grid then we have to remove these points from the shadow storage.
 */
void PrewaveletGridGenerator::consolidateShadow() {
  GridStorage* temp = new GridStorage(storage->dim());

  for (size_t i = 0; i < shadowstorage->size(); i++) {
    temp->insert(*shadowstorage->get(i));
  }

  shadowstorage->emptyStorage();

  for (size_t i = 0; i < temp->size(); i++) {
    if (storage->find(temp->get(i)) == storage->end()) {
      shadowstorage->insert(*temp->get(i));
    }
  }

  delete temp;
}

int main() {
  // create a two-dimensional piecewise bilinear grid
  size_t dim = 2;
  Grid* grid = Grid::createLinearGrid(dim);
  GridStorage* gridStorage = grid->getStorage();
  std::cout << "dimensionality:                   " << gridStorage->dim() << std::endl;

  // create regular grid, level 3
  size_t level = 3;
  GridGenerator* gridGen = grid->createGridGenerator();
  gridGen->regular(level);
  std::cout << "number of initial grid points:    " << gridStorage->size() << std::endl;

  // create coefficient vector
  DataVector alpha(gridStorage->size());
  alpha.setAll(0.0);
  std::cout << "length of alpha vector:           " << alpha.getSize() << std::endl;

  GridIndex* gp;

  // refine adaptively 5 times
  for (int step = 0; step < 5; step++) {
    // set function values in alpha
    for (size_t i = 0; i < gridStorage->size(); i++) {
      gp = gridStorage->get(i);
      alpha[i] = f(gp->getCoord(0), gp->getCoord(1));
    }

    // hierarchize
    SGPP::op_factory::createOperationHierarchisation(*grid)->doHierarchisation(
      alpha);

    // refine a single grid point each time
    gridGen->refine(new SurplusRefinementFunctor(&alpha, 1));
    std::cout << "refinement step " << step + 1 << ", new grid size: " << alpha.getSize()
         << std::endl;

    // extend alpha vector (new entries uninitialized)
    alpha.resize(gridStorage->size());
  }

  delete grid;
}

 virtual typename GridStorageAccessInterface<TBASE>::iterator begin()
 {
     return AccessIteratorImpl<T, TBASE, std::vector<T> >(storage->begin());
 }

 virtual typename GridStorageAccessInterface<TBASE>::const_iterator end() const
 {
     return ConstAccessIteratorImpl<T, TBASE, std::vector<T> >(storage->end());
 }

void testHierarchisationDehierarchisation(SGPP::base::Grid* grid, size_t level,
                                          SGPP::float_t (*func)(DataVector&),
                                          SGPP::float_t tolerance = 0.0, bool naiveOp = false,
                                          bool doStretch = false) {
  GridGenerator* gridGen = grid->createGridGenerator();
  gridGen->regular(level);
  GridStorage* gridStore = grid->getStorage();
  size_t dim = gridStore->dim();
  Stretching* stretch = NULL;

  if (doStretch) stretch = gridStore->getStretching();

  DataVector node_values = DataVector(gridStore->size());
  DataVector coords = DataVector(dim);

  for (size_t n = 0; n < gridStore->size(); n++) {
    if (doStretch) {
      gridStore->get(n)->getCoordsStretching(coords, *stretch);
    } else {
      gridStore->get(n)->getCoords(coords);
    }

    node_values[n] = func(coords);
  }

  DataVector alpha = DataVector(node_values);
  OperationHierarchisation* hierarchisation =
      SGPP::op_factory::createOperationHierarchisation(*grid);
  hierarchisation->doHierarchisation(alpha);

  if (naiveOp == true) {
    OperationNaiveEval* op = SGPP::op_factory::createOperationNaiveEval(*grid);

    for (size_t n = 0; n < gridStore->size(); n++) {
      if (doStretch) {
        gridStore->get(n)->getCoordsStretching(coords, *stretch);
      } else {
        gridStore->get(n)->getCoords(coords);
      }

      SGPP::float_t eval = op->eval(alpha, coords);
      BOOST_CHECK_CLOSE(eval, node_values[n], tolerance);
    }
  } else {
    OperationEval* op = SGPP::op_factory::createOperationEval(*grid);

    for (size_t n = 0; n < gridStore->size(); n++) {
      if (doStretch) {
        gridStore->get(n)->getCoordsStretching(coords, *stretch);
      } else {
        gridStore->get(n)->getCoords(coords);
      }

      SGPP::float_t eval = op->eval(alpha, coords);
      BOOST_CHECK_CLOSE(eval, node_values[n], tolerance);
    }
  }

  DataVector node_values_back = DataVector(alpha);
  hierarchisation->doDehierarchisation(node_values_back);

  for (size_t n = 0; n < gridStore->size(); n++) {
    BOOST_CHECK_CLOSE(node_values_back[n], node_values[n], tolerance);
  }
}

 /**
  * Create a new sweep object with a copied functor
  *
  * @param functor the functor that is executed on the grid
  * @param storage the storage that contains the grid points
  */
 sweep(FUNC& functor, GridStorage* storage) :
   functor(functor), storage(storage),
   algoDims(storage->getAlgorithmicDimensions()),
   numAlgoDims_(storage->getAlgorithmicDimensions().size()) {
 }

  /**
   * Descends on all dimensions beside dim_sweep. Class functor for dim_sweep.
   * Boundaries are not regarded
   *
   * @param source coefficients of the sparse grid
   * @param result coefficients of the function computed by sweep
   * @param index current grid position
   * @param dim_list list of dimensions, that should be handled
   * @param dim_rem number of remaining dims
   * @param dim_sweep static dimension, in this dimension the functor is executed
  */
  void sweep_rec(DataMatrix& source, DataMatrix& result, grid_iterator& index,
                 std::vector<size_t>& dim_list, size_t dim_rem,
                 size_t dim_sweep) {
    functor(source, result, index, dim_sweep);

    // dimension recursion unrolled
    for (size_t d = 0; d < dim_rem; d++) {
      size_t current_dim = dim_list[d];

      if (index.hint()) {
        continue;
      }

      index.leftChild(current_dim);

      if (!storage->end(index.seq())) {
        sweep_rec(source, result, index, dim_list, d + 1, dim_sweep);
      }

      index.stepRight(current_dim);

      if (!storage->end(index.seq())) {
        sweep_rec(source, result, index, dim_list, d + 1, dim_sweep);
      }

      index.up(current_dim);
    }
  }

  /**
   * Descends on all dimensions beside dim_sweep. Class functor for dim_sweep
   * Boundaries are not regarded
   *
   * @param source a DataVector containing the source coefficients of the grid points
   * @param result a DataVector containing the result coefficients of the grid points
   * @param dim_sweep the dimension in which the functor is executed
   */
  void sweep1D(DataVector& source, DataVector& result, size_t dim_sweep) {
    // generate a list of all dimension (-dim_sweep)
    // from dimension recursion unrolling
    std::vector<size_t> dim_list;

    for (size_t i = 0; i < storage->dim(); i++) {
      if (i != dim_sweep) {
        dim_list.push_back(i);
      }
    }

    grid_iterator index(storage);

    sweep_rec(source, result, index, dim_list, storage->dim() - 1, dim_sweep);
  }

 virtual void clear()
 {
     storage->clear();
 };

 virtual const Vector2ui &getNumCells() const
 {
     return storage->getNumCells();
 };

 virtual TBASE& at(size_t x, size_t y)
 {
     return static_cast<TBASE &>(storage->at(x, y));
 }

 virtual void moveBy(Index idx)
 {
     storage->moveBy(idx);
 }

 virtual void resize(Vector2ui newSize)
 {
     storage->resize(newSize);
 };