Context Navigation

← Previous Change
Next Change →

connection.py

Timestamp:

08/12/10 22:19:31 (22 months ago)

Author:

thesamovar

Message:

Reorganised connection.py into new subpackage connections

For the moment, the old connection.py can still be used by changing the _usenew variable in your copy of connection.py

Files:

: 1 modified

trunk/brian/connection.py (modified) (2 diffs)

Legend:

: Unmodified
: Added
: Removed

trunk/brian/connection.py

r2081	r2117
1		# ----------------------------------------------------------------------------------
2		# Copyright ENS, INRIA, CNRS
3		# Contributors: Romain Brette (brette@di.ens.fr) and Dan Goodman (goodman@di.ens.fr)
4		#
5		# Brian is a computer program whose purpose is to simulate models
6		# of biological neural networks.
7		#
8		# This software is governed by the CeCILL license under French law and
9		# abiding by the rules of distribution of free software. You can use,
10		# modify and/ or redistribute the software under the terms of the CeCILL
11		# license as circulated by CEA, CNRS and INRIA at the following URL
12		# "http://www.cecill.info".
13		#
14		# As a counterpart to the access to the source code and rights to copy,
15		# modify and redistribute granted by the license, users are provided only
16		# with a limited warranty and the software's author, the holder of the
17		# economic rights, and the successive licensors have only limited
18		# liability.
19		#
20		# In this respect, the user's attention is drawn to the risks associated
21		# with loading, using, modifying and/or developing or reproducing the
22		# software by the user in light of its specific status of free software,
23		# that may mean that it is complicated to manipulate, and that also
24		# therefore means that it is reserved for developers and experienced
25		# professionals having in-depth computer knowledge. Users are therefore
26		# encouraged to load and test the software's suitability as regards their
27		# requirements in conditions enabling the security of their systems and/or
28		# data to be ensured and, more generally, to use and operate it in the
29		# same conditions as regards security.
30		#
31		# The fact that you are presently reading this means that you have had
32		# knowledge of the CeCILL license and that you accept its terms.
33		# ----------------------------------------------------------------------------------
34		#
35		__all__ = [
36		'Connection', 'IdentityConnection', 'MultiConnection',
37		'DelayConnection',
38		'random_matrix_fixed_column',
39		'random_matrix',
40		'ConstructionMatrix', 'SparseConstructionMatrix', 'DenseConstructionMatrix', 'DynamicConstructionMatrix',
41		'ConnectionMatrix', 'SparseConnectionMatrix', 'DenseConnectionMatrix', 'DynamicConnectionMatrix',
42		'ConnectionVector', 'SparseConnectionVector', 'DenseConnectionVector',
43		]
	1	_usenew = True
44	2
45		import copy
46		from itertools import izip
47		import itertools
48		from random import sample
49		import bisect
50		from units import *
51		import types
52		import magic
53		from log import *
54		from numpy import *
55		from scipy import sparse, stats, rand, weave, linalg
56		import scipy
57		import scipy.sparse
58		import numpy
59		from numpy.random import binomial, exponential
60		import random as pyrandom
61		from scipy import random as scirandom
62		from utils.approximatecomparisons import is_within_absolute_tolerance
63		from globalprefs import *
64		from base import *
65		from stdunits import ms
66		from operator import isSequenceType
67
68		use_sparse_matrix = 'scipy_patch' # values are own, own_scipy, scipy, scipy_patch
69		if use_sparse_matrix == 'own':
70		from utils import sparse_matrix as sparse
71		elif use_sparse_matrix == 'own_scipy':
72		from utils import sparse # Brian's version of scipy sparse matrix library
73		elif use_sparse_matrix == 'scipy_patch':
74		from utils import sparse_patch as sparse
75		# We should do this:
76		# from network import network_operation
77		# but instead we do it at the bottom of the module (see comments there for explanation)
78
79		effective_zero = 1e-40
80
81		def random_row_func(N, p, weight=1., initseed=None):
82		'''
83		Returns a random connectivity ``row_func`` for use with :class:`UserComputedConnectionMatrix`
84
85		Gives equivalent output to the :meth:`Connection.connect_random` method.
86
87		Arguments:
88
89		``N``
90		The number of target neurons.
91		``p``
92		The probability of a synapse.
93		``weight``
94		The connection weight (must be a single value).
95		``initseed``
96		The initial seed value (for reproducible results).
97		'''
98		if initseed is None:
99		initseed = pyrandom.randint(100000, 1000000) # replace this
100		cur_row = numpy.zeros(N)
101		myrange = numpy.arange(N, dtype=int)
102
103		def row_func(i):
104		pyrandom.seed(initseed + int(i))
105		scirandom.seed(initseed + int(i))
106		k = scirandom.binomial(N, p, 1)[0]
107		cur_row[:] = 0.0
108		cur_row[pyrandom.sample(myrange, k)] = weight
109		return cur_row
110
111		return row_func
112
113		colon_slice = slice(None, None, None)
114
115		def todense(x):
116		if hasattr(x, 'todense'):
117		return x.todense()
118		return array(x, copy=False)
119
120
121		class ConnectionVector(object):
122		'''
123		Base class for connection vectors, just used for defining the interface
124
125		ConnectionVector objects are returned by ConnectionMatrix objects when
126		they retrieve rows or columns. At the moment, there are two choices,
127		sparse or dense.
128
129		This class has no real function at the moment.
130		'''
131		def todense(self):
132		return NotImplemented
133
134		def tosparse(self):
135		return NotImplemented
136
137
138		class DenseConnectionVector(ConnectionVector, numpy.ndarray):
139		'''
140		Just a numpy array.
141		'''
142		def __new__(subtype, arr):
143		return numpy.array(arr, copy=False).view(subtype)
144
145		def todense(self):
146		return self
147
148		def tosparse(self):
149		return SparseConnectionVector(len(self), self.nonzero(), self)
150
151
152		class SparseConnectionVector(ConnectionVector, numpy.ndarray):
153		'''
154		Sparse vector class
155
156		A sparse vector is typically a row or column of a sparse matrix. This
157		class can be treated in many cases as if it were just a vector without
158		worrying about the fact that it is sparse. For example, if you write
159		``2*v`` it will evaluate to a new sparse vector. There is one aspect
160		of the semantics which is potentially confusing. In a binary operation
161		with a dense vector such as ``sv+dv`` where ``sv`` is sparse and ``dv``
162		is dense, the result will be a sparse vector with zeros where ``sv``
163		has zeros, the potentially nonzero elements of ``dv`` where ``sv`` has
164		no entry will be simply ignored. It is for this reason that it is a
165		``SparseConnectionVector`` and not a general ``SparseVector``, because
166		these semantics make sense for rows and columns of connection matrices
167		but not in general.
168
169		Implementation details:
170
171		The underlying numpy array contains the values, the attribute ``n`` is
172		the length of the sparse vector, and ``ind`` is an array of the indices
173		of the nonzero elements.
174		'''
175		def __new__(subtype, n, ind, data):
176		x = numpy.array(data, copy=False).view(subtype)
177		x.n = n
178		x.ind = ind
179		return x
180
181		def __array_finalize__(self, orig):
182		# the array is passed through this function after standard numpy operations,
183		# this ensures that the indices are kept from the original array. This makes,
184		# for example, sin(x) do the right thing for x a sparse vector.
185		try:
186		self.ind = orig.ind
187		self.n = orig.n
188		except AttributeError:
189		pass
190		return self
191
192		def todense(self):
193		x = zeros(self.n)
194		x[self.ind] = self
195		return x
196
197		def tosparse(self):
198		return self
199		# This is a list of the binary operations that numpy arrays support.
200		modifymeths = ['__add__', '__and__',
201		'__div__', '__divmod__', '__eq__',
202		'__floordiv__', '__ge__', '__gt__', '__iadd__', '__iand__', '__idiv__',
203		'__ifloordiv__', '__ilshift__', '__imod__', '__imul__',
204		'__ior__', '__ipow__', '__irshift__', '__isub__', '__itruediv__',
205		'__ixor__', '__le__', '__lshift__', '__lt__', '__mod__', '__mul__',
206		'__ne__', '__or__', '__pow__', '__radd__', '__rand__', '__rdiv__',
207		'__rdivmod__', '__rfloordiv__', '__rlshift__',
208		'__rmod__', '__rmul__', '__ror__', '__rpow__', '__rrshift__', '__rshift__',
209		'__rsub__', '__rtruediv__', '__rxor__', '__sub__', '__truediv__', '__xor__']
210		# This template function (where __add__ is replaced by any of the methods above) implements
211		# the semantics described in this class' docstring when operating with a dense vector.
212		template = '''
213		def __add__(self, other):
214		if isinstance(other, SparseConnectionVector):
215		# Note that removing this check is potentially dangerous, but only in weird circumstances would it cause
216		# any problems, and leaving it in causes problems for STDP with DelayConnection (because the indices are
217		# not the same, but should be presumed to be equal).
218		#if other.ind is not self.ind:
219		# raise TypeError('__add__(SparseConnectionVector, SparseConnectionVector) only defined if indices are the same')
220		return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), asarray(other)))
221		if isinstance(other, numpy.ndarray):
222		return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), other[self.ind]))
223		return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), other))
224		'''.strip()
225		# this substitutes any of the method names in the modifymeths list for __add__ in the template
226		# above and then executes them, i.e. adding them as methods to the class. When the behaviour is
227		# stable, this can be replaced by the explicit definitions but it may as well be left as it is for
228		# the moment (it's slower at import time, but not at run time, and errors are more difficult to
229		# catch when it is done like this).
230		for m in modifymeths:
231		s = template.replace('__add__', m)
232		exec s
233		del modifymeths, template
234
235
236		class ConstructionMatrix(object):
237		'''
238		Base class for construction matrices
239
240		A construction matrix is used to initialise and build connection matrices.
241		A ``ConstructionMatrix`` class has to implement a method
242		``connection_matrix(args, *kwds)`` which returns a :class:`ConnectionMatrix`
243		object of the appropriate type.
244		'''
245		def connection_matrix(self, args, *kwds):
246		return NotImplemented
247
248
249		class DenseConstructionMatrix(ConstructionMatrix, numpy.ndarray):
250		'''
251		Dense construction matrix. Essentially just numpy.ndarray.
252
253		The ``connection_matrix`` method returns a :class:`DenseConnectionMatrix`
254		object.
255
256		The ``__setitem__`` method is overloaded so that you can set values with
257		a sparse matrix.
258		'''
259		def __init__(self, val, **kwds):
260		self[:] = 0
261		self.init_kwds = kwds
262
263		def connection_matrix(self, **additional_kwds):
264		self.init_kwds.update(additional_kwds)
265		kwds = self.init_kwds
266		return DenseConnectionMatrix(self, **kwds)
267
268		def __setitem__(self, index, W):
269		# Make it work for sparse matrices
270		#if isinstance(W,sparse.spmatrix):
271		if isinstance(W, sparse.spmatrix):
272		ndarray.__setitem__(self, index, W.todense())
273		else:
274		ndarray.__setitem__(self, index, W)
275
276		def todense(self):
277		return asarray(self)
278
279		# set this to True using the sparse library packaged with Brian, from scipy 0.7.1
280		if use_sparse_matrix == 'own' or use_sparse_matrix == 'own_scipy':
281		oldscipy = True
	3	if _usenew:
	4	from connections import *
282	5	else:
283		olscipy = scipy.__version__.startswith('0.6.') or scipy.__version__.startswith('0.7.1')
284
285		if use_sparse_matrix == 'own' or use_sparse_matrix == 'scipy_patch':
286		SparseMatrix = sparse.lil_matrix
287		else:
288		class SparseMatrix(sparse.lil_matrix):
289		'''
290		Used as the base for sparse construction matrix classes, essentially just scipy's lil_matrix.
291
292		The scipy lil_matrix class allows you to specify slices in ``__setitem__`` but the
293		performance is cripplingly slow. This class has a faster implementation.
294		'''
295		# Unfortunately we still need to implement this because although scipy 0.7.0
296		# now supports X[a:b,c:d] for sparse X it is unbelievably slow (shabby code
297		# on their part).
	6	__all__ = [
	7	'Connection', 'IdentityConnection', 'MultiConnection',
	8	'DelayConnection',
	9	'random_matrix_fixed_column',
	10	'random_matrix',
	11	'ConstructionMatrix', 'SparseConstructionMatrix', 'DenseConstructionMatrix', 'DynamicConstructionMatrix',
	12	'ConnectionMatrix', 'SparseConnectionMatrix', 'DenseConnectionMatrix', 'DynamicConnectionMatrix',
	13	'ConnectionVector', 'SparseConnectionVector', 'DenseConnectionVector',
	14	]
	15
	16	import copy
	17	from itertools import izip
	18	import itertools
	19	from random import sample
	20	import bisect
	21	from units import *
	22	import types
	23	import magic
	24	from log import *
	25	from numpy import *
	26	from scipy import sparse, stats, rand, weave, linalg
	27	import scipy
	28	import scipy.sparse
	29	import numpy
	30	from numpy.random import binomial, exponential
	31	import random as pyrandom
	32	from scipy import random as scirandom
	33	from utils.approximatecomparisons import is_within_absolute_tolerance
	34	from globalprefs import *
	35	from base import *
	36	from stdunits import ms
	37	from operator import isSequenceType
	38
	39	use_sparse_matrix = 'scipy_patch' # values are own, own_scipy, scipy, scipy_patch
	40	if use_sparse_matrix == 'own':
	41	from utils import sparse_matrix as sparse
	42	elif use_sparse_matrix == 'own_scipy':
	43	from utils import sparse # Brian's version of scipy sparse matrix library
	44	elif use_sparse_matrix == 'scipy_patch':
	45	from utils import sparse_patch as sparse
	46	# We should do this:
	47	# from network import network_operation
	48	# but instead we do it at the bottom of the module (see comments there for explanation)
	49
	50	effective_zero = 1e-40
	51
	52	def random_row_func(N, p, weight=1., initseed=None):
	53	'''
	54	Returns a random connectivity ``row_func`` for use with :class:`UserComputedConnectionMatrix`
	55
	56	Gives equivalent output to the :meth:`Connection.connect_random` method.
	57
	58	Arguments:
	59
	60	``N``
	61	The number of target neurons.
	62	``p``
	63	The probability of a synapse.
	64	``weight``
	65	The connection weight (must be a single value).
	66	``initseed``
	67	The initial seed value (for reproducible results).
	68	'''
	69	if initseed is None:
	70	initseed = pyrandom.randint(100000, 1000000) # replace this
	71	cur_row = numpy.zeros(N)
	72	myrange = numpy.arange(N, dtype=int)
	73
	74	def row_func(i):
	75	pyrandom.seed(initseed + int(i))
	76	scirandom.seed(initseed + int(i))
	77	k = scirandom.binomial(N, p, 1)[0]
	78	cur_row[:] = 0.0
	79	cur_row[pyrandom.sample(myrange, k)] = weight
	80	return cur_row
	81
	82	return row_func
	83
	84	colon_slice = slice(None, None, None)
	85
	86	def todense(x):
	87	if hasattr(x, 'todense'):
	88	return x.todense()
	89	return array(x, copy=False)
	90
	91
	92	class ConnectionVector(object):
	93	'''
	94	Base class for connection vectors, just used for defining the interface
	95
	96	ConnectionVector objects are returned by ConnectionMatrix objects when
	97	they retrieve rows or columns. At the moment, there are two choices,
	98	sparse or dense.
	99
	100	This class has no real function at the moment.
	101	'''
	102	def todense(self):
	103	return NotImplemented
	104
	105	def tosparse(self):
	106	return NotImplemented
	107
	108
	109	class DenseConnectionVector(ConnectionVector, numpy.ndarray):
	110	'''
	111	Just a numpy array.
	112	'''
	113	def __new__(subtype, arr):
	114	return numpy.array(arr, copy=False).view(subtype)
	115
	116	def todense(self):
	117	return self
	118
	119	def tosparse(self):
	120	return SparseConnectionVector(len(self), self.nonzero(), self)
	121
	122
	123	class SparseConnectionVector(ConnectionVector, numpy.ndarray):
	124	'''
	125	Sparse vector class
	126
	127	A sparse vector is typically a row or column of a sparse matrix. This
	128	class can be treated in many cases as if it were just a vector without
	129	worrying about the fact that it is sparse. For example, if you write
	130	``2*v`` it will evaluate to a new sparse vector. There is one aspect
	131	of the semantics which is potentially confusing. In a binary operation
	132	with a dense vector such as ``sv+dv`` where ``sv`` is sparse and ``dv``
	133	is dense, the result will be a sparse vector with zeros where ``sv``
	134	has zeros, the potentially nonzero elements of ``dv`` where ``sv`` has
	135	no entry will be simply ignored. It is for this reason that it is a
	136	``SparseConnectionVector`` and not a general ``SparseVector``, because
	137	these semantics make sense for rows and columns of connection matrices
	138	but not in general.
	139
	140	Implementation details:
	141
	142	The underlying numpy array contains the values, the attribute ``n`` is
	143	the length of the sparse vector, and ``ind`` is an array of the indices
	144	of the nonzero elements.
	145	'''
	146	def __new__(subtype, n, ind, data):
	147	x = numpy.array(data, copy=False).view(subtype)
	148	x.n = n
	149	x.ind = ind
	150	return x
	151
	152	def __array_finalize__(self, orig):
	153	# the array is passed through this function after standard numpy operations,
	154	# this ensures that the indices are kept from the original array. This makes,
	155	# for example, sin(x) do the right thing for x a sparse vector.
	156	try:
	157	self.ind = orig.ind
	158	self.n = orig.n
	159	except AttributeError:
	160	pass
	161	return self
	162
	163	def todense(self):
	164	x = zeros(self.n)
	165	x[self.ind] = self
	166	return x
	167
	168	def tosparse(self):
	169	return self
	170	# This is a list of the binary operations that numpy arrays support.
	171	modifymeths = ['__add__', '__and__',
	172	'__div__', '__divmod__', '__eq__',
	173	'__floordiv__', '__ge__', '__gt__', '__iadd__', '__iand__', '__idiv__',
	174	'__ifloordiv__', '__ilshift__', '__imod__', '__imul__',
	175	'__ior__', '__ipow__', '__irshift__', '__isub__', '__itruediv__',
	176	'__ixor__', '__le__', '__lshift__', '__lt__', '__mod__', '__mul__',
	177	'__ne__', '__or__', '__pow__', '__radd__', '__rand__', '__rdiv__',
	178	'__rdivmod__', '__rfloordiv__', '__rlshift__',
	179	'__rmod__', '__rmul__', '__ror__', '__rpow__', '__rrshift__', '__rshift__',
	180	'__rsub__', '__rtruediv__', '__rxor__', '__sub__', '__truediv__', '__xor__']
	181	# This template function (where __add__ is replaced by any of the methods above) implements
	182	# the semantics described in this class' docstring when operating with a dense vector.
	183	template = '''
	184	def __add__(self, other):
	185	if isinstance(other, SparseConnectionVector):
	186	# Note that removing this check is potentially dangerous, but only in weird circumstances would it cause
	187	# any problems, and leaving it in causes problems for STDP with DelayConnection (because the indices are
	188	# not the same, but should be presumed to be equal).
	189	#if other.ind is not self.ind:
	190	# raise TypeError('__add__(SparseConnectionVector, SparseConnectionVector) only defined if indices are the same')
	191	return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), asarray(other)))
	192	if isinstance(other, numpy.ndarray):
	193	return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), other[self.ind]))
	194	return SparseConnectionVector(self.n, self.ind, numpy.ndarray.__add__(asarray(self), other))
	195	'''.strip()
	196	# this substitutes any of the method names in the modifymeths list for __add__ in the template
	197	# above and then executes them, i.e. adding them as methods to the class. When the behaviour is
	198	# stable, this can be replaced by the explicit definitions but it may as well be left as it is for
	199	# the moment (it's slower at import time, but not at run time, and errors are more difficult to
	200	# catch when it is done like this).
	201	for m in modifymeths:
	202	s = template.replace('__add__', m)
	203	exec s
	204	del modifymeths, template
	205
	206
	207	class ConstructionMatrix(object):
	208	'''
	209	Base class for construction matrices
	210
	211	A construction matrix is used to initialise and build connection matrices.
	212	A ``ConstructionMatrix`` class has to implement a method
	213	``connection_matrix(args, *kwds)`` which returns a :class:`ConnectionMatrix`
	214	object of the appropriate type.
	215	'''
	216	def connection_matrix(self, args, *kwds):
	217	return NotImplemented
	218
	219
	220	class DenseConstructionMatrix(ConstructionMatrix, numpy.ndarray):
	221	'''
	222	Dense construction matrix. Essentially just numpy.ndarray.
	223
	224	The ``connection_matrix`` method returns a :class:`DenseConnectionMatrix`
	225	object.
	226
	227	The ``__setitem__`` method is overloaded so that you can set values with
	228	a sparse matrix.
	229	'''
	230	def __init__(self, val, **kwds):
	231	self[:] = 0
	232	self.init_kwds = kwds
	233
	234	def connection_matrix(self, **additional_kwds):
	235	self.init_kwds.update(additional_kwds)
	236	kwds = self.init_kwds
	237	return DenseConnectionMatrix(self, **kwds)
	238
298	239	def __setitem__(self, index, W):
299		"""
300		Speed-up if x is a sparse matrix.
301		TODO: checks (first remove the data).
302		"""
303		try:
304		i, j = index
305		except (ValueError, TypeError):
306		raise IndexError, "invalid index"
307
308		if isinstance(i, slice) and isinstance(j, slice) and\
309		(i.step is None) and (j.step is None) and\
310		(isinstance(W, sparse.spmatrix) or isinstance(W, numpy.ndarray)):
311		rows = self.rows[i]
312		datas = self.data[i]
313		j0 = j.start
314		if isinstance(W, sparse.lil_matrix):
315		for row, data, rowW, dataW in izip(rows, datas, W.rows, W.data):
316		jj = bisect.bisect(row, j0) # Find the insertion point
317		row[jj:jj] = [j0 + k for k in rowW]
318		data[jj:jj] = dataW
319		elif isinstance(W, ndarray):
320		nq = W.shape[1]
321		for row, data, rowW in izip(rows, datas, W):
322		jj = bisect.bisect(row, j0) # Find the insertion point
323		row[jj:jj] = range(j0, j0 + nq)
324		data[jj:jj] = rowW
325		elif oldscipy and isinstance(i, int) and isinstance(j, (list, tuple, numpy.ndarray)):
326		row = dict(izip(self.rows[i], self.data[i]))
	240	# Make it work for sparse matrices
	241	#if isinstance(W,sparse.spmatrix):
	242	if isinstance(W, sparse.spmatrix):
	243	ndarray.__setitem__(self, index, W.todense())
	244	else:
	245	ndarray.__setitem__(self, index, W)
	246
	247	def todense(self):
	248	return asarray(self)
	249
	250	# set this to True using the sparse library packaged with Brian, from scipy 0.7.1
	251	if use_sparse_matrix == 'own' or use_sparse_matrix == 'own_scipy':
	252	oldscipy = True
	253	else:
	254	olscipy = scipy.__version__.startswith('0.6.') or scipy.__version__.startswith('0.7.1')
	255
	256	if use_sparse_matrix == 'own' or use_sparse_matrix == 'scipy_patch':
	257	SparseMatrix = sparse.lil_matrix
	258	else:
	259	class SparseMatrix(sparse.lil_matrix):
	260	'''
	261	Used as the base for sparse construction matrix classes, essentially just scipy's lil_matrix.
	262
	263	The scipy lil_matrix class allows you to specify slices in ``__setitem__`` but the
	264	performance is cripplingly slow. This class has a faster implementation.
	265	'''
	266	# Unfortunately we still need to implement this because although scipy 0.7.0
	267	# now supports X[a:b,c:d] for sparse X it is unbelievably slow (shabby code
	268	# on their part).
	269	def __setitem__(self, index, W):
	270	"""
	271	Speed-up if x is a sparse matrix.
	272	TODO: checks (first remove the data).
	273	"""
327	274	try:
328		row.update(dict(izip(j, W)))
329		except TypeError:
330		row.update(dict(izip(j, itertools.repeat(W))))
331		items = row.items()
332		items.sort()
333		row, data = izip(*items)
334		self.rows[i] = list(row)
335		self.data[i] = list(data)
336		elif isinstance(i, slice) and isinstance(j, int) and isSequenceType(W):
337		# This corrects a bug in scipy sparse matrix as of version 0.7.0, but
338		# it is not efficient!
339		for w, k in izip(W, xrange(*i.indices(self.shape[0]))):
340		sparse.lil_matrix.__setitem__(self, (k, j), w)
341		else:
342		sparse.lil_matrix.__setitem__(self, index, W)
343
344
345		class SparseConstructionMatrix(ConstructionMatrix, SparseMatrix):
346		'''
347		SparseConstructionMatrix is converted to SparseConnectionMatrix.
348		'''
349		def __init__(self, arg, **kwds):
350		SparseMatrix.__init__(self, arg)
351		self.init_kwds = kwds
352
353		def connection_matrix(self, **additional_kwds):
354		self.init_kwds.update(additional_kwds)
355		return SparseConnectionMatrix(self, **self.init_kwds)
356
357
358		class DynamicConstructionMatrix(ConstructionMatrix, SparseMatrix):
359		'''
360		DynamicConstructionMatrix is converted to DynamicConnectionMatrix.
361		'''
362		def __init__(self, arg, **kwds):
363		SparseMatrix.__init__(self, arg)
364		self.init_kwds = kwds
365
366		def connection_matrix(self, **additional_kwds):
367		self.init_kwds.update(additional_kwds)
368		return DynamicConnectionMatrix(self, **self.init_kwds)
369
370		# this is used to look up str->class conversions for structure=... keyword
371		construction_matrix_register = {
372		'dense':DenseConstructionMatrix,
373		'sparse':SparseConstructionMatrix,
374		'dynamic':DynamicConstructionMatrix,
375		}
376
377
378		class ConnectionMatrix(object):
379		'''
380		Base class for connection matrix objects
381
382		Connection matrix objects support a subset of the following methods:
383
384		``get_row(i)``, ``get_col(i)``
385		Returns row/col ``i`` as a :class:`DenseConnectionVector` or
386		:class:`SparseConnectionVector` as appropriate for the class.
387		``set_row(i, val)``, ``set_col(i, val)``
388		Sets row/col with an array, :class:`DenseConnectionVector` or
389		:class:`SparseConnectionVector` (if supported).
390		``get_element(i, j)``, ``set_element(i, j, val)``
391		Gets or sets a single value.
392		``get_rows(rows)``
393		Returns a list of rows, should be implemented without Python
394		function calls for efficiency if possible.
395		``get_cols(cols)``
396		Returns a list of cols, should be implemented without Python
397		function calls for efficiency if possible.
398		``insert(i,j,x)``, ``remove(i,j)``
399		For sparse connection matrices which support it, insert a new
400		entry or remove an existing one.
401		``getnnz()``
402		Return the number of nonzero entries.
403		``todense()``
404		Return the matrix as a dense array.
405
406		The ``__getitem__`` and ``__setitem__`` methods are implemented by
407		default, and automatically select the appropriate methods from the
408		above in the cases where the item to be got or set is of the form
409		``:``, ``i,:``, ``:,j`` or ``i,j``.
410		'''
411		# methods to be implemented by subclass
412		def get_row(self, i):
413		return NotImplemented
414
415		def get_col(self, i):
416		return NotImplemented
417
418		def set_row(self, i, x):
419		return NotImplemented
420
421		def set_col(self, i, x):
422		return NotImplemented
423
424		def set_element(self, i, j, x):
425		return NotImplemented
426
427		def get_element(self, i, j):
428		return NotImplemented
429
430		def get_rows(self, rows):
431		return [self.get_row(i) for i in rows]
432
433		def get_cols(self, cols):
434		return [self.get_col(i) for i in cols]
435
436		def insert(self, i, j, x):
437		return NotImplemented
438
439		def remove(self, i, j):
440		return NotImplemented
441
442		def getnnz(self):
443		return NotImplemented
444
445		def todense(self):
446		return array([todense(r) for r in self])
447		# we support the following indexing schemes:
448		# - s[:]
449		# - s[i,:]
450		# - s[:,i]
451		# - s[i,j]
452
453		def __getitem__(self, item):
454		if isinstance(item, tuple) and isinstance(item[0], int) and item[1] == colon_slice:
455		return self.get_row(item[0])
456		if isinstance(item, slice):
	275	i, j = index
	276	except (ValueError, TypeError):
	277	raise IndexError, "invalid index"
	278
	279	if isinstance(i, slice) and isinstance(j, slice) and\
	280	(i.step is None) and (j.step is None) and\
	281	(isinstance(W, sparse.spmatrix) or isinstance(W, numpy.ndarray)):
	282	rows = self.rows[i]
	283	datas = self.data[i]
	284	j0 = j.start
	285	if isinstance(W, sparse.lil_matrix):
	286	for row, data, rowW, dataW in izip(rows, datas, W.rows, W.data):
	287	jj = bisect.bisect(row, j0) # Find the insertion point
	288	row[jj:jj] = [j0 + k for k in rowW]
	289	data[jj:jj] = dataW
	290	elif isinstance(W, ndarray):
	291	nq = W.shape[1]
	292	for row, data, rowW in izip(rows, datas, W):
	293	jj = bisect.bisect(row, j0) # Find the insertion point
	294	row[jj:jj] = range(j0, j0 + nq)
	295	data[jj:jj] = rowW
	296	elif oldscipy and isinstance(i, int) and isinstance(j, (list, tuple, numpy.ndarray)):
	297	row = dict(izip(self.rows[i], self.data[i]))
	298	try:
	299	row.update(dict(izip(j, W)))
	300	except TypeError:
	301	row.update(dict(izip(j, itertools.repeat(W))))
	302	items = row.items()
	303	items.sort()
	304	row, data = izip(*items)
	305	self.rows[i] = list(row)
	306	self.data[i] = list(data)
	307	elif isinstance(i, slice) and isinstance(j, int) and isSequenceType(W):
	308	# This corrects a bug in scipy sparse matrix as of version 0.7.0, but
	309	# it is not efficient!
	310	for w, k in izip(W, xrange(*i.indices(self.shape[0]))):
	311	sparse.lil_matrix.__setitem__(self, (k, j), w)
	312	else:
	313	sparse.lil_matrix.__setitem__(self, index, W)
	314
	315
	316	class SparseConstructionMatrix(ConstructionMatrix, SparseMatrix):
	317	'''
	318	SparseConstructionMatrix is converted to SparseConnectionMatrix.
	319	'''
	320	def __init__(self, arg, **kwds):
	321	SparseMatrix.__init__(self, arg)
	322	self.init_kwds = kwds
	323
	324	def connection_matrix(self, **additional_kwds):
	325	self.init_kwds.update(additional_kwds)
	326	return SparseConnectionMatrix(self, **self.init_kwds)
	327
	328
	329	class DynamicConstructionMatrix(ConstructionMatrix, SparseMatrix):
	330	'''
	331	DynamicConstructionMatrix is converted to DynamicConnectionMatrix.
	332	'''
	333	def __init__(self, arg, **kwds):
	334	SparseMatrix.__init__(self, arg)
	335	self.init_kwds = kwds
	336
	337	def connection_matrix(self, **additional_kwds):
	338	self.init_kwds.update(additional_kwds)
	339	return DynamicConnectionMatrix(self, **self.init_kwds)
	340
	341	# this is used to look up str->class conversions for structure=... keyword
	342	construction_matrix_register = {
	343	'dense':DenseConstructionMatrix,
	344	'sparse':SparseConstructionMatrix,
	345	'dynamic':DynamicConstructionMatrix,
	346	}
	347
	348
	349	class ConnectionMatrix(object):
	350	'''
	351	Base class for connection matrix objects
	352
	353	Connection matrix objects support a subset of the following methods:
	354
	355	``get_row(i)``, ``get_col(i)``
	356	Returns row/col ``i`` as a :class:`DenseConnectionVector` or
	357	:class:`SparseConnectionVector` as appropriate for the class.
	358	``set_row(i, val)``, ``set_col(i, val)``
	359	Sets row/col with an array, :class:`DenseConnectionVector` or
	360	:class:`SparseConnectionVector` (if supported).
	361	``get_element(i, j)``, ``set_element(i, j, val)``
	362	Gets or sets a single value.
	363	``get_rows(rows)``
	364	Returns a list of rows, should be implemented without Python
	365	function calls for efficiency if possible.
	366	``get_cols(cols)``
	367	Returns a list of cols, should be implemented without Python
	368	function calls for efficiency if possible.
	369	``insert(i,j,x)``, ``remove(i,j)``
	370	For sparse connection matrices which support it, insert a new
	371	entry or remove an existing one.
	372	``getnnz()``
	373	Return the number of nonzero entries.
	374	``todense()``
	375	Return the matrix as a dense array.
	376
	377	The ``__getitem__`` and ``__setitem__`` methods are implemented by
	378	default, and automatically select the appropriate methods from the
	379	above in the cases where the item to be got or set is of the form
	380	``:``, ``i,:``, ``:,j`` or ``i,j``.
	381	'''
	382	# methods to be implemented by subclass
	383	def get_row(self, i):
	384	return NotImplemented
	385
	386	def get_col(self, i):
	387	return NotImplemented
	388
	389	def set_row(self, i, x):
	390	return NotImplemented
	391
	392	def set_col(self, i, x):
	393	return NotImplemented
	394
	395	def set_element(self, i, j, x):
	396	return NotImplemented
	397
	398	def get_element(self, i, j):
	399	return NotImplemented
	400
	401	def get_rows(self, rows):
	402	return [self.get_row(i) for i in rows]
	403
	404	def get_cols(self, cols):
	405	return [self.get_col(i) for i in cols]
	406
	407	def insert(self, i, j, x):
	408	return NotImplemented
	409
	410	def remove(self, i, j):
	411	return NotImplemented
	412
	413	def getnnz(self):
	414	return NotImplemented
	415
	416	def todense(self):
	417	return array([todense(r) for r in self])
	418	# we support the following indexing schemes:
	419	# - s[:]
	420	# - s[i,:]
	421	# - s[:,i]
	422	# - s[i,j]
	423
	424	def __getitem__(self, item):
	425	if isinstance(item, tuple) and isinstance(item[0], int) and item[1] == colon_slice:
	426	return self.get_row(item[0])
	427	if isinstance(item, slice):
	428	if item == colon_slice:
	429	return self
	430	else:
	431	raise ValueError(str(item) + ' not supported.')
	432	if isinstance(item, int):
	433	return self.get_row(item)
	434	if isinstance(item, tuple):
	435	if len(item) != 2:
	436	raise TypeError('Only 2D indexing supported.')
	437	item_i, item_j = item
	438	if isinstance(item_i, int) and isinstance(item_j, slice):
	439	if item_j == colon_slice:
	440	return self.get_row(item_i)
	441	raise ValueError('Only ":" indexing supported.')
	442	if isinstance(item_i, slice) and isinstance(item_j, int):
	443	if item_i == colon_slice:
	444	return self.get_col(item_j)
	445	raise ValueError('Only ":" indexing supported.')
	446	if isinstance(item_i, int) and isinstance(item_j, int):
	447	return self.get_element(item_i, item_j)
	448	raise TypeError('Only (i,:), (:,j), (i,j) indexing supported.')
	449	raise TypeError('Can only get items of type slice or tuple')
	450
	451	def __setitem__(self, item, value):
	452	if isinstance(item, tuple) and isinstance(item[0], int) and item[1] == colon_slice:
	453	return self.set_row(item[0], value)
	454	if isinstance(item, slice):
	455	raise ValueError(str(item) + ' not supported.')
	456	if isinstance(item, int):
	457	return self.set_row(item, value)
	458	if isinstance(item, tuple):
	459	if len(item) != 2:
	460	raise TypeError('Only 2D indexing supported.')
	461	item_i, item_j = item
	462	if isinstance(item_i, int) and isinstance(item_j, slice):
	463	if item_j == colon_slice:
	464	return self.set_row(item_i, value)
	465	raise ValueError('Only ":" indexing supported.')
	466	if isinstance(item_i, slice) and isinstance(item_j, int):
	467	if item_i == colon_slice:
	468	return self.set_col(item_j, value)
	469	raise ValueError('Only ":" indexing supported.')
	470	if isinstance(item_i, int) and isinstance(item_j, int):
	471	return self.set_element(item_i, item_j, value)
	472	raise TypeError('Only (i,:), (:,j), (i,j) indexing supported.')
	473	raise TypeError('Can only set items of type slice or tuple')
	474
	475
	476	class DenseConnectionMatrix(ConnectionMatrix, numpy.ndarray):
	477	'''
	478	Dense connection matrix
	479
	480	See documentation for :class:`ConnectionMatrix` for details on
	481	connection matrix types.
	482
	483	This matrix implements a dense connection matrix. It is just
	484	a numpy array. The ``get_row`` and ``get_col`` methods return
	485	:class:`DenseConnectionVector`` objects.
	486	'''
	487	def __new__(subtype, data, **kwds):
	488	if 'copy' not in kwds:
	489	kwds = dict(kwds.iteritems())
	490	kwds['copy'] = False
	491	return numpy.array(data, **kwds).view(subtype)
	492
	493	def __init__(self, val, **kwds):
	494	# precompute rows and cols for fast returns by get_rows etc.
	495	self.rows = [DenseConnectionVector(numpy.ndarray.__getitem__(self, i)) for i in xrange(val.shape[0])]
	496	self.cols = [DenseConnectionVector(numpy.ndarray.__getitem__(self, (slice(None), i))) for i in xrange(val.shape[1])]
	497
	498	def get_rows(self, rows):
	499	return [self.rows[i] for i in rows]
	500
	501	def get_cols(self, cols):
	502	return [self.cols[i] for i in cols]
	503
	504	def get_row(self, i):
	505	return self.rows[i]
	506
	507	def get_col(self, i):
	508	return self.cols[i]
	509
	510	def set_row(self, i, x):
	511	numpy.ndarray.__setitem__(self, i, todense(x))
	512
	513	def set_col(self, i, x):
	514	numpy.ndarray.__setitem__(self, (colon_slice, i), todense(x))
	515	#self[:, i] = todense(x)
	516
	517	def get_element(self, i, j):
	518	return numpy.ndarray.__getitem__(self, (i, j))
	519	#return self[i,j]
	520
	521	def set_element(self, i, j, val):
	522	numpy.ndarray.__setitem__(self, (i, j), val)
	523	#self[i,j] = val
	524	insert = set_element
	525
	526	def remove(self, i, j):
	527	numpy.ndarray.__setitem__(self, (i, j), 0)
	528	#self[i, j] = 0
	529
	530
	531	class SparseConnectionMatrix(ConnectionMatrix):
	532	'''
	533	Sparse connection matrix
	534
	535	See documentation for :class:`ConnectionMatrix` for details on
	536	connection matrix types.
	537
	538	This class implements a sparse matrix with a fixed number of nonzero
	539	entries. Row access is very fast, and if the ``column_access`` keyword
	540	is ``True`` then column access is also supported (but is not as fast
	541	as row access).
	542
	543	The matrix should be initialised with a scipy sparse matrix.
	544
	545	The ``get_row`` and ``get_col`` methods return
	546	:class:`SparseConnectionVector` objects. In addition to the
	547	usual slicing operations supported, ``M[:]=val`` is supported, where
	548	``val`` must be a scalar or an array of length ``nnz``.
	549
	550	Implementation details:
	551
	552	The values are stored in an array ``alldata`` of length ``nnz`` (number
	553	of nonzero entries). The slice ``alldata[rowind[i]:rowind[i+1]]`` gives
	554	the values for row ``i``. These slices are stored in the list ``rowdata``
	555	so that ``rowdata[i]`` is the data for row ``i``. The array ``rowj[i]``
	556	gives the corresponding column ``j`` indices. For row access, the
	557	memory requirements are 12 bytes per entry (8 bytes for the float value,
	558	and 4 bytes for the column indices). The array ``allj`` of length ``nnz``
	559	gives the column ``j`` coordinates for each element in ``alldata`` (the
	560	elements of ``rowj`` are slices of this array so no extra memory is
	561	used).
	562
	563	If column access is being used, then in addition to the above there are
	564	lists ``coli`` and ``coldataindices``. For column ``j``, the array
	565	``coli[j]`` gives the row indices for the data values in column ``j``,
	566	while ``coldataindices[j]`` gives the indices in the array ``alldata``
	567	for the values in column ``j``. Column access therefore involves a
	568	copy operation rather than a slice operation. Column access increases
	569	the memory requirements to 20 bytes per entry (4 extra bytes for the
	570	row indices and 4 extra bytes for the data indices).
	571	'''
	572	def __init__(self, val, column_access=True, **kwds):
	573	self._useaccel = get_global_preference('useweave')
	574	self._cpp_compiler = get_global_preference('weavecompiler')
	575	self._extra_compile_args = ['-O3']
	576	if self._cpp_compiler == 'gcc':
	577	self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
	578	self.nnz = nnz = val.getnnz()# nnz stands for number of nonzero entries
	579	alldata = numpy.zeros(nnz)
	580	if column_access:
	581	colind = numpy.zeros(val.shape[1] + 1, dtype=int)
	582	allj = numpy.zeros(nnz, dtype=int)
	583	rowind = numpy.zeros(val.shape[0] + 1, dtype=int)
	584	rowdata = []
	585	rowj = []
	586	if column_access:
	587	coli = []
	588	coldataindices = []
	589	i = 0 # i points to the current index in the alldata array as we go through row by row
	590	for c in xrange(val.shape[0]):
	591	# extra the row values and column indices of row c of the initialising matrix
	592	# this works for any of the scipy sparse matrix formats
	593	if isinstance(val, sparse.lil_matrix):
	594	r = val.rows[c]
	595	d = val.data[c]
	596	else:
	597	sr = val[c, :]
	598	sr = sr.tolil()
	599	r = sr.rows[0]
	600	d = sr.data[0]
	601	# copy the values into the alldata array, the indices into the allj array, and
	602	# so forth
	603	rowind[c] = i
	604	alldata[i:i + len(d)] = d
	605	allj[i:i + len(r)] = r
	606	rowdata.append(alldata[i:i + len(d)])
	607	rowj.append(allj[i:i + len(r)])
	608	i = i + len(r)
	609	rowind[val.shape[0]] = i
	610	if column_access:
	611	# counts the number of nonzero elements in each column
	612	counts = zeros(val.shape[1], dtype=int)
	613	if len(allj):
	614	bincounts = numpy.bincount(allj)
	615	else:
	616	bincounts = numpy.array([], dtype=int)
	617	counts[:len(bincounts)] = bincounts # ensure that counts is the right length
	618	# two algorithms depending on whether weave is available
	619	if self._useaccel:
	620	# this algorithm just goes through one by one adding each
	621	# element to the appropriate bin whose sizes we have
	622	# precomputed. alldi will contain all the data indices
	623	# in blocks alldi[s[i]:s[i+1]] of length counts[i], and
	624	# curcdi[i] is the current offset into each block. s is
	625	# therefore just the cumulative sum of counts.
	626	curcdi = numpy.zeros(val.shape[1], dtype=int)
	627	allcoldataindices = numpy.zeros(nnz, dtype=int)
	628	colind[:] = numpy.hstack(([0], cumsum(counts)))
	629	colalli = numpy.zeros(nnz, dtype=int)
	630	numrows = val.shape[0]
	631	code = '''
	632	int i = 0;
	633	for(int k=0;k<nnz;k++)
	634	{
	635	while(k>=rowind[i+1]) i++;
	636	int j = allj[k];
	637	allcoldataindices[colind[j]+curcdi[j]] = k;
	638	colalli[colind[j]+curcdi[j]] = i;
	639	curcdi[j]++;
	640	}
	641	'''
	642	weave.inline(code, ['nnz', 'allj', 'allcoldataindices',
	643	'rowind', 'numrows',
	644	'curcdi', 'colind', 'colalli'],
	645	compiler=self._cpp_compiler,
	646	extra_compile_args=self._extra_compile_args)
	647	# now store the blocks of allcoldataindices in coldataindices and update coli too
	648	for i in xrange(len(colind) - 1):
	649	D = allcoldataindices[colind[i]:colind[i + 1]]
	650	I = colalli[colind[i]:colind[i + 1]]
	651	coldataindices.append(D)
	652	coli.append(I)
	653	else:
	654	# now allj[a] will be the columns in order, so that
	655	# the first counts[0] elements of allj[a] will be 0,
	656	# or in other words the first counts[0] elements of a
	657	# will be the data indices of the elements (i,j) with j==0
	658	# mergesort is necessary because we want the relative ordering
	659	# of the elements of a within a block to be maintained
	660	allcoldataindices = a = argsort(allj, kind='mergesort')
	661	# this defines colind so that a[colind[i]:colind[i+1]] are the data
	662	# indices where j==i
	663	colind[:] = numpy.hstack(([0], cumsum(counts)))
	664	# this computes the row index of each entry by first generating
	665	# expanded_row_indices which gives the corresponding row index
	666	# for each entry enumerated row-by-row, and then using the
	667	# array allcoldataindices to index this array to convert into
	668	# the corresponding row index for each entry enumerated
	669	# col-by-col.
	670	if len(a):
	671	expanded_row_indices = empty(len(a), dtype=int)
	672	for k, (i, j) in enumerate(zip(rowind[:-1], rowind[1:])):
	673	expanded_row_indices[i:j] = k
	674	colalli = expanded_row_indices[a]
	675	else:
	676	colalli = numpy.zeros(nnz, dtype=int)
	677	# in this loop, I are the data indices where j==i
	678	# and alli[I} are the corresponding i coordinates
	679	for i in xrange(len(colind) - 1):
	680	D = a[colind[i]:colind[i + 1]]
	681	I = colalli[colind[i]:colind[i + 1]]
	682	coldataindices.append(D)
	683	coli.append(I)
	684
	685	self.alldata = alldata
	686	self.rowdata = rowdata
	687	self.allj = allj
	688	self.rowj = rowj
	689	self.rowind = rowind
	690	self.shape = val.shape
	691	self.column_access = column_access
	692	if column_access:
	693	self.colalli = colalli
	694	self.coli = coli
	695	self.coldataindices = coldataindices
	696	self.allcoldataindices = allcoldataindices
	697	self.colind = colind
	698	self.rows = [SparseConnectionVector(self.shape[1], self.rowj[i], self.rowdata[i]) for i in xrange(self.shape[0])]
	699
	700	def getnnz(self):
	701	return self.nnz
	702
	703	def get_element(self, i, j):
	704	n = searchsorted(self.rowj[i], j)
	705	if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
	706	return 0
	707	return self.rowdata[i][n]
	708
	709	def set_element(self, i, j, x):
	710	n = searchsorted(self.rowj[i], j)
	711	if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
	712	raise ValueError('Insertion of new elements not supported for SparseConnectionMatrix.')
	713	self.rowdata[i][n] = x
	714
	715	def get_row(self, i):
	716	return self.rows[i]
	717
	718	def get_rows(self, rows):
	719	return [self.rows[i] for i in rows]
	720
	721	def get_col(self, j):
	722	if self.column_access:
	723	return SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataindices[j]])
	724	else:
	725	raise TypeError('No column access.')
	726
	727	def get_cols(self, cols):
	728	if self.column_access:
	729	return [SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataindices[j]]) for j in cols]
	730	else:
	731	raise TypeError('No column access.')
	732
	733	def set_row(self, i, val):
	734	if isinstance(val, SparseConnectionVector):
	735	if val.ind is not self.rowj[i]:
	736	if not (val.ind == self.rowj[i]).all():
	737	raise ValueError('Sparse row setting must use same indices.')
	738	self.rowdata[i][:] = val
	739	else:
	740	if isinstance(val, numpy.ndarray):
	741	val = asarray(val)
	742	self.rowdata[i][:] = val[self.rowj[i]]
	743	else:
	744	self.rowdata[i][:] = val
	745
	746	def set_col(self, j, val):
	747	if self.column_access:
	748	if isinstance(val, SparseConnectionVector):
	749	if val.ind is not self.coli[j]:
	750	if not (val.ind == self.coli[j]).all():
	751	raise ValueError('Sparse col setting must use same indices.')
	752	self.alldata[self.coldataindices[j]] = val
	753	else:
	754	if isinstance(val, numpy.ndarray):
	755	val = asarray(val)
	756	self.alldata[self.coldataindices[j]] = val[self.coli[j]]
	757	else:
	758	self.alldata[self.coldataindices[j]] = val
	759	else:
	760	raise TypeError('No column access.')
	761
	762	def __setitem__(self, item, value):
457	763	if item == colon_slice:
458		return self
459		else:
460		raise ValueError(str(item) + ' not supported.')
461		if isinstance(item, int):
462		return self.get_row(item)
463		if isinstance(item, tuple):
464		if len(item) != 2:
465		raise TypeError('Only 2D indexing supported.')
466		item_i, item_j = item
467		if isinstance(item_i, int) and isinstance(item_j, slice):
468		if item_j == colon_slice:
469		return self.get_row(item_i)
470		raise ValueError('Only ":" indexing supported.')
471		if isinstance(item_i, slice) and isinstance(item_j, int):
472		if item_i == colon_slice:
473		return self.get_col(item_j)
474		raise ValueError('Only ":" indexing supported.')
475		if isinstance(item_i, int) and isinstance(item_j, int):
476		return self.get_element(item_i, item_j)
477		raise TypeError('Only (i,:), (:,j), (i,j) indexing supported.')
478		raise TypeError('Can only get items of type slice or tuple')
479
480		def __setitem__(self, item, value):
481		if isinstance(item, tuple) and isinstance(item[0], int) and item[1] == colon_slice:
482		return self.set_row(item[0], value)
483		if isinstance(item, slice):
484		raise ValueError(str(item) + ' not supported.')
485		if isinstance(item, int):
486		return self.set_row(item, value)
487		if isinstance(item, tuple):
488		if len(item) != 2:
489		raise TypeError('Only 2D indexing supported.')
490		item_i, item_j = item
491		if isinstance(item_i, int) and isinstance(item_j, slice):
492		if item_j == colon_slice:
493		return self.set_row(item_i, value)
494		raise ValueError('Only ":" indexing supported.')
495		if isinstance(item_i, slice) and isinstance(item_j, int):
496		if item_i == colon_slice:
497		return self.set_col(item_j, value)
498		raise ValueError('Only ":" indexing supported.')
499		if isinstance(item_i, int) and isinstance(item_j, int):
500		return self.set_element(item_i, item_j, value)
501		raise TypeError('Only (i,:), (:,j), (i,j) indexing supported.')
502		raise TypeError('Can only set items of type slice or tuple')
503
504
505		class DenseConnectionMatrix(ConnectionMatrix, numpy.ndarray):
506		'''
507		Dense connection matrix
508
509		See documentation for :class:`ConnectionMatrix` for details on
510		connection matrix types.
511
512		This matrix implements a dense connection matrix. It is just
513		a numpy array. The ``get_row`` and ``get_col`` methods return
514		:class:`DenseConnectionVector`` objects.
515		'''
516		def __new__(subtype, data, **kwds):
517		if 'copy' not in kwds:
518		kwds = dict(kwds.iteritems())
519		kwds['copy'] = False
520		return numpy.array(data, **kwds).view(subtype)
521
522		def __init__(self, val, **kwds):
523		# precompute rows and cols for fast returns by get_rows etc.
524		self.rows = [DenseConnectionVector(numpy.ndarray.__getitem__(self, i)) for i in xrange(val.shape[0])]
525		self.cols = [DenseConnectionVector(numpy.ndarray.__getitem__(self, (slice(None), i))) for i in xrange(val.shape[1])]
526
527		def get_rows(self, rows):
528		return [self.rows[i] for i in rows]
529
530		def get_cols(self, cols):
531		return [self.cols[i] for i in cols]
532
533		def get_row(self, i):
534		return self.rows[i]
535
536		def get_col(self, i):
537		return self.cols[i]
538
539		def set_row(self, i, x):
540		numpy.ndarray.__setitem__(self, i, todense(x))
541
542		def set_col(self, i, x):
543		numpy.ndarray.__setitem__(self, (colon_slice, i), todense(x))
544		#self[:, i] = todense(x)
545
546		def get_element(self, i, j):
547		return numpy.ndarray.__getitem__(self, (i, j))
548		#return self[i,j]
549
550		def set_element(self, i, j, val):
551		numpy.ndarray.__setitem__(self, (i, j), val)
552		#self[i,j] = val
553		insert = set_element
554
555		def remove(self, i, j):
556		numpy.ndarray.__setitem__(self, (i, j), 0)
557		#self[i, j] = 0
558
559
560		class SparseConnectionMatrix(ConnectionMatrix):
561		'''
562		Sparse connection matrix
563
564		See documentation for :class:`ConnectionMatrix` for details on
565		connection matrix types.
566
567		This class implements a sparse matrix with a fixed number of nonzero
568		entries. Row access is very fast, and if the ``column_access`` keyword
569		is ``True`` then column access is also supported (but is not as fast
570		as row access).
571
572		The matrix should be initialised with a scipy sparse matrix.
573
574		The ``get_row`` and ``get_col`` methods return
575		:class:`SparseConnectionVector` objects. In addition to the
576		usual slicing operations supported, ``M[:]=val`` is supported, where
577		``val`` must be a scalar or an array of length ``nnz``.
578
579		Implementation details:
580
581		The values are stored in an array ``alldata`` of length ``nnz`` (number
582		of nonzero entries). The slice ``alldata[rowind[i]:rowind[i+1]]`` gives
583		the values for row ``i``. These slices are stored in the list ``rowdata``
584		so that ``rowdata[i]`` is the data for row ``i``. The array ``rowj[i]``
585		gives the corresponding column ``j`` indices. For row access, the
586		memory requirements are 12 bytes per entry (8 bytes for the float value,
587		and 4 bytes for the column indices). The array ``allj`` of length ``nnz``
588		gives the column ``j`` coordinates for each element in ``alldata`` (the
589		elements of ``rowj`` are slices of this array so no extra memory is
590		used).
591
592		If column access is being used, then in addition to the above there are
593		lists ``coli`` and ``coldataindices``. For column ``j``, the array
594		``coli[j]`` gives the row indices for the data values in column ``j``,
595		while ``coldataindices[j]`` gives the indices in the array ``alldata``
596		for the values in column ``j``. Column access therefore involves a
597		copy operation rather than a slice operation. Column access increases
598		the memory requirements to 20 bytes per entry (4 extra bytes for the
599		row indices and 4 extra bytes for the data indices).
600		'''
601		def __init__(self, val, column_access=True, **kwds):
602		self._useaccel = get_global_preference('useweave')
603		self._cpp_compiler = get_global_preference('weavecompiler')
604		self._extra_compile_args = ['-O3']
605		if self._cpp_compiler == 'gcc':
606		self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
607		self.nnz = nnz = val.getnnz()# nnz stands for number of nonzero entries
608		alldata = numpy.zeros(nnz)
609		if column_access:
610		colind = numpy.zeros(val.shape[1] + 1, dtype=int)
611		allj = numpy.zeros(nnz, dtype=int)
612		rowind = numpy.zeros(val.shape[0] + 1, dtype=int)
613		rowdata = []
614		rowj = []
615		if column_access:
616		coli = []
617		coldataindices = []
618		i = 0 # i points to the current index in the alldata array as we go through row by row
619		for c in xrange(val.shape[0]):
620		# extra the row values and column indices of row c of the initialising matrix
621		# this works for any of the scipy sparse matrix formats
622		if isinstance(val, sparse.lil_matrix):
623		r = val.rows[c]
624		d = val.data[c]
625		else:
626		sr = val[c, :]
627		sr = sr.tolil()
628		r = sr.rows[0]
629		d = sr.data[0]
630		# copy the values into the alldata array, the indices into the allj array, and
631		# so forth
632		rowind[c] = i
633		alldata[i:i + len(d)] = d
634		allj[i:i + len(r)] = r
635		rowdata.append(alldata[i:i + len(d)])
636		rowj.append(allj[i:i + len(r)])
637		i = i + len(r)
638		rowind[val.shape[0]] = i
639		if column_access:
	764	self.alldata[:] = value
	765	else:
	766	ConnectionMatrix.__setitem__(self, item, value)
	767
	768
	769	class DynamicConnectionMatrix(ConnectionMatrix):
	770	'''
	771	Dynamic (sparse) connection matrix
	772
	773	See documentation for :class:`ConnectionMatrix` for details on
	774	connection matrix types.
	775
	776	This class implements a sparse matrix with a variable number of nonzero
	777	entries. Row access and column access are provided, but are not as fast
	778	as for :class:`SparseConnectionMatrix`.
	779
	780	The matrix should be initialised with a scipy sparse matrix.
	781
	782	The ``get_row`` and ``get_col`` methods return
	783	:class:`SparseConnectionVector` objects. In addition to the
	784	usual slicing operations supported, ``M[:]=val`` is supported, where
	785	``val`` must be a scalar or an array of length ``nnz``.
	786
	787	Implementation details
	788
	789	The values are stored in an array ``alldata`` of length ``nnzmax`` (maximum
	790	number of nonzero entries). This is a dynamic array, see:
	791
	792	http://en.wikipedia.org/wiki/Dynamic_array
	793
	794	You can set the resizing constant with the argument ``dynamic_array_const``.
	795	Normally the default value 2 is fine but if memory is a worry it could be
	796	made smaller.
	797
	798	Rows and column point in to this data array, and the list ``rowj`` consists
	799	of an array of column indices for each row, with ``coli`` containing arrays
	800	of row indices for each column. Similarly, ``rowdataind`` and ``coldataind``
	801	consist of arrays of pointers to the indices in the ``alldata`` array.
	802	'''
	803	def __init__(self, val, nnzmax=None, dynamic_array_const=2, **kwds):
	804	self.shape = val.shape
	805	self.dynamic_array_const = dynamic_array_const
	806	if nnzmax is None or nnzmax < val.getnnz():
	807	nnzmax = val.getnnz()
	808	self.nnzmax = nnzmax
	809	self.nnz = val.getnnz()
	810	self.alldata = numpy.zeros(nnzmax)
	811	self.unusedinds = range(self.nnz, self.nnzmax)
	812	i = 0
	813	self.rowj = []
	814	self.rowdataind = []
	815	alli = zeros(self.nnz, dtype=int)
	816	allj = zeros(self.nnz, dtype=int)
	817	for c in xrange(val.shape[0]):
	818	# extra the row values and column indices of row c of the initialising matrix
	819	# this works for any of the scipy sparse matrix formats
	820	if isinstance(val, sparse.lil_matrix):
	821	r = val.rows[c]
	822	d = val.data[c]
	823	else:
	824	sr = val[c, :]
	825	sr = sr.tolil()
	826	r = sr.rows[0]
	827	d = sr.data[0]
	828	self.alldata[i:i + len(d)] = d
	829	self.rowj.append(array(r, dtype=int))
	830	self.rowdataind.append(arange(i, i + len(d)))
	831	allj[i:i + len(d)] = r
	832	alli[i:i + len(d)] = c
	833	i += len(d)
	834	# now update the coli and coldataind variables
	835	self.coli = []
	836	self.coldataind = []
640	837	# counts the number of nonzero elements in each column
641		counts = zeros(val.shape[1], dtype=int)
642		if len(allj):
643		bincounts = numpy.bincount(allj)
644		else:
645		bincounts = numpy.array([], dtype=int)
646		counts[:len(bincounts)] = bincounts # ensure that counts is the right length
647		# two algorithms depending on whether weave is available
648		if self._useaccel:
649		# this algorithm just goes through one by one adding each
650		# element to the appropriate bin whose sizes we have
651		# precomputed. alldi will contain all the data indices
652		# in blocks alldi[s[i]:s[i+1]] of length counts[i], and
653		# curcdi[i] is the current offset into each block. s is
654		# therefore just the cumulative sum of counts.
655		curcdi = numpy.zeros(val.shape[1], dtype=int)
656		allcoldataindices = numpy.zeros(nnz, dtype=int)
657		colind[:] = numpy.hstack(([0], cumsum(counts)))
658		colalli = numpy.zeros(nnz, dtype=int)
659		numrows = val.shape[0]
660		code = '''
661		int i = 0;
662		for(int k=0;k<nnz;k++)
663		{
664		while(k>=rowind[i+1]) i++;
665		int j = allj[k];
666		allcoldataindices[colind[j]+curcdi[j]] = k;
667		colalli[colind[j]+curcdi[j]] = i;
668		curcdi[j]++;
669		}
670		'''
671		weave.inline(code, ['nnz', 'allj', 'allcoldataindices',
672		'rowind', 'numrows',
673		'curcdi', 'colind', 'colalli'],
674		compiler=self._cpp_compiler,
675		extra_compile_args=self._extra_compile_args)
676		# now store the blocks of allcoldataindices in coldataindices and update coli too
677		for i in xrange(len(colind) - 1):
678		D = allcoldataindices[colind[i]:colind[i + 1]]
679		I = colalli[colind[i]:colind[i + 1]]
680		coldataindices.append(D)
681		coli.append(I)
682		else:
683		# now allj[a] will be the columns in order, so that
684		# the first counts[0] elements of allj[a] will be 0,
685		# or in other words the first counts[0] elements of a
686		# will be the data indices of the elements (i,j) with j==0
687		# mergesort is necessary because we want the relative ordering
688		# of the elements of a within a block to be maintained
689		allcoldataindices = a = argsort(allj, kind='mergesort')
690		# this defines colind so that a[colind[i]:colind[i+1]] are the data
691		# indices where j==i
692		colind[:] = numpy.hstack(([0], cumsum(counts)))
693		# this computes the row index of each entry by first generating
694		# expanded_row_indices which gives the corresponding row index
695		# for each entry enumerated row-by-row, and then using the
696		# array allcoldataindices to index this array to convert into
697		# the corresponding row index for each entry enumerated
698		# col-by-col.
699		if len(a):
700		expanded_row_indices = empty(len(a), dtype=int)
701		for k, (i, j) in enumerate(zip(rowind[:-1], rowind[1:])):
702		expanded_row_indices[i:j] = k
703		colalli = expanded_row_indices[a]
	838	if numpy.__version__ >= '1.3.0':
	839	counts = numpy.histogram(allj, numpy.arange(val.shape[1] + 1, dtype=int))[0]
	840	else:
	841	counts = numpy.histogram(allj, numpy.arange(val.shape[1] + 1, dtype=int), new=True)[0]
	842	# now we have to go through one by one unfortunately, and so we keep curcdi, the
	843	# current column data index for each column
	844	curcdi = numpy.zeros(val.shape[1], dtype=int)
	845	# initialise the memory for the column data indices
	846	for j in xrange(val.shape[1]):
	847	self.coldataind.append(numpy.zeros(counts[j], dtype=int))
	848	# one by one for every element, update the dataindices and curcdi data pointers
	849	for i, j in enumerate(allj):
	850	self.coldataind[j][curcdi[j]] = i
	851	curcdi[j] += 1
	852	for j in xrange(val.shape[1]):
	853	self.coli.append(alli[self.coldataind[j]])
	854
	855	def getnnz(self):
	856	return self.nnz
	857
	858	def insert(self, i, j, x):
	859	n = searchsorted(self.rowj[i], j)
	860	if n < len(self.rowj[i]) and self.rowj[i][n] == j:
	861	self.alldata[self.rowdataind[i][n]] = x
	862	return
	863	m = searchsorted(self.coli[j], i)
	864	if self.nnz == self.nnzmax:
	865	# reallocate memory using a dynamic array structure (amortized O(1) cost for append)
	866	newnnzmax = int(self.nnzmax * self.dynamic_array_const)
	867	if newnnzmax <= self.nnzmax:
	868	newnnzmax += 1
	869	if newnnzmax > self.shape[0] * self.shape[1]:
	870	newnnzmax = self.shape[0] * self.shape[1]
	871	self.alldata = hstack((self.alldata, numpy.zeros(newnnzmax - self.nnzmax, dtype=self.alldata.dtype)))
	872	self.unusedinds.extend(range(self.nnz, newnnzmax))
	873	self.nnzmax = newnnzmax
	874	newind = self.unusedinds.pop(-1)
	875	self.alldata[newind] = x
	876	self.nnz += 1
	877	# update row
	878	newrowj = numpy.zeros(len(self.rowj[i]) + 1, dtype=int)
	879	newrowj[:n] = self.rowj[i][:n]
	880	newrowj[n] = j
	881	newrowj[n + 1:] = self.rowj[i][n:]
	882	self.rowj[i] = newrowj
	883	newrowdataind = numpy.zeros(len(self.rowdataind[i]) + 1, dtype=int)
	884	newrowdataind[:n] = self.rowdataind[i][:n]
	885	newrowdataind[n] = newind
	886	newrowdataind[n + 1:] = self.rowdataind[i][n:]
	887	self.rowdataind[i] = newrowdataind
	888	# update col
	889	newcoli = numpy.zeros(len(self.coli[j]) + 1, dtype=int)
	890	newcoli[:m] = self.coli[j][:m]
	891	newcoli[m] = i
	892	newcoli[m + 1:] = self.coli[j][m:]
	893	self.coli[j] = newcoli
	894	newcoldataind = numpy.zeros(len(self.coldataind[j]) + 1, dtype=int)
	895	newcoldataind[:m] = self.coldataind[j][:m]
	896	newcoldataind[m] = newind
	897	newcoldataind[m + 1:] = self.coldataind[j][m:]
	898	self.coldataind[j] = newcoldataind
	899
	900	def remove(self, i, j):
	901	n = searchsorted(self.rowj[i], j)
	902	if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
	903	raise ValueError('No element to remove at position ' + str(i, j))
	904	oldind = self.rowdataind[i][n]
	905	self.unusedinds.append(oldind)
	906	self.nnz -= 1
	907	m = searchsorted(self.coli[j], i)
	908	# update row
	909	newrowj = numpy.zeros(len(self.rowj[i]) - 1, dtype=int)
	910	newrowj[:n] = self.rowj[i][:n]
	911	newrowj[n:] = self.rowj[i][n + 1:]
	912	self.rowj[i] = newrowj
	913	newrowdataind = numpy.zeros(len(self.rowdataind[i]) - 1, dtype=int)
	914	newrowdataind[:n] = self.rowdataind[i][:n]
	915	newrowdataind[n:] = self.rowdataind[i][n + 1:]
	916	self.rowdataind[i] = newrowdataind
	917	# update col
	918	newcoli = numpy.zeros(len(self.coli[j]) - 1, dtype=int)
	919	newcoli[:m] = self.coli[j][:m]
	920	newcoli[m:] = self.coli[j][m + 1:]
	921	self.coli[j] = newcoli
	922	newcoldataind = numpy.zeros(len(self.coldataind[j]) - 1, dtype=int)
	923	newcoldataind[:m] = self.coldataind[j][:m]
	924	newcoldataind[m:] = self.coldataind[j][m + 1:]
	925	self.coldataind[j] = newcoldataind
	926
	927	def get_element(self, i, j):
	928	n = searchsorted(self.rowj[i], j)
	929	if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
	930	return 0
	931	return self.alldata[self.rowdataind[i][n]]
	932
	933	set_element = insert
	934
	935	def get_row(self, i):
	936	return SparseConnectionVector(self.shape[1], self.rowj[i], self.alldata[self.rowdataind[i]])
	937
	938	def get_rows(self, rows):
	939	return [SparseConnectionVector(self.shape[1], self.rowj[i], self.alldata[self.rowdataind[i]]) for i in rows]
	940
	941	def get_col(self, j):
	942	return SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataind[j]])
	943
	944	def get_cols(self, cols):
	945	return [SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataind[j]]) for j in cols]
	946
	947	def set_row(self, i, val):
	948	if isinstance(val, SparseConnectionVector):
	949	if val.ind is not self.rowj[i]:
	950	if not (val.ind == self.rowj[i]).all():
	951	raise ValueError('Sparse row setting must use same indices.')
	952	self.alldata[self.rowdataind[i]] = val
	953	else:
	954	if isinstance(val, numpy.ndarray):
	955	val = asarray(val)
	956	self.alldata[self.rowdataind[i]] = val[self.rowj[i]]
704	957	else:
705		colalli = numpy.zeros(nnz, dtype=int)
706		# in this loop, I are the data indices where j==i
707		# and alli[I} are the corresponding i coordinates
708		for i in xrange(len(colind) - 1):
709		D = a[colind[i]:colind[i + 1]]
710		I = colalli[colind[i]:colind[i + 1]]
711		coldataindices.append(D)
712		coli.append(I)
713
714		self.alldata = alldata
715		self.rowdata = rowdata
716		self.allj = allj
717		self.rowj = rowj
718		self.rowind = rowind
719		self.shape = val.shape
720		self.column_access = column_access
721		if column_access:
722		self.colalli = colalli
723		self.coli = coli
724		self.coldataindices = coldataindices
725		self.allcoldataindices = allcoldataindices
726		self.colind = colind
727		self.rows = [SparseConnectionVector(self.shape[1], self.rowj[i], self.rowdata[i]) for i in xrange(self.shape[0])]
728
729		def getnnz(self):
730		return self.nnz
731
732		def get_element(self, i, j):
733		n = searchsorted(self.rowj[i], j)
734		if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
735		return 0
736		return self.rowdata[i][n]
737
738		def set_element(self, i, j, x):
739		n = searchsorted(self.rowj[i], j)
740		if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
741		raise ValueError('Insertion of new elements not supported for SparseConnectionMatrix.')
742		self.rowdata[i][n] = x
743
744		def get_row(self, i):
745		return self.rows[i]
746
747		def get_rows(self, rows):
748		return [self.rows[i] for i in rows]
749
750		def get_col(self, j):
751		if self.column_access:
752		return SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataindices[j]])
753		else:
754		raise TypeError('No column access.')
755
756		def get_cols(self, cols):
757		if self.column_access:
758		return [SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataindices[j]]) for j in cols]
759		else:
760		raise TypeError('No column access.')
761
762		def set_row(self, i, val):
763		if isinstance(val, SparseConnectionVector):
764		if val.ind is not self.rowj[i]:
765		if not (val.ind == self.rowj[i]).all():
766		raise ValueError('Sparse row setting must use same indices.')
767		self.rowdata[i][:] = val
768		else:
769		if isinstance(val, numpy.ndarray):
770		val = asarray(val)
771		self.rowdata[i][:] = val[self.rowj[i]]
772		else:
773		self.rowdata[i][:] = val
774
775		def set_col(self, j, val):
776		if self.column_access:
	958	self.alldata[self.rowdataind[i]] = val
	959
	960	def set_col(self, j, val):
777	961	if isinstance(val, SparseConnectionVector):
778	962	if val.ind is not self.coli[j]:
779	963	if not (val.ind == self.coli[j]).all():
780		raise ValueError('Sparse ~~col~~ setting must use same indices.')
781		self.alldata[self.coldataind~~ices~~[j]] = val
	964	raise ValueError('Sparse row setting must use same indices.')
	965	self.alldata[self.coldataind[j]] = val
782	966	else:
783	967	if isinstance(val, numpy.ndarray):
784	968	val = asarray(val)
785		self.alldata[self.coldataind~~ices~~[j]] = val[self.coli[j]]
	969	self.alldata[self.coldataind[j]] = val[self.coli[j]]
786	970	else:
787		self.alldata[self.coldataindices[j]] = val
788		else:
789		raise TypeError('No column access.')
790
791		def __setitem__(self, item, value):
792		if item == colon_slice:
793		self.alldata[:] = value
794		else:
795		ConnectionMatrix.__setitem__(self, item, value)
796
797
798		class DynamicConnectionMatrix(ConnectionMatrix):
799		'''
800		Dynamic (sparse) connection matrix
801
802		See documentation for :class:`ConnectionMatrix` for details on
803		connection matrix types.
804
805		This class implements a sparse matrix with a variable number of nonzero
806		entries. Row access and column access are provided, but are not as fast
807		as for :class:`SparseConnectionMatrix`.
808
809		The matrix should be initialised with a scipy sparse matrix.
810
811		The ``get_row`` and ``get_col`` methods return
812		:class:`SparseConnectionVector` objects. In addition to the
813		usual slicing operations supported, ``M[:]=val`` is supported, where
814		``val`` must be a scalar or an array of length ``nnz``.
815
816		Implementation details
817
818		The values are stored in an array ``alldata`` of length ``nnzmax`` (maximum
819		number of nonzero entries). This is a dynamic array, see:
820
821		http://en.wikipedia.org/wiki/Dynamic_array
822
823		You can set the resizing constant with the argument ``dynamic_array_const``.
824		Normally the default value 2 is fine but if memory is a worry it could be
825		made smaller.
826
827		Rows and column point in to this data array, and the list ``rowj`` consists
828		of an array of column indices for each row, with ``coli`` containing arrays
829		of row indices for each column. Similarly, ``rowdataind`` and ``coldataind``
830		consist of arrays of pointers to the indices in the ``alldata`` array.
831		'''
832		def __init__(self, val, nnzmax=None, dynamic_array_const=2, **kwds):
833		self.shape = val.shape
834		self.dynamic_array_const = dynamic_array_const
835		if nnzmax is None or nnzmax < val.getnnz():
836		nnzmax = val.getnnz()
837		self.nnzmax = nnzmax
838		self.nnz = val.getnnz()
839		self.alldata = numpy.zeros(nnzmax)
840		self.unusedinds = range(self.nnz, self.nnzmax)
841		i = 0
842		self.rowj = []
843		self.rowdataind = []
844		alli = zeros(self.nnz, dtype=int)
845		allj = zeros(self.nnz, dtype=int)
846		for c in xrange(val.shape[0]):
847		# extra the row values and column indices of row c of the initialising matrix
848		# this works for any of the scipy sparse matrix formats
849		if isinstance(val, sparse.lil_matrix):
850		r = val.rows[c]
851		d = val.data[c]
852		else:
853		sr = val[c, :]
854		sr = sr.tolil()
855		r = sr.rows[0]
856		d = sr.data[0]
857		self.alldata[i:i + len(d)] = d
858		self.rowj.append(array(r, dtype=int))
859		self.rowdataind.append(arange(i, i + len(d)))
860		allj[i:i + len(d)] = r
861		alli[i:i + len(d)] = c
862		i += len(d)
863		# now update the coli and coldataind variables
864		self.coli = []
865		self.coldataind = []
866		# counts the number of nonzero elements in each column
867		if numpy.__version__ >= '1.3.0':
868		counts = numpy.histogram(allj, numpy.arange(val.shape[1] + 1, dtype=int))[0]
869		else:
870		counts = numpy.histogram(allj, numpy.arange(val.shape[1] + 1, dtype=int), new=True)[0]
871		# now we have to go through one by one unfortunately, and so we keep curcdi, the
872		# current column data index for each column
873		curcdi = numpy.zeros(val.shape[1], dtype=int)
874		# initialise the memory for the column data indices
875		for j in xrange(val.shape[1]):
876		self.coldataind.append(numpy.zeros(counts[j], dtype=int))
877		# one by one for every element, update the dataindices and curcdi data pointers
878		for i, j in enumerate(allj):
879		self.coldataind[j][curcdi[j]] = i
880		curcdi[j] += 1
881		for j in xrange(val.shape[1]):
882		self.coli.append(alli[self.coldataind[j]])
883
884		def getnnz(self):
885		return self.nnz
886
887		def insert(self, i, j, x):
888		n = searchsorted(self.rowj[i], j)
889		if n < len(self.rowj[i]) and self.rowj[i][n] == j:
890		self.alldata[self.rowdataind[i][n]] = x
891		return
892		m = searchsorted(self.coli[j], i)
893		if self.nnz == self.nnzmax:
894		# reallocate memory using a dynamic array structure (amortized O(1) cost for append)
895		newnnzmax = int(self.nnzmax * self.dynamic_array_const)
896		if newnnzmax <= self.nnzmax:
897		newnnzmax += 1
898		if newnnzmax > self.shape[0] * self.shape[1]:
899		newnnzmax = self.shape[0] * self.shape[1]
900		self.alldata = hstack((self.alldata, numpy.zeros(newnnzmax - self.nnzmax, dtype=self.alldata.dtype)))
901		self.unusedinds.extend(range(self.nnz, newnnzmax))
902		self.nnzmax = newnnzmax
903		newind = self.unusedinds.pop(-1)
904		self.alldata[newind] = x
905		self.nnz += 1
906		# update row
907		newrowj = numpy.zeros(len(self.rowj[i]) + 1, dtype=int)
908		newrowj[:n] = self.rowj[i][:n]
909		newrowj[n] = j
910		newrowj[n + 1:] = self.rowj[i][n:]
911		self.rowj[i] = newrowj
912		newrowdataind = numpy.zeros(len(self.rowdataind[i]) + 1, dtype=int)
913		newrowdataind[:n] = self.rowdataind[i][:n]
914		newrowdataind[n] = newind
915		newrowdataind[n + 1:] = self.rowdataind[i][n:]
916		self.rowdataind[i] = newrowdataind
917		# update col
918		newcoli = numpy.zeros(len(self.coli[j]) + 1, dtype=int)
919		newcoli[:m] = self.coli[j][:m]
920		newcoli[m] = i
921		newcoli[m + 1:] = self.coli[j][m:]
922		self.coli[j] = newcoli
923		newcoldataind = numpy.zeros(len(self.coldataind[j]) + 1, dtype=int)
924		newcoldataind[:m] = self.coldataind[j][:m]
925		newcoldataind[m] = newind
926		newcoldataind[m + 1:] = self.coldataind[j][m:]
927		self.coldataind[j] = newcoldataind
928
929		def remove(self, i, j):
930		n = searchsorted(self.rowj[i], j)
931		if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
932		raise ValueError('No element to remove at position ' + str(i, j))
933		oldind = self.rowdataind[i][n]
934		self.unusedinds.append(oldind)
935		self.nnz -= 1
936		m = searchsorted(self.coli[j], i)
937		# update row
938		newrowj = numpy.zeros(len(self.rowj[i]) - 1, dtype=int)
939		newrowj[:n] = self.rowj[i][:n]
940		newrowj[n:] = self.rowj[i][n + 1:]
941		self.rowj[i] = newrowj
942		newrowdataind = numpy.zeros(len(self.rowdataind[i]) - 1, dtype=int)
943		newrowdataind[:n] = self.rowdataind[i][:n]
944		newrowdataind[n:] = self.rowdataind[i][n + 1:]
945		self.rowdataind[i] = newrowdataind
946		# update col
947		newcoli = numpy.zeros(len(self.coli[j]) - 1, dtype=int)
948		newcoli[:m] = self.coli[j][:m]
949		newcoli[m:] = self.coli[j][m + 1:]
950		self.coli[j] = newcoli
951		newcoldataind = numpy.zeros(len(self.coldataind[j]) - 1, dtype=int)
952		newcoldataind[:m] = self.coldataind[j][:m]
953		newcoldataind[m:] = self.coldataind[j][m + 1:]
954		self.coldataind[j] = newcoldataind
955
956		def get_element(self, i, j):
957		n = searchsorted(self.rowj[i], j)
958		if n >= len(self.rowj[i]) or self.rowj[i][n] != j:
959		return 0
960		return self.alldata[self.rowdataind[i][n]]
961
962		set_element = insert
963
964		def get_row(self, i):
965		return SparseConnectionVector(self.shape[1], self.rowj[i], self.alldata[self.rowdataind[i]])
966
967		def get_rows(self, rows):
968		return [SparseConnectionVector(self.shape[1], self.rowj[i], self.alldata[self.rowdataind[i]]) for i in rows]
969
970		def get_col(self, j):
971		return SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataind[j]])
972
973		def get_cols(self, cols):
974		return [SparseConnectionVector(self.shape[0], self.coli[j], self.alldata[self.coldataind[j]]) for j in cols]
975
976		def set_row(self, i, val):
977		if isinstance(val, SparseConnectionVector):
978		if val.ind is not self.rowj[i]:
979		if not (val.ind == self.rowj[i]).all():
980		raise ValueError('Sparse row setting must use same indices.')
981		self.alldata[self.rowdataind[i]] = val
982		else:
983		if isinstance(val, numpy.ndarray):
984		val = asarray(val)
985		self.alldata[self.rowdataind[i]] = val[self.rowj[i]]
986		else:
987		self.alldata[self.rowdataind[i]] = val
988
989		def set_col(self, j, val):
990		if isinstance(val, SparseConnectionVector):
991		if val.ind is not self.coli[j]:
992		if not (val.ind == self.coli[j]).all():
993		raise ValueError('Sparse row setting must use same indices.')
994		self.alldata[self.coldataind[j]] = val
995		else:
996		if isinstance(val, numpy.ndarray):
997		val = asarray(val)
998		self.alldata[self.coldataind[j]] = val[self.coli[j]]
999		else:
1000		self.alldata[self.coldataind[j]] = val
1001
1002		def __setitem__(self, item, value):
1003		if item == colon_slice:
1004		self.alldata[:self.nnz] = value
1005		else:
1006		ConnectionMatrix.__setitem__(self, item, value)
1007
1008
1009		class Connection(magic.InstanceTracker, ObjectContainer):
1010		'''
1011		Mechanism for propagating spikes from one group to another
1012
1013		A Connection object declares that when spikes in a source
1014		group are generated, certain neurons in the target group
1015		should have a value added to specific states. See
1016		Tutorial 2: Connections to understand this better.
1017
1018		With arguments:
1019
1020		``source``
1021		The group from which spikes will be propagated.
1022		``target``
1023		The group to which spikes will be propagated.
1024		``state``
1025		The state variable name or number that spikes will be
1026		propagated to in the target group.
1027		``delay``
1028		The delay between a spike being generated at the source
1029		and received at the target. Depending on the type of ``delay``
1030		it has different effects. If ``delay`` is a scalar value, then
1031		the connection will be initialised with all neurons having
1032		that delay. For very long delays, this may raise an error. If
1033		``delay=True`` then the connection will be initialised as a
1034		:class:`DelayConnection`, allowing heterogeneous delays (a
1035		different delay for each synapse). ``delay`` can also be a
1036		pair ``(min,max)`` or a function of one or two variables, in
1037		both cases it will be initialised as a :class:`DelayConnection`,
1038		see the documentation for that class for details. Note that in
1039		these cases, initialisation of delays will only have the
1040		intended effect if used with the ``weight`` and ``sparseness``
1041		arguments below.
1042		``max_delay``
1043		If you are using a connection with heterogeneous delays, specify
1044		this to set the maximum allowed delay (smaller values use less
1045		memory). The default is 5ms.
1046		``modulation``
1047		The state variable name from the source group that scales
1048		the synaptic weights (for short-term synaptic plasticity).
1049		``structure``
1050		Data structure: ``sparse`` (default), ``dense`` or
1051		``dynamic``. See below for more information on structures.
1052		``weight``
1053		If specified, the connection matrix will be initialised with
1054		values specified by ``weight``, which can be any of the values
1055		allowed in the methods `connect*`` below.
1056		``sparseness``
1057		If ``weight`` is specified and ``sparseness`` is not, a full
1058		connection is assumed, otherwise random connectivity with this
1059		level of sparseness is assumed.
1060
1061		Methods
1062
1063		``connect_random(P,Q,p[,weight=1[,fixed=False[,seed=None]]])``
1064		Connects each neuron in ``P`` to each neuron in ``Q`` with independent
1065		probability ``p`` and weight ``weight`` (this is the amount that
1066		gets added to the target state variable). If ``fixed`` is True, then
1067		the number of presynaptic neurons per neuron is constant. If ``seed``
1068		is given, it is used as the seed to the random number generators, for
1069		exactly repeatable results.
1070		``connect_full(P,Q[,weight=1])``
1071		Connect every neuron in ``P`` to every neuron in ``Q`` with the given
1072		weight.
1073		``connect_one_to_one(P,Q)``
1074		If ``P`` and ``Q`` have the same number of neurons then neuron ``i``
1075		in ``P`` will be connected to neuron ``i`` in ``Q`` with weight 1.
1076		``connect(P,Q,W)``
1077		You can specify a matrix of weights directly (can be in any format
1078		recognised by NumPy). Note that due to internal implementation details,
1079		passing a full matrix rather than a sparse one may slow down your code
1080		(because zeros will be propagated as well as nonzero values).
1081		WARNING: No unit checking is done at the moment.
1082
1083		Additionally, you can directly access the matrix of weights by writing::
1084
1085		C = Connection(P,Q)
1086		print C[i,j]
1087		C[i,j] = ...
1088
1089		Where here ``i`` is the source neuron and ``j`` is the target neuron.
1090		Note: if ``C[i,j]`` should be zero, it is more efficient not to write
1091		``C[i,j]=0``, if you write this then when neuron ``i`` fires all the
1092		targets will have the value 0 added to them rather than just the
1093		nonzero ones.
1094		WARNING: No unit checking is currently done if you use this method.
1095		Take care to set the right units.
1096
1097		Connection matrix structures
1098
1099		Brian currently features three types of connection matrix structures,
1100		each of which is suited for different situations. Brian has two stages
1101		of connection matrix. The first is the construction stage, used for
1102		building a weight matrix. This stage is optimised for the construction
1103		of matrices, with lots of features, but would be slow for runtime
1104		behaviour. Consequently, the second stage is the connection stage,
1105		used when Brian is being run. The connection stage is optimised for
1106		run time behaviour, but many features which are useful for construction
1107		are absent (e.g. the ability to add or remove synapses). Conversion
1108		between construction and connection stages is done by the
1109		``compress()`` method of :class:`Connection` which is called
1110		automatically when it is used for the first time.
1111
1112		The structures are:
1113
1114		``dense``
1115		A dense matrix. Allows runtime modification of all values. If
1116		connectivity is close to being dense this is probably the most
1117		efficient, but in most cases it is less efficient. In addition,
1118		a dense connection matrix will often do the wrong thing if
1119		using STDP. Because a synapse will be considered to exist but
1120		with weight 0, STDP will be able to create new synapses where
1121		there were previously none. Memory requirements are ``8NM``
1122		bytes where ``(N,M)`` are the dimensions. (A ``double`` float
1123		value uses 8 bytes.)
1124		``sparse``
1125		A sparse matrix. See :class:`SparseConnectionMatrix` for
1126		details on implementation. This class features very fast row
1127		access, and slower column access if the ``column_access=True``
1128		keyword is specified (making it suitable for learning
1129		algorithms such as STDP which require this). Memory
1130		requirements are 12 bytes per nonzero entry for row access
1131		only, or 20 bytes per nonzero entry if column access is
1132		specified. Synapses cannot be created or deleted at runtime
1133		with this class (although weights can be set to zero).
1134		``dynamic``
1135		A sparse matrix which allows runtime insertion and removal
1136		of synapses. See :class:`DynamicConnectionMatrix` for
1137		implementation details. This class features row and column
1138		access. The row access is slower than for ``sparse`` so this
1139		class should only be used when insertion and removal of
1140		synapses is crucial. Memory requirements are 24 bytes per
1141		nonzero entry. However, note that more memory than this
1142		may be required because memory is allocated using a
1143		dynamic array which grows by doubling its size when it runs
1144		out. If you know the maximum number of nonzero entries you will
1145		have in advance, specify the ``nnzmax`` keyword to set the
1146		initial size of the array.
1147
1148		Advanced information
1149
1150		The following methods are also defined and used internally, if you are
1151		writing your own derived connection class you need to understand what
1152		these do.
1153
1154		``propagate(spikes)``
1155		Action to take when source neurons with indices in ``spikes``
1156		fired.
1157		``do_propagate()``
1158		The method called by the :class:`Network` ``update()`` step,
1159		typically just propagates the spikes obtained by calling
1160		the ``get_spikes`` method of the ``source`` :class:`NeuronGroup`.
1161		'''
1162		#@check_units(delay=second)
1163		def __init__(self, source, target, state=0, delay=0 * msecond, modulation=None,
1164		structure='sparse', weight=None, sparseness=None, max_delay=5 * ms, **kwds):
1165		if not isinstance(delay, float):
1166		if delay is True:
1167		delay = None # this instructs us to use DelayConnection, but not initialise any delays
1168		self.__class__ = DelayConnection
1169		self.__init__(source, target, state=state, modulation=modulation, structure=structure,
1170		weight=weight, sparseness=sparseness, delay=delay, max_delay=max_delay, **kwds)
1171		return
1172		self.source = source # pointer to source group
1173		self.target = target # pointer to target group
1174		if isinstance(state, str): # named state variable
1175		self.nstate = target.get_var_index(state)
1176		else:
1177		self.nstate = state # target state index
1178		if isinstance(modulation, str): # named state variable
1179		self._nstate_mod = source.get_var_index(modulation)
1180		else:
1181		self._nstate_mod = modulation # source state index
1182		if isinstance(structure, str):
1183		structure = construction_matrix_register[structure]
1184		self.W = structure((len(source), len(target)), **kwds)
1185		self.iscompressed = False # True if compress() has been called
1186		source.set_max_delay(delay)
1187		if not isinstance(self, DelayConnection):
1188		self.delay = int(delay / source.clock.dt) # Synaptic delay in time bins
1189		self._useaccel = get_global_preference('useweave')
1190		self._cpp_compiler = get_global_preference('weavecompiler')
1191		self._extra_compile_args = ['-O3']
1192		if self._cpp_compiler == 'gcc':
1193		self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
1194		self._keyword_based_init(weight=weight, sparseness=sparseness)
1195
1196		def _keyword_based_init(self, weight=None, sparseness=None, **kwds):
1197		# Initialisation of weights
1198		# TODO: check consistency of weight and sparseness
1199		# TODO: select dense or sparse according to sparseness
1200		if weight is not None or sparseness is not None:
1201		if weight is None:
1202		weight = 1.0
1203		if sparseness is None:
1204		sparseness = 1
1205		if isinstance(weight, sparse.lil_matrix) or isinstance(weight, ndarray):
1206		self.connect(W=weight)
1207		elif sparseness == 1:
1208		self.connect_full(weight=weight)
1209		else:
1210		self.connect_random(weight=weight, p=sparseness)
1211
1212		def propagate(self, spikes):
1213		if not self.iscompressed:
1214		self.compress()
1215		if len(spikes):
1216		# Target state variable
1217		sv = self.target._S[self.nstate]
1218		# If specified, modulation state variable
1219		if self._nstate_mod is not None:
1220		sv_pre = self.source._S[self._nstate_mod]
1221		# Get the rows of the connection matrix, each row will be either a
1222		# DenseConnectionVector or a SparseConnectionVector.
1223		rows = self.W.get_rows(spikes)
1224		if not self._useaccel: # Pure Python version is easier to understand, but slower than C++ version below
1225		if isinstance(rows[0], SparseConnectionVector):
1226		if self._nstate_mod is None:
1227		# Rows stored as sparse vectors without modulation
1228		for row in rows:
1229		sv[row.ind] += row
	971	self.alldata[self.coldataind[j]] = val
	972
	973	def __setitem__(self, item, value):
	974	if item == colon_slice:
	975	self.alldata[:self.nnz] = value
	976	else:
	977	ConnectionMatrix.__setitem__(self, item, value)
	978
	979
	980	class Connection(magic.InstanceTracker, ObjectContainer):
	981	'''
	982	Mechanism for propagating spikes from one group to another
	983
	984	A Connection object declares that when spikes in a source
	985	group are generated, certain neurons in the target group
	986	should have a value added to specific states. See
	987	Tutorial 2: Connections to understand this better.
	988
	989	With arguments:
	990
	991	``source``
	992	The group from which spikes will be propagated.
	993	``target``
	994	The group to which spikes will be propagated.
	995	``state``
	996	The state variable name or number that spikes will be
	997	propagated to in the target group.
	998	``delay``
	999	The delay between a spike being generated at the source
	1000	and received at the target. Depending on the type of ``delay``
	1001	it has different effects. If ``delay`` is a scalar value, then
	1002	the connection will be initialised with all neurons having
	1003	that delay. For very long delays, this may raise an error. If
	1004	``delay=True`` then the connection will be initialised as a
	1005	:class:`DelayConnection`, allowing heterogeneous delays (a
	1006	different delay for each synapse). ``delay`` can also be a
	1007	pair ``(min,max)`` or a function of one or two variables, in
	1008	both cases it will be initialised as a :class:`DelayConnection`,
	1009	see the documentation for that class for details. Note that in
	1010	these cases, initialisation of delays will only have the
	1011	intended effect if used with the ``weight`` and ``sparseness``
	1012	arguments below.
	1013	``max_delay``
	1014	If you are using a connection with heterogeneous delays, specify
	1015	this to set the maximum allowed delay (smaller values use less
	1016	memory). The default is 5ms.
	1017	``modulation``
	1018	The state variable name from the source group that scales
	1019	the synaptic weights (for short-term synaptic plasticity).
	1020	``structure``
	1021	Data structure: ``sparse`` (default), ``dense`` or
	1022	``dynamic``. See below for more information on structures.
	1023	``weight``
	1024	If specified, the connection matrix will be initialised with
	1025	values specified by ``weight``, which can be any of the values
	1026	allowed in the methods `connect*`` below.
	1027	``sparseness``
	1028	If ``weight`` is specified and ``sparseness`` is not, a full
	1029	connection is assumed, otherwise random connectivity with this
	1030	level of sparseness is assumed.
	1031
	1032	Methods
	1033
	1034	``connect_random(P,Q,p[,weight=1[,fixed=False[,seed=None]]])``
	1035	Connects each neuron in ``P`` to each neuron in ``Q`` with independent
	1036	probability ``p`` and weight ``weight`` (this is the amount that
	1037	gets added to the target state variable). If ``fixed`` is True, then
	1038	the number of presynaptic neurons per neuron is constant. If ``seed``
	1039	is given, it is used as the seed to the random number generators, for
	1040	exactly repeatable results.
	1041	``connect_full(P,Q[,weight=1])``
	1042	Connect every neuron in ``P`` to every neuron in ``Q`` with the given
	1043	weight.
	1044	``connect_one_to_one(P,Q)``
	1045	If ``P`` and ``Q`` have the same number of neurons then neuron ``i``
	1046	in ``P`` will be connected to neuron ``i`` in ``Q`` with weight 1.
	1047	``connect(P,Q,W)``
	1048	You can specify a matrix of weights directly (can be in any format
	1049	recognised by NumPy). Note that due to internal implementation details,
	1050	passing a full matrix rather than a sparse one may slow down your code
	1051	(because zeros will be propagated as well as nonzero values).
	1052	WARNING: No unit checking is done at the moment.
	1053
	1054	Additionally, you can directly access the matrix of weights by writing::
	1055
	1056	C = Connection(P,Q)
	1057	print C[i,j]
	1058	C[i,j] = ...
	1059
	1060	Where here ``i`` is the source neuron and ``j`` is the target neuron.
	1061	Note: if ``C[i,j]`` should be zero, it is more efficient not to write
	1062	``C[i,j]=0``, if you write this then when neuron ``i`` fires all the
	1063	targets will have the value 0 added to them rather than just the
	1064	nonzero ones.
	1065	WARNING: No unit checking is currently done if you use this method.
	1066	Take care to set the right units.
	1067
	1068	Connection matrix structures
	1069
	1070	Brian currently features three types of connection matrix structures,
	1071	each of which is suited for different situations. Brian has two stages
	1072	of connection matrix. The first is the construction stage, used for
	1073	building a weight matrix. This stage is optimised for the construction
	1074	of matrices, with lots of features, but would be slow for runtime
	1075	behaviour. Consequently, the second stage is the connection stage,
	1076	used when Brian is being run. The connection stage is optimised for
	1077	run time behaviour, but many features which are useful for construction
	1078	are absent (e.g. the ability to add or remove synapses). Conversion
	1079	between construction and connection stages is done by the
	1080	``compress()`` method of :class:`Connection` which is called
	1081	automatically when it is used for the first time.
	1082
	1083	The structures are:
	1084
	1085	``dense``
	1086	A dense matrix. Allows runtime modification of all values. If
	1087	connectivity is close to being dense this is probably the most
	1088	efficient, but in most cases it is less efficient. In addition,
	1089	a dense connection matrix will often do the wrong thing if
	1090	using STDP. Because a synapse will be considered to exist but
	1091	with weight 0, STDP will be able to create new synapses where
	1092	there were previously none. Memory requirements are ``8NM``
	1093	bytes where ``(N,M)`` are the dimensions. (A ``double`` float
	1094	value uses 8 bytes.)
	1095	``sparse``
	1096	A sparse matrix. See :class:`SparseConnectionMatrix` for
	1097	details on implementation. This class features very fast row
	1098	access, and slower column access if the ``column_access=True``
	1099	keyword is specified (making it suitable for learning
	1100	algorithms such as STDP which require this). Memory
	1101	requirements are 12 bytes per nonzero entry for row access
	1102	only, or 20 bytes per nonzero entry if column access is
	1103	specified. Synapses cannot be created or deleted at runtime
	1104	with this class (although weights can be set to zero).
	1105	``dynamic``
	1106	A sparse matrix which allows runtime insertion and removal
	1107	of synapses. See :class:`DynamicConnectionMatrix` for
	1108	implementation details. This class features row and column
	1109	access. The row access is slower than for ``sparse`` so this
	1110	class should only be used when insertion and removal of
	1111	synapses is crucial. Memory requirements are 24 bytes per
	1112	nonzero entry. However, note that more memory than this
	1113	may be required because memory is allocated using a
	1114	dynamic array which grows by doubling its size when it runs
	1115	out. If you know the maximum number of nonzero entries you will
	1116	have in advance, specify the ``nnzmax`` keyword to set the
	1117	initial size of the array.
	1118
	1119	Advanced information
	1120
	1121	The following methods are also defined and used internally, if you are
	1122	writing your own derived connection class you need to understand what
	1123	these do.
	1124
	1125	``propagate(spikes)``
	1126	Action to take when source neurons with indices in ``spikes``
	1127	fired.
	1128	``do_propagate()``
	1129	The method called by the :class:`Network` ``update()`` step,
	1130	typically just propagates the spikes obtained by calling
	1131	the ``get_spikes`` method of the ``source`` :class:`NeuronGroup`.
	1132	'''
	1133	#@check_units(delay=second)
	1134	def __init__(self, source, target, state=0, delay=0 * msecond, modulation=None,
	1135	structure='sparse', weight=None, sparseness=None, max_delay=5 * ms, **kwds):
	1136	if not isinstance(delay, float):
	1137	if delay is True:
	1138	delay = None # this instructs us to use DelayConnection, but not initialise any delays
	1139	self.__class__ = DelayConnection
	1140	self.__init__(source, target, state=state, modulation=modulation, structure=structure,
	1141	weight=weight, sparseness=sparseness, delay=delay, max_delay=max_delay, **kwds)
	1142	return
	1143	self.source = source # pointer to source group
	1144	self.target = target # pointer to target group
	1145	if isinstance(state, str): # named state variable
	1146	self.nstate = target.get_var_index(state)
	1147	else:
	1148	self.nstate = state # target state index
	1149	if isinstance(modulation, str): # named state variable
	1150	self._nstate_mod = source.get_var_index(modulation)
	1151	else:
	1152	self._nstate_mod = modulation # source state index
	1153	if isinstance(structure, str):
	1154	structure = construction_matrix_register[structure]
	1155	self.W = structure((len(source), len(target)), **kwds)
	1156	self.iscompressed = False # True if compress() has been called
	1157	source.set_max_delay(delay)
	1158	if not isinstance(self, DelayConnection):
	1159	self.delay = int(delay / source.clock.dt) # Synaptic delay in time bins
	1160	self._useaccel = get_global_preference('useweave')
	1161	self._cpp_compiler = get_global_preference('weavecompiler')
	1162	self._extra_compile_args = ['-O3']
	1163	if self._cpp_compiler == 'gcc':
	1164	self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
	1165	self._keyword_based_init(weight=weight, sparseness=sparseness)
	1166
	1167	def _keyword_based_init(self, weight=None, sparseness=None, **kwds):
	1168	# Initialisation of weights
	1169	# TODO: check consistency of weight and sparseness
	1170	# TODO: select dense or sparse according to sparseness
	1171	if weight is not None or sparseness is not None:
	1172	if weight is None:
	1173	weight = 1.0
	1174	if sparseness is None:
	1175	sparseness = 1
	1176	if isinstance(weight, sparse.lil_matrix) or isinstance(weight, ndarray):
	1177	self.connect(W=weight)
	1178	elif sparseness == 1:
	1179	self.connect_full(weight=weight)
	1180	else:
	1181	self.connect_random(weight=weight, p=sparseness)
	1182
	1183	def propagate(self, spikes):
	1184	if not self.iscompressed:
	1185	self.compress()
	1186	if len(spikes):
	1187	# Target state variable
	1188	sv = self.target._S[self.nstate]
	1189	# If specified, modulation state variable
	1190	if self._nstate_mod is not None:
	1191	sv_pre = self.source._S[self._nstate_mod]
	1192	# Get the rows of the connection matrix, each row will be either a
	1193	# DenseConnectionVector or a SparseConnectionVector.
	1194	rows = self.W.get_rows(spikes)
	1195	if not self._useaccel: # Pure Python version is easier to understand, but slower than C++ version below
	1196	if isinstance(rows[0], SparseConnectionVector):
	1197	if self._nstate_mod is None:
	1198	# Rows stored as sparse vectors without modulation
	1199	for row in rows:
	1200	sv[row.ind] += row
	1201	else:
	1202	# Rows stored as sparse vectors with modulation
	1203	for i, row in izip(spikes, rows):
	1204	# note we call the numpy __mul__ directly because row is
	1205	# a SparseConnectionVector with different mul semantics
	1206	sv[row.ind] += numpy.ndarray.__mul__(row, sv_pre[i])
1230	1207	else:
1231		# Rows stored as sparse vectors with modulation
1232		for i, row in izip(spikes, rows):
1233		# note we call the numpy __mul__ directly because row is
1234		# a SparseConnectionVector with different mul semantics
1235		sv[row.ind] += numpy.ndarray.__mul__(row, sv_pre[i])
	1208	if self._nstate_mod is None:
	1209	# Rows stored as dense vectors without modulation
	1210	for row in rows:
	1211	sv += row
	1212	else:
	1213	# Rows stored as dense vectors with modulation
	1214	for i, row in izip(spikes, rows):
	1215	sv += numpy.ndarray.__mul__(row, sv_pre[i])
	1216	else: # C++ accelerated code, does the same as the code above but faster and less pretty
	1217	if isinstance(rows[0], SparseConnectionVector):
	1218	if self._nstate_mod is None:
	1219	nspikes = len(spikes)
	1220	rowinds = [r.ind for r in rows]
	1221	datas = rows
	1222	code = """
	1223	for(int j=0;j<nspikes;j++)
	1224	{
	1225	PyObject* _rowind = rowinds[j];
	1226	PyArrayObject* _row = convert_to_numpy(_rowind, "row");
	1227	conversion_numpy_check_type(_row, PyArray_LONG, "row");
	1228	conversion_numpy_check_size(_row, 1, "row");
	1229	//blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
	1230	long* row = (long*)_row->data;
	1231	PyObject* _datasj = datas[j];
	1232	PyArrayObject* _data = convert_to_numpy(_datasj, "data");
	1233	conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
	1234	conversion_numpy_check_size(_data, 1, "data");
	1235	//blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
	1236	double* data = (double*)_data->data;
	1237	//int m = row.numElements();
	1238	int m = _row->dimensions[0];
	1239	for(int k=0;k<m;k++)
	1240	//sv(row(k)) += data(k);
	1241	sv[row[k]] += data[k];
	1242	Py_DECREF(_rowind);
	1243	Py_DECREF(_datasj);
	1244	}
	1245	"""
	1246	weave.inline(code, ['sv', 'rowinds', 'datas', 'spikes', 'nspikes'],
	1247	compiler=self._cpp_compiler,
	1248	#type_converters=weave.converters.blitz,
	1249	extra_compile_args=self._extra_compile_args)
	1250	else:
	1251	nspikes = len(spikes)
	1252	rowinds = [r.ind for r in rows]
	1253	datas = rows
	1254	code = """
	1255	for(int j=0;j<nspikes;j++)
	1256	{
	1257	PyObject* _rowind = rowinds[j];
	1258	PyArrayObject* _row = convert_to_numpy(_rowind, "row");
	1259	conversion_numpy_check_type(_row, PyArray_LONG, "row");
	1260	conversion_numpy_check_size(_row, 1, "row");
	1261	//blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
	1262	long* row = (long*)_row->data;
	1263	PyObject* _datasj = datas[j];
	1264	PyArrayObject* _data = convert_to_numpy(_datasj, "data");
	1265	conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
	1266	conversion_numpy_check_size(_data, 1, "data");
	1267	//blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
	1268	double* data = (double*)_data->data;
	1269	//int m = row.numElements();
	1270	int m = _row->dimensions[0];
	1271	//double mod = sv_pre(spikes(j));
	1272	double mod = sv_pre[spikes[j]];
	1273	for(int k=0;k<m;k++)
	1274	//sv(row(k)) += data(k)*mod;
	1275	sv[row[k]] += data[k]*mod;
	1276	Py_DECREF(_rowind);
	1277	Py_DECREF(_datasj);
	1278	}
	1279	"""
	1280	weave.inline(code, ['sv', 'sv_pre', 'rowinds', 'datas', 'spikes', 'nspikes'],
	1281	compiler=self._cpp_compiler,
	1282	#type_converters=weave.converters.blitz,
	1283	extra_compile_args=self._extra_compile_args)
	1284	else:
	1285	if self._nstate_mod is None:
	1286	if not isinstance(spikes, numpy.ndarray):
	1287	spikes = array(spikes, dtype=int)
	1288	nspikes = len(spikes)
	1289	N = len(sv)
	1290	code = """
	1291	for(int j=0;j<nspikes;j++)
	1292	{
	1293	PyObject* _rowsj = rows[j];
	1294	PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
	1295	conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
	1296	conversion_numpy_check_size(_row, 1, "row");
	1297	//blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
	1298	double row = (double )_row->data;
	1299	for(int k=0;k<N;k++)
	1300	//sv(k) += row(k);
	1301	sv[k] += row[k];
	1302	Py_DECREF(_rowsj);
	1303	}
	1304	"""
	1305	weave.inline(code, ['sv', 'spikes', 'nspikes', 'N', 'rows'],
	1306	compiler=self._cpp_compiler,
	1307	#type_converters=weave.converters.blitz,
	1308	extra_compile_args=self._extra_compile_args)
	1309	else:
	1310	if not isinstance(spikes, numpy.ndarray):
	1311	spikes = array(spikes, dtype=int)
	1312	nspikes = len(spikes)
	1313	N = len(sv)
	1314	code = """
	1315	for(int j=0;j<nspikes;j++)
	1316	{
	1317	PyObject* _rowsj = rows[j];
	1318	PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
	1319	conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
	1320	conversion_numpy_check_size(_row, 1, "row");
	1321	//blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
	1322	//double mod = sv_pre(spikes(j));
	1323	double row = (double )_row->data;
	1324	double mod = sv_pre[spikes[j]];
	1325	for(int k=0;k<N;k++)
	1326	sv[k] += row[k]*mod;
	1327	//sv(k) += row(k)*mod;
	1328	Py_DECREF(_rowsj);
	1329	}
	1330	"""
	1331	weave.inline(code, ['sv', 'sv_pre', 'spikes', 'nspikes', 'N', 'rows'],
	1332	compiler=self._cpp_compiler,
	1333	#type_converters=weave.converters.blitz,
	1334	extra_compile_args=self._extra_compile_args)
	1335
	1336	def compress(self):
	1337	if not self.iscompressed:
	1338	self.W = self.W.connection_matrix()
	1339	self.iscompressed = True
	1340
	1341	def reinit(self):
	1342	'''
	1343	Resets the variables.
	1344	'''
	1345	pass
	1346
	1347	def do_propagate(self):
	1348	self.propagate(self.source.get_spikes(self.delay))
	1349
	1350	def origin(self, P, Q):
	1351	'''
	1352	Returns the starting coordinate of the given groups in
	1353	the connection matrix W.
	1354	'''
	1355	return (P._origin - self.source._origin, Q._origin - self.target._origin)
	1356
	1357	# TODO: rewrite all the connection functions to work row by row for memory and time efficiency
	1358
	1359	# TODO: change this
	1360	def connect(self, source=None, target=None, W=None):
	1361	'''
	1362	Connects (sub)groups P and Q with the weight matrix W (any type).
	1363	Internally: inserts W as a submatrix.
	1364	TODO: checks if the submatrix has already been specified.
	1365	'''
	1366	P = source or self.source
	1367	Q = target or self.target
	1368	i0, j0 = self.origin(P, Q)
	1369	self.W[i0:i0 + len(P), j0:j0 + len(Q)] = W
	1370
	1371	def connect_random(self, source=None, target=None, p=1., weight=1., fixed=False, seed=None, sparseness=None):
	1372	'''
	1373	Connects the neurons in group P to neurons in group Q with probability p,
	1374	with given weight (default 1).
	1375	The weight can be a quantity or a function of i (in P) and j (in Q).
	1376	If ``fixed`` is True, then the number of presynaptic neurons per neuron is constant.
	1377	'''
	1378	P = source or self.source
	1379	Q = target or self.target
	1380	if sparseness is not None: p = sparseness # synonym
	1381	if seed is not None:
	1382	numpy.random.seed(seed) # numpy's random number seed
	1383	pyrandom.seed(seed) # Python's random number seed
	1384	if fixed:
	1385	random_matrix_function = random_matrix_fixed_column
	1386	else:
	1387	random_matrix_function = random_matrix
	1388
	1389	if callable(weight):
	1390	# Check units
	1391	try:
	1392	if weight.func_code.co_argcount == 2:
	1393	weight(0, 0) + Q._S0[self.nstate]
	1394	else:
	1395	weight() + Q._S0[self.nstate]
	1396	except DimensionMismatchError, inst:
	1397	raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
	1398	self.connect(P, Q, random_matrix_function(len(P), len(Q), p, value=weight))
	1399	else:
	1400	# Check units
	1401	try:
	1402	weight + Q._S0[self.nstate]
	1403	except DimensionMismatchError, inst:
	1404	raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
	1405	self.connect(P, Q, random_matrix_function(len(P), len(Q), p, value=float(weight)))
	1406
	1407	def connect_full(self, source=None, target=None, weight=1.):
	1408	'''
	1409	Connects the neurons in group P to all neurons in group Q,
	1410	with given weight (default 1).
	1411	The weight can be a quantity or a function of i (in P) and j (in Q).
	1412	'''
	1413	P = source or self.source
	1414	Q = target or self.target
	1415	# TODO: check units
	1416	if callable(weight):
	1417	# Check units
	1418	try:
	1419	weight(0, 0) + Q._S0[self.nstate]
	1420	except DimensionMismatchError, inst:
	1421	raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
	1422	W = zeros((len(P), len(Q)))
	1423	try:
	1424	x = weight(0, 1. * arange(0, len(Q)))
	1425	failed = False
	1426	if array(x).size != len(Q):
	1427	failed = True
	1428	except:
	1429	failed = True
	1430	if failed: # vector-based not possible
	1431	log_debug('connections', 'Cannot build the connection matrix by rows')
	1432	for i in range(len(P)):
	1433	for j in range(len(Q)):
	1434	w = float(weight(i, j))
	1435	#if not is_within_absolute_tolerance(w,0.,effective_zero): # for sparse matrices
	1436	W[i, j] = w
1236	1437	else:
1237		if self._nstate_mod is None:
1238		# Rows stored as dense vectors without modulation
1239		for row in rows:
1240		sv += row
1241		else:
1242		# Rows stored as dense vectors with modulation
1243		for i, row in izip(spikes, rows):
1244		sv += numpy.ndarray.__mul__(row, sv_pre[i])
1245		else: # C++ accelerated code, does the same as the code above but faster and less pretty
1246		if isinstance(rows[0], SparseConnectionVector):
1247		if self._nstate_mod is None:
1248		nspikes = len(spikes)
1249		rowinds = [r.ind for r in rows]
1250		datas = rows
1251		code = """
1252		for(int j=0;j<nspikes;j++)
1253		{
1254		PyObject* _rowind = rowinds[j];
1255		PyArrayObject* _row = convert_to_numpy(_rowind, "row");
1256		conversion_numpy_check_type(_row, PyArray_LONG, "row");
1257		conversion_numpy_check_size(_row, 1, "row");
1258		//blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
1259		long* row = (long*)_row->data;
1260		PyObject* _datasj = datas[j];
1261		PyArrayObject* _data = convert_to_numpy(_datasj, "data");
1262		conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
1263		conversion_numpy_check_size(_data, 1, "data");
1264		//blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
1265		double* data = (double*)_data->data;
1266		//int m = row.numElements();
1267		int m = _row->dimensions[0];
1268		for(int k=0;k<m;k++)
1269		//sv(row(k)) += data(k);
1270		sv[row[k]] += data[k];
1271		Py_DECREF(_rowind);
1272		Py_DECREF(_datasj);
1273		}
1274		"""
1275		weave.inline(code, ['sv', 'rowinds', 'datas', 'spikes', 'nspikes'],
1276		compiler=self._cpp_compiler,
1277		#type_converters=weave.converters.blitz,
1278		extra_compile_args=self._extra_compile_args)
1279		else:
1280		nspikes = len(spikes)
1281		rowinds = [r.ind for r in rows]
1282		datas = rows
1283		code = """
1284		for(int j=0;j<nspikes;j++)
1285		{
1286		PyObject* _rowind = rowinds[j];
1287		PyArrayObject* _row = convert_to_numpy(_rowind, "row");
1288		conversion_numpy_check_type(_row, PyArray_LONG, "row");
1289		conversion_numpy_check_size(_row, 1, "row");
1290		//blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
1291		long* row = (long*)_row->data;
1292		PyObject* _datasj = datas[j];
1293		PyArrayObject* _data = convert_to_numpy(_datasj, "data");
1294		conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
1295		conversion_numpy_check_size(_data, 1, "data");
1296		//blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
1297		double* data = (double*)_data->data;
1298		//int m = row.numElements();
1299		int m = _row->dimensions[0];
1300		//double mod = sv_pre(spikes(j));
1301		double mod = sv_pre[spikes[j]];
1302		for(int k=0;k<m;k++)
1303		//sv(row(k)) += data(k)*mod;
1304		sv[row[k]] += data[k]*mod;
1305		Py_DECREF(_rowind);
1306		Py_DECREF(_datasj);
1307		}
1308		"""
1309		weave.inline(code, ['sv', 'sv_pre', 'rowinds', 'datas', 'spikes', 'nspikes'],
1310		compiler=self._cpp_compiler,
1311		#type_converters=weave.converters.blitz,
1312		extra_compile_args=self._extra_compile_args)
1313		else:
1314		if self._nstate_mod is None:
1315		if not isinstance(spikes, numpy.ndarray):
1316		spikes = array(spikes, dtype=int)
1317		nspikes = len(spikes)
1318		N = len(sv)
1319		code = """
1320		for(int j=0;j<nspikes;j++)
1321		{
1322		PyObject* _rowsj = rows[j];
1323		PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
1324		conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
1325		conversion_numpy_check_size(_row, 1, "row");
1326		//blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
1327		double row = (double )_row->data;
1328		for(int k=0;k<N;k++)
1329		//sv(k) += row(k);
1330		sv[k] += row[k];
1331		Py_DECREF(_rowsj);
1332		}
1333		"""
1334		weave.inline(code, ['sv', 'spikes', 'nspikes', 'N', 'rows'],
1335		compiler=self._cpp_compiler,
1336		#type_converters=weave.converters.blitz,
1337		extra_compile_args=self._extra_compile_args)
1338		else:
1339		if not isinstance(spikes, numpy.ndarray):
1340		spikes = array(spikes, dtype=int)
1341		nspikes = len(spikes)
1342		N = len(sv)
1343		code = """
1344		for(int j=0;j<nspikes;j++)
1345		{
1346		PyObject* _rowsj = rows[j];
1347		PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
1348		conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
1349		conversion_numpy_check_size(_row, 1, "row");
1350		//blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
1351		//double mod = sv_pre(spikes(j));
1352		double row = (double )_row->data;
1353		double mod = sv_pre[spikes[j]];
1354		for(int k=0;k<N;k++)
1355		sv[k] += row[k]*mod;
1356		//sv(k) += row(k)*mod;
1357		Py_DECREF(_rowsj);
1358		}
1359		"""
1360		weave.inline(code, ['sv', 'sv_pre', 'spikes', 'nspikes', 'N', 'rows'],
1361		compiler=self._cpp_compiler,
1362		#type_converters=weave.converters.blitz,
1363		extra_compile_args=self._extra_compile_args)
1364
1365		def compress(self):
1366		if not self.iscompressed:
1367		self.W = self.W.connection_matrix()
1368		self.iscompressed = True
1369
1370		def reinit(self):
1371		'''
1372		Resets the variables.
1373		'''
1374		pass
1375
1376		def do_propagate(self):
1377		self.propagate(self.source.get_spikes(self.delay))
1378
1379		def origin(self, P, Q):
1380		'''
1381		Returns the starting coordinate of the given groups in
1382		the connection matrix W.
1383		'''
1384		return (P._origin - self.source._origin, Q._origin - self.target._origin)
1385
1386		# TODO: rewrite all the connection functions to work row by row for memory and time efficiency
1387
1388		# TODO: change this
1389		def connect(self, source=None, target=None, W=None):
1390		'''
1391		Connects (sub)groups P and Q with the weight matrix W (any type).
1392		Internally: inserts W as a submatrix.
1393		TODO: checks if the submatrix has already been specified.
1394		'''
1395		P = source or self.source
1396		Q = target or self.target
1397		i0, j0 = self.origin(P, Q)
1398		self.W[i0:i0 + len(P), j0:j0 + len(Q)] = W
1399
1400		def connect_random(self, source=None, target=None, p=1., weight=1., fixed=False, seed=None, sparseness=None):
1401		'''
1402		Connects the neurons in group P to neurons in group Q with probability p,
1403		with given weight (default 1).
1404		The weight can be a quantity or a function of i (in P) and j (in Q).
1405		If ``fixed`` is True, then the number of presynaptic neurons per neuron is constant.
1406		'''
1407		P = source or self.source
1408		Q = target or self.target
1409		if sparseness is not None: p = sparseness # synonym
1410		if seed is not None:
1411		numpy.random.seed(seed) # numpy's random number seed
1412		pyrandom.seed(seed) # Python's random number seed
1413		if fixed:
1414		random_matrix_function = random_matrix_fixed_column
1415		else:
1416		random_matrix_function = random_matrix
1417
1418		if callable(weight):
1419		# Check units
1420		try:
1421		if weight.func_code.co_argcount == 2:
1422		weight(0, 0) + Q._S0[self.nstate]
1423		else:
1424		weight() + Q._S0[self.nstate]
1425		except DimensionMismatchError, inst:
1426		raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
1427		self.connect(P, Q, random_matrix_function(len(P), len(Q), p, value=weight))
1428		else:
1429		# Check units
1430		try:
1431		weight + Q._S0[self.nstate]
1432		except DimensionMismatchError, inst:
1433		raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
1434		self.connect(P, Q, random_matrix_function(len(P), len(Q), p, value=float(weight)))
1435
1436		def connect_full(self, source=None, target=None, weight=1.):
1437		'''
1438		Connects the neurons in group P to all neurons in group Q,
1439		with given weight (default 1).
1440		The weight can be a quantity or a function of i (in P) and j (in Q).
1441		'''
1442		P = source or self.source
1443		Q = target or self.target
1444		# TODO: check units
1445		if callable(weight):
1446		# Check units
1447		try:
1448		weight(0, 0) + Q._S0[self.nstate]
1449		except DimensionMismatchError, inst:
1450		raise DimensionMismatchError("Incorrects unit for the synaptic weights.", *inst._dims)
1451		W = zeros((len(P), len(Q)))
1452		try:
1453		x = weight(0, 1. * arange(0, len(Q)))
1454		failed = False
1455		if array(x).size != len(Q):
1456		failed = True
1457		except:
1458		failed = True
1459		if failed: # vector-based not possible
1460		log_debug('connections', 'Cannot build the connection matrix by rows')
1461		for i in range(len(P)):
1462		for j in range(len(Q)):
1463		w = float(weight(i, j))
1464		#if not is_within_absolute_tolerance(w,0.,effective_zero): # for sparse matrices
1465		W[i, j] = w
1466		else:
1467		for i in range(len(P)): # build W row by row
1468		#Below: for sparse matrices (?)
1469		#w = weight(i,1.*arange(0,len(Q)))
1470		#I = (abs(w)>effective_zero).nonzero()[0]
1471		#print w, I, w[I]
1472		#W[i,I] = w[I]
1473		W[i, :] = weight(i, 1. * arange(0, len(Q)))
1474		self.connect(P, Q, W)
1475		else:
1476		try:
1477		weight + Q._S0[self.nstate]
1478		except DimensionMismatchError, inst:
1479		raise DimensionMismatchError("Incorrect unit for the synaptic weights.", *inst._dims)
1480		self.connect(P, Q, float(weight) * ones((len(P), len(Q))))
1481
1482		def connect_one_to_one(self, source=None, target=None, weight=1):
1483		'''
1484		Connects source[i] to target[i] with weights 1 (or weight).
1485		'''
1486		P = source or self.source
1487		Q = target or self.target
1488		if (len(P) != len(Q)):
1489		raise AttributeError, 'The connected (sub)groups must have the same size.'
1490		# TODO: unit checking
1491		self.connect(P, Q, float(weight) * eye_lil_matrix(len(P)))
1492
1493		def __getitem__(self, i):
1494		return self.W.__getitem__(i)
1495
1496		def __setitem__(self, i, x):
1497		# TODO: unit checking
1498		self.W.__setitem__(i, x)
1499
1500
1501		class DelayConnection(Connection):
1502		'''
1503		Connection which implements heterogeneous postsynaptic delays
1504
1505		Initialised as for a :class:`Connection`, but with the additional
1506		keyword:
1507
1508		``max_delay``
1509		Specifies the maximum delay time for any
1510		neuron. Note, the smaller you make this the less memory will be
1511		used.
1512
1513		Overrides the following attribute of :class:`Connection`:
1514
1515		.. attribute:: delay
1516
1517		A matrix of delays. This array can be changed during a run,
1518		but at no point should it be greater than ``max_delay``.
1519
1520		In addition, the methods ``connect``, ``connect_random``, ``connect_full``,
1521		and ``connect_one_to_one`` have a new keyword ``delay=...`` for setting the
1522		initial values of the delays, where ``delay`` can be one of:
1523
1524		* A float, all delays will be set to this value
1525		* A pair (min, max), delays will be uniform between these two
1526		values.
1527		* A function of no arguments, will be called for each nonzero
1528		entry in the weight matrix.
1529		* A function of two argument ``(i,j)`` will be called for each
1530		nonzero entry in the weight matrix.
1531		* A matrix of an appropriate type (e.g. ndarray or lil_matrix).
1532
1533		Finally, there is a method:
1534
1535		``set_delays(source, target, delay)``
1536		Where ``delay`` must be of one of the types above.
1537
1538		Notes
1539
1540		This class implements post-synaptic delays. This means that the spike is
1541		propagated immediately from the presynaptic neuron with the synaptic
1542		weight at the time of the spike, but arrives at the postsynaptic neuron
1543		with the given delay. At the moment, Brian only provides support for
1544		presynaptic delays if they are homogeneous, using the ``delay`` keyword
1545		of a standard ``Connection``.
1546
1547		Implementation
1548
1549		:class:`DelayConnection` stores an array of size ``(n,m)`` where
1550		``n`` is ``max_delay/dt`` for ``dt`` of the target :class:`NeuronGroup`'s clock,
1551		and ``m`` is the number of neurons in the target. This array can potentially
1552		be quite large. Each row in this array represents the array that should be
1553		added to the target state variable at some particular future time. Which
1554		row corresponds to which time is tracked using a circular indexing scheme.
1555
1556		When a spike from neuron ``i`` in the source is encountered, the delay time
1557		of neuron ``i`` is looked up, the row corresponding to the current time
1558		plus that delay time is found using the circular indexing scheme, and then
1559		the spike is propagated to that row as for a standard connection (although
1560		this won't be propagated to the target until a later time).
1561
1562		Warning
1563
1564		If you are using a dynamic connection matrix, it is your responsibility to
1565		ensure that the nonzero entries of the weight matrix and the delay matrix
1566		exactly coincide. This is not an issue for sparse or dense matrices.
1567		'''
1568
1569		def __init__(self, source, target, state=0, modulation=None,
1570		structure='sparse',
1571		weight=None, sparseness=None, delay=None,
1572		max_delay=5 * msecond, **kwds):
1573		Connection.__init__(self, source, target, state=state, modulation=modulation,
1574		structure=structure, weight=weight, sparseness=sparseness, **kwds)
1575		self._max_delay = int(max_delay / target.clock.dt) + 1
1576		source.set_max_delay(max_delay)
1577		# Each row of the following array stores the cumulative effect of spikes at some
1578		# particular time, defined by a circular indexing scheme. The _cur_delay_ind attribute
1579		# stores the row corresponding to the current time, so that _cur_delay_ind+1 corresponds
1580		# to that time + target.clock.dt, and so on. When _cur_delay_ind reaches _max_delay it
1581		# resets to zero.
1582		self._delayedreaction = numpy.zeros((self._max_delay, len(target)))
1583		# vector of delay times, can be changed during a run
1584		if isinstance(structure, str):
1585		structure = construction_matrix_register[structure]
1586		self.delayvec = structure((len(source), len(target)), **kwds)
1587		self._cur_delay_ind = 0
1588		# this network operation is added to the Network object via the contained_objects
1589		# protocol (see the line after the function definition). The standard Connection.propagate
1590		# function propagates spikes to _delayedreaction rather than the target, and this
1591		# function which is called after the usual propagations propagates that data from
1592		# _delayedreaction to the target. It only needs to be called each target.clock update.
1593		@network_operation(clock=target.clock, when='after_connections')
1594		def delayed_propagate():
1595		# propagate from _delayedreaction -> target group
1596		target._S[self.nstate] += self._delayedreaction[self._cur_delay_ind, :]
1597		# reset the current row of _delayedreaction
1598		self._delayedreaction[self._cur_delay_ind, :] = 0.0
1599		# increase the index for the circular indexing scheme
1600		self._cur_delay_ind = (self._cur_delay_ind + 1) % self._max_delay
1601		self.contained_objects = [delayed_propagate]
1602		# this is just used to convert delayvec's which are in ms to integers, precalculating it makes it faster
1603		self._invtargetdt = 1 / self.target.clock._dt
1604		self._useaccel = get_global_preference('useweave')
1605		self._cpp_compiler = get_global_preference('weavecompiler')
1606		self._extra_compile_args = ['-O3']
1607		if self._cpp_compiler == 'gcc':
1608		self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
1609		if delay is not None:
1610		self.set_delays(delay=delay)
1611
1612		def propagate(self, spikes):
1613		if not self.iscompressed:
1614		self.compress()
1615		if len(spikes):
1616		# Target state variable
1617		dr = self._delayedreaction
1618		# If specified, modulation state variable
1619		if self._nstate_mod is not None:
1620		sv_pre = self.source._S[self._nstate_mod]
1621		# Get the rows of the connection matrix, each row will be either a
1622		# DenseConnectionVector or a SparseConnectionVector.
1623		rows = self.W.get_rows(spikes)
1624		dvecrows = self.delayvec.get_rows(spikes)
1625		if not self._useaccel: # Pure Python version is easier to understand, but slower than C++ version below
1626		if isinstance(rows[0], SparseConnectionVector):
1627		if self._nstate_mod is None:
1628		# Rows stored as sparse vectors without modulation
1629		for row, dvecrow in izip(rows, dvecrows):
1630		if not len(row.ind) == len(dvecrow.ind):
1631		raise RuntimeError('Weight and delay matrices must be kept in synchrony for sparse matrices.')
1632		drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
1633		dr[drind, dvecrow.ind] += row
1634		else:
1635		# Rows stored as sparse vectors with modulation
1636		for i, row, dvecrow in izip(spikes, rows, dvecrows):
1637		if not len(row.ind) == len(dvecrow.ind):
1638		raise RuntimeError('Weight and delay matrices must be kept in synchrony for sparse matrices.')
1639		drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
1640		# note we call the numpy __mul__ directly because row is
1641		# a SparseConnectionVector with different mul semantics
1642		dr[drind, dvecrow.ind] += numpy.ndarray.__mul__(row, sv_pre[i])
1643		else:
1644		if self._nstate_mod is None:
1645		# Rows stored as dense vectors without modulation
1646		drjind = numpy.arange(len(self.target), dtype=int)
1647		for row, dvecrow in izip(rows, dvecrows):
1648		drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
1649		dr[drind, drjind[:len(drind)]] += row
1650		else:
1651		# Rows stored as dense vectors with modulation
1652		drjind = numpy.arange(len(self.target), dtype=int)
1653		for i, row, dvecrow in izip(spikes, rows, dvecrows):
1654		drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
1655		dr[drind, drjind[:len(drind)]] += numpy.ndarray.__mul__(row, sv_pre[i])
1656		else: # C++ accelerated code, does the same as the code above but faster and less pretty
1657		if isinstance(rows[0], SparseConnectionVector):
1658		if self._nstate_mod is None:
1659		nspikes = len(spikes)
1660		rowinds = [r.ind for r in rows]
1661		datas = rows
1662		cdi = self._cur_delay_ind
1663		idt = self._invtargetdt
1664		md = self._max_delay
1665		code = """
1666		for(int j=0;j<nspikes;j++)
1667		{
1668		PyObject* _rowind = rowinds[j];
1669		PyArrayObject* _row = convert_to_numpy(_rowind, "row");
1670		conversion_numpy_check_type(_row, PyArray_LONG, "row");
1671		conversion_numpy_check_size(_row, 1, "row");
1672		blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
1673		PyObject* _datasj = datas[j];
1674		PyArrayObject* _data = convert_to_numpy(_datasj, "data");
1675		conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
1676		conversion_numpy_check_size(_data, 1, "data");
1677		blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
1678		PyObject* _dvecrowsj = dvecrows[j];
1679		PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
1680		conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
1681		conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
1682		blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
1683		int m = row.numElements();
1684		for(int k=0;k<m;k++)
1685		{
1686		dr((cdi+(int)(idt*dvecrow(k)))%md, (int)row(k)) += data(k);
1687		}
1688		Py_DECREF(_rowind);
1689		Py_DECREF(_datasj);
1690		Py_DECREF(_dvecrowsj);
1691		}
1692		"""
1693		weave.inline(code, ['rowinds', 'datas', 'spikes', 'nspikes',
1694		'dvecrows', 'dr', 'cdi', 'idt', 'md'],
1695		compiler=self._cpp_compiler,
1696		type_converters=weave.converters.blitz,
1697		extra_compile_args=self._extra_compile_args)
1698		else:
1699		nspikes = len(spikes)
1700		rowinds = [r.ind for r in rows]
1701		datas = rows
1702		cdi = self._cur_delay_ind
1703		idt = self._invtargetdt
1704		md = self._max_delay
1705		code = """
1706		for(int j=0;j<nspikes;j++)
1707		{
1708		PyObject* _rowind = rowinds[j];
1709		PyArrayObject* _row = convert_to_numpy(_rowind, "row");
1710		conversion_numpy_check_type(_row, PyArray_LONG, "row");
1711		conversion_numpy_check_size(_row, 1, "row");
1712		blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
1713		PyObject* _datasj = datas[j];
1714		PyArrayObject* _data = convert_to_numpy(_datasj, "data");
1715		conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
1716		conversion_numpy_check_size(_data, 1, "data");
1717		blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
1718		PyObject* _dvecrowsj = dvecrows[j];
1719		PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
1720		conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
1721		conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
1722		blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
1723		int m = row.numElements();
1724		double mod = sv_pre(spikes(j));
1725		for(int k=0;k<m;k++)
1726		{
1727		dr((cdi+(int)(idtdvecrow(k)))%md, (int)row(k)) += data(k)mod;
1728		}
1729		Py_DECREF(_rowind);
1730		Py_DECREF(_datasj);
1731		Py_DECREF(_dvecrowsj);
1732		}
1733		"""
1734		weave.inline(code, ['sv_pre', 'rowinds', 'datas', 'spikes', 'nspikes',
1735		'dvecrows', 'dr', 'cdi', 'idt', 'md'],
1736		compiler=self._cpp_compiler,
1737		type_converters=weave.converters.blitz,
1738		extra_compile_args=self._extra_compile_args)
1739		else:
1740		if self._nstate_mod is None:
1741		if not isinstance(spikes, numpy.ndarray):
1742		spikes = array(spikes, dtype=int)
1743		nspikes = len(spikes)
1744		N = len(self.target)
1745		cdi = self._cur_delay_ind
1746		idt = self._invtargetdt
1747		md = self._max_delay
1748		code = """
1749		for(int j=0;j<nspikes;j++)
1750		{
1751		PyObject* _rowsj = rows[j];
1752		PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
1753		conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
1754		conversion_numpy_check_size(_row, 1, "row");
1755		blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
1756		PyObject* _dvecrowsj = dvecrows[j];
1757		PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
1758		conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
1759		conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
1760		blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
1761		for(int k=0;k<N;k++)
1762		dr((cdi+(int)(idt*dvecrow(k)))%md, k) += row(k);
1763		Py_DECREF(_rowsj);
1764		Py_DECREF(_dvecrowsj);
1765		}
1766		"""
1767		weave.inline(code, ['spikes', 'nspikes', 'N', 'rows',
1768		'dr', 'cdi', 'idt', 'md', 'dvecrows'],
1769		compiler=self._cpp_compiler,
1770		type_converters=weave.converters.blitz,
1771		extra_compile_args=self._extra_compile_args)
1772		else:
1773		if not isinstance(spikes, numpy.ndarray):
1774		spikes = array(spikes, dtype=int)
1775		nspikes = len(spikes)
1776		N = len(self.target)
1777		cdi = self._cur_delay_ind
1778		idt = self._invtargetdt
1779		md = self._max_delay
1780		code = """
1781		for(int j=0;j<nspikes;j++)
1782		{
1783		PyObject* _rowsj = rows[j];
1784		PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
1785		conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
1786		conversion_numpy_check_size(_row, 1, "row");
1787		blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
1788		PyObject* _dvecrowsj = dvecrows[j];
1789		PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
1790		conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
1791		conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
1792		blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
1793		double mod = sv_pre(spikes(j));
1794		for(int k=0;k<N;k++)
1795		dr((cdi+(int)(idtdvecrow(k)))%md, k) += row(k)mod;
1796		Py_DECREF(_rowsj);
1797		Py_DECREF(_dvecrowsj);
1798		}
1799		"""
1800		weave.inline(code, ['sv_pre', 'spikes', 'nspikes', 'N', 'rows',
1801		'dr', 'cdi', 'idt', 'md', 'dvecrows'],
1802		compiler=self._cpp_compiler,
1803		type_converters=weave.converters.blitz,
1804		extra_compile_args=self._extra_compile_args)
1805
1806		def do_propagate(self):
1807		self.propagate(self.source.get_spikes(0))
1808
1809		def _set_delay_property(self, val):
1810		self.delayvec[:] = val
1811
1812		delay = property(fget=lambda self:self.delayvec, fset=_set_delay_property)
1813
1814		def compress(self):
1815		if not self.iscompressed:
1816		# We want delayvec to have nonzero entries at the same places as
1817		# W does, so we use W to initialise the compressed version of
1818		# delayvec, and then copy the values from the old delayvec to
1819		# the new compressed one, allowing delayvec and W to not have
1820		# to be perfectly intersected at the initialisation stage. If
1821		# the structure is dynamic, it will be the user's
1822		# responsibility to update them in sequence
1823		delayvec = self.delayvec
1824		# The delayvec[i,:] operation for sparse.lil_matrix format
1825		# is VERY slow, but the CSR format is fine.
1826		if isinstance(delayvec, sparse.lil_matrix):
1827		delayvec = delayvec.tocsr()
1828		self.delayvec = self.W.connection_matrix(copy=True)
1829		for i in xrange(self.W.shape[0]):
1830		self.delayvec.set_row(i, array(todense(delayvec[i, :]), copy=False).flatten())
1831		Connection.compress(self)
1832
1833		def set_delays(self, source=None, target=None, delay=None):
1834		'''
1835		Set the delays corresponding to the weight matrix
1836
1837		``delay`` must be one of:
1838
	1438	for i in range(len(P)): # build W row by row
	1439	#Below: for sparse matrices (?)
	1440	#w = weight(i,1.*arange(0,len(Q)))
	1441	#I = (abs(w)>effective_zero).nonzero()[0]
	1442	#print w, I, w[I]
	1443	#W[i,I] = w[I]
	1444	W[i, :] = weight(i, 1. * arange(0, len(Q)))
	1445	self.connect(P, Q, W)
	1446	else:
	1447	try:
	1448	weight + Q._S0[self.nstate]
	1449	except DimensionMismatchError, inst:
	1450	raise DimensionMismatchError("Incorrect unit for the synaptic weights.", *inst._dims)
	1451	self.connect(P, Q, float(weight) * ones((len(P), len(Q))))
	1452
	1453	def connect_one_to_one(self, source=None, target=None, weight=1):
	1454	'''
	1455	Connects source[i] to target[i] with weights 1 (or weight).
	1456	'''
	1457	P = source or self.source
	1458	Q = target or self.target
	1459	if (len(P) != len(Q)):
	1460	raise AttributeError, 'The connected (sub)groups must have the same size.'
	1461	# TODO: unit checking
	1462	self.connect(P, Q, float(weight) * eye_lil_matrix(len(P)))
	1463
	1464	def __getitem__(self, i):
	1465	return self.W.__getitem__(i)
	1466
	1467	def __setitem__(self, i, x):
	1468	# TODO: unit checking
	1469	self.W.__setitem__(i, x)
	1470
	1471
	1472	class DelayConnection(Connection):
	1473	'''
	1474	Connection which implements heterogeneous postsynaptic delays
	1475
	1476	Initialised as for a :class:`Connection`, but with the additional
	1477	keyword:
	1478
	1479	``max_delay``
	1480	Specifies the maximum delay time for any
	1481	neuron. Note, the smaller you make this the less memory will be
	1482	used.
	1483
	1484	Overrides the following attribute of :class:`Connection`:
	1485
	1486	.. attribute:: delay
	1487
	1488	A matrix of delays. This array can be changed during a run,
	1489	but at no point should it be greater than ``max_delay``.
	1490
	1491	In addition, the methods ``connect``, ``connect_random``, ``connect_full``,
	1492	and ``connect_one_to_one`` have a new keyword ``delay=...`` for setting the
	1493	initial values of the delays, where ``delay`` can be one of:
	1494
1839	1495	* A float, all delays will be set to this value
1840	1496	* A pair (min, max), delays will be uniform between these two
…	…
1845	1501	nonzero entry in the weight matrix.
1846	1502	* A matrix of an appropriate type (e.g. ndarray or lil_matrix).
1847		'''
1848		if delay is None:
1849		return
1850		W = self.W
1851		P = source or self.source
1852		Q = target or self.target
1853		i0, j0 = self.origin(P, Q)
1854		i1 = i0 + len(P)
1855		j1 = j0 + len(Q)
1856		if isinstance(W, sparse.lil_matrix):
1857		def getrow(i):
1858		inds = array(W.rows[i], dtype=int)
1859		inds = inds[logical_and(inds >= j0, inds < j1)]
1860		return inds, len(inds)
1861		else:
1862		def getrow(i):
1863		inds = (W[i, j0:j1] != 0).nonzero()[0] + j0
1864		return inds, len(inds)
1865		#return slice(j0, j1), j1-j0
1866		if isinstance(delay, (float, int)):
1867		for i in xrange(i0, i1):
1868		inds, L = getrow(i)
1869		self.delayvec[i, inds] = delay
1870		elif isinstance(delay, (tuple, list)) and len(delay) == 2:
1871		delaymin, delaymax = delay
1872		for i in xrange(i0, i1):
1873		inds, L = getrow(i)
1874		rowdelay = rand(L) * (delaymax - delaymin) + delaymin
1875		self.delayvec[i, inds] = rowdelay
1876		elif callable(delay) and delay.func_code.co_argcount == 0:
1877		for i in xrange(i0, i1):
1878		inds, L = getrow(i)
1879		rowdelay = [delay() for _ in xrange(L)]
1880		self.delayvec[i, inds] = rowdelay
1881		elif callable(delay) and delay.func_code.co_argcount == 2:
1882		for i in xrange(i0, i1):
1883		inds, L = getrow(i)
1884		if isinstance(inds, slice):
1885		inds = numpy.arange(inds.start, inds.stop)
1886		self.delayvec[i, inds] = delay(i - i0, inds - j0)
1887		else:
1888		#raise TypeError('delays must be float, pair or function of 0 or 2 arguments')
1889		self.delayvec[i0:i1, j0:j1] = delay # probably won't work, but then it will raise an error
1890
1891		def connect(self, source=None, target=None, W=None, delay=None):
1892		Connection.connect(self, source=source, target=target, W=W)
1893		if delay is not None:
1894		self.set_delays(source, target, delay)
1895
1896		def connect_random(self, source=None, target=None, p=1.0, weight=1.0,
1897		fixed=False, seed=None, sparseness=None, delay=None):
1898		Connection.connect_random(self, source=source, target=target, p=p,
1899		weight=weight, fixed=fixed, seed=seed,
1900		sparseness=sparseness)
1901		if delay is not None:
1902		self.set_delays(source, target, delay)
1903
1904		def connect_full(self, source=None, target=None, weight=1.0, delay=None):
1905		Connection.connect_full(self, source=source, target=target, weight=weight)
1906		if delay is not None:
1907		self.set_delays(source, target, delay)
1908
1909		def connect_one_to_one(self, source=None, target=None, weight=1.0, delay=None):
1910		Connection.connect_one_to_one(self, source=source, target=target,
1911		weight=weight)
1912		if delay is not None:
1913		self.set_delays(source, target, delay)
1914
1915
1916		class IdentityConnection(Connection):
1917		'''
1918		A :class:`Connection` between two groups of the same size, where neuron ``i`` in the
1919		source group is connected to neuron ``i`` in the target group.
1920
1921		Initialised with arguments:
1922
1923		``source``, ``target``
1924		The source and target :class:`NeuronGroup` objects.
1925		``state``
1926		The target state variable.
1927		``weight``
1928		The weight of the synapse, must be a scalar.
1929		``delay``
1930		Only homogeneous delays are allowed.
1931
1932		The benefit of this class is that it has no storage requirements and is optimised for
1933		this special case.
1934		'''
1935		@check_units(delay=second)
1936		def __init__(self, source, target, state=0, weight=1, delay=0 * msecond):
1937		if (len(source) != len(target)):
1938		raise AttributeError, 'The connected (sub)groups must have the same size.'
1939		self.source = source # pointer to source group
1940		self.target = target # pointer to target group
1941		if type(state) == types.StringType: # named state variable
1942		self.nstate = target.get_var_index(state)
1943		else:
1944		self.nstate = state # target state index
1945		self.W = float(weight) # weight
1946		source.set_max_delay(delay)
1947		self.delay = int(delay / source.clock.dt) # Synaptic delay in time bins
1948
1949		def propagate(self, spikes):
1950		'''
1951		Propagates the spikes to the target.
1952		'''
1953		self.target._S[self.nstate, spikes] += self.W
1954
1955		def compress(self):
1956		pass
1957
1958
1959		class MultiConnection(Connection):
1960		'''
1961		A hub for multiple connections with a common source group.
1962		'''
1963		def __init__(self, source, connections=[]):
1964		self.source = source
1965		self.connections = connections
1966		self.iscompressed = False
1967		self.delay = connections[0].delay
1968
1969		def propagate(self, spikes):
1970		'''
1971		Propagates the spikes to the targets.
1972		'''
1973		for C in self.connections:
1974		C.propagate(spikes)
1975
1976		def compress(self):
1977		if not self.iscompressed:
	1503
	1504	Finally, there is a method:
	1505
	1506	``set_delays(source, target, delay)``
	1507	Where ``delay`` must be of one of the types above.
	1508
	1509	Notes
	1510
	1511	This class implements post-synaptic delays. This means that the spike is
	1512	propagated immediately from the presynaptic neuron with the synaptic
	1513	weight at the time of the spike, but arrives at the postsynaptic neuron
	1514	with the given delay. At the moment, Brian only provides support for
	1515	presynaptic delays if they are homogeneous, using the ``delay`` keyword
	1516	of a standard ``Connection``.
	1517
	1518	Implementation
	1519
	1520	:class:`DelayConnection` stores an array of size ``(n,m)`` where
	1521	``n`` is ``max_delay/dt`` for ``dt`` of the target :class:`NeuronGroup`'s clock,
	1522	and ``m`` is the number of neurons in the target. This array can potentially
	1523	be quite large. Each row in this array represents the array that should be
	1524	added to the target state variable at some particular future time. Which
	1525	row corresponds to which time is tracked using a circular indexing scheme.
	1526
	1527	When a spike from neuron ``i`` in the source is encountered, the delay time
	1528	of neuron ``i`` is looked up, the row corresponding to the current time
	1529	plus that delay time is found using the circular indexing scheme, and then
	1530	the spike is propagated to that row as for a standard connection (although
	1531	this won't be propagated to the target until a later time).
	1532
	1533	Warning
	1534
	1535	If you are using a dynamic connection matrix, it is your responsibility to
	1536	ensure that the nonzero entries of the weight matrix and the delay matrix
	1537	exactly coincide. This is not an issue for sparse or dense matrices.
	1538	'''
	1539
	1540	def __init__(self, source, target, state=0, modulation=None,
	1541	structure='sparse',
	1542	weight=None, sparseness=None, delay=None,
	1543	max_delay=5 * msecond, **kwds):
	1544	Connection.__init__(self, source, target, state=state, modulation=modulation,
	1545	structure=structure, weight=weight, sparseness=sparseness, **kwds)
	1546	self._max_delay = int(max_delay / target.clock.dt) + 1
	1547	source.set_max_delay(max_delay)
	1548	# Each row of the following array stores the cumulative effect of spikes at some
	1549	# particular time, defined by a circular indexing scheme. The _cur_delay_ind attribute
	1550	# stores the row corresponding to the current time, so that _cur_delay_ind+1 corresponds
	1551	# to that time + target.clock.dt, and so on. When _cur_delay_ind reaches _max_delay it
	1552	# resets to zero.
	1553	self._delayedreaction = numpy.zeros((self._max_delay, len(target)))
	1554	# vector of delay times, can be changed during a run
	1555	if isinstance(structure, str):
	1556	structure = construction_matrix_register[structure]
	1557	self.delayvec = structure((len(source), len(target)), **kwds)
	1558	self._cur_delay_ind = 0
	1559	# this network operation is added to the Network object via the contained_objects
	1560	# protocol (see the line after the function definition). The standard Connection.propagate
	1561	# function propagates spikes to _delayedreaction rather than the target, and this
	1562	# function which is called after the usual propagations propagates that data from
	1563	# _delayedreaction to the target. It only needs to be called each target.clock update.
	1564	@network_operation(clock=target.clock, when='after_connections')
	1565	def delayed_propagate():
	1566	# propagate from _delayedreaction -> target group
	1567	target._S[self.nstate] += self._delayedreaction[self._cur_delay_ind, :]
	1568	# reset the current row of _delayedreaction
	1569	self._delayedreaction[self._cur_delay_ind, :] = 0.0
	1570	# increase the index for the circular indexing scheme
	1571	self._cur_delay_ind = (self._cur_delay_ind + 1) % self._max_delay
	1572	self.contained_objects = [delayed_propagate]
	1573	# this is just used to convert delayvec's which are in ms to integers, precalculating it makes it faster
	1574	self._invtargetdt = 1 / self.target.clock._dt
	1575	self._useaccel = get_global_preference('useweave')
	1576	self._cpp_compiler = get_global_preference('weavecompiler')
	1577	self._extra_compile_args = ['-O3']
	1578	if self._cpp_compiler == 'gcc':
	1579	self._extra_compile_args += get_global_preference('gcc_options') # ['-march=native', '-ffast-math']
	1580	if delay is not None:
	1581	self.set_delays(delay=delay)
	1582
	1583	def propagate(self, spikes):
	1584	if not self.iscompressed:
	1585	self.compress()
	1586	if len(spikes):
	1587	# Target state variable
	1588	dr = self._delayedreaction
	1589	# If specified, modulation state variable
	1590	if self._nstate_mod is not None:
	1591	sv_pre = self.source._S[self._nstate_mod]
	1592	# Get the rows of the connection matrix, each row will be either a
	1593	# DenseConnectionVector or a SparseConnectionVector.
	1594	rows = self.W.get_rows(spikes)
	1595	dvecrows = self.delayvec.get_rows(spikes)
	1596	if not self._useaccel: # Pure Python version is easier to understand, but slower than C++ version below
	1597	if isinstance(rows[0], SparseConnectionVector):
	1598	if self._nstate_mod is None:
	1599	# Rows stored as sparse vectors without modulation
	1600	for row, dvecrow in izip(rows, dvecrows):
	1601	if not len(row.ind) == len(dvecrow.ind):
	1602	raise RuntimeError('Weight and delay matrices must be kept in synchrony for sparse matrices.')
	1603	drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
	1604	dr[drind, dvecrow.ind] += row
	1605	else:
	1606	# Rows stored as sparse vectors with modulation
	1607	for i, row, dvecrow in izip(spikes, rows, dvecrows):
	1608	if not len(row.ind) == len(dvecrow.ind):
	1609	raise RuntimeError('Weight and delay matrices must be kept in synchrony for sparse matrices.')
	1610	drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
	1611	# note we call the numpy __mul__ directly because row is
	1612	# a SparseConnectionVector with different mul semantics
	1613	dr[drind, dvecrow.ind] += numpy.ndarray.__mul__(row, sv_pre[i])
	1614	else:
	1615	if self._nstate_mod is None:
	1616	# Rows stored as dense vectors without modulation
	1617	drjind = numpy.arange(len(self.target), dtype=int)
	1618	for row, dvecrow in izip(rows, dvecrows):
	1619	drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
	1620	dr[drind, drjind[:len(drind)]] += row
	1621	else:
	1622	# Rows stored as dense vectors with modulation
	1623	drjind = numpy.arange(len(self.target), dtype=int)
	1624	for i, row, dvecrow in izip(spikes, rows, dvecrows):
	1625	drind = (self._cur_delay_ind + numpy.array(self._invtargetdt * dvecrow, dtype=int)) % self._max_delay
	1626	dr[drind, drjind[:len(drind)]] += numpy.ndarray.__mul__(row, sv_pre[i])
	1627	else: # C++ accelerated code, does the same as the code above but faster and less pretty
	1628	if isinstance(rows[0], SparseConnectionVector):
	1629	if self._nstate_mod is None:
	1630	nspikes = len(spikes)
	1631	rowinds = [r.ind for r in rows]
	1632	datas = rows
	1633	cdi = self._cur_delay_ind
	1634	idt = self._invtargetdt
	1635	md = self._max_delay
	1636	code = """
	1637	for(int j=0;j<nspikes;j++)
	1638	{
	1639	PyObject* _rowind = rowinds[j];
	1640	PyArrayObject* _row = convert_to_numpy(_rowind, "row");
	1641	conversion_numpy_check_type(_row, PyArray_LONG, "row");
	1642	conversion_numpy_check_size(_row, 1, "row");
	1643	blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
	1644	PyObject* _datasj = datas[j];
	1645	PyArrayObject* _data = convert_to_numpy(_datasj, "data");
	1646	conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
	1647	conversion_numpy_check_size(_data, 1, "data");
	1648	blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
	1649	PyObject* _dvecrowsj = dvecrows[j];
	1650	PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
	1651	conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
	1652	conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
	1653	blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
	1654	int m = row.numElements();
	1655	for(int k=0;k<m;k++)
	1656	{
	1657	dr((cdi+(int)(idt*dvecrow(k)))%md, (int)row(k)) += data(k);
	1658	}
	1659	Py_DECREF(_rowind);
	1660	Py_DECREF(_datasj);
	1661	Py_DECREF(_dvecrowsj);
	1662	}
	1663	"""
	1664	weave.inline(code, ['rowinds', 'datas', 'spikes', 'nspikes',
	1665	'dvecrows', 'dr', 'cdi', 'idt', 'md'],
	1666	compiler=self._cpp_compiler,
	1667	type_converters=weave.converters.blitz,
	1668	extra_compile_args=self._extra_compile_args)
	1669	else:
	1670	nspikes = len(spikes)
	1671	rowinds = [r.ind for r in rows]
	1672	datas = rows
	1673	cdi = self._cur_delay_ind
	1674	idt = self._invtargetdt
	1675	md = self._max_delay
	1676	code = """
	1677	for(int j=0;j<nspikes;j++)
	1678	{
	1679	PyObject* _rowind = rowinds[j];
	1680	PyArrayObject* _row = convert_to_numpy(_rowind, "row");
	1681	conversion_numpy_check_type(_row, PyArray_LONG, "row");
	1682	conversion_numpy_check_size(_row, 1, "row");
	1683	blitz::Array<long,1> row = convert_to_blitz<long,1>(_row,"row");
	1684	PyObject* _datasj = datas[j];
	1685	PyArrayObject* _data = convert_to_numpy(_datasj, "data");
	1686	conversion_numpy_check_type(_data, PyArray_DOUBLE, "data");
	1687	conversion_numpy_check_size(_data, 1, "data");
	1688	blitz::Array<double,1> data = convert_to_blitz<double,1>(_data,"data");
	1689	PyObject* _dvecrowsj = dvecrows[j];
	1690	PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
	1691	conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
	1692	conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
	1693	blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
	1694	int m = row.numElements();
	1695	double mod = sv_pre(spikes(j));
	1696	for(int k=0;k<m;k++)
	1697	{
	1698	dr((cdi+(int)(idtdvecrow(k)))%md, (int)row(k)) += data(k)mod;
	1699	}
	1700	Py_DECREF(_rowind);
	1701	Py_DECREF(_datasj);
	1702	Py_DECREF(_dvecrowsj);
	1703	}
	1704	"""
	1705	weave.inline(code, ['sv_pre', 'rowinds', 'datas', 'spikes', 'nspikes',
	1706	'dvecrows', 'dr', 'cdi', 'idt', 'md'],
	1707	compiler=self._cpp_compiler,
	1708	type_converters=weave.converters.blitz,
	1709	extra_compile_args=self._extra_compile_args)
	1710	else:
	1711	if self._nstate_mod is None:
	1712	if not isinstance(spikes, numpy.ndarray):
	1713	spikes = array(spikes, dtype=int)
	1714	nspikes = len(spikes)
	1715	N = len(self.target)
	1716	cdi = self._cur_delay_ind
	1717	idt = self._invtargetdt
	1718	md = self._max_delay
	1719	code = """
	1720	for(int j=0;j<nspikes;j++)
	1721	{
	1722	PyObject* _rowsj = rows[j];
	1723	PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
	1724	conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
	1725	conversion_numpy_check_size(_row, 1, "row");
	1726	blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
	1727	PyObject* _dvecrowsj = dvecrows[j];
	1728	PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
	1729	conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
	1730	conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
	1731	blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
	1732	for(int k=0;k<N;k++)
	1733	dr((cdi+(int)(idt*dvecrow(k)))%md, k) += row(k);
	1734	Py_DECREF(_rowsj);
	1735	Py_DECREF(_dvecrowsj);
	1736	}
	1737	"""
	1738	weave.inline(code, ['spikes', 'nspikes', 'N', 'rows',
	1739	'dr', 'cdi', 'idt', 'md', 'dvecrows'],
	1740	compiler=self._cpp_compiler,
	1741	type_converters=weave.converters.blitz,
	1742	extra_compile_args=self._extra_compile_args)
	1743	else:
	1744	if not isinstance(spikes, numpy.ndarray):
	1745	spikes = array(spikes, dtype=int)
	1746	nspikes = len(spikes)
	1747	N = len(self.target)
	1748	cdi = self._cur_delay_ind
	1749	idt = self._invtargetdt
	1750	md = self._max_delay
	1751	code = """
	1752	for(int j=0;j<nspikes;j++)
	1753	{
	1754	PyObject* _rowsj = rows[j];
	1755	PyArrayObject* _row = convert_to_numpy(_rowsj, "row");
	1756	conversion_numpy_check_type(_row, PyArray_DOUBLE, "row");
	1757	conversion_numpy_check_size(_row, 1, "row");
	1758	blitz::Array<double,1> row = convert_to_blitz<double,1>(_row,"row");
	1759	PyObject* _dvecrowsj = dvecrows[j];
	1760	PyArrayObject* _dvecrow = convert_to_numpy(_dvecrowsj, "dvecrow");
	1761	conversion_numpy_check_type(_dvecrow, PyArray_DOUBLE, "dvecrow");
	1762	conversion_numpy_check_size(_dvecrow, 1, "dvecrow");
	1763	blitz::Array<double,1> dvecrow = convert_to_blitz<double,1>(_dvecrow,"dvecrow");
	1764	double mod = sv_pre(spikes(j));
	1765	for(int k=0;k<N;k++)
	1766	dr((cdi+(int)(idtdvecrow(k)))%md, k) += row(k)mod;
	1767	Py_DECREF(_rowsj);
	1768	Py_DECREF(_dvecrowsj);
	1769	}
	1770	"""
	1771	weave.inline(code, ['sv_pre', 'spikes', 'nspikes', 'N', 'rows',
	1772	'dr', 'cdi', 'idt', 'md', 'dvecrows'],
	1773	compiler=self._cpp_compiler,
	1774	type_converters=weave.converters.blitz,
	1775	extra_compile_args=self._extra_compile_args)
	1776
	1777	def do_propagate(self):
	1778	self.propagate(self.source.get_spikes(0))
	1779
	1780	def _set_delay_property(self, val):
	1781	self.delayvec[:] = val
	1782
	1783	delay = property(fget=lambda self:self.delayvec, fset=_set_delay_property)
	1784
	1785	def compress(self):
	1786	if not self.iscompressed:
	1787	# We want delayvec to have nonzero entries at the same places as
	1788	# W does, so we use W to initialise the compressed version of
	1789	# delayvec, and then copy the values from the old delayvec to
	1790	# the new compressed one, allowing delayvec and W to not have
	1791	# to be perfectly intersected at the initialisation stage. If
	1792	# the structure is dynamic, it will be the user's
	1793	# responsibility to update them in sequence
	1794	delayvec = self.delayvec
	1795	# The delayvec[i,:] operation for sparse.lil_matrix format
	1796	# is VERY slow, but the CSR format is fine.
	1797	if isinstance(delayvec, sparse.lil_matrix):
	1798	delayvec = delayvec.tocsr()
	1799	self.delayvec = self.W.connection_matrix(copy=True)
	1800	for i in xrange(self.W.shape[0]):
	1801	self.delayvec.set_row(i, array(todense(delayvec[i, :]), copy=False).flatten())
	1802	Connection.compress(self)
	1803
	1804	def set_delays(self, source=None, target=None, delay=None):
	1805	'''
	1806	Set the delays corresponding to the weight matrix
	1807
	1808	``delay`` must be one of:
	1809
	1810	* A float, all delays will be set to this value
	1811	* A pair (min, max), delays will be uniform between these two
	1812	values.
	1813	* A function of no arguments, will be called for each nonzero
	1814	entry in the weight matrix.
	1815	* A function of two argument ``(i,j)`` will be called for each
	1816	nonzero entry in the weight matrix.
	1817	* A matrix of an appropriate type (e.g. ndarray or lil_matrix).
	1818	'''
	1819	if delay is None:
	1820	return
	1821	W = self.W
	1822	P = source or self.source
	1823	Q = target or self.target
	1824	i0, j0 = self.origin(P, Q)
	1825	i1 = i0 + len(P)
	1826	j1 = j0 + len(Q)
	1827	if isinstance(W, sparse.lil_matrix):
	1828	def getrow(i):
	1829	inds = array(W.rows[i], dtype=int)
	1830	inds = inds[logical_and(inds >= j0, inds < j1)]
	1831	return inds, len(inds)
	1832	else:
	1833	def getrow(i):
	1834	inds = (W[i, j0:j1] != 0).nonzero()[0] + j0
	1835	return inds, len(inds)
	1836	#return slice(j0, j1), j1-j0
	1837	if isinstance(delay, (float, int)):
	1838	for i in xrange(i0, i1):
	1839	inds, L = getrow(i)
	1840	self.delayvec[i, inds] = delay
	1841	elif isinstance(delay, (tuple, list)) and len(delay) == 2:
	1842	delaymin, delaymax = delay
	1843	for i in xrange(i0, i1):
	1844	inds, L = getrow(i)
	1845	rowdelay = rand(L) * (delaymax - delaymin) + delaymin
	1846	self.delayvec[i, inds] = rowdelay
	1847	elif callable(delay) and delay.func_code.co_argcount == 0:
	1848	for i in xrange(i0, i1):
	1849	inds, L = getrow(i)
	1850	rowdelay = [delay() for _ in xrange(L)]
	1851	self.delayvec[i, inds] = rowdelay
	1852	elif callable(delay) and delay.func_code.co_argcount == 2:
	1853	for i in xrange(i0, i1):
	1854	inds, L = getrow(i)
	1855	if isinstance(inds, slice):
	1856	inds = numpy.arange(inds.start, inds.stop)
	1857	self.delayvec[i, inds] = delay(i - i0, inds - j0)
	1858	else:
	1859	#raise TypeError('delays must be float, pair or function of 0 or 2 arguments')
	1860	self.delayvec[i0:i1, j0:j1] = delay # probably won't work, but then it will raise an error
	1861
	1862	def connect(self, source=None, target=None, W=None, delay=None):
	1863	Connection.connect(self, source=source, target=target, W=W)
	1864	if delay is not None:
	1865	self.set_delays(source, target, delay)
	1866
	1867	def connect_random(self, source=None, target=None, p=1.0, weight=1.0,
	1868	fixed=False, seed=None, sparseness=None, delay=None):
	1869	Connection.connect_random(self, source=source, target=target, p=p,
	1870	weight=weight, fixed=fixed, seed=seed,
	1871	sparseness=sparseness)
	1872	if delay is not None:
	1873	self.set_delays(source, target, delay)
	1874
	1875	def connect_full(self, source=None, target=None, weight=1.0, delay=None):
	1876	Connection.connect_full(self, source=source, target=target, weight=weight)
	1877	if delay is not None:
	1878	self.set_delays(source, target, delay)
	1879
	1880	def connect_one_to_one(self, source=None, target=None, weight=1.0, delay=None):
	1881	Connection.connect_one_to_one(self, source=source, target=target,
	1882	weight=weight)
	1883	if delay is not None:
	1884	self.set_delays(source, target, delay)
	1885
	1886
	1887	class IdentityConnection(Connection):
	1888	'''
	1889	A :class:`Connection` between two groups of the same size, where neuron ``i`` in the
	1890	source group is connected to neuron ``i`` in the target group.
	1891
	1892	Initialised with arguments:
	1893
	1894	``source``, ``target``
	1895	The source and target :class:`NeuronGroup` objects.
	1896	``state``
	1897	The target state variable.
	1898	``weight``
	1899	The weight of the synapse, must be a scalar.
	1900	``delay``
	1901	Only homogeneous delays are allowed.
	1902
	1903	The benefit of this class is that it has no storage requirements and is optimised for
	1904	this special case.
	1905	'''
	1906	@check_units(delay=second)
	1907	def __init__(self, source, target, state=0, weight=1, delay=0 * msecond):
	1908	if (len(source) != len(target)):
	1909	raise AttributeError, 'The connected (sub)groups must have the same size.'
	1910	self.source = source # pointer to source group
	1911	self.target = target # pointer to target group
	1912	if type(state) == types.StringType: # named state variable
	1913	self.nstate = target.get_var_index(state)
	1914	else:
	1915	self.nstate = state # target state index
	1916	self.W = float(weight) # weight
	1917	source.set_max_delay(delay)
	1918	self.delay = int(delay / source.clock.dt) # Synaptic delay in time bins
	1919
	1920	def propagate(self, spikes):
	1921	'''
	1922	Propagates the spikes to the target.
	1923	'''
	1924	self.target._S[self.nstate, spikes] += self.W
	1925
	1926	def compress(self):
	1927	pass
	1928
	1929
	1930	class MultiConnection(Connection):
	1931	'''
	1932	A hub for multiple connections with a common source group.
	1933	'''
	1934	def __init__(self, source, connections=[]):
	1935	self.source = source
	1936	self.connections = connections
	1937	self.iscompressed = False
	1938	self.delay = connections[0].delay
	1939
	1940	def propagate(self, spikes):
	1941	'''
	1942	Propagates the spikes to the targets.
	1943	'''
1978	1944	for C in self.connections:
1979		C.compress()
1980		self.iscompressed = True
1981
1982
1983		# Generation of matrices
1984		def random_matrix(n, m, p, value=1.):
1985		'''
1986		Generates a sparse random matrix with size (n,m).
1987		Entries are 1 (or optionnally value) with probability p.
1988		If value is a function, then that function is called for each
1989		non zero element as value() or value(i,j).
1990		'''
1991		# TODO:
1992		# Simplify (by using valuef)
1993		W = sparse.lil_matrix((n, m))
1994		if callable(value) and callable(p):
1995		if value.func_code.co_argcount == 0:
1996		valuef = lambda i, j:[value() for _ in j] # value function
1997		elif value.func_code.co_argcount == 2:
1998		try:
1999		failed = (array(value(0, arange(m))).size != m)
2000		except:
2001		failed = True
2002		if failed: # vector-based not possible
2003		log_debug('connections', 'Cannot build the connection matrix by rows')
2004		valuef = lambda i, j:[value(i, k) for k in j]
2005		else:
2006		valuef = value
2007		else:
2008		raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"
2009
2010		if p.func_code.co_argcount == 2:
2011		# Check if p(i,j) is vectorisable
2012		try:
2013		failed = (array(p(0, arange(m))).size != m)
2014		except:
2015		failed = True
2016		if failed: # vector-based not possible
2017		log_debug('connections', 'Cannot build the connection matrix by rows')
	1945	C.propagate(spikes)
	1946
	1947	def compress(self):
	1948	if not self.iscompressed:
	1949	for C in self.connections:
	1950	C.compress()
	1951	self.iscompressed = True
	1952
	1953
	1954	# Generation of matrices
	1955	def random_matrix(n, m, p, value=1.):
	1956	'''
	1957	Generates a sparse random matrix with size (n,m).
	1958	Entries are 1 (or optionnally value) with probability p.
	1959	If value is a function, then that function is called for each
	1960	non zero element as value() or value(i,j).
	1961	'''
	1962	# TODO:
	1963	# Simplify (by using valuef)
	1964	W = sparse.lil_matrix((n, m))
	1965	if callable(value) and callable(p):
	1966	if value.func_code.co_argcount == 0:
	1967	valuef = lambda i, j:[value() for _ in j] # value function
	1968	elif value.func_code.co_argcount == 2:
	1969	try:
	1970	failed = (array(value(0, arange(m))).size != m)
	1971	except:
	1972	failed = True
	1973	if failed: # vector-based not possible
	1974	log_debug('connections', 'Cannot build the connection matrix by rows')
	1975	valuef = lambda i, j:[value(i, k) for k in j]
	1976	else:
	1977	valuef = value
	1978	else:
	1979	raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"
	1980
	1981	if p.func_code.co_argcount == 2:
	1982	# Check if p(i,j) is vectorisable
	1983	try:
	1984	failed = (array(p(0, arange(m))).size != m)
	1985	except:
	1986	failed = True
	1987	if failed: # vector-based not possible
	1988	log_debug('connections', 'Cannot build the connection matrix by rows')
	1989	for i in xrange(n):
	1990	W.rows[i] = [j for j in range(m) if rand() < p(i, j)]
	1991	W.data[i] = list(valuef(i, array(W.rows[i])))
	1992	else: # vector-based possible
	1993	for i in xrange(n):
	1994	W.rows[i] = list((rand(m) < p(i, arange(m))).nonzero()[0])
	1995	W.data[i] = list(valuef(i, array(W.rows[i])))
	1996	elif p.func_code.co_argcount == 0:
2018	1997	for i in xrange(n):
2019		W.rows[i] = [j for j in range(m) if rand() < p(~~i, j~~)]
	1998	W.rows[i] = [j for j in range(m) if rand() < p()]
2020	1999	W.data[i] = list(valuef(i, array(W.rows[i])))
2021		else: # vector-based possible
2022		for i in xrange(n):
2023		W.rows[i] = list((rand(m) < p(i, arange(m))).nonzero()[0])
2024		W.data[i] = list(valuef(i, array(W.rows[i])))
2025		elif p.func_code.co_argcount == 0:
2026		for i in xrange(n):
2027		W.rows[i] = [j for j in range(m) if rand() < p()]
2028		W.data[i] = list(valuef(i, array(W.rows[i])))
2029		else:
2030		raise AttributeError, "Bad number of arguments in p function (should be 2)"
2031		elif callable(value):
2032		if value.func_code.co_argcount == 0: # TODO: should work with partial objects
2033		for i in xrange(n):
2034		k = random.binomial(m, p, 1)[0]
2035		W.rows[i] = sample(xrange(m), k)
2036		W.rows[i].sort()
2037		W.data[i] = [value() for _ in xrange(k)]
2038		elif value.func_code.co_argcount == 2:
2039		try:
2040		failed = (array(value(0, arange(m))).size != m)
2041		except:
2042		failed = True
2043		if failed: # vector-based not possible
2044		log_debug('connections', 'Cannot build the connection matrix by rows')
	2000	else:
	2001	raise AttributeError, "Bad number of arguments in p function (should be 2)"
	2002	elif callable(value):
	2003	if value.func_code.co_argcount == 0: # TODO: should work with partial objects
2045	2004	for i in xrange(n):
2046	2005	k = random.binomial(m, p, 1)[0]
2047	2006	W.rows[i] = sample(xrange(m), k)
2048	2007	W.rows[i].sort()
	2008	W.data[i] = [value() for _ in xrange(k)]
	2009	elif value.func_code.co_argcount == 2:
	2010	try:
	2011	failed = (array(value(0, arange(m))).size != m)
	2012	except:
	2013	failed = True
	2014	if failed: # vector-based not possible
	2015	log_debug('connections', 'Cannot build the connection matrix by rows')
	2016	for i in xrange(n):
	2017	k = random.binomial(m, p, 1)[0]
	2018	W.rows[i] = sample(xrange(m), k)
	2019	W.rows[i].sort()
	2020	W.data[i] = [value(i, j) for j in W.rows[i]]
	2021	else:
	2022	for i in xrange(n):
	2023	k = random.binomial(m, p, 1)[0]
	2024	W.rows[i] = sample(xrange(m), k)
	2025	W.rows[i].sort()
	2026	W.data[i] = list(value(i, array(W.rows[i])))
	2027	else:
	2028	raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"
	2029	elif callable(p):
	2030	if p.func_code.co_argcount == 2:
	2031	# Check if p(i,j) is vectorisable
	2032	try:
	2033	failed = (array(p(0, arange(m))).size != m)
	2034	except:
	2035	failed = True
	2036	if failed: # vector-based not possible
	2037	log_debug('connections', 'Cannot build the connection matrix by rows')
	2038	for i in xrange(n):
	2039	W.rows[i] = [j for j in range(m) if rand() < p(i, j)]
	2040	W.data[i] = [value] * len(W.rows[i])
	2041	else: # vector-based possible
	2042	for i in xrange(n):
	2043	W.rows[i] = list((rand(m) < p(i, arange(m))).nonzero()[0])
	2044	W.data[i] = [value] * len(W.rows[i])
	2045	elif p.func_code.co_argcount == 0:
	2046	for i in xrange(n):
	2047	W.rows[i] = [j for j in range(m) if rand() < p()]
	2048	W.data[i] = [value] * len(W.rows[i])
	2049	else:
	2050	raise AttributeError, "Bad number of arguments in p function (should be 2)"
	2051	else:
	2052	for i in xrange(n):
	2053	k = random.binomial(m, p, 1)[0]
	2054	# Not significantly faster to generate all random numbers in one pass
	2055	# N.B.: the sample method is implemented in Python and it is not in Scipy
	2056	W.rows[i] = sample(xrange(m), k)
	2057	W.rows[i].sort()
	2058	W.data[i] = [value] * k
	2059
	2060	return W
	2061
	2062	def random_matrix_fixed_column(n, m, p, value=1.):
	2063	'''
	2064	Generates a sparse random matrix with size (n,m).
	2065	Entries are 1 (or optionnally value) with probability p.
	2066	The number of non-zero entries by per column is fixed: (int)(p*n)
	2067	If value is a function, then that function is called for each
	2068	non zero element as value() or value(i,j).
	2069	'''
	2070	W = sparse.lil_matrix((n, m))
	2071	k = (int)(p * n)
	2072	for j in xrange(m):
	2073	# N.B.: the sample method is implemented in Python and it is not in Scipy
	2074	for i in sample(xrange(n), k):
	2075	W.rows[i].append(j)
	2076
	2077	if callable(value):
	2078	if value.func_code.co_argcount == 0:
	2079	for i in xrange(n):
	2080	W.data[i] = [value() for _ in xrange(len(W.rows[i]))]
	2081	elif value.func_code.co_argcount == 2:
	2082	for i in xrange(n):
2049	2083	W.data[i] = [value(i, j) for j in W.rows[i]]
2050	2084	else:
2051		for i in xrange(n):
2052		k = random.binomial(m, p, 1)[0]
2053		W.rows[i] = sample(xrange(m), k)
2054		W.rows[i].sort()
2055		W.data[i] = list(value(i, array(W.rows[i])))
	2085	raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"
2056	2086	else:
2057		~~raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"~~
2058		~~elif callable(p):~~
2059		~~if p.func_code.co_argcount == 2:~~
2060		~~# Check if p(i,j) is vectorisable~~
2061		~~try:~~
2062		~~failed = (array(p(0, arange(m))).size != m)~~
2063		~~except:~~
2064		~~failed = True~~
2065		~~if failed: # vector-based not possible~~
2066		~~log_debug('connections', 'Cannot build the connection matrix by rows')~~
2067		~~for i in xrange(n):~~
2068		~~W.rows[i] = [j for j in range(m) if rand() < p(i, j)]~~
2069		~~W.data[i] = [value] * len(W.rows[i])~~
2070		~~else: # vector-based possible~~
2071		~~for i in xrange(n):~~
2072		~~W.rows[i] = list((rand(m) < p(i, arange(m))).nonzero()[0])~~
2073		~~W.data[i] = [value] * len(W.rows[i])~~
2074		~~elif p.func_code.co_argcount == 0:~~
2075	2087	for i in xrange(n):
2076		~~W.rows[i] = [j for j in range(m) if rand() < p()]~~
2077	2088	W.data[i] = [value] * len(W.rows[i])
2078		else:
2079		raise AttributeError, "Bad number of arguments in p function (should be 2)"
2080		else:
2081		for i in xrange(n):
2082		k = random.binomial(m, p, 1)[0]
2083		# Not significantly faster to generate all random numbers in one pass
2084		# N.B.: the sample method is implemented in Python and it is not in Scipy
2085		W.rows[i] = sample(xrange(m), k)
2086		W.rows[i].sort()
2087		W.data[i] = [value] * k
	2089
	2090	return W
	2091
	2092	def eye_lil_matrix(n):
	2093	'''
	2094	Returns the identity matrix of size n as a lil_matrix
	2095	(sparse matrix).
	2096	'''
	2097	M = sparse.lil_matrix((n, n))
	2098	M.setdiag([1.] * n)
	2099	return M
	2100
	2101	# this is put down here because it is needed in the code above but the order of imports is messed up
	2102	# if it is included at the top
	2103	from network import network_operation
2088	2104
2089		~~return W~~
2090
2091		~~def random_matrix_fixed_column(n, m, p, value=1.):~~
2092		~~'''~~
2093		~~Generates a sparse random matrix with size (n,m).~~
2094		~~Entries are 1 (or optionnally value) with probability p.~~
2095		~~The number of non-zero entries by per column is fixed: (int)(p*n)~~
2096		~~If value is a function, then that function is called for each~~
2097		~~non zero element as value() or value(i,j).~~
2098		~~'''~~
2099		~~W = sparse.lil_matrix((n, m))~~
2100		~~k = (int)(p * n)~~
2101		~~for j in xrange(m):~~
2102		~~# N.B.: the sample method is implemented in Python and it is not in Scipy~~
2103		~~for i in sample(xrange(n), k):~~
2104		~~W.rows[i].append(j)~~
2105
2106		~~if callable(value):~~
2107		~~if value.func_code.co_argcount == 0:~~
2108		~~for i in xrange(n):~~
2109		~~W.data[i] = [value() for _ in xrange(len(W.rows[i]))]~~
2110		~~elif value.func_code.co_argcount == 2:~~
2111		~~for i in xrange(n):~~
2112		~~W.data[i] = [value(i, j) for j in W.rows[i]]~~
2113		~~else:~~
2114		~~raise AttributeError, "Bad number of arguments in value function (should be 0 or 2)"~~
2115		~~else:~~
2116		~~for i in xrange(n):~~
2117		~~W.data[i] = [value] * len(W.rows[i])~~
2118
2119		~~return W~~
2120
2121		~~def eye_lil_matrix(n):~~
2122		~~'''~~
2123		~~Returns the identity matrix of size n as a lil_matrix~~
2124		~~(sparse matrix).~~
2125		~~'''~~
2126		~~M = sparse.lil_matrix((n, n))~~
2127		~~M.setdiag([1.] * n)~~
2128		~~return M~~
2129
2130		~~# this is put down here because it is needed in the code above but the order of imports is messed up~~
2131		~~# if it is included at the top~~
2132		~~from network import network_operation~~

Context Navigation

Changeset 2117 for trunk/brian/connection.py

Legend:

trunk/brian/connection.py

Download in other formats: