free_layer.h

Go to the documentation of this file.

/*
 *  free_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 FREE_LAYER_H
#define FREE_LAYER_H

#include <limits>
#include <sstream>
#include <algorithm>
#include "layer.h"
#include "topology_names.h"
#include "dictutils.h"
#include "ntree_impl.h"

namespace nest
{

  template <int D>
  class FreeLayer: public Layer<D>
  {
  public:
    Position<D> get_position(index sind) const;
    void set_status(const DictionaryDatum&);
    void get_status(DictionaryDatum&) const;

  protected:

    template <class Ins>
    void communicate_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);

    std::vector<Position<D> > positions_;

    class NodePositionData {
    public:
      index get_gid() const { return gid_; }
      Position<D> get_position() const { return Position<D>(pos_); }
      bool operator<(const NodePositionData& other) const
        { return gid_<other.gid_; }
      bool operator==(const NodePositionData& other) const
        { return gid_==other.gid_; }

    private:
      double_t gid_;
      double_t pos_[D];
    };

  };

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

    // Read positions from dictionary
    if(d->known(names::positions)) {
      TokenArray pos = getValue<TokenArray>(d, names::positions);
      if(this->global_size()/this->depth_ != pos.size()) {
        std::stringstream expected;
        std::stringstream got;
        expected << "position array with length " << this->global_size()/this->depth_;
        got << "position array with length" << pos.size();
        throw TypeMismatch(expected.str(), got.str());
      }

      positions_.clear();
      positions_.reserve(this->local_size());

      const index nodes_per_depth = this->global_size()/this->depth_;
      const index first_lid = this->nodes_[0]->get_lid();

      for(vector<Node*>::iterator i = this->local_begin(); i != this->local_end(); ++i) {

        // Nodes are grouped by depth. When lid % nodes_per_depth ==
        // first_lid, we have "wrapped around", and do not need to gather
        // more positions.
        if ( ((*i)->get_lid() != first_lid) &&
             ((*i)->get_lid() % nodes_per_depth == first_lid ) ) {
          break;
        }

        Position<D> point =
          getValue<std::vector<double_t> >(pos[(*i)->get_lid() % nodes_per_depth]);

        if (not ((point >= this->lower_left_) and (point < this->lower_left_+this->extent_)))
          throw BadProperty("Node position outside of layer");

        positions_.push_back(point);

      }

    }
  }

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

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

    TokenArray points;
    for(typename std::vector<Position<D> >::const_iterator it = positions_.begin();
        it != positions_.end(); ++it) {
      points.push_back(it->getToken());
    }
    def2<TokenArray, ArrayDatum>(topology_dict, names::positions, points);

  }

  template <int D>
  Position<D> FreeLayer<D>::get_position(index sind) const
  {
    // If sind > positions_.size(), we must have "wrapped around" when
    // storing positions, so we may simply mod with the size
    return positions_[sind % positions_.size()];
  }

  template <int D>
  template <class Ins>
  void FreeLayer<D>::communicate_positions_(Ins iter, const Selector& filter)
  {
    assert(this->nodes_.size() >= positions_.size());
    
    // This array will be filled with GID,pos_x,pos_y[,pos_z] for local nodes:
    std::vector<double_t> local_gid_pos;
    std::vector<Node*>::const_iterator nodes_begin;
    std::vector<Node*>::const_iterator nodes_end;

    // Nodes in the subnet are grouped by depth, so to select by depth, we
    // just adjust the begin and end pointers:
    if (filter.select_depth()) {
      local_gid_pos.reserve((D+1)*(this->nodes_.size()/this->depth_ + 1));
      nodes_begin = this->local_begin(filter.depth);
      nodes_end = this->local_end(filter.depth);
    } else {
      local_gid_pos.reserve((D+1)*this->nodes_.size());
      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;

      // Push GID into array to communicate
      local_gid_pos.push_back((*node_it)->get_gid());
      // Push coordinates one by one
      for(int j=0;j<D;++j)
        local_gid_pos.push_back(positions_[(*node_it)->get_subnet_index() % positions_.size()][j]);
    }

    // This array will be filled with GID,pos_x,pos_y[,pos_z] for global nodes:
    std::vector<double_t> global_gid_pos;
    std::vector<int> displacements;
    Communicator::communicate(local_gid_pos,global_gid_pos,displacements);

    // To avoid copying the vector one extra time in order to sort, we
    // sneakishly use reinterpret_cast
    NodePositionData * pos_ptr;
    NodePositionData * pos_end;
    pos_ptr = reinterpret_cast<NodePositionData*>(&global_gid_pos[0]);
    pos_end = pos_ptr + global_gid_pos.size()/(D+1);

    // Get rid of any multiple entries
    std::sort(pos_ptr, pos_end);
    pos_end = std::unique(pos_ptr, pos_end);

    // Unpack GIDs and coordinates
    for(;pos_ptr < pos_end; pos_ptr++) {
      *iter++ = std::pair<Position<D>,index>(pos_ptr->get_position(), pos_ptr->get_gid());
    }

  }

  template <int D>
  void FreeLayer<D>::insert_global_positions_ntree_(Ntree<D,index> & tree, const Selector& filter)
  {

    communicate_positions_(std::inserter(tree, tree.end()), filter);

  }

  template <int D>
  void FreeLayer<D>::insert_local_positions_ntree_(Ntree<D,index> & tree, const Selector& filter)
  {
    assert(this->nodes_.size() >= positions_.size());
    
    std::vector<Node*>::const_iterator nodes_begin;
    std::vector<Node*>::const_iterator nodes_end;

    // Nodes in the subnet are grouped by depth, so to select by depth, we
    // just adjust the begin and end pointers:
    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>(positions_[(*node_it)->get_subnet_index() % positions_.size()],(*node_it)->get_gid()));
    }

  }

  // Helper function to compare GIDs used for sorting (Position,GID) pairs
  template<int D>
  static bool gid_less(const std::pair<Position<D>,index>& a, const std::pair<Position<D>,index>& b)
  {
    return a.second < b.second;
  }

  template <int D>
  void FreeLayer<D>::insert_global_positions_vector_(std::vector<std::pair<Position<D>,index> > & vec, const Selector& filter)
  {

    communicate_positions_(std::back_inserter(vec), filter);

    // Sort vector to ensure consistent results
    std::sort(vec.begin(),vec.end(),gid_less<D>);

  }

} // namespace nest

#endif