Thanks for an interesting question. My review starts with your main interest -- speeding up what you've already got. It then goes to discuss a number of other topics you may be interested in.
Speedup
To figure out what's slow in your code, you've got to profile it. Python make this easy with the line_profiler
package.
I used this tool on your code and found that the np.concatenate
was responsible for 88% of the run time.
After thinking about it, this didn't surprise me. The way you build ee_i
is very inefficient. Keep in mind the output of np.concatenate()
is a new array. So every loop iteration requires instantiating a new NumPy array and copying all the existing ee_i
data into it.
A better way is to pre-allocate the ee_i
array with something like this:
ee_j = np.repeat(np.arange(num_neurons),
num_cnxns)
ee_i = np.zeros_like(ee_j, dtype = 'int')
### REST OF CODE ###
Then, instead of copying all of the existing ee_i
data into a new array every time you go through the loop, just update the relevant positions of ee_i
using NumPy slice notation.
ee_i[start_idx:stop_idx] = ##<code that supplies new values>##
This sped up the code by >5-fold for me.
Other issues
Why encoders = np.random.choice(np.arange(N_E), (p, int(f * N_E)), False)
instead of just arbitrarily specificying that the first f * N_E / p
neurons are part of the first population of encoders, the next, from f * N_E /p + 1
to 2 * f * N_E /p
, the second population of encoders, etc. Since the number and ordering neurons is (presumably) arbitrary, it shouldn't matter if you number them in a way that's convenient for you.
It's tough to understand your variable names. What is f
? What is c
? It would be nice to have a docstring that explained at least this much.
A different (faster) approach
Conceptually, the connections between neurons can be thought of as a square matrix of boolean values. If there is a connection from neuron i to neuron j, then this matrix will have a True
at position (i, j)
. This is a much different representation than the array (i.e. list) of sources and targets you are trying to generate. I find it much more intuitive.
The downside of the matrix representation is that many of the entries will be 0, and it requires creating a very large matrix, which is memory (and sometimes speed) inefficient. Thus, the lists of sources and targets you are trying to create are much more efficient.
However sparse matrices offer a way to get the intuitive advantages of the matrix representation with the speed and memory properties of the "linked list" representation you are shooting for. In particular SciPy's COO sparse matrix format has an underlying representation that is nearly exactly the same as your list of target nodes and source nodes (i.e. a list of the column indices and a list of the row indices with nonzero values).
Here's a version of that approach:
from scipy import sparse
import numpy as np
def rand_col_sparse_bool_arr(m, n, c):
"""
Returns a sparse boolean matrix (COO format) of size m by n with a total of c nonzero entries per column
:param: m int number of rows in output matrix
:param: n int number of columns in output matrix; every column has c nonzero entries
:param: c int number of nonzero entries in each column of output matrix
"""
if c > m:
raise(ValueError, 'c must be less than m.')
data = np.repeat([True], repeats = c*n)
row_idxs = np.random.randint(m, size = c*n)
col_idxs = np.repeat(np.arange(n), repeats = c)
# return the sparse matrix
return sparse.coo_matrix((data, (row_idxs, col_idxs)),
dtype = bool,
shape = (m, n))
This first function is a generalized function that creates random boolean matrices with non-zero elements in every column at random row locations. It can be used to create all the "sub-matrices" for each population of selective/encoding neurons and also for the other neurons.
I tried to make a function that has the population you want and I think I did so. But you should check it carefully to make sure the population structure is really want you want. I may have gotten the rows and columns (sources and sinks) backwards, for example. Even if there are errors, I think this code serves to illustrate the approach:
def generate_network(n_e = 8000, p = 5, f = 0.1, c = 0.2):
# Generate a matrix of connections for a neural network as described in at Code Review
# http://codereview.stackexchange.com/questions/129645/assembling-edges-of-a-graph
# unpack parameters
num_selec_cnxns_per_pop = int(c * f * n_e)
num_selec_cnxns = int(num_selec_cnxns_per_pop * p)
num_nonselec_cnxns = int(c * (1 - f*p) * n_e)
num_cnxns = num_selec_cnxns + num_nonselec_cnxns
num_encoders_per_pop = int(f * n_e / p)
num_encoders = int(num_encoders_per_pop * p)
num_nonselec = n_e - num_encoders
# make the encoder matrices
encoder_connectivity_mat = sparse.hstack(
[rand_col_sparse_bool_arr(n_e,
num_encoders_per_pop,
num_selec_cnxns_per_pop) for _ in xrange(p)])
# make the nonselective matrix
nonselective_connectivity_mat = rand_col_sparse_bool_arr(n_e,
num_nonselec,
num_nonselec_cnxns)
# stack the matrices together
connectivity_mat = sparse.hstack([encoder_connectivity_mat,
nonselective_connectivity_mat])
# check the size of the output
assert connectivity_mat.shape == (n_e, n_e)
# derive and return connectivity in requested format
ee_i = connectivity_mat.col
ee_j = connectivity_mat.row
return ee_i, ee_j
This approach relies on a predefined ordering of the different nodes. The first rows are the first population of selective neurons, etc., and the final rows are the non-selective neurons.
This approach took about 0.25 seconds for your parameters on my machine, compared to more than five minutes for your original code.