Plotting of arbitrary functions

Support
1

Hi!

I want to plot some functions that I defined for my models.

I did the following but there is probably a better way:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
from brian2 import run, NeuronGroup, StateMonitor, defaultclock
from brian2 import second

group = NeuronGroup(1, "x = (t/second - 100)**2 : 1", method='euler')
statemon = StateMonitor(group, ['x'], record=True)
defaultclock.dt = 0.1*second
run(200*second)

plt.plot(statemon.t/second-100, statemon.x[0])
plt.xlabel('x')
plt.ylabel('x*x')
plt.savefig("brian_plotting.png")

brian_plotting

Can you please point me to the correct way to plot equations used in Brian2?

Thanks in advance,

Sebastian

2

Hi. I agree that this is something that we should support in a more convenient way… Here’s what I usually do: I set the function argument (often the membrane potential v) to the full range of values in a NeuronGroup, and then directly plot the variable that is described by an equation. This has the advantage of not having to change the equations and is also not actually running a simulation (the expression gets directly evaluated by numpy). It is still not a perfect solution, since you usually have to create a NeuronGroup just for this purpose (if the existing NeuronGroup does not happen to have exactly the size you want for the resolution). In your example, x would probably be a differential equation, so it would look something like this:

import matplotlib.pyplot as plt
import numpy as np
from brian2 import NeuronGroup

group = NeuronGroup(100,  # the resolution I want
                    """dx/dt = -x / (10*ms) : 1  # this would be v typically
                       y = x**2 : 1  # the function I'm interested in""",
                    method='euler')
group.x = np.linspace(-100, 100, len(group))

plt.plot(group.x, group.y)
plt.xlabel('x')
plt.ylabel('x*x')
plt.savefig("brian_plotting.png")

This will give you the same plot.

2 Likes
3

Hi!

Thanks for that alternative. Without running the simulation it seems that variable are not replaced, e.g., say that y = p*x**2 : 1 with p a python variable. That would work in a full run but not pure numpy. Is there are workaround?

Thanks,

Sebastian

4

Ah, no, it does work in principle but there is a surprisingly complex issue behind this problem that I’ve never managed to address in any satisfying way. Without going too much into detail, what group.y returns is not just an array with the evaluated values, but a VariableView. This adds functionality, e.g. you can use Brian’s string-based indexing, say group.y['y > 10']. The problem is that when this is finally asked to return all the values in the depths of matplotlib’s plotting code, it does not know where to find the value for p any more… Long story short, there’s an easy solution (and in some situations, e.g. when you try print(group.y) we directly mention this workaround): use group.y[:] (to say, “evaluate me now”) instead of group.y in the plot command.
Alternatively, you could provide p to the NeuronGroup via its namespace argument (i.e. namespace={'p': p}) instead of relying on Brian’s “magic” to find it in the surrounding code.

1 Like
5

Thank you! Is it possible to disable all magic?

1 Like
6

If you use the namespace arguments, it will not look up things in the surrounding Python namespace. If you do this on the level of the NeuronGroup (which is necessary if you use the above “trick” to plot functions), then using namespace={} for the run-call makes sure that nothing is looked up elsewhere. See Namespaces — Brian 2 2.4.2 documentation for more details.

7

Thanks. I’ve put together a short doc for it: How to plot model functions by schmitts · Pull Request #1308 · brian-team/brian2 · GitHub

8

Thanks, that’s great! :+1:

9

@mstimberg There seems to be some overhead though, when not using the [:] approach. In the example below my user defined function gets called 4 times. Any idea why?

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

from brian2 import check_units, NeuronGroup

@check_units(x=1, result=1)
def heaviside(x):
    print("heaviside called", x)
    return (x > 0).astype(int)

group = NeuronGroup(10,
                    """dc/dt = -c/(10*ms) : 1
                    f = heaviside(c - 0.3) : 1
                    """,
                    method='euler',
                    namespace={'heaviside' : heaviside})

group.c = np.linspace(-1, 1, len(group))

#plt.plot(group.c, group.f[:]) # heaviside is called once
plt.plot(group.c, group.f) # heaviside is called 4 times
10

This is due to group.f being a VariableView. As the name suggests, this is not an array of values, but a view on the variable (your equation f is also a Variable in Brian’s terminology here). As soon as you use this view as if it were an array (indexing it, calculating with it, converting it with something like np.array, …), it will actually evaluate the expression. So if you pass this object into a function like plt.plot, it is possible that it is asked for its values several times, triggering an evaluation each time. We cannot really avoid this without introducing some error-prone caching.
All that said, I think most of calls by matplotlib are acutally related to calling group.f.shape, which internally is basically implemented as (group.f[:]).shape. We should be able to handle this more elegantly.

1 Like