connection_creator_impl.h

Go to the documentation of this file.

/*
 *  connection_creator_impl.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 CONNECTION_CREATOR_IMPL_H
#define CONNECTION_CREATOR_IMPL_H
 
#include "connection_creator.h"
 
// C++ includes:
#include <vector>
 
// Includes from nestkernel:
#include "kernel_manager.h"
 
namespace nest
{
template < int D >
void
ConnectionCreator::connect( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  switch ( type_ )
  {
  case Pairwise_bernoulli_on_source:
 
    pairwise_bernoulli_on_source_( source, source_nc, target, target_nc );
    break;
 
  case Fixed_indegree:
 
    fixed_indegree_( source, source_nc, target, target_nc );
    break;
 
  case Fixed_outdegree:
 
    fixed_outdegree_( source, source_nc, target, target_nc );
    break;
 
  case Pairwise_bernoulli_on_target:
 
    pairwise_bernoulli_on_target_( source, source_nc, target, target_nc );
    break;
 
  case Pairwise_poisson:
 
    pairwise_poisson_( source, source_nc, target, target_nc );
    break;
 
  default:
    throw BadProperty( "Unknown connection type." );
  }
}
 
template < typename Iterator, int D >
void
ConnectionCreator::connect_to_target_( Iterator from,
  Iterator to,
  Node* tgt_ptr,
  const Position< D >& tgt_pos,
  size_t tgt_thread,
  const Layer< D >& source )
{
  RngPtr rng = get_vp_specific_rng( tgt_thread );
 
  // We create a source pos vector here that can be updated with the
  // source position. This is done to avoid creating and destroying
  // unnecessarily many vectors.
  std::vector< double > source_pos( D );
  const std::vector< double > target_pos = tgt_pos.get_vector();
 
  for ( Iterator iter = from; iter != to; ++iter )
  {
    if ( not allow_autapses_ and ( iter->second == tgt_ptr->get_node_id() ) )
    {
      continue;
    }
    iter->first.get_vector( source_pos );
 
    if ( not kernel_ or rng->drand() < kernel_->value( rng, source_pos, target_pos, source, tgt_ptr ) )
    {
      for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
      {
        kernel().connection_manager.connect( iter->second,
          tgt_ptr,
          tgt_thread,
          synapse_model_[ indx ],
          param_dicts_[ indx ][ tgt_thread ],
          delay_[ indx ]->value( rng, source_pos, target_pos, source, tgt_ptr ),
          weight_[ indx ]->value( rng, source_pos, target_pos, source, tgt_ptr ) );
      }
    }
  }
}
 
template < typename Iterator, int D >
void
ConnectionCreator::connect_to_target_poisson_( Iterator from,
  Iterator to,
  Node* tgt_ptr,
  const Position< D >& tgt_pos,
  size_t tgt_thread,
  const Layer< D >& source )
{
  RngPtr rng = get_vp_specific_rng( tgt_thread );
 
  // We create a source pos vector here that can be updated with the
  // source position. This is done to avoid creating and destroying
  // unnecessarily many vectors.
  std::vector< double > source_pos( D );
  const std::vector< double > target_pos = tgt_pos.get_vector();
  poisson_distribution poi_dist;
 
  for ( Iterator iter = from; iter != to; ++iter )
  {
    if ( not allow_autapses_ and ( iter->second == tgt_ptr->get_node_id() ) )
    {
      continue;
    }
    iter->first.get_vector( source_pos );
 
    // Sample number of connections that are to be established
    poisson_distribution::param_type param( kernel_->value( rng, source_pos, target_pos, source, tgt_ptr ) );
    const unsigned long num_conns = poi_dist( rng, param );
    for ( unsigned long conn_counter = 0; conn_counter < num_conns; ++conn_counter )
    {
      for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
      {
        kernel().connection_manager.connect( iter->second,
          tgt_ptr,
          tgt_thread,
          synapse_model_[ indx ],
          param_dicts_[ indx ][ tgt_thread ],
          delay_[ indx ]->value( rng, source_pos, target_pos, source, tgt_ptr ),
          weight_[ indx ]->value( rng, source_pos, target_pos, source, tgt_ptr ) );
      }
    }
  }
}
 
template < int D >
ConnectionCreator::PoolWrapper_< D >::PoolWrapper_()
  : masked_layer_( 0 )
  , positions_( 0 )
{
}
 
template < int D >
ConnectionCreator::PoolWrapper_< D >::~PoolWrapper_()
{
  if ( masked_layer_ )
  {
    delete masked_layer_;
  }
}
 
template < int D >
void
ConnectionCreator::PoolWrapper_< D >::define( MaskedLayer< D >* ml )
{
  assert( masked_layer_ == 0 );
  assert( positions_ == 0 );
  assert( ml != 0 );
  masked_layer_ = ml;
}
 
template < int D >
void
ConnectionCreator::PoolWrapper_< D >::define( std::vector< std::pair< Position< D >, size_t > >* pos )
{
  assert( masked_layer_ == 0 );
  assert( positions_ == 0 );
  assert( pos != 0 );
  positions_ = pos;
}
 
template < int D >
typename Ntree< D, size_t >::masked_iterator
ConnectionCreator::PoolWrapper_< D >::masked_begin( const Position< D >& pos ) const
{
  return masked_layer_->begin( pos );
}
 
template < int D >
typename Ntree< D, size_t >::masked_iterator
ConnectionCreator::PoolWrapper_< D >::masked_end() const
{
  return masked_layer_->end();
}
 
template < int D >
typename std::vector< std::pair< Position< D >, size_t > >::iterator
ConnectionCreator::PoolWrapper_< D >::begin() const
{
  return positions_->begin();
}
 
template < int D >
typename std::vector< std::pair< Position< D >, size_t > >::iterator
ConnectionCreator::PoolWrapper_< D >::end() const
{
  return positions_->end();
}
 
 
template < int D >
void
ConnectionCreator::pairwise_bernoulli_on_source_( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  // Connect using pairwise Bernoulli drawing source nodes (target driven)
  // For each local target node:
  //  1. Apply Mask to source layer
  //  2. For each source node: Compute probability, draw random number, make
  //     connection conditionally
 
  // retrieve global positions, either for masked or unmasked pool
  PoolWrapper_< D > pool;
  if ( mask_.get() )  // MaskedLayer will be freed by PoolWrapper d'tor
  {
    pool.define( new MaskedLayer< D >( source, mask_, allow_oversized_, source_nc ) );
  }
  else
  {
    pool.define( source.get_global_positions_vector( source_nc ) );
  }
 
  std::vector< std::exception_ptr > exceptions_raised_( kernel().vp_manager.get_num_threads() );
 
#pragma omp parallel
  {
    const int tid = kernel().vp_manager.get_thread_id();
    try
    {
      NodeCollection::const_iterator target_begin = target_nc->begin();
      NodeCollection::const_iterator target_end = target_nc->end();
 
      for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
      {
        Node* const tgt = kernel().node_manager.get_node_or_proxy( ( *tgt_it ).node_id, tid );
 
        if ( not tgt->is_proxy() )
        {
          const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
          if ( mask_.get() )
          {
            connect_to_target_( pool.masked_begin( target_pos ), pool.masked_end(), tgt, target_pos, tid, source );
          }
          else
          {
            connect_to_target_( pool.begin(), pool.end(), tgt, target_pos, tid, source );
          }
        }
      }  // for target_begin
    }
    catch ( ... )
    {
      // Capture the current exception object and create an std::exception_ptr
      exceptions_raised_.at( tid ) = std::current_exception();
    }
  }  // omp parallel
 
  // check if any exceptions have been raised
  for ( auto eptr : exceptions_raised_ )
  {
    if ( eptr )
    {
      std::rethrow_exception( eptr );
    }
  }
}
 
 
template < int D >
void
ConnectionCreator::pairwise_bernoulli_on_target_( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  // Connecting using pairwise Bernoulli drawing target nodes (source driven)
  // It is actually implemented as pairwise Bernoulli on source nodes,
  // but with displacements computed in the target layer. The Mask has been
  // reversed so that it can be applied to the source instead of the target.
  // For each local target node:
  //  1. Apply (Converse)Mask to source layer
  //  2. For each source node: Compute probability, draw random number, make
  //     connection conditionally
 
  PoolWrapper_< D > pool;
  if ( mask_.get() )  // MaskedLayer will be freed by PoolWrapper d'tor
  {
    // By supplying the target layer to the MaskedLayer constructor, the
    // mask is mirrored so it may be applied to the source layer instead
    pool.define( new MaskedLayer< D >( source, mask_, allow_oversized_, target, source_nc ) );
  }
  else
  {
    pool.define( source.get_global_positions_vector( source_nc ) );
  }
 
  std::vector< std::exception_ptr > exceptions_raised_( kernel().vp_manager.get_num_threads() );
 
  // We only need to check the first in the NodeCollection
  Node* const first_in_tgt = kernel().node_manager.get_node_or_proxy( target_nc->operator[]( 0 ) );
  if ( not first_in_tgt->has_proxies() )
  {
    throw IllegalConnection( "Spatial Connect with pairwise_bernoulli to devices is not possible." );
  }
 
#pragma omp parallel
  {
    const int tid = kernel().vp_manager.get_thread_id();
    try
    {
      NodeCollection::const_iterator target_begin = target_nc->thread_local_begin();
      NodeCollection::const_iterator target_end = target_nc->end();
 
      for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
      {
        Node* const tgt = kernel().node_manager.get_node_or_proxy( ( *tgt_it ).node_id, tid );
 
        assert( not tgt->is_proxy() );
 
        const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
        if ( mask_.get() )
        {
          // We do the same as in the target driven case, except that we calculate displacements in the target layer.
          // We therefore send in target as last parameter.
          connect_to_target_( pool.masked_begin( target_pos ), pool.masked_end(), tgt, target_pos, tid, target );
        }
        else
        {
          // We do the same as in the target driven case, except that we calculate displacements in the target layer.
          // We therefore send in target as last parameter.
          connect_to_target_( pool.begin(), pool.end(), tgt, target_pos, tid, target );
        }
 
      }  // end for
    }
    catch ( ... )
    {
      // Capture the current exception object and create an std::exception_ptr
      exceptions_raised_.at( tid ) = std::current_exception();
    }
  }  // omp parallel
 
  // check if any exceptions have been raised
  for ( auto eptr : exceptions_raised_ )
  {
    if ( eptr )
    {
      std::rethrow_exception( eptr );
    }
  }
}
 
 
template < int D >
void
ConnectionCreator::pairwise_poisson_( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  // Connect using pairwise Poisson drawing source nodes (target driven)
  // For each local target node:
  //  1. Apply Mask to source layer
  //  2. For each source node: Compute probability, draw random number, make
  //     connection conditionally
 
  // retrieve global positions, either for masked or unmasked pool
  PoolWrapper_< D > pool;
  if ( mask_.get() )  // MaskedLayer will be freed by PoolWrapper d'tor
  {
    pool.define( new MaskedLayer< D >( source, mask_, allow_oversized_, source_nc ) );
  }
  else
  {
    pool.define( source.get_global_positions_vector( source_nc ) );
  }
 
  std::vector< std::exception_ptr > exceptions_raised_( kernel().vp_manager.get_num_threads() );
 
#pragma omp parallel
  {
    const int thread_id = kernel().vp_manager.get_thread_id();
    try
    {
      NodeCollection::const_iterator target_begin = target_nc->begin();
      NodeCollection::const_iterator target_end = target_nc->end();
 
      for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
      {
        Node* const tgt = kernel().node_manager.get_node_or_proxy( ( *tgt_it ).node_id, thread_id );
 
        if ( not tgt->is_proxy() )
        {
          const Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
          if ( mask_.get() )
          {
            connect_to_target_poisson_(
              pool.masked_begin( target_pos ), pool.masked_end(), tgt, target_pos, thread_id, source );
          }
          else
          {
            connect_to_target_poisson_( pool.begin(), pool.end(), tgt, target_pos, thread_id, source );
          }
        }
      }  // for target_begin
    }
    catch ( ... )
    {
      exceptions_raised_.at( thread_id ) = std::current_exception();
    }
  }  // omp parallel
 
  // check if any exceptions have been raised
  for ( auto eptr : exceptions_raised_ )
  {
    if ( eptr )
    {
      std::rethrow_exception( eptr );
    }
  }
}
 
 
template < int D >
void
ConnectionCreator::fixed_indegree_( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  // fixed_indegree connections (fixed fan in)
  //
  // For each local target node:
  // 1. Apply Mask to source layer
  // 2. Compute connection probability for each source position
  // 3. Draw source nodes and make connections
 
  // We only need to check the first in the NodeCollection
  Node* const first_in_tgt = kernel().node_manager.get_node_or_proxy( target_nc->operator[]( 0 ) );
  if ( not first_in_tgt->has_proxies() )
  {
    throw IllegalConnection( "Spatial Connect with fixed_indegree to devices is not possible." );
  }
 
  NodeCollection::const_iterator target_begin = target_nc->rank_local_begin();
  NodeCollection::const_iterator target_end = target_nc->end();
 
  // protect against connecting to devices without proxies
  // we need to do this before creating the first connection to leave
  // the network untouched if any target does not have proxies
  for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
  {
    Node* const tgt = kernel().node_manager.get_node_or_proxy( ( *tgt_it ).node_id );
 
    assert( not tgt->is_proxy() );
  }
 
  if ( mask_.get() )
  {
    MaskedLayer< D > masked_source( source, mask_, allow_oversized_, source_nc );
    const auto masked_source_end = masked_source.end();
 
    std::vector< std::pair< Position< D >, size_t > > positions;
 
    for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
    {
      size_t target_id = ( *tgt_it ).node_id;
      Node* const tgt = kernel().node_manager.get_node_or_proxy( target_id );
 
      size_t target_thread = tgt->get_thread();
      RngPtr rng = get_vp_specific_rng( target_thread );
      Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
      // We create a source pos vector here that can be updated with the
      // source position. This is done to avoid creating and destroying
      // unnecessarily many vectors.
      std::vector< double > source_pos_vector( D );
      const std::vector< double > target_pos_vector = target_pos.get_vector();
 
      unsigned long target_number_connections =
        std::round( number_of_connections_->value( rng, source_pos_vector, target_pos_vector, source, tgt ) );
 
      // Get (position,node ID) pairs for sources inside mask
      positions.resize( std::distance( masked_source.begin( target_pos ), masked_source_end ) );
      std::copy( masked_source.begin( target_pos ), masked_source_end, positions.begin() );
 
      // We will select `number_of_connections_` sources within the mask.
      // If there is no kernel, we can just draw uniform random numbers,
      // but with a kernel we have to set up a probability distribution
      // function using a discrete_distribution.
      if ( kernel_ )
      {
 
        std::vector< double > probabilities;
        probabilities.reserve( positions.size() );
 
        // Collect probabilities for the sources
        for ( typename std::vector< std::pair< Position< D >, size_t > >::iterator iter = positions.begin();
          iter != positions.end();
          ++iter )
        {
          iter->first.get_vector( source_pos_vector );
          probabilities.push_back( kernel_->value( rng, source_pos_vector, target_pos_vector, source, tgt ) );
        }
 
        if ( positions.empty()
          or ( not allow_autapses_ and ( positions.size() == 1 ) and positions[ 0 ].second == target_id )
          or ( not allow_multapses_ and ( positions.size() < target_number_connections ) ) )
        {
          std::string msg = String::compose( "Global target ID %1: Not enough sources found inside mask", target_id );
          throw KernelException( msg.c_str() );
        }
 
        // A discrete_distribution draws random integers with a non-uniform
        // distribution.
        discrete_distribution lottery;
        const discrete_distribution::param_type param( probabilities.begin(), probabilities.end() );
        lottery.param( param );
 
        // If multapses are not allowed, we must keep track of which
        // sources have been selected already.
        std::vector< bool > is_selected( positions.size() );
 
        // Draw `target_number_connections` sources
        while ( target_number_connections > 0 )
        {
          size_t random_id = lottery( rng );
          if ( not allow_multapses_ and is_selected[ random_id ] )
          {
            continue;
          }
 
          size_t source_id = positions[ random_id ].second;
          if ( not allow_autapses_ and source_id == target_id )
          {
            continue;
          }
          positions[ random_id ].first.get_vector( source_pos_vector );
          for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
          {
            const double w = weight_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            const double d = delay_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            kernel().connection_manager.connect(
              source_id, tgt, target_thread, synapse_model_[ indx ], param_dicts_[ indx ][ target_thread ], d, w );
          }
 
          is_selected[ random_id ] = true;
          --target_number_connections;
        }
      }
      else
      {
 
        // no kernel
 
        if ( positions.empty()
          or ( not allow_autapses_ and ( positions.size() == 1 ) and positions[ 0 ].second == target_id )
          or ( not allow_multapses_ and ( positions.size() < target_number_connections ) ) )
        {
          std::string msg = String::compose( "Global target ID %1: Not enough sources found inside mask", target_id );
          throw KernelException( msg.c_str() );
        }
 
        // If multapses are not allowed, we must keep track of which
        // sources have been selected already.
        std::vector< bool > is_selected( positions.size() );
 
        // Draw `target_number_connections` sources
        while ( target_number_connections > 0 )
        {
          const size_t random_id = rng->ulrand( positions.size() );
          if ( not allow_multapses_ and is_selected[ random_id ] )
          {
            continue;
          }
          positions[ random_id ].first.get_vector( source_pos_vector );
          const size_t source_id = positions[ random_id ].second;
          for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
          {
            const double w = weight_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            const double d = delay_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            kernel().connection_manager.connect(
              source_id, tgt, target_thread, synapse_model_[ indx ], param_dicts_[ indx ][ target_thread ], d, w );
          }
 
          is_selected[ random_id ] = true;
          --target_number_connections;
        }
      }
    }
  }
  else
  {
    // no mask
 
    // Get (position,node ID) pairs for all nodes in source layer
    std::vector< std::pair< Position< D >, size_t > >* positions = source.get_global_positions_vector( source_nc );
 
    for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
    {
      size_t target_id = ( *tgt_it ).node_id;
      Node* const tgt = kernel().node_manager.get_node_or_proxy( target_id );
      size_t target_thread = tgt->get_thread();
      RngPtr rng = get_vp_specific_rng( target_thread );
      Position< D > target_pos = target.get_position( ( *tgt_it ).nc_index );
 
      unsigned long target_number_connections = std::round( number_of_connections_->value( rng, tgt ) );
 
      std::vector< double > source_pos_vector( D );
      const std::vector< double > target_pos_vector = target_pos.get_vector();
 
      if ( ( positions->size() == 0 )
        or ( not allow_autapses_ and ( positions->size() == 1 ) and ( ( *positions )[ 0 ].second == target_id ) )
        or ( not allow_multapses_ and ( positions->size() < target_number_connections ) ) )
      {
        std::string msg = String::compose( "Global target ID %1: Not enough sources found", target_id );
        throw KernelException( msg.c_str() );
      }
 
      // We will select `target_number_connections` sources within the mask.
      // If there is no kernel, we can just draw uniform random numbers,
      // but with a kernel we have to set up a probability distribution
      // function using a discrete_distribution.
      if ( kernel_ )
      {
 
        std::vector< double > probabilities;
        probabilities.reserve( positions->size() );
 
        // Collect probabilities for the sources
        for ( typename std::vector< std::pair< Position< D >, size_t > >::iterator iter = positions->begin();
          iter != positions->end();
          ++iter )
        {
          iter->first.get_vector( source_pos_vector );
          probabilities.push_back( kernel_->value( rng, source_pos_vector, target_pos_vector, source, tgt ) );
        }
 
        // A discrete_distribution draws random integers with a non-uniform
        // distribution.
        discrete_distribution lottery;
        const discrete_distribution::param_type param( probabilities.begin(), probabilities.end() );
        lottery.param( param );
 
        // If multapses are not allowed, we must keep track of which
        // sources have been selected already.
        std::vector< bool > is_selected( positions->size() );
 
        // Draw `target_number_connections` sources
        while ( target_number_connections > 0 )
        {
          const size_t random_id = lottery( rng );
          if ( not allow_multapses_ and is_selected[ random_id ] )
          {
            continue;
          }
 
          const size_t source_id = ( *positions )[ random_id ].second;
          if ( not allow_autapses_ and source_id == target_id )
          {
            continue;
          }
 
          ( *positions )[ random_id ].first.get_vector( source_pos_vector );
          for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
          {
            const double w = weight_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            const double d = delay_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            kernel().connection_manager.connect(
              source_id, tgt, target_thread, synapse_model_[ indx ], param_dicts_[ indx ][ target_thread ], d, w );
          }
 
          is_selected[ random_id ] = true;
          --target_number_connections;
        }
      }
      else
      {
 
        // no kernel
 
        // If multapses are not allowed, we must keep track of which
        // sources have been selected already.
        std::vector< bool > is_selected( positions->size() );
 
        // Draw `target_number_connections` sources
        while ( target_number_connections > 0 )
        {
          const size_t random_id = rng->ulrand( positions->size() );
          if ( not allow_multapses_ and is_selected[ random_id ] )
          {
            continue;
          }
 
          const size_t source_id = ( *positions )[ random_id ].second;
          if ( not allow_autapses_ and source_id == target_id )
          {
            continue;
          }
 
          ( *positions )[ random_id ].first.get_vector( source_pos_vector );
          for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
          {
            const double w = weight_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            const double d = delay_[ indx ]->value( rng, source_pos_vector, target_pos_vector, source, tgt );
            kernel().connection_manager.connect(
              source_id, tgt, target_thread, synapse_model_[ indx ], param_dicts_[ indx ][ target_thread ], d, w );
          }
 
          is_selected[ random_id ] = true;
          --target_number_connections;
        }
      }
    }
  }
}
 
 
template < int D >
void
ConnectionCreator::fixed_outdegree_( Layer< D >& source,
  NodeCollectionPTR source_nc,
  Layer< D >& target,
  NodeCollectionPTR target_nc )
{
  // protect against connecting to devices without proxies
  // we need to do this before creating the first connection to leave
  // the network untouched if any target does not have proxies
 
  // We only need to check the first in the NodeCollection
  Node* const first_in_tgt = kernel().node_manager.get_node_or_proxy( target_nc->operator[]( 0 ) );
  if ( not first_in_tgt->has_proxies() )
  {
    throw IllegalConnection( "Spatial Connect with fixed_outdegree to devices is not possible." );
  }
 
  NodeCollection::const_iterator target_begin = target_nc->rank_local_begin();
  NodeCollection::const_iterator target_end = target_nc->end();
 
  for ( NodeCollection::const_iterator tgt_it = target_begin; tgt_it < target_end; ++tgt_it )
  {
    Node* const tgt = kernel().node_manager.get_node_or_proxy( ( *tgt_it ).node_id );
 
    assert( not tgt->is_proxy() );
  }
 
  // Fixed_outdegree connections (fixed fan out)
  //
  // For each (global) source: (All connections made on all mpi procs)
  // 1. Apply mask to global targets
  // 2. If using kernel: Compute connection probability for each global target
  // 3. Draw connections to make using global rng
 
  MaskedLayer< D > masked_target( target, mask_, allow_oversized_, target_nc );
  const auto masked_target_end = masked_target.end();
 
  // We create a target positions vector here that can be updated with the
  // position and node ID pairs. This is done to avoid creating and destroying
  // unnecessarily many vectors.
  std::vector< std::pair< Position< D >, size_t > > target_pos_node_id_pairs;
  std::vector< std::pair< Position< D >, size_t > > source_pos_node_id_pairs =
    *source.get_global_positions_vector( source_nc );
 
  for ( const auto& source_pos_node_id_pair : source_pos_node_id_pairs )
  {
    const Position< D > source_pos = source_pos_node_id_pair.first;
    const size_t source_id = source_pos_node_id_pair.second;
    const auto src = kernel().node_manager.get_node_or_proxy( source_id );
    const std::vector< double > source_pos_vector = source_pos.get_vector();
 
    // We create a target pos vector here that can be updated with the
    // target position. This is done to avoid creating and destroying
    // unnecessarily many vectors.
    std::vector< double > target_pos_vector( D );
    std::vector< double > probabilities;
 
    // Find potential targets and probabilities
    RngPtr grng = get_rank_synced_rng();
    target_pos_node_id_pairs.resize( std::distance( masked_target.begin( source_pos ), masked_target_end ) );
    std::copy( masked_target.begin( source_pos ), masked_target_end, target_pos_node_id_pairs.begin() );
 
    probabilities.reserve( target_pos_node_id_pairs.size() );
    if ( kernel_ )
    {
      for ( const auto& target_pos_node_id_pair : target_pos_node_id_pairs )
      {
        // TODO: Why is probability calculated in source layer, but weight and delay in target layer?
        target_pos_node_id_pair.first.get_vector( target_pos_vector );
        const auto tgt = kernel().node_manager.get_node_or_proxy( target_pos_node_id_pair.second );
        probabilities.push_back( kernel_->value( grng, source_pos_vector, target_pos_vector, source, tgt ) );
      }
    }
    else
    {
      probabilities.resize( target_pos_node_id_pairs.size(), 1.0 );
    }
 
    unsigned long number_of_connections = std::round( number_of_connections_->value( grng, src ) );
 
    if ( target_pos_node_id_pairs.empty()
      or ( not allow_multapses_ and target_pos_node_id_pairs.size() < number_of_connections ) )
    {
      std::string msg = String::compose( "Global source ID %1: Not enough targets found", source_id );
      throw KernelException( msg.c_str() );
    }
 
    // Draw targets.  A discrete_distribution draws random integers with a
    // non-uniform distribution.
    discrete_distribution lottery;
    const discrete_distribution::param_type param( probabilities.begin(), probabilities.end() );
    lottery.param( param );
 
    // If multapses are not allowed, we must keep track of which
    // targets have been selected already.
    std::vector< bool > is_selected( target_pos_node_id_pairs.size() );
 
    // Draw `number_of_connections` targets
    while ( number_of_connections > 0 )
    {
      const size_t random_id = lottery( get_rank_synced_rng() );
      if ( not allow_multapses_ and is_selected[ random_id ] )
      {
        continue;
      }
      const size_t target_id = target_pos_node_id_pairs[ random_id ].second;
      if ( not allow_autapses_ and source_id == target_id )
      {
        continue;
      }
 
      is_selected[ random_id ] = true;
 
      target_pos_node_id_pairs[ random_id ].first.get_vector( target_pos_vector );
 
      std::vector< double > rng_weight_vec;
      std::vector< double > rng_delay_vec;
      for ( size_t indx = 0; indx < weight_.size(); ++indx )
      {
        const auto tgt = kernel().node_manager.get_node_or_proxy( target_pos_node_id_pairs[ indx ].second );
        rng_weight_vec.push_back( weight_[ indx ]->value( grng, source_pos_vector, target_pos_vector, target, tgt ) );
        rng_delay_vec.push_back( delay_[ indx ]->value( grng, source_pos_vector, target_pos_vector, target, tgt ) );
      }
 
      // Each VP has now decided to create this connection and drawn any random parameter values
      // required for it. Each VP thus counts the connection as created, but only the VP hosting the
      // target neuron actually creates the connection.
      --number_of_connections;
      if ( not kernel().node_manager.is_local_node_id( target_id ) )
      {
        continue;
      }
 
      Node* target_ptr = kernel().node_manager.get_node_or_proxy( target_id );
      const size_t target_thread = target_ptr->get_thread();
 
      for ( size_t indx = 0; indx < synapse_model_.size(); ++indx )
      {
        kernel().connection_manager.connect( source_id,
          target_ptr,
          target_thread,
          synapse_model_[ indx ],
          param_dicts_[ indx ][ target_thread ],
          rng_delay_vec[ indx ],
          rng_weight_vec[ indx ] );
      }
    }
  }
}
 
}  // namespace nest
 
#endif