grid_layer.h

Go to the documentation of this file.

/*
 *  grid_layer.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#ifndef GRID_LAYER_H
#define GRID_LAYER_H

#include "layer.h"

namespace nest
{

  template <int D>
  class GridLayer: public Layer<D>
  {
  public:
    typedef Position<D> key_type;
    typedef index mapped_type;
    typedef std::pair<Position<D>,index> value_type;
    typedef value_type& reference;
    typedef const value_type& const_reference;

    class masked_iterator {
    public:
      masked_iterator(const GridLayer<D> & layer) :
        layer_(layer), node_(), depth_(-1)
        {}

      masked_iterator(const GridLayer<D> & layer, const Mask<D> &mask, const Position<D> &anchor, const Selector & filter);

      value_type operator*();

      masked_iterator & operator++();

      masked_iterator operator++(int)
        {
          masked_iterator tmp = *this;
          ++*this;
          return tmp;
        }

      bool operator==(const masked_iterator &other) const
        { return (other.layer_.get_gid()==layer_.get_gid()) && (other.node_==node_) && (other.depth_==depth_); }
      bool operator!=(const masked_iterator &other) const
        { return (other.layer_.get_gid()!=layer_.get_gid()) || (other.node_!=node_) || (other.depth_!=depth_); }

    protected:

      const GridLayer<D> & layer_;
      int_t layer_size_;
      const Mask<D> * mask_;
      Position<D> anchor_;
      Selector filter_;
      MultiIndex<D> node_;
      int_t depth_;
    };

    GridLayer():
      Layer<D>()
      {}

    GridLayer(const GridLayer &l):
      Layer<D>(l),
      dims_(l.dims_)
      {}

    Position<D> get_position(index sind) const;

    Position<D> lid_to_position(index lid) const;

    index gridpos_to_lid(Position<D,int_t> pos) const;

    Position<D> gridpos_to_position(Position<D,int_t> gridpos) const;

    std::vector<index> get_nodes(Position<D,int_t> pos);

    using Layer<D>::get_global_positions_vector;

    std::vector<std::pair<Position<D>,index> > get_global_positions_vector(Selector filter, const AbstractMask& mask, const Position<D>& anchor, bool allow_oversized);

    masked_iterator masked_begin(const Mask<D> &mask, const Position<D> &anchor, const Selector & filter);
    masked_iterator masked_end();

    Position<D,index> get_dims() const;

    void set_status(const DictionaryDatum &d);
    void get_status(DictionaryDatum &d) const;

  protected:
    Position<D,index> dims_;   

    template<class Ins>
    void insert_global_positions_(Ins iter, const Selector& filter);
    void insert_global_positions_ntree_(Ntree<D,index> & tree, const Selector& filter);
    void insert_global_positions_vector_(std::vector<std::pair<Position<D>,index> > & vec, const Selector& filter);
    void insert_local_positions_ntree_(Ntree<D,index> & tree, const Selector& filter);
  };

  template <int D>
  Position<D,index> GridLayer<D>::get_dims() const
  {
    return dims_;
  }

  template <int D>
  void GridLayer<D>::set_status(const DictionaryDatum &d)
  {
    Position<D,index> new_dims = dims_;
    updateValue<long_t>(d, names::columns, new_dims[0]);
    if (D>=2) updateValue<long_t>(d, names::rows, new_dims[1]);
    if (D>=3) updateValue<long_t>(d, names::layers, new_dims[2]);

    index new_size = this->depth_;
    for(int i=0;i<D;++i) {
      new_size *= new_dims[i];
    }

    if (new_size != this->global_size()) {
      throw BadProperty("Total size of layer must be unchanged.");
    }

    this->dims_ = new_dims;

    Layer<D>::set_status(d);
  }

  template <int D>
  void GridLayer<D>::get_status(DictionaryDatum &d) const
  {
    Layer<D>::get_status(d);

    DictionaryDatum topology_dict = getValue<DictionaryDatum>((*d)[names::topology]);

    (*topology_dict)[names::columns] = dims_[0];
    if (D>=2) (*topology_dict)[names::rows] = dims_[1];
    if (D>=3) (*topology_dict)[names::layers] = dims_[2];
  }

  template <int D>
  Position<D> GridLayer<D>::lid_to_position(index lid) const
  {
    lid %= this->global_size()/this->depth_;
    Position<D,int_t> gridpos;
    for(int i=D-1;i>0;--i) {
      gridpos[i] = lid % dims_[i];
      lid = lid / dims_[i];
    }
    assert(lid < dims_[0]);
    gridpos[0] = lid;
    return gridpos_to_position(gridpos);
  }

  template <int D>
  Position<D> GridLayer<D>::gridpos_to_position(Position<D,int_t> gridpos) const
  {
    // grid layer uses "matrix convention", i.e. reversed y axis
    Position<D> ext = this->extent_;
    Position<D> upper_left = this->lower_left_;
    if (D>1) {
      upper_left[1] += ext[1];
      ext[1] = -ext[1];
    }
    return upper_left + ext/dims_ * gridpos + ext/dims_ * 0.5;
  }
    
  template <int D>
  Position<D> GridLayer<D>::get_position(index sind) const
  {
    return lid_to_position(this->nodes_[sind]->get_lid());
  }

  template <int D>
  index GridLayer<D>::gridpos_to_lid(Position<D,int_t> pos) const
  {
    index lid = 0;

    // In case of periodic boundaries, allow grid positions outside layer
    for(int i=0;i<D;++i) {
      if (this->periodic_[i]) {
        pos[i] %= int(dims_[i]);
        if (pos[i]<0) pos[i] += dims_[i];
      }
    }

    for(int i=0;i<D;++i) {
      lid *= dims_[i];
      lid += pos[i];
    }

    return lid;
  }

  template <int D>
  std::vector<index> GridLayer<D>::get_nodes(Position<D,int_t> pos)
  {
    std::vector<index> gids;
    index lid = gridpos_to_lid(pos);
    index layer_size = this->global_size()/this->depth_;

    for(int d=0;d<this->depth_;++d) {
      gids.push_back(this->gids_[lid + d*layer_size]);
    }

    return gids;
  }

  template <int D>
  void GridLayer<D>::insert_local_positions_ntree_(Ntree<D,index> & tree, const Selector& filter)
  {
    std::vector<Node*>::const_iterator nodes_begin;
    std::vector<Node*>::const_iterator nodes_end;

    if (filter.select_depth()) {
      nodes_begin = this->local_begin(filter.depth);
      nodes_end = this->local_end(filter.depth);
    } else {
      nodes_begin = this->local_begin();
      nodes_end = this->local_end();
    }

    for(std::vector<Node*>::const_iterator node_it = nodes_begin; node_it != nodes_end; ++node_it) {

      if (filter.select_model() && ((*node_it)->get_model_id() != filter.model))
        continue;

      tree.insert(std::pair<Position<D>,index>(lid_to_position((*node_it)->get_lid()), (*node_it)->get_gid()));
    }
  }

  template <int D>
  template <class Ins>
  void GridLayer<D>::insert_global_positions_(Ins iter, const Selector& filter)
  {
    index i = 0;
    index lid_end = this->gids_.size();

    if (filter.select_depth()) {
      const index nodes_per_layer = this->gids_.size() / this->depth_;
      i = nodes_per_layer * filter.depth;
      lid_end = nodes_per_layer * (filter.depth+1);
      if ((i >= this->gids_.size()) or (lid_end > this->gids_.size()))
        throw BadProperty("Selected depth out of range");
    }

    Multirange::iterator gi = this->gids_.begin();
    for (index j=0;j<i;++j)  // Advance iterator to first gid at selected depth
      ++gi;

    for(; (gi != this->gids_.end()) && (i < lid_end); ++gi, ++i) {

        if (filter.select_model() && ((int)this->net_->get_model_id_of_gid(*gi) != filter.model))
        continue;

      *iter++ = std::pair<Position<D>,index>(lid_to_position(i), *gi);
    }
  }

  template <int D>
  void GridLayer<D>::insert_global_positions_ntree_(Ntree<D,index> & tree, const Selector& filter)
  {
    insert_global_positions_(std::inserter(tree, tree.end()), filter);
  }

  template <int D>
  void GridLayer<D>::insert_global_positions_vector_(std::vector<std::pair<Position<D>,index> > & vec, const Selector& filter)
  {
    insert_global_positions_(std::back_inserter(vec), filter);
  }

  template <int D>
  inline
  typename GridLayer<D>::masked_iterator GridLayer<D>::masked_begin(const Mask<D> &mask, const Position<D> &anchor, const Selector & filter)
  {
    return masked_iterator(*this, mask, anchor, filter);
  }

  template <int D>
  inline
  typename GridLayer<D>::masked_iterator GridLayer<D>::masked_end()
  {
    return masked_iterator(*this);
  }

  template <int D>
  GridLayer<D>::masked_iterator::masked_iterator(const GridLayer<D> & layer, const Mask<D> &mask, const Position<D> &anchor, const Selector & filter) :
    layer_(layer), mask_(&mask), anchor_(anchor), filter_(filter)
  {
    layer_size_ = layer.global_size()/layer.depth_;

    Position<D,int> ll;
    Position<D,int> ur;
    Box<D> bbox=mask.get_bbox();
    bbox.lower_left += anchor;
    bbox.upper_right += anchor;
    for(int i=0;i<D;++i) {
      if (layer.periodic_[i]) {
        ll[i] = ceil((bbox.lower_left[i] - layer.lower_left_[i])*layer_.dims_[i]/layer.extent_[i] - 0.5);
        ur[i] = round((bbox.upper_right[i] - layer.lower_left_[i])*layer_.dims_[i]/layer.extent_[i]);
      } else {
        ll[i] = std::min(index(std::max(ceil((bbox.lower_left[i] - layer.lower_left_[i])*layer_.dims_[i]/layer.extent_[i] - 0.5), 0.0)), layer.dims_[i]);
        ur[i] = std::min(index(std::max(round((bbox.upper_right[i] - layer.lower_left_[i])*layer_.dims_[i]/layer.extent_[i]), 0.0)), layer.dims_[i]);
      }
    }
    if (D>1) {
      // grid layer uses "matrix convention", i.e. reversed y axis
      int tmp = ll[1];
      ll[1] = layer.dims_[1] - ur[1];
      ur[1] = layer.dims_[1] - tmp;
    }

    node_ = MultiIndex<D>(ll,ur);

    if (filter_.select_depth()) 
      depth_ = filter_.depth;
    else
      depth_ = 0;

    if ((not mask_->inside(layer_.gridpos_to_position(node_)-anchor_)) or
        (filter_.select_model() && (layer_.net_->get_model_id_of_gid(layer_.gids_[depth_*layer_size_]) != index(filter_.model))))
      ++(*this);
  
  }

  template <int D>
  inline
  std::pair<Position<D>,index> GridLayer<D>::masked_iterator::operator*()
  {
    assert(depth_>=0);
    return std::pair<Position<D>,index>(layer_.gridpos_to_position(node_),
                                        layer_.gids_[layer_.gridpos_to_lid(node_) + depth_*layer_size_]);
  }

  template <int D>
  typename GridLayer<D>::masked_iterator & GridLayer<D>::masked_iterator::operator++()
  {
    if (depth_ == -1) return *this;  // Invalid (end) iterator

    if (not filter_.select_depth()) {
      depth_++;
      if (depth_ >= layer_.depth_) {
        depth_ = 0;
      } else {
        if (filter_.select_model() && (layer_.net_->get_model_id_of_gid(layer_.gids_[depth_*layer_size_]) != index(filter_.model)))
          return operator++();
        else
          return *this;
      }
    }

    do {
      ++node_;

      if (node_ == node_.get_upper_right()) {
        // Mark as invalid
        depth_ = -1;
        node_ = MultiIndex<D>();
        return *this;
      }

    } while (not mask_->inside(layer_.gridpos_to_position(node_)-anchor_));

    if (filter_.select_model() && (layer_.net_->get_model_id_of_gid(layer_.gids_[depth_*layer_size_]) != index(filter_.model)))
      return operator++();

    return *this;
  }

  template <int D>
  std::vector<std::pair<Position<D>,index> > GridLayer<D>::get_global_positions_vector(Selector filter, const AbstractMask& mask, const Position<D>& anchor, bool)
  {
    std::vector<std::pair<Position<D>,index> > positions;

    const Mask<D>& mask_d = dynamic_cast<const Mask<D>&>(mask);
    for(typename GridLayer<D>::masked_iterator mi = masked_begin(mask_d,anchor,filter); mi != masked_end(); ++mi) {
      positions.push_back(*mi);
    }

    return positions;
  }

} // namespace nest

#endif