Merge pull request #2387 from ackurth/gabor_distribution

Add rectified Gabor distribution for nest.spatial

doc/htmldoc/networks/scripts/connections.py

@@ -408,6 +408,26 @@ def kernel_fig(fig, loc, cdict, kern=None):
 # { end #}
 kernel_fig(fig, 235, conndict)
 
+# { conn42d #}
+conndict = {
+    "rule": "pairwise_bernoulli",
+    "p": nest.spatial_distributions.gaussian2D(nest.spatial.distance.x, nest.spatial.distance.y, std_x=1.0, std_y=3.0),
+    "mask": {"circular": {"radius": 4.0}},
+}
+# { end #}
+# { conn4gab #}
+conndict = {
+    "rule": "pairwise_bernoulli",
+    "p": nest.spatial_distributions.gabor(
+        nest.spatial.target_pos.x - nest.spatial.source_pos.x,
+        nest.spatial.target_pos.y - nest.spatial.source_pos.y,
+        theta=np.pi / 4,
+        gamma=0.7,
+    ),
+}
+# { end #}
+kernel_fig(fig, 236, conndict)
+
 plt.savefig("../user_manual_figures/conn4.png", bbox_inches="tight")
 
 # -----------------------------------------------

doc/htmldoc/networks/spatially_structured_networks.rst

@@ -914,6 +914,16 @@ parameters drawing values from random distributions.
   |                                              | | rho              |    {2(1-\rho^2)}}                                    |
   |                                              |                    |                                                      |
   +----------------------------------------------+--------------------+------------------------------------------------------+
+  |                                              |                    | .. math::                                            |
+  |                                              | | x,               |                                                      |
+  |                                              | | y,               |    p(x) = \big[\cos(360^{\circ}                      |
+  |                                              | | theta,           |    \frac{y^{\prime}}{\lambda}                        |
+  | ``nest.spatial_distributions.gabor()``       | | gamma,           |    + \psi)\big]^{+} e^{-\frac{                       |
+  |                                              | | std,             |    \gamma^{2}x^{\prime 2}+y^{\prime 2}}{             |
+  |                                              | | lam,             |    2\text{std}^{2}}}                                 |
+  |                                              | | psi              |    \\ x^{\prime} = x\cos\theta + y\sin\theta         |
+  |                                              |                    |    \\ y^{\prime} = -x\sin\theta + y\cos\theta        |
+  +----------------------------------------------+--------------------+------------------------------------------------------+
   |                                              |                    | .. math:: p(x) = \frac{x^{\kappa-1}e^{-\frac{x}      |
   | ``nest.spatial_distributions.gamma()``       | | x,               |     {\theta}}}{\theta^\kappa\Gamma(\kappa)}          |
   |                                              | | kappa            |                                                      |
@@ -976,15 +986,24 @@ Cut-off Gaussian
     :end-before: # { end #}
 
 2D Gaussian
-   We conclude with an example using a two-dimensional Gaussian
-   distribution, i.e., a Gaussian with different widths in :math:`x`- and
-   :math:`y`- directions. This probability depends on displacement, not
-   only on distance:
+   Here we use a two-dimensional Gaussian distribution, i.e., a Gaussian with
+   different widths in :math:`x`- and :math:`y`- directions. This probability
+   depends on displacement, not only on distance:
 
 .. literalinclude:: scripts/connections.py
     :start-after: # { conn42d #}
     :end-before: # { end #}
 
+Rectified Gabor Function
+   We conclude with an example of a rectified Gabor distribution, i.e., a
+   two-dimensional Gaussian distribution modulated with a spatial oscillation
+   perpendicular to :math:`\theta`. This probability depends on the
+   displacement along the coordinates axes, not the distance:
+
+.. literalinclude:: spatially_structured_networks/scripts/connections.py
+    :start-after: #{ conn4gab #}
+    :end-before: #{ end #}
+
 Note that for pool layers with periodic boundary conditions, NEST
 always uses the shortest possible displacement vector from driver to
 pool neuron as ``nest.spatial.distance``.

nestkernel/nestmodule.cpp

@@ -2192,6 +2192,7 @@ NestModule::init( SLIInterpreter* i )
   register_parameter< Gaussian2DParameter >( "gaussian2d" );
   register_parameter< GammaParameter >( "gamma" );
   register_parameter< ExpDistParameter >( "exp_distribution" );
+  register_parameter< GaborParameter >( "gabor" );
 
 #ifdef HAVE_LIBNEUROSIM
   i->createcommand( "CGParse", &cgparse_sfunction );

nestkernel/parameter.cpp

@@ -351,6 +351,50 @@ Gaussian2DParameter::value( RngPtr rng,
 }
 
 
+GaborParameter::GaborParameter( const DictionaryDatum& d )
+  : Parameter( true )
+  , px_( getValue< ParameterDatum >( d, "x" ) )
+  , py_( getValue< ParameterDatum >( d, "y" ) )
+  , cos_( std::cos( getValue< double >( d, "theta" ) * numerics::pi / 180. ) )
+  , sin_( std::sin( getValue< double >( d, "theta" ) * numerics::pi / 180. ) )
+  , gamma_( getValue< double >( d, "gamma" ) )
+  , inv_two_std2_( 1.0 / ( 2 * getValue< double >( d, "std" ) * getValue< double >( d, "std" ) ) )
+  , lambda_( getValue< double >( d, "lam" ) )
+  , psi_( getValue< double >( d, "psi" ) )
+{
+  const auto gamma = getValue< double >( d, "gamma" );
+  const auto std = getValue< double >( d, "std" );
+  if ( std <= 0 )
+  {
+    throw BadProperty( "std > 0 required for gabor function parameter, got std=" + std::to_string( std ) );
+  }
+  if ( gamma <= 0 )
+  {
+    throw BadProperty( "gamma > 0 required for gabor function parameter, got gamma=" + std::to_string( gamma ) );
+  }
+}
+
+double
+GaborParameter::value( RngPtr rng,
+  const std::vector< double >& source_pos,
+  const std::vector< double >& target_pos,
+  const AbstractLayer& layer,
+  Node* node )
+{
+  const auto dx = px_->value( rng, source_pos, target_pos, layer, node );
+  const auto dy = py_->value( rng, source_pos, target_pos, layer, node );
+  const auto dx_prime = dx * cos_ + dy * sin_;
+  const auto dy_prime = -dx * sin_ + dy * cos_;
+  const auto gabor_exp =
+    std::exp( -gamma_ * gamma_ * dx_prime * dx_prime * inv_two_std2_ - dy_prime * dy_prime * inv_two_std2_ );
+  const auto gabor_cos_plus =
+    std::max( std::cos( 2 * numerics::pi * dy_prime / lambda_ + psi_ * numerics::pi / 180. ), 0. );
+  const auto gabor_res = gabor_exp * gabor_cos_plus;
+
+  return gabor_res;
+}
+
+
 GammaParameter::GammaParameter( const DictionaryDatum& d )
   : Parameter( true )
   , p_( getValue< ParameterDatum >( d, "x" ) )

nestkernel/parameter.h

@@ -1366,6 +1366,58 @@ class Gaussian2DParameter : public Parameter
   const double xy_term_const_;
 };
 
+class GaborParameter : public Parameter
+{
+public:
+  using Parameter::value;
+
+  /**
+   * Construct the parameter from a dictionary of arguments.
+   */
+  GaborParameter( const DictionaryDatum& d );
+
+  /**
+   * Copy constructor.
+   */
+  GaborParameter( const GaborParameter& p )
+    : Parameter( p )
+    , px_( p.px_ )
+    , py_( p.py_ )
+    , cos_( p.cos_ )
+    , sin_( p.sin_ )
+    , gamma_( p.gamma_ )
+    , inv_two_std2_( p.inv_two_std2_ )
+    , lambda_( p.lambda_ )
+    , psi_( p.psi_ )
+  {
+  }
+
+  /**
+   * @returns the value of the parameter.
+   */
+  double
+  value( RngPtr, Node* ) override
+  {
+    throw BadParameterValue( "Gabor parameter can only be used when connecting." );
+  }
+
+  double value( RngPtr rng,
+    const std::vector< double >& source_pos,
+    const std::vector< double >& target_pos,
+    const AbstractLayer& layer,
+    Node* node ) override;
+
+protected:
+  std::shared_ptr< Parameter > const px_;
+  std::shared_ptr< Parameter > const py_;
+  const double cos_;
+  const double sin_;
+  const double gamma_;
+  const double inv_two_std2_;
+  const double lambda_;
+  const double psi_;
+};
+
 
 /**
  * Parameter class representing a gamma distribution applied on a parameter.

pynest/nest/spatial_distributions/hl_api_spatial_distributions.py

@@ -21,12 +21,7 @@
 
 from ..lib.hl_api_types import CreateParameter
 
-__all__ = [
-    "exponential",
-    "gaussian",
-    "gaussian2D",
-    "gamma",
-]
+__all__ = ["exponential", "gaussian", "gaussian2D", "gamma", "gabor"]
 
 
 def exponential(x, beta=1.0):
@@ -143,3 +138,49 @@ def gamma(x, kappa=1.0, theta=1.0):
         Object yielding values drawn from the distribution.
     """
     return CreateParameter("gamma", {"x": x, "kappa": kappa, "theta": theta})
+
+
+def gabor(x, y, theta=0.0, gamma=1.0, std=1.0, lam=1.0, psi=0.0):
+    r"""
+    Applies a rectified Gabor function on two Parameters, representing displacements in the x and y direction.
+
+    .. math::
+        p(x) = \big[\cos(360^{\circ} \frac{y^{\prime}}{\lambda} + \psi)\big]^{+}
+               e^{-\frac{\gamma^{2}x^{\prime 2}+y^{\prime 2}}{2\text{std}^{2}}}
+        \\ x^{\prime} = x\cos\theta + y\sin\theta
+        \\ y^{\prime} = -x\sin\theta + y\cos\theta
+
+    Parameters
+    ----------
+    x : Parameter
+        Input Parameter for the x-direction.
+    y : Parameter
+        Input Parameter for the y-direction.
+    theta : float, optional
+        Orientation of profile in degree. Default is 0.0.
+    gamma : float, optional
+        Spatial aspect ratio. Controls ratio of major and minor axis. Default is 1.0.
+    std : float, optional
+        Standard deviation of the distribution. Default is 1.0.
+    lam : float, optional
+        Wavelength of spatial oscillations. Default is 1.0.
+    psi : float, optional
+        Phase of spatial oscillations in degree. Default is 0.0
+
+    Returns
+    -------
+    Parameter:
+        Object yielding values drawn from the distribution.
+    """
+    return CreateParameter(
+        "gabor",
+        {
+            "x": x,
+            "y": y,
+            "theta": theta,
+            "gamma": gamma,
+            "std": std,
+            "lam": lam,
+            "psi": psi,
+        },
+    )

testsuite/pytests/test_spatial/test_spatial_distributions.py

@@ -122,6 +122,7 @@ def __init__(self, seed, dim, L, N, spatial_distribution, distribution_params=No
             "gaussian": self._gauss,
             "gaussian2d": self._gauss2d,
             "gamma": self._gamma,
+            "gabor": self._gabor,
         }
         self._distribution = spatial_distributions[spatial_distribution]
 
@@ -140,6 +141,7 @@ def __init__(self, seed, dim, L, N, spatial_distribution, distribution_params=No
                 "c": 0.0,
             },
             "gamma": {"kappa": 3.0, "theta": self._L / 4.0},
+            "gabor": {"std": self._L, "gamma": 0.7, "lam": self._L, "theta": 360.0 / 8.0, "psi": -360.0 / 6.0},
         }
         self._params = default_params[spatial_distribution]
         if distribution_params is not None:
@@ -184,6 +186,16 @@ def __init__(self, seed, dim, L, N, spatial_distribution, distribution_params=No
             distribution = nest.spatial_distributions.gamma(
                 nest.spatial.distance, kappa=self._params["kappa"], theta=self._params["theta"]
             )
+        elif spatial_distribution == "gabor":
+            distribution = nest.spatial_distributions.gabor(
+                nest.spatial.distance,
+                nest.spatial.distance,
+                std=self._params["std"],
+                gamma=self._params["gamma"],
+                lam=self._params["lam"],
+                theta=self._params["theta"],
+                psi=self._params["psi"],
+            )
 
         self._conndict = {"rule": "pairwise_bernoulli", "p": distribution, "mask": maskdict}
 
@@ -256,6 +268,14 @@ def _gamma(self, D):
         """Gamma spatial distribution"""
         return D**self.kappa_m_1 / self.gamma_kappa_mul_theta_pow_kappa * math.exp(-D / self._params["theta"])
 
+    def _gabor(self, D):
+        """Rectified Gabor distribution"""
+        x_prime = D * np.cos(self._params["theta"] * np.pi / 180) + D * np.sin(self._params["theta"] * np.pi / 180)
+        y_prime = -D * np.sin(self._params["theta"] * np.pi / 180) + D * np.cos(self._params["theta"] * np.pi / 180)
+        gabor_cos = np.cos(2 * np.pi * y_prime / self._params["lam"] + self._params["psi"] * np.pi / 180)
+        gabor_exp = np.exp(-(self._params["gamma"] ** 2 * x_prime**2 + y_prime**2) / (2 * self._params["std"] ** 2))
+        return np.maximum(gabor_cos, 0) * gabor_exp
+
     def _create_distance_data(self):
         self._reset(self._seed)
         self._build()
@@ -615,6 +635,14 @@ def test_gamma(self):
         self.assertGreater(p_ks, P_MIN, "{} failed KS-test".format(distribution))
         self.assertGreater(p_Z, P_MIN, "{} failed Z-test".format(distribution))
 
+    def test_gabor(self):
+        distribution = "gabor"
+        test = SpatialTester(seed=SEED, dim=2, L=1.0, N=10000, spatial_distribution=distribution)
+        _, p_ks = test.ks_test()
+        _, p_Z = test.z_test()
+        self.assertGreater(p_ks, P_MIN, "{} failed KS-test".format(distribution))
+        self.assertGreater(p_Z, P_MIN, "{} failed Z-test".format(distribution))
+
 
 class TestSpatial2DOBC(unittest.TestCase):
     """
@@ -669,6 +697,14 @@ def test_gamma(self):
         self.assertGreater(p_ks, P_MIN, "{} failed KS-test".format(distribution))
         self.assertGreater(p_Z, P_MIN, "{} failed Z-test".format(distribution))
 
+    def test_gabor(self):
+        distribution = "gabor"
+        test = SpatialTester(seed=SEED, dim=2, L=1.0, N=10000, spatial_distribution=distribution, open_bc=True)
+        _, p_ks = test.ks_test()
+        _, p_Z = test.z_test()
+        self.assertGreater(p_ks, P_MIN, "{} failed KS-test".format(distribution))
+        self.assertGreater(p_Z, P_MIN, "{} failed Z-test".format(distribution))
+
 
 class TestSpatial3D(unittest.TestCase):
     """