brian-team / brian2 Goto Github PK
View Code? Open in Web Editor NEWBrian is a free, open source simulator for spiking neural networks.
Home Page: http://briansimulator.org
License: Other
Brian is a free, open source simulator for spiking neural networks.
Home Page: http://briansimulator.org
License: Other
Consistently deal with user-provided functions, make them work with code generation. Unit-checking for functions will be supported via the check_units
decorator (that we'll make mandatory), quoting from an older issue:
So a function that returns time-dependent voltage would use
@check_units(t=second, result=volt)
and a simple square function would use@check_units(x=None, result=lambda x_unit: x_unit**2)
.
One general question: If we make the code-generation language a preference and keep using the same language across objects, do we make it unnecessarily difficult to use user-provided functions in synapse creation or for setting state variables? The scenario I'm thinking of is using C++ code-generation, but using simple Python functions in the expressions for synapse creation or state variable setting. With a "one language" approach, one would be forced to provide C++ code for these functions.
As we discussed before, we don't want to recreate the (very sparsely populated) model library from brian, e.g. the Izhikevich
, HH
, or alpha_current
functions. Rather we want to have a webpage (or maybe wikipage) listing these. Ideally, all should be in a standard format, e.g. have a reference and the mathematical equations (autogenerated from the code :) ?). The code itself should then be in a way that is easily copy/pasteable. This is different from full example script it would really only focus on the neural or synaptic model. We should probably agree on some conventions before hand, e.g. always call the membrane potential v
or V
, etc.
This would also be a nice way for users who are not willing/able/confident enough to work on the Brian code itself to contribute to Brian.
For functions, users will need to specify the input/output relationships with an @check_units
decorator. In some cases, we might be able to automatically parse simple functions and determine the units automatically.
After fixing #38 I thought a bit about the current refractoriness system and what kind of refractoriness we can implement with it. The docstring of NeuronGroup states about the refractory
variable: "This value can be modified in the reset code. (TODO: more modifiability)" -- @thesamovar do you remember what you meant by this?
From my experience, there are three types of absolute refractoriness models (and combinations of them)
(active)
flagis_active
variable: threshold='(v > v_t) * is_active'
pre
code of a Synapses, e.g.: pre='v += w*is_active
(is_active
is is_active_post
).Anything missing here that we don't yet support? Regardless of this, I wonder whether it would not make more sense to replace refractory_until
by lastspike
, i.e. do
is_active = (t - lastspike) < refractory
instead of
is_active = t > refractory_until
On the other hand, we'd need some dummy value for the initial lastspike, that would make things a bit ugly -- so probably it's better to stay with refractory_until
(one can always do refractory_until - refractory
to get the last spike, anyway). There would be also some subtle issues, e.g. the semantics of changing the refractory
variable in the reset code changes.
Two changes to Network:
net.objects['my_obj']
, net['my_obj']
or net.my_obj
(obviously I like the last one best).Like this, users can construct and move around Network objects without having to pass lists or dicts of all their objects.
e.g. when writing Quantity
you should get a link to the reference docs for Quantity, but often you don't. Similarly, if you have a class X with a method f then writing f
in the docs for X should give you a link to that method.
The existing string tools get_identifiers
and analyse_identifiers
should be given more rigorous version using ast parsing. We can still keep these functions though, for example for analysing C++ code.
Previously we have been using sympy for translation from Python syntax to C++, but this doesn't allow for Java and has some restrictions. It is also not difficult to implement this ourselves, allowing us more control over the process and also making it easier to do things like Java for the BrianDROID project.
This module should ideally be a relatively independent piece of code with enough options that it can do all the things we'd like it to. Here are some requirements, we might need more:
Language
, or even more simply just a function for each language with the same signature.Anything else?
If the above is enough, a nice minimal interface would be just to pass the expression as a string, and a dict of variable/function name replacements. I don't see any obvious reason to complicate it beyond this.
Currently, this line (https://github.com/brian-team/brian2/blob/master/brian2/tests/test_network.py#L286):
assert_raises(MagicError, lambda: run(1*ms, level=3))
fails to raise the expected MagicError
when tested with Python 3.3 -- it's working as expected under Python 2.x
We should be able to declare state variables as being of boolean type by writing bool in the units slot.
Some of the code in units assumes there is a numpy.float128
object which doesn't exist on all platforms. This should be checked for. It crashes on Windows 7, 32 bit with Python 2.7 and numpy 1.6.2 and 1.7.1.
I stumbled across a quite nasty bug that led to wrong results when using refractoriness via the (active)
flag -- it took me quite a while to realize what the problem is... We set the is_active
variable via is_active = (t >= refractory_until)
, the result is an array of boolean values. Unfortunately, there's one major difference between a boolean value in Python and in numpy: The minus sign works as a negation in numpy whereas it just treats the boolean value as an integer in Python:
>>> -True
-1
>>> -False
0
>>> print -np.array([True, False])
[False True]
A nasty consequence of this behaviour is that symbol rearrangements done by sympy might change the value:
-is_active*(x - 1) # -is_active will be logically inverted
is_active*(1 - x) # is_active will be treated as a number
This is a Python-specific problem and we can circumvent it by explicitly multiplying is_active
by 1
when it is calculated the first time. We should keep this issue in mind whenever we deal with boolean values.
This is a follow-up to the work done in #6:
This is of course a big chunk of work but I think we should rather start it sooner than later because it is a core component and probably a lot of unforeseen problems will pop up during its implementation.
I had a quick look through the code in Brian 1 and will use this issue to jot down some observations:
SynapticEquations
class (it was only introduce to allow for the (event-driven)
keyword, support for _pre
and _post
variables can be easily implemented using the specifiers, no need for implementing them via static equations, etc.SynapticVariable
class, for example, is in fact very similar to a specifier in the current system (this also entails that Synapses
can use the standard __getattr__
mechanism of Group
).SpikeQueue
fits into the current system: It could quite naturally be handled in a Devices
system (see discussion in #23 ) because it is about storing data. On the other hand, we probably want to have a C version available in the standard runtime version for performance reasons, this would rather suggest to connect it to Language
(i.e. automatically use the C version if we use C code generation).This will be necessary for Python 3 because weave is not being ported to Python 3 (very sad, I know).
At least on Windows 32, Python 2.6, if a program raises an exception with Brian 2 installed, the exception is silently discarded and the program terminates.
To allow for reproducible simulations, we'll need a mechanism to reproduce random number generation. While setting the seed is sufficient in many cases, this is not guaranteed to give the same results across platforms, numpy versions etc. The problem gets even trickier when considering multiple processes or GPU code where we do not want the results to depend on the number of processes, for example.
OK a thought on how to efficiently implement synapse creation. The problem to solve is that if you just compute cond(i,j)
for all pairs (i, j) it is O(N^2)
and in many cases this is highly inefficient. We should list some of the cases and propose solutions.
The two most important that I can think of are:
rand()<p
should be replaced with binomial distribution sampling, and also compounds like rand()<p and cond(i,j)
should be replaced by binomial distribution sampling followed by condition evaluation only on that subset.i==j
should be replaced with one-to-oneSo there are various possible solutions here - it would be nice to come up with something quite general that can adapt and become more efficient as we add special cases, but that doesn't require us to rewrite all the templates each time.
The most obvious implementation of i==j
is just to explicitly detect it and convert it to connect_one_to_one
and require that devices implement this connectivity function. Similarly for rand()<p
and connect_random
. I don't really like this approach because it's not very general.
Another option is to have a set of flags that are passed to codegen and the template can make use of that extra information. For example, if the condition was rand()<p and f(i,j)
we could pass to codegen sparseness=p
and cond=f(i,j)
. The code would then look something like:
for i in xrange(_num_source_neurons):
if sparseness<1:
j = binomial_sampling_code_here
else:
j = arange(_num_target_neurons)
_cond = f(i, j)
...
A more general and flexible mechanism would be to handle all conditions of the form a and b and c
, parse them into a set [a, b, c]
, decode this set a bit and then pass that to codegen. For example, a
might be rand()<p
, b
might be abs(i-j)<=width
and c
might be i!=j
. This would be parsed into:
a: ('rand()<p', 'binomial', p)
b: ('abs(i-j)<=width', 'abswidth', width)
c: ('i!=j', None, None)
The parsed form is (original_expr, specialisation_name, argument)
. The codegen template could then look through this list, check which ones it supported, implement them directly if they are supported and if not fall back on the default mechanism.
It might be the case though that it's only worth doing one specialisation, because you won't improve the order of the computation by doing multiple specialisations assuming that each one returns an array of candidate indices. In fact, you might even slow it down. So another option is that each specialisation has an estimate of the number of candidate synapses j it will generate for any given i, and we simply choose the best one. For example if there was rand()<0.5 and i==j
then rand()<0.5
has an estimated 0.5*_num_target_neurons
for each i, whereas i==j
has precisely 1, and in most cases 1<0.5*_num_target_neurons
so this will be selected.
Any thoughts on the above?
We decided on the preferences system detailed in:
https://github.com/brian-team/brian2/blob/master/docs_sphinx/developer/preferences.rst
The _cond
variable is given a double data type in codegen, it should be bool. This used to be handled by OutputVariable
specifier.
And integrate with the system described in #11
See the notes about several issues surrounding linear and "multi-linear" equations in the document refractoriness and linearity
This issue supersedes #23 (although read the discussion thread for lots of interesting background).
We agreed to create a new framework for standalone code generation which requires a redesign of all the Brian objects. The idea is to have a new object Device
that handles allocation of memory and computational work. Core Brian objects like NeuronGroup
will handle parsing arguments and deciding what memory to allocate and what computational work to do, and request this from the Device
object. In runtime mode, this memory will actually be allocated and code compiled. In standalone mode, code will only be generated to do this. Like this, the core Brian objects will be the same for both standalone and runtime modes, only the Device
object will change. We will have a RuntimeDevice
and StandaloneDevice
for this.
Soon, I'll write a more detailed document outlining the system.
Some issues that immediately jump out are:
NeuronGroup
is pretty much worked out, but how will Synapses
work with this framework?CodeObject
or will we do something different for these objects?How should we proceed in implementing this since it requires major changes to everything in Brian. I suggest we try to work out an initial design for Device
at least (an API say), then we will need a concentrated burst of work to rework the existing code to work with that.
At this point, there is a big danger of different branches of the repository going too far out of sync, so I think we need to coordinate when and how we do this.
To not clutter #42 too much, here's something of peripheral importance that occurred to me: We could easily provide access to the number of incoming and outgoing synapses per source/target neuron. That would be somewhat useful for normalizing a weight attribute (but in this case one can get the information easily by summing along an axis, etc.) but if it is provided as a specifier, then one can access it in code as well. That would provide enormous memory benefits in some simple models with weights that are constant up to a normalization.
Example how I imagine it would work:
w = 1
S = Synapses(source, target, pre='v += w / _incoming_synapses`)
Alternatively, one could implement this probably using the synaptic_sum
mechanism and access to the number of synapses. e.g.:
G = NeuronGroup(N, '''dv/dt = (-v + I)/ tau : 1
dI/dt = -I / tau : 1
normalization : 1''')
S = Synapses(source, target, pre='I += 1`)
G.normalization = synaptic_sum(S, 'num_synapses')
Integrate with the system described in #11
It would be nice to have slightly more sophisticated progress reporting and notification when the simulation is complete.
For progress reporting:
ProgressReporter
almost has this, but it's not recursive.For notification, the idea would be to let the user know a long-running simulation has finished which doesn't require them to monitor a console. For example, a pop up notification on the system tray in Windows, or an email.
There might be standard apps that you can use for notification - so perhaps it's not worth adding it to Brian?
In the boolean_parameter
branch, I made some adjustments to handle the case of subexpressions that reference other subexpressions. Codegen does not handle this case correctly, though. I tried a model like this:
G = NeuronGroup(1, '''
dv/dt = 100*Hz : 1 (unless-refractory)
time_since_spike = t - lastspike : second
ref_subexpression = time_since_spike < 5*ms : bool''',
threshold='v>1', reset='v=0',
refractory='ref_subexpression')
This now works in principle, but the generated code is wrong, it looks like:
...
ref_subexpression = time_since_spike < 5 * ms
time_since_spike = t - lastspike
not_refractory = (not_refractory) + (logical_not(ref_subexpression))
...
But ref_subexpression
should be updated after time_since_spike
. We already do the topological sorting in the Equations
object, but this information is not propagated to codegen (I'm also not sure that, say, attaching a simple number to the Subexpression
specifier, derived from the position in the Equations
object is enough. This might become a mess when we are including subexpressions from other groups in Synapses
). Any good idea? An easy solution would be to substitute all subexpressions, but I remember you arguing against this a while ago :)
Two features for preferences that would be nice would be:
Some issues are:
We should port Brian's mechanism for automatically running benchmarks and come up with a couple of benchmarks, both micro-benchmarks and realistic model simulations (e.g. the classical CUBA, COBA sims). This should include comparing the different code generation targets and maybe comparisons to as a reference.
It would also be good to include memory usage in addition to speed, but I'm not sure how feasible that is.
This issue is not a high priority, though, first we need to make the simulations run, then we can worry about performance :)
See title and issue #49.
[This is a nice-to-have and not high priority but a nice opportunity for external contributions]
Neo (http://neuralensemble.org/neo/) is becoming the standard for the representation of electrophysiological data in the Python/neuroscience community. Recently, the neurotools project (http://neuralensemble.org/NeuroTools/) has been restarted with the goal of bundling existing analysis functions, using neo as the data representation format. I think it would be useful if Brian results (e.g. from a StateMonitor
and SpikeMonitor
) could be easily converted into neo objects to make use of these tools.
neo objects are basically numpy arrays with additional metadata (such as physical units), so generating them shouldn't be too difficult. The conversion would have to include a conversion from Brian's unit system to the unit system provided by quantities (http://pythonhosted.org/quantities/) but that should be quite straightforward.
I guess having this would make the tools.statistic package in brian1 obsolete?
For the standalone framework, if we stick with something like the current system of overriding methods with standalone-generating methods, it would be useful to have an automatic way to take a function argument specifier and turn it into a dictionary. So something like, if you have a function defined as f(x, y, z=1) then the user could write f(1,2,3) or f(y=2,x=1,z=3) or any combination. It would be nice if we could automatically construct a dictionary of the keyword arguments by name, removing the arguments all together, so that they could be treated in the same way. This should be relatively simple to implement but may not be necessary if we switch to the new Device
way of implementing standalone.
Supporting Python 3 is not a high priority at the moment but if it is easy to do, there's no reason not to support it. Possibly, using 2to3 as part of the setup process is already enough:
http://python3porting.com/2to3.html#running-2to3-on-install
If yes, do this and add python 3 to the testing targets.
There are also a couple of uses of e.g. dict.items()
that should be dict.iteritems()
, etc.
The example I-F_curve_LIF
example behaves wrongly now that the refractoriness branch is merged into master.
We could use something like https://github.com/ActiveState/appdirs to get more adequate directory names, in particular for the user-specific configuration file. ~/.brian/...
works but is not the recommended directory on any platform (for Windows one wouldn't use a dot in the name, on Linux it should rather be in .config/brian
, etc.). Similarly, the temporary directory might not be the ideal place for log files (but we should make this configurable anyway).
In Brian 2, we want to have specific stateupdaters for for equations containing xi
variables. We also want to allow more than one noise variable. To make this all work, we need the following things to be done:
Expression.split_stochastic
so that it works with several random variables. Maybe this should move to Equations
as the same variable can be used in several expressions?The old Brian Connection
class had multiple internal data structures possible (dense, sparse, dynamic). I wonder if Synapses could do with something similar? In particular, there might be some data structures that could be more efficient in certain cases?
Pros: possible efficiency gains
Cons: complicates implementation
Thoughts?
We should have methods to save and load network state that are more portable than pickle, for example for use in standalone C++ code. Two options occur to me:
Or we can do both, but then we have effectively three file format (pickle, internal, HDF5) to support, which seems like a bad idea.
We need to make clear which strings are interpreted by sympy and which are not. Most of the time using sympy and executing the string in Python/C gives the same result, but there are important differences, e.g. integer fractions (see #38 ).
In principle we want to use sympy only for "mathematics" (differential equations, state updater descriptions), not for "programming" (reset statements, pre/post codes). On the other hand, using a sympy expression allows for using the CCode printer, which replaces some function calls and **
by pow
. But we also already have a mechanism in place that deals with functions directly.
Proposed procedure:
parse_to_sympy
and sympy_to_str
in brian2.codegen.parsing
.Function
system in code generation (e.g. in C we would rather use #define abs fabs
instead of using sympy to replace abs
calls by fabs
)The only thing I'm not clear about is the translation about **
-- maybe we simply require the use of a power
function? Note that this is not as bad as it sounds, it would not affect equations (there we could even support ^
) but only abstract code where this will not be used often, anyway.
Currently, the same sympy symbols are used for parsing of the state updaters in ExplicitStateUpdater
and the equations. This has the unintended effect that the identifiers f
and g
in equations is interpreted as functions and not as a variables.
More generally, identifiers in the namespace should be parsed according to their role (e.g. if f
refers to a Function
specifier, it should be interpreted as a function, if it refers to a variable, it should be interpreted as a symbol).
We should investigate templating packages and either use one or create a Brian specific one that we use consistently throughout Brian.
This issue supersedes #22.
A new mechanism for doing Synapses and Monitor using CodeObject. We allow resizeable arrays and if an index is accessed that is outside the range of the array, it is extended. In addition, we have conditional indexing, saying that we iterate only over the values where a certain condition holds.
The workaround possiblity becomes unnecessary after fixes in sphinx, see:
https://bitbucket.org/birkenfeld/sphinx/issue/904/autodocs-does-not-work-well-with-instance
When specifying a name, you should be able to write something like myname*
and this is interpreted as: if there is no object myname
then name it myname
, otherwise name it myname_2
, myname_3
, ... until you find one that isn't taken. This is already done with the default name for a class, but it would be nice to make it more widely accessible.
We want to allow to have parameters change between runs, the easiest way to support this is to have codeobjects (state updater, reset, threshold) get regenerated on every run. This is also more robust then to only do changes in the namespace, as a state updater might have replaced constant values with literals already.
I think the cleanest way to do this is to replace the current prepare
method in Network and BrianObject with pre_run
and post_run
. In contrast to the prepare
method, which was called at somewhat unspecified times and could be called multiple times, pre_run
and post_run
should be called exactly once before and after a run.
In NeuronGroup, pre_run
would take care of building the code objects.
I just started work on this in the language_rewrite
branch and realised that it's not entirely trivial because the current ad-hoc templating system has support for some quite nice features like getting the indents right, etc. I wonder if we actually need to rewrite it after all? I think given the deadlines, it would be better not to unless there are specific missing features that we want to see. Are there? Another reason would be to avoid locking ourselves into an API that will make our life more difficult later on, but perhaps this can be put off until after the developer preview / 2.0alpha version?
This will make it integrate better with the Sphinx build system.
Thinking about our "automatic model documentation" and looking a bit into NineML I realized that it would be a nice thing if we could attach textual descriptions to variable names (including parameters and static equations [did we agree on a new name for them? Subexpressions?]). Then I realized that we don't need a new syntax for that, we could simply parse the comments instead of discarding them. I'm thinking along the lines of:
>>> G = NeuronGroup(100, '''
dv/dt = -v / tau : 1 # membrane potential
v_t : 1 # threshold
''', threshold='v > v_t', reset='v=0')
>>> print G.description['v']
membrane potential
>>> print G.description['v_t']
threshold
or maybe, use something like G.specifiers['v'].description
instead -- I don't think this feature has to be particularly accessible, it's mostly for automatic processing.
The only problem I'm seeing is what to do with comments like:
G = NeuronGroup(100, '''
dv/dt = (g_L*(E_L - v) + # leak current
g_e*(E_e - v) + # excitatory input
g_i*(E_i - v)) / tau : 1 # inhibitory input
v_t : 1 # threshold
''')
but maybe one should use static equations/subexpressions in this case, anyway?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.