[alexfikl]

This has two performance improvements:

1. Uses Horner to evaluate some of the cubic polynomials. As far as I can tell, this has a very negligible effect (likely because -O3 does it anyway).
2. Passes anchors around by pointer instead of by value. This seems to have a noticeable effect for larger systems, because it avoids a bunch of copies.

I've tested the changes on the examples/kuramoto_network.py example (using hyperfine to do full timing) and using this benchmark (to time separate parts). The results can be seen below.

This was run on my machine (AMD Ryzen AI 9 HX 370) and Python 3.14 using GCC and SymEngine 0.14.0. The huge amount of time is spent in the compilation step, but not sure we can do much there (between symengine and gcc doing work). The improvement I was hoping for is in the run part, which seems pretty nice 😁

What do you think? Also checked that there are no new leaks reported by valgrind.

Before

> hyperfine 'python -O examples/kuramoto_network.py' 
Benchmark 1: python -O examples/kuramoto_network.py
  Time (mean ± σ):      2.774 s ±  0.036 s    [User: 2.765 s, System: 0.121 s]
  Range (min … max):    2.725 s …  2.827 s    10 runs

> python -O benchmark_kuramoto.py --reps 10 --n 200    
Results (omp=off)  (n=200, 10 reps, mean ± stdev)
  setup        1.1 ms  ±  1.3 ms
  compile      4.806 s  ±  228.4 ms
  transient    2.5 ms  ±  1.8 ms
  run          201.5 ms  ±  84.6 ms
  total        5.011 s  ±  184.9 ms

After

> hyperfine 'python -O examples/kuramoto_network.py'
Benchmark 1: python examples/kuramoto_network.py
  Time (mean ± σ):      1.679 s ±  0.067 s    [User: 1.568 s, System: 0.109 s]
  Range (min … max):    1.626 s …  1.858 s    10 runs

> python -O benchmark_kuramoto.py --reps 10 --n 200    
Results (n=200, 10 reps, mean ± stdev)
  setup        1.1 ms  ±  1.3 ms
  compile      4.679 s  ±  178.8 ms
  transient    1.8 ms  ±  1.4 ms
  run          123.2 ms  ±  52.1 ms
  total        4.805 s  ±  150.8 ms

I'm not sure why the hyperfine run was so much faster after, since it really should be dominated by the compilation part.

Wilson-Cowan

I've also tested this on my own code (a similarly sized Wilson-Cowan system with fewer delays) and it seems to improve the situation quite a lot. When running a Lyapunov exponent calculation, I'm seeing

# before 
eps = 0.000 tau = 1.610 behavior 'Periodic' (wall 50.14251s)
# after
eps = 0.000 tau = 1.610 behavior 'Periodic' (wall 5.90889s)

This is a bit dramatic and I'm not quite sure why it improves it so much, but it seems repeatable with the changes on this branch. I can make an smaller script for it, if you think it's useful.

[Wrzlprmft]

1. Uses Horner to evaluate some of the cubic polynomials. As far as I can tell, this has a very negligible effect (likely because -O3 does it anyway).

If I am not mistaken, this is not a beneficial change.

The Horner scheme is best if you have the coefficients of the polynomial in its canonical form. But that’s not what we have here. Instead we have states and derivative at certain points, from which you first need to compute the coefficients of the polynomial, which itself takes more time than evaluating the polynomial itself. The reason why this evaluation of the polynomial is so unusual is that it was the best I could find for this case.

The best way to compare performances here is to look at the total number of additions (including subtractions) and multiplications and here my code beats yours. For example for getting the state, I have seven additions and multiplications each, while you have nine additions and eight multiplications. I am not super sure whether multiplications by 2 are somewhat cheaper, so you might redeem a bit of performance here, but I don’t see it being enough.

Now, where you might save some runtime in evaluation is if you store the coefficients of the polynomial for a pair of anchors that will be re-used later. However, I would estimate that it won’t make up for the runtime lost by administrating this storage, let alone the pains to implement something like this.

2. Passes anchors around by pointer instead of by value. This seems to have a noticeable effect for larger systems, because it avoids a bunch of copies.

I haven’t tried to understand the benchmarks in detail yet and haven’t checked all the changes, but I expect this to be a purely beneficial change. At worst, it doesn’t do anything.

I guess that I didn’t implement it like this from the start, because I hoped the compiler would take care of this, as anchors are rather lightweight structs. Seems like I was wrong about that.

[alexfikl]

If I am not mistaken, this is not a beneficial change.

Yeah, agreed, I'll remove that part. I haven't managed to measure any impact that didn't look like noise there.

I guess that I didn’t implement it like this from the start, because I hoped the compiler would take care of this, as anchors are rather lightweight structs. Seems like I was wrong about that.

IIUC, the speedup from this should only show up for larger systems due to the state and diff members in anchor

typedef struct anchor
{
	double time;
	double state[{{n}}];
	double diff[{{n}}];
	struct anchor * next;
	struct anchor * previous;
	...
} anchor;

I've very briefly tried it with the mackey_glass.py example and there's pretty much no change at all there. I'm guessing I'm seeing an even bigger improvement with my Lyapunov run because the system is even larger there.