diff --git a/doc/htmldoc/networks/scripts/connections.py b/doc/htmldoc/networks/scripts/connections.py index 5663636a7b..b420a30cb1 100644 --- a/doc/htmldoc/networks/scripts/connections.py +++ b/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") # ----------------------------------------------- diff --git a/doc/htmldoc/networks/spatially_structured_networks.rst b/doc/htmldoc/networks/spatially_structured_networks.rst index e3416dfc67..01cb0bdaa4 100644 --- a/doc/htmldoc/networks/spatially_structured_networks.rst +++ b/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``. diff --git a/nestkernel/nestmodule.cpp b/nestkernel/nestmodule.cpp index 6665868e67..927790abc1 100644 --- a/nestkernel/nestmodule.cpp +++ b/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 ); diff --git a/nestkernel/parameter.cpp b/nestkernel/parameter.cpp index 4c2d1137ac..837c43c8de 100644 --- a/nestkernel/parameter.cpp +++ b/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" ) ) diff --git a/nestkernel/parameter.h b/nestkernel/parameter.h index 1cbee8ce83..dfae3cfcb9 100644 --- a/nestkernel/parameter.h +++ b/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. diff --git a/pynest/nest/spatial_distributions/hl_api_spatial_distributions.py b/pynest/nest/spatial_distributions/hl_api_spatial_distributions.py index a0512993d6..bd46da15c6 100644 --- a/pynest/nest/spatial_distributions/hl_api_spatial_distributions.py +++ b/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, + }, + ) diff --git a/testsuite/pytests/test_spatial/test_spatial_distributions.py b/testsuite/pytests/test_spatial/test_spatial_distributions.py index d4bcb775ed..89ff1c2955 100644 --- a/testsuite/pytests/test_spatial/test_spatial_distributions.py +++ b/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): """