/*
* arr.sli
*
* 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/>.
*
*/
/*
A library of SLI-routines operating on arrays
Authors:
R Kupper
M-O Gewaltig
*/
/arr ($Revision: 9952 $) provide
/arr namespace
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::Sum - Return sum of the elements in array.
Synopsis: [array] Sum -> number
Description:
Returns the sum of all elements in the array. The array is expected
to contain only numbers, and is flattened before the computation.
Parameters:
[array]: array of numbers
Examples:
[1 2 3 4] Sum -> 10
[1 2 3 4.0] Sum -> 10.0
Author: R Kupper
FirstVersion: 17-sep-2007
Remarks:
The array must contain no other elements than numbers.
The return type (integer or double) depends on the type of elements.
Availability: library "arr"
SeeAlso: arr::Product, arr::Mean, arr::SDev, arr::Var
*/
/Sum[/arraytype]
{
Flatten arrayload 1 sub {add} repeat
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::Product - Return product of the elements in array.
Synopsis: [array] Product -> number
Description:
Returns the product of all elements in the array. The array is expected
to contain only numbers, and is flattened before the computation.
Parameters:
[array]: array of numbers
Examples:
[1 2 3 4] Product -> 24
[1 2 3 4.0] Product -> 24.0
Author: R Kupper
FirstVersion: 17-sep-2007
Remarks:
The array must contain no other elements than numbers.
The return type (integer or double) depends on the type of elements.
Availability: library "arr"
SeeAlso: arr::Sum, arr::Mean, arr::SDev, arr::Var
*/
/Product[/arraytype]
{
Flatten arrayload 1 sub {mul} repeat
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::Mean - Return mean of the elements in array.
Synopsis: [array] Mean -> number
Description:
Returns the sum of all elements in the array. The array is expected
to contain only numbers, and is flattened before the computation.
Parameters:
[array]: array of numbers
Examples:
[1 2 3 4] Mean -> 2.5
Author: R Kupper
FirstVersion: 17-sep-2007
Remarks:
The array must contain no other elements than numbers.
The return type is always double.
Availability: library "arr"
SeeAlso: arr::Product, arr::Sum, arr::SDev, arr::Var
*/
/Mean[/arraytype]
{
Flatten size exch
arrayload 1 sub {add} repeat
cvd
exch div
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::Var - Return variance of the elements in array.
Synopsis: [array] Var -> number
Description:
Returns the variance of all elements in the array. The array is expected
to contain only numbers, and is flattened before the computation.
Parameters:
[array]: array of numbers
Examples:
[0 1 2] Var -> 1.0
[1 2 3 4] Var -> 1.666667
Author: R Kupper
FirstVersion: 17-sep-2007
Remarks:
The array must contain no other elements than numbers.
The return type is always double.
Availability: library "arr"
SeeAlso: arr::Product, arr::Sum, arr::Mean, arr::SDev
*/
/Var[/arraytype]
{
Flatten
size exch
dup Mean sub
{sqr} Map Sum
exch 1 sub div
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::SDev - Return standard deviation of the elements in array.
Synopsis: [array] SDev -> number
Description:
Returns the standard deviation of all elements in the array. The array is expected
to contain only numbers, and is flattened before the computation.
Parameters:
[array]: array of numbers
Examples:
[1 2 3 4] SDev ->
Author: R Kupper
FirstVersion: 17-sep-2007
Remarks:
The array must contain no other elements than numbers.
The return type is always double.
Availability: library "arr"
SeeAlso: arr::Product, arr::Sum, arr::Mean, arr::Var
*/
/SDev[/arraytype]
{
Var sqrt
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::Reform - Reform the dimensions of a hyperrectengular array
Synopsis: [array] [dim1 dim2 ... dimn] Reform -> [reformedArray]
Description:
"Reform" changes the dimensionality i.e. the nesting of an array. The
result is a hyperrectangular array of the given dimensions. The order
of elements is not changed.
"Reform" relates to "Dimensions" in the following way:
Given a hyperrectangular array [a] of "Dimensions" [d], these
statements are true:
a d Reform -> a
a Flatten d Reform -> a
a dv Reform d Reform -> a, for any valid dimension vector [dv]
Diagnostics:
If the new dimensionality is not consistent with the numer of elements
in the array, /RangeCheck is raised.
Availability:
"Namespace"-dictionary "arr".
Author: Ruediger Kupper
FirstVersion: 10.3.2003
Remarks:
The new dimensionality must not change the number of elements in the
whole array, that is, dim1*...*dimn must equal the number of the
flattened source array.
References:
This routine is inspired by IDL's REFORM() routine.
SeeAlso: Dimensions, TensorRank, Flatten
*/
/Reform [/arraytype /arraytype]
{
%stack: [array] [dim1 dim2 ... dimn]
% first do consistency check:
dup Times
%stack: [array] [dim1 dim2 ... dimn] product
2 pick Flatten size
%stack: [array] [dim1 dim2 ... dimn] product [flatarray] size
rolld neq
{
%stack: [array] [dim1 dim2 ... dimn] [flatarray]
M_ERROR (Reform) (New dimensions must not change the number of elements in the array.) message
pop
%stack: [array] [dim1 dim2 ... dimn]
/Reform /RangeCheck raiseerror
} if
%stack: [array] [dim1 dim2 ... dimn] [flatarray]
rolld pop exch
%stack: [flatarray] [dim1 dim2 ... dimn]
% having passed this check, we can discard the first (redundant) dimension:
Rest
%stack: [flatarray] [dim2 ... dimn]
reverse {Partition} Fold
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::getdeep - get an element from a nested container (array, dict, ...).
Synopsis:
(general form:)
nested_container index_array getdeep -> element
(specific forms:)
[nested_array] [i1 i2 ... in] getdeep -> element
{nested_procedure} [i1 i2 ... in] getdeep -> element
{nested_litprocedure} [i1 i2 ... in] getdeep -> element
<<nested_dictionary>> [/n1 /n2 ... /nn] getdeep -> element
Description:
"Getdeep" resolves to a repeated call to "get". That is, it retrieves
a single element from a nested array. The call iterates as deep into
the nested structure as is indicated by the number of elements in the
index array.
In short: The call "a [i1 i2 ... in] getdeep" is identical to
"a i1 get i2 get ... in get".
Parameters:
"Getdeep" can be used on all containers that "get" can be used on.
Since strings cannot be nested, it does not make much sense for
strings, though.
The first argument is a (probably nested) container.
The second argument is a flat index array. For all array type containers,
this must be an array of integers. For nested dictionaries, this must
be an array of literals.
The call iterates as deep into
the nested structure as is indicated by the number of elements in the
index array.
Examples:
[[1 2 3] [4 5 6] [7 8 9]] [0 2] getdeep -> 3
[[1 2 3] [4 5 6] [7 8 9]] [1] getdeep -> [4 5 6]
<< /a << /b 1 /c 2 >> >> [/a /c] getdeep -> 2
Diagnostics:
The number of elements in the index array must not exceed the number
of nested levels in the container. Otherwise, the nested call to get will
raise /ArgumentType.
The elements in the index array must be suited to index the container.
Otherwise, the nested call to get will raise /ArgumentType.
For array containers, the index at position i must be smaller than the
number of elements at level i in the container. Otherwise, the nested
call to get will raise /RangeCheck.
Author: R Kupper
FirstVersion: 19-jun-2006
Remarks:
There currently is no "putdeep" command, but this functionality is
provided by "put". See discussion thread on the nest developer list
from 19-jun-2006.
Availability: library "arr"
SeeAlso: put, get, Dimensions
*/
/getdeep[/anytype /arraytype]
{
{get} forall
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::EdgeTruncate - Truncate 2-d array indices at array edges
Synopsis: [2d-indices] height width EdgeTruncate -> [truncated-2d-indices]
Description:
This function iterates through the given array indices and checks if
they lie inside the bounds [0,height) and [0,width), respectively.
Indices are modified according to the following rules:
1. If both indices lie inside [0,height) and [0,width), respectively,
they are left untouched.
2. If either the row index lies outside [0,height), or the column
index lies outside [0,width), the respective index is replaced by "false".
Note that by NEST convention, for index pairs, the first index denotes
the row, and the second index denotes the column.
Dimensions of the index array are preserved.
Diagnostics:
The index array is expected to be a (nested) array of integer values
only. Code will break otherwise.
Availability:
"Namespace"-dictionary "arr".
Author: Ruediger Kupper
FirstVersion: 17.3.2003
Remarks:
The index array is expected to be a (nested) array of integer values only.
SeeAlso: arr::EdgeWrap, arr::EdgeClip, area2
*/
/EdgeTruncate [/arraytype /integertype /integertype]
{
%stack: [2d-indices] height width
2 arraystore
%stack: [2d-indices] [height width]
% we store array dimensions before operation and restore it afterwards:
exch dup Dimensions rollu
%stack: [Dimensions] [height width] [2d-indices]
Flatten 2 Partition
%stack: [Dimensions] [height width] [[y1 x1] [y2 x2] ... [yn xn]]
{
%entry stack: [height width] [yi xi]
1 pick 2 arraystore
%stack: [height width] [ [yi xi] [height width] ]
{
%entry stack: index maxvalue
1 pick rollu
%stack: index index maxvalue
geq_ii
{
pop
false
}
{
%stack: index
dup 0 lt_ii
{
pop
false
} if
} ifelse
%exit stack: truncated_index
} MapThread
%exit stack : [height width] [yi_truncated xi_truncated]
} Map
%stack: [Dimensions] [height width] [[y1_truncated x1_truncated] [y2_truncated x2_truncated] ... [yn_truncated xn_truncated]]
rolld Reform
%stack: [height width] [truncated-2d-indices-in-right-dimensions]
exch pop
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::EdgeClip - Clip 2-d array indices at array edges
Synopsis: [2d-indices] height width EdgeClip -> [clipped-2d-indices]
Description:
This function iterates through the given array indices and checks if
they lie inside the specified range [0,height) and [0,width),
respectively.
Indices are modified according to the following rules:
1. If both indices lie inside [0,height) and [0,width), respectively,
they are left untouched.
2. If the row index lies below 0, this index is replaced by 0.
3. If the row index lies above "height", this index is replaced by "height".
4. If the column index lies below 0, this index is replaced by 0.
5. If the column index lies above "width", this index is replaced by "width".
Thus, the indices are effectively hard clipped to the array bounds;
that is, to the range [0,height), [0,width) respectively.
Note that by NEST convention, for index pairs, the first index denotes
the row, and the second index denotes the column.
Dimensions of the index array are preserved.
Diagnostics:
The index array is expected to be a (nested) array of integer values
only. Code will break otherwise.
Availability:
"Namespace"-dictionary "arr".
Author: Ruediger Kupper
FirstVersion: 17.3.2003
Remarks:
The index array is expected to be a (nested) array of integer values only.
SeeAlso: arr::EdgeWrap, arr::EdgeTruncate, area2
*/
/EdgeClip [/arraytype /integertype /integertype]
{
%stack: [2d-indices] height width
2 arraystore
%stack: [2d-indices] [height width]
{1 sub_ii} Map
%stack: [2d-indices] [height-1 width-1]
% we store array dimensions before operation and restore it afterwards:
exch dup Dimensions rollu
%stack: [Dimensions] [height-1 width-1] [2d-indices]
Flatten 2 Partition
%stack: [Dimensions] [height-1 width-1] [[y1 x1] [y2 x2] ... [yn xn]]
{
%entry stack: [height-1 width-1] [yi xi]
1 pick 2 arraystore
%stack: [height-1 width-1] [ [yi xi] [height-1 width-1] ]
{
%entry stack: index maxvalue
min_i_i
0 max_i_i
%exit stack: clipped_index
} MapThread
%exit stack : [height-1 width-1] [yi_clipped xi_clipped]
} Map
%stack: [Dimensions] [height-1 width-1] [[y1_clipped x1_clipped] [y2_clipped x2_clipped] ... [yn_clipped xn_clipped]]
rolld Reform
%stack: [height-1 width-1] [clipped-2d-indices-in-right-dimensions]
exch pop
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::EdgeWrap - Wrap 2-d array indices around edges (toriodal)
Synopsis: [2d-indices] height width EdgeWrap -> [wrapped-2d-indices]
Description:
This function iterates through the given array indices and checks if
they lie inside the specified range [0,height) and [0,width),
respectively.
Indices are modified according to the following rules:
1. If both indices lie inside [0,height) and [0,width), respectively,
they are left untouched.
2. If the row index lies outside [0,height) it is cyclicly wrapped
around. That is, a suitable multiple of "height" is added or
substracted, that makes the index fall inside [0,height).
3. If the column index lies outside [0,width) it is cyclicly wrapped
around. That is, a suitable multiple of "width" is added or
substracted, that makes the index fall inside [0,width).
Thus, the indices are effectively wrapped around the array edges;
that is, they are mapped onto a torus of dimensions height,width.
Note that by NEST convention, for index pairs, the first index denotes
the row, and the second index denotes the column.
Dimensions of the index array are preserved.
Diagnostics:
The index array is expected to be a (nested) array of integer values
only. Code will break otherwise.
Availability:
"Namespace"-dictionary "arr".
Author: Ruediger Kupper
FirstVersion: 14.3.2003 (Einstein's birthday)
Remarks:
The index array is expected to be a (nested) array of integer values only.
SeeAlso: arr::IndexWrap, arr::EdgeTruncate, arr::EdgeClip, area2
*/
/EdgeWrap [/arraytype /integertype /integertype]
{
%stack: [2d-indices] height width
2 arraystore
%stack: [2d-indices] [height width]
% we store array dimensions before operation and restore it afterwards:
exch dup Dimensions rollu
%stack: [Dimensions] [height width] [2d-indices]
Flatten 2 Partition
%stack: [Dimensions] [height width] [[y1 x1] [y2 x2] ... [yn xn]]
{
%entry stack: [height width] [yi xi]
1 pick 2 arraystore
%stack: [height width] [ [yi xi] [height width] ]
{
%entry stack: index maxvalue
IndexWrap_i_i
%exit stack: wrapped_index
} MapThread
%exit stack : [height width] [yi_wrapped xi_wrapped]
} Map
%stack: [Dimensions] [height width] [[y1_wrapped x1_wrapped] [y2_wrapped x2_wrapped] ... [yn_wrapped xn_wrapped]]
rolld Reform
%stack: [height width] [wrapped-2d-indices-in-right-dimensions]
exch pop
} bind def
%----------------------------------------------------------------------------
/* BeginDocumentation
Name: arr::IndexWrap - project a cyclic index value onto interval [0,N).
Synopsis: index N CyclicValue -> normindex
Description:
"IndexWrap" projects a cyclic integer index in the range (-oo,oo),
of periodicy N, onto its norm interval [0,N).
This function can be used to "wrap around" array indices in order to
index an array
Alternatives: Function IndexWrap_i_i (undocumented) -> behaviour and
synopsis are the same, except that no warnings or error messages are
thrown.
Parameters:
In : index: integer value in (-oo,oo).
N: Peroidicity of the cyclic index.
"index" is projected on the half-open interval [0,N).
N must be positive (and different from 0).
Out: The cyclic equivalent of the given index, regarding period N.
Diagnostics:
N must be positive (and different from 0). If N <= 0, /RangeCheck is raised.
Note that the variant IndexWrap_i_i does not do this check for
efficiency, and will break or yield invalid results in this case.
Examples:
-6 3 IndexWrap -> 0
-5 3 IndexWrap -> 1
-4 3 IndexWrap -> 2
-3 3 IndexWrap -> 0
-2 3 IndexWrap -> 1
-1 3 IndexWrap -> 2
0 3 IndexWrap -> 0
1 3 IndexWrap -> 1
2 3 IndexWrap -> 2
3 3 IndexWrap -> 0
4 3 IndexWrap -> 1
5 3 IndexWrap -> 2
6 3 IndexWrap -> 0
Availability:
"Namespace"-dictionary "arr".
Author: Ruediger Kupper
FirstVersion: 14.2.2003 (Einstein's birthday)
Remarks:
This function behaves like Mathematica's Mod function (which is
different from the mathematical definition of MOD).
SeeAlso: arr::EdgeWrap, mod, CyclicValue
*/
/IndexWrap_i_i
{
%stack: index N
% index might be negative, so we add a suitable multiple of N:
% posindex = index + (abs(index)/N + 1)*N with integer div.
exch 1 pick 1 pick abs_i 1 pick
%stack: N index N abs(index) N
div_ii 1 add_ii mul_ii add_ii
%stack: N posindex
% now that we know everything is positive, we can use mod:
exch_ mod
} bind def
/IndexWrap [/integertype /integertype]
{
dup_ 0 gt_ii
{
IndexWrap_i_i
}
{
M_ERROR (IndexWrap) (Period N must be positive (and different from 0).) message
/IndexWrap /RangeCheck raiseerror
}ifelse
} bind def
/* BeginDocumentation
Name: arr::GaborPatch - Return a two-dimensional array with Gabor function.
Synopsis:
nrows ncols GaborPatch -> [[] .. []]
Parameters:
nrows: number of rows of the result matrix.
ncols: number of columns of the result matrix.
Description:
This function returns a matrix of nrows by ncols, with the amplitudes
of a Gabor function, computed over the range of [x_min,x_max] by [y_min,y_max].
These and other parameters can be set by changing the function's options.
An orientation of 0 RAD results in a vertically oriented Gabor.
An orientation of Pi/2 results in an horizontally oriented Gabor.
Angles are measured counter clockwise from the positive x-axis.
The implementation follows the description given by
N. Petkov and R. Kruizinga Biol.Cybern. 76, 83-97 (1997).
Note that GaborPatch automatically performs a coordinate transformation
from the mathematical x-y plane to the matrix row-column system.
Options:
Options which determine the argument range of the Gabor patch.
x_min doubletype - smallest x coordinate value [-2Pi]
x_max doubletype - largest x coordinate value. [ 2Pi]
y_min doubletype - smallest y coordinate value. [-2Pi]
y_max doubletype - largest y coordinate value. [ 2Pi]
Options which determine the Gabor function:
lambda doubletype - Wavelength of the Gabor in RAD. [ 2Pi]
phase doubletype - Phase of the Gabor in RAD. [ 0.0]
orientation doubletype - Rotation of the Gabor in RAD [ 0.0]
sigma doubletype - Width of the Gaussian envelope. [ Pi]
gamma doubletype - Spatial aspect ratio. [ 1.0]
The options correspond to the status dictionary entries of the
gabor_device. Thus, it is possible to use the option dictionary
of the GaborPatch to set the properties of a gabor device.
Author: Marc-Oliver Gewaltig
References: Petkov N and Kruizinga P: Biol. Cybern. 76, 83-96 (1997)
SeeAlso: gabor_
*/
/GaborPatch
[/integertype /integertype]
{
/arr::GaborPatch GetOptions
begin
x_min x_max
y_min y_max
lambda orientation phase sigma gamma
end
gabor_
} def
%% initialization values correspond to those of the NEST
%% gabor_generator device (see file synod2/nest/gabor.cpp).
/arr::GaborPatch
<<
/x_min -2.0 Pi mul
/x_max 2.0 Pi mul
/y_min -2.0 Pi mul
/y_max 2.0 Pi mul
/lambda 2.0 Pi mul
/orientation 0.0
/phase 0.0
/sigma Pi
/gamma 1.0
>> Options
/* BeginDocumentation
Name: arr::GaussPatch - Return a two-dimensional array with Gauss function.
Synopsis:
nrows ncols GaussPatch -> [[] .. []]
Parameters:
nrows: number of rows of the result matrix.
ncols: number of columns of the result matrix.
Description:
This function returns a matrix of nrows by ncols, with the amplitudes
of a two dimensional Gauss function, computed over the range of
[x_min,x_max] by [y_min,y_max].
These and other parameters can be set by changing the function's options.
If the aspect ratio gamma < 1, an orientation of 0 RAD results in
a vertically oriented eliptic Gauss .
Angles are measured counter clockwise from the positive x-axis.
Note that GaussPatch automatically performs a coordinate transformation
from the mathematical x-y plane to the matrix row-column system.
Options:
Options which determine the argument range of the Gauss patch.
x_min doubletype - smallest x coordinate value [-3.0]
x_max doubletype - largest x coordinate value. [ 3.0]
y_min doubletype - smallest y coordinate value. [-3.0]
y_max doubletype - largest y coordinate value. [ 3.0]
Options which determine the Gaussian function:
orientation doubletype - Rotation of the Gaussian in RAD [ 0.0]
sigma doubletype - Width of the Gaussian. [ 1.0]
gamma doubletype - Spatial aspect ratio. [ 1.0]
Normalization options:
Normalize booltype - Normalize the Gaussian to [false]
integral 1.0 (before sampling
to output array)
NormalizeSum booltype - After sampling the Gaussian, [false]
normalize the output array
to have total sum of 1.0
Author: Marc-Oliver Gewaltig
SeeAlso: gauss2d_
*/
/GaussPatch
[/integertype /integertype]
{
/arr::GaussPatch GetOptions
begin
x_min x_max
y_min y_max
orientation sigma gamma
gauss2d_
Normalize
{%normalize the integral by a factor (sigma_x*sqrt(2Pi))*(sigma_y*sqrt(2Pi))
% = sigma_x*sigma_y*2Pi
% where sigma_x=sigma, sigma_y=sigma/gamma
2 Pi mul sigma mul sigma gamma div mul div
} if
NormalizeSum
{%normalize to total sum of array elements
dup Flatten Plus div
} if
end
} def
/arr::GaussPatch
<<
/x_min -3.0
/x_max 3.0
/y_min -3.0
/y_max 3.0
/orientation 0.0
/sigma 1.0
/gamma 1.0
/Normalize false
/NormalizeSum false
>> Options
end % namespace arr
%% NOTE: There must be a carriage return after the last line, i.e., here: