Make sure simulations match Kylie's ~exactly~ well enough

[MichaelClerx]

We want to recreate Kylie's IKr work in Web Lab (like Aidan has done) but using Pints as a back-end. In order to test Pints - before it becomes possible to do so in the web lab - we need to be able to recreate the simulation predictions exactly. Currently there's some small differences. Will need to fix this.

[jonc125]

Or explain how the differences arise - solvers written in different languages, slightly different ODE solver algorithms, slightly different compiler optimisation, etc. all can add up!

[MichaelClerx]

I created some more accurate matlab files of Kylie's stuff, and now getting this for the difference between her protocol and the myokit one: All in the order of 1e-13, so that seems good enough to me. Here's a zoom on the sine wave part: That looks like truncation error (or maybe underflow?) to me!

This is the difference between the currents (kylie's matlab one versus my myokit one) There is still a big-ish error at the point where the sine waves stop. The protocol doesn't have any particular differences at this point, so it must be down to (how we use) the solver.

In the region before the end of the sine waves, the maximum error is 1.5e-6, so pretty close!

[MichaelClerx]

This is the difference in currents again, using abstol=reltol=1e-10 in both Matlab and Myokit

I think this means it's good enough! @jonc125 @mirams ?

To summarise:

• The first step in the protocol is 0.1ms too long (clearly explained in Kylie's m files)
• I had a small Myokit bug where the order of logging / protocol changing events wasn't always consistent (because I had a >= somewhere that should have been a >) resulting in sharp errors at some of the voltage transitions
• At this point, the errors got small enough that the default matlab output didn't have enough decimals, so I modified Kylie's code to add -double and get full precision decimals
• There was still an error after the sine-wave region. At this point the protocols were 'identical', differing only due to some form of numerical/truncation error, so the difference had to lie in the solver. Since we use the same solver, the difference had to be in the way we use it.
• Reducing the tolerances reduced the error in simulated currents to O(10^-8). This means that the differences we see are due to small differences in code, and not because of some fundamental difference in the mathematics of the problem.

[MichaelClerx]

When I evaluate the log likelihoods I get:

Myokit signal comparison:
  Log-likelihood: -1.51042757394234696e+06
Matlab signal comparison:
  Log-likelihood: -1.51042762220851122e+06

So that seems pretty good!

The value reported by Kylie's matlab code is slighty different though:

My value for Kylie's signal:
  -1.51042762220851122e+06
The value in her stored matlab files
  -1.51063175578740e+006

Not really sure how to find the source of this mismatch without digging deeply into Kylie's code. Don't think it's very important though (probably some truncation in storing the value in a file somewhere) so will accept the above as good enough!

[mirams]

Yep looks as if it will be good enough. Next step is probably just doing the CMA-ES and checking it gets there too, and then repeating the MCMC and comparing with both Kylie and Aidan's posteriors.

[jonc125]

If you've changed the tolerances then I'd expect a different value from her Matlab code?

But yes, this is probably fine!

[mirams]

Good detective work by the way!

[MichaelClerx]

Thanks!