run_once or run_at method

[thesamovar]

We have a run_regularly method and network_operation, but people often want to have something run just once and it's a bit of a fiddle at the moment. How about we have a Group.run_once or to set something to happen just a few times Group.run_at. We could also similarly have a version of network_operation that only runs once. There are workarounds now, but it seems to be a pretty common use case.

[TheSquake]

Hello,

I am interested in working on this 😃 . I saw brian2 when I was working on my dissertation, where I was making my own rather simple simulator. I am trying to get into computational neuroscience masters, and I think contributing to this project would be a great way to learn.
I've been playing around with brian2 now and I think I have found a pretty simple way to make that work.

import brian2 as brian
import matplotlib as mpl


# this could be a special code runner added to brian2.groups,
#   or added to the regular CodeRunner
class MyCodeRunner(brian.CodeRunner):

    def __init__(self, activate_at, *args, **kwargs):
        self.activate_at = activate_at
        brian.CodeRunner.__init__(self, *args, **kwargs)

    def run(self):
        if (self.clock.t in self.activate_at):
            super().run()


# this would be added as a method to the Group class, with some tweaks
def run_at(group: brian.Group, times_to_run: list[int], code):
    runner = MyCodeRunner(
        times_to_run, G, 'stateupdate', code=code,
        name="regular_update_test", dt=None,
        when="start", order=0,
        codeobj_class=None
    )

    group.contained_objects.append(runner)


tau = 10*brian.ms
eqs = '''
dv/dt = (1-v)/tau : 1
'''

N = brian.Network()
G = brian.NeuronGroup(1, eqs)
M = brian.StateMonitor(G, 'v', record=True)

run_at(G, [1*tau, 7*tau], "v = v + 5")

N.add(G)
N.add(M)

N.run(100*brian.ms)

and this produces the following result

let me know if it's fine to work on this, and if I am going in the right direction with my solution 😁

[mstimberg]

Hi @TheSquake . Many thanks for your interest in working on this, it would be much appreciated 😄

Also thanks for your proof of concept. Your approach would work (at least in the context of the default "runtime mode" – for standalone mode, no Python code will be executed during the run). I think it is not the optimal approach, though: the problem is that it keeps a clock that makes the object execute the run method every time step, before checking whether it should run at the respective time step. For a long list of times, this check will not be very efficient, but maybe more importantly, it might also do unexpected things due to floating point issues. E.g, here's a minimal version of your approach that appears to work:

times = [0.5*ms]
for timestep in range(10):
    t = timestep * 0.1*ms
    if t in times:
        print(f"Executing at t={t!s}")

(Prints "Executing at t=0.5 ms")
But for times = [0.3*ms] it will not print anything. This is due to floating point issues (0.3 ≠ 3 * 0.1) and could be dealt with, but things can become messy quite quickly 😬 Either way, if we want to make it work with standalone mode (and we do 😉 ), we will have to use a different approach.

I think a more promising and not that complicated approach would be to introduce a new type of clock, an EventClock. This would also come with the advantage that it could be used for other things that aren't run_at operations, e.g. I could image it being useful for a StateMonitor to record at pre-defined times only. Let me explain this a bit (in a separate comment)…

[mstimberg]

Here's a simplified version of what a network does during a run to support multiple clocks and multiple runs. This code does not actually simulate anything, it just demonstrates the decisions of how to advance time and how to decide which objects should be simulated next:

class ClockNetwork:
    def __init__(self, clocks):
        self._clocks = clocks
        self.t = 0 * second

    def _nextclocks(self):
        clocks_times_dt = [(c, c.t, c.dt) for c in self._clocks]
        minclock, min_time, minclock_dt = min(clocks_times_dt, key=lambda k: k[1])
        curclocks = {
            clock
            for clock, time, dt in clocks_times_dt
            if (
                time == min_time
                or abs(time - min_time) / min(minclock_dt, dt) < Clock.epsilon_dt
            )
        }
        return minclock, curclocks

    def run(self, runtime):
        print(f"Starting a simulation run at t={self.t!s} for {runtime!s}")
        for clock in self._clocks:
            clock.set_interval(self.t, self.t + runtime)

        clock, curclocks = self._nextclocks()
        running = clock.timestep < clock._i_end
        while running:
            # Run objects for current clock
            print(
                f"t={clock.t[:]!s}, running objects for: {', '.join(sorted(c.name for c in curclocks))}"
            )
            # Advance clocks
            for c in curclocks:
                # Since these attributes are read-only for normal users, we have to
                # update them via the variables object directly
                c.variables["timestep"].set_value(c.timestep + 1)
                c.variables["t"].set_value(c.timestep * c.dt)

            clock, curclocks = self._nextclocks()
            running = clock.timestep < clock._i_end
        print("End of simulation")
        self.t = clock.t[:]

For a single clock, all this is unnecessarily complex, but if you have some objects that should be executed every 0.1ms, and some that should be executed every 0.2ms, it will do the expected thing:

>>> net = ClockNetwork(
        [Clock(dt=0.1 * ms, name="fast"), Clock(dt=0.2 * ms, name="slow")]
    )
>>> net.run(1 * ms)
Starting a simulation run at t=0. s for 1. ms
t=0. s, running objects for: fast, slow
t=100. us, running objects for: fast
t=200. us, running objects for: fast, slow
t=300. us, running objects for: fast
t=0.4 ms, running objects for: fast, slow
t=0.5 ms, running objects for: fast
t=0.6 ms, running objects for: fast, slow
t=0.7 ms, running objects for: fast
t=0.8 ms, running objects for: fast, slow
t=0.9 ms, running objects for: fast
End of simulation

The situation above could also be handled with something that is closer to your approach, i.e. by determining the faster clock and advancing it, and checking each timestep whether the other clock should be considered as well. Brian's approach is more general, though, it also works if the dt values of the two clocks are not multiples of each other:

>>> net = ClockNetwork(
        [Clock(dt=0.1 * ms, name="fast"), Clock(dt=0.11 * ms, name="less_fast")]
    )
>>> net.run(.5 * ms)
Starting a simulation run at t=0. s for 0.5 ms
t=0. s, running objects for: fast, less_fast
t=100. us, running objects for: fast
t=110. us, running objects for: less_fast
t=200. us, running objects for: fast
t=220. us, running objects for: less_fast
t=300. us, running objects for: fast
t=0.33 ms, running objects for: less_fast
t=0.4 ms, running objects for: fast
t=0.44 ms, running objects for: less_fast
End of simulation

Now my idea of an EventClock would be that it would be handled in the _nextclock comparisons in the same way as the other clocks, and it would also internally advance an integer "timestep". In contrast to a normal clock, though, it would not calculate its time t as t = timestep * dt, but instead as t = event_times[timestep].

Does that make sense to you and would you be interested in working on that? I can give you more details on the specific changes we'd need to support all that (and, obviously, don't hesitate to ask if anything is unclear). Oh, and don't worry about the standalone mode, I can take care of that part later.

[TheSquake]

Hey @mstimberg

Thank you for this comment, yes I think making this a new clock would be much better. I think I understand the idea and while making my first solution I was wondering if this should be a new clock. Of course I am still interested in doing this 😀 I have started working on that today, will comment again once I get version of that working. I guess my only question is whether that clock should be a subclass of the regular Clock and then what should be the dt of that clock since it doesn't get fired regularly. That is used in nextclocks function and maybe in some other places in the actual run function of a Network, or other parts of the codebase (sorry didn't check this yet)

[mstimberg]

I think the cleanest solution would be to have a new BaseClock with everything that is common across all clocks, and make Clock and EventClock inherit from it. As you said, a dt does not really make sense for EventClock, so it would be a bit of a hassle to have EventClock inherit directly from Clock. The main things that need to be adapted are indeed in the Network – basically everything that assumes that a clock has a dt needs to go and be replaced by logic in the Clock/EventClock classes themselves. Currently, the Networks's run function advances the clocks by increasing their timestep and calculating the new t, but this should instead be done by a BaseClock.advance() method. For the comparisons in nextclock, they could instead be handled by implementing the comparison operators in the clock (__lt__, etc.) and a function that checks for a clock whether another time is "close enough" to be considered equals, so that the Network does not have to order clocks itself.

PS: There is also SchedulingSummary that makes direct references to clock attributes, but this is a low priority and will be simple to fix.

[TheSquake]

yep seems good, I was wondering what about making EventClock the BaseClock maybe under a slightly different name for the regular Clock, then for the Clock it would go over the event_times which would simply be the subsequent times of that clock (if start time is 0 end time 10 and dt 1 then it would go over [0,1,2,3,4,5,6,7,8,9,10]). It's like those event_times would be the hours on a real clock, t would be the hand, and advance would be the turning of that hand. Does that seem like a reasonable idea ? Would that perhaps require too many changes ? Or maybe that would create a sequence that is too large for some simulations ?

[mstimberg]

Sorry for the very late reply – end-of-year illness and end-of-year holidays got into the way… I see your point, it would probably make sense to unify things in this way. It's probably slightly less efficient (not even sure), but advancing the clock is negligible compared to everything that happens within the time step. Actually, now that I am thinking about it, we should probably simply wrap this into some kind of on-the-fly evaluation (creating a huge array for long simulations is probably a bit of a waste). Something along these lines:

class ClockArray:
  def __init__(self, clock):
    self.clock = clock
  def __getitem__(self, item):
    return item * self.clock.dt

The standard clock would then use self.times = ClockArray(self) as the "array" that stores the times, whereas the EventClock would use self.times = event_times – both would advance their timestep attribute and would set self.t = self.times[timestep] after each advance.
Does that make sense?