Port key matlab script with self-contained figure to python

[slarson]

The matlab script here can generate the key voltage clamp figure of the Boyle & Cohen paper. This can be used as a standard to measure the rest of the implementation against very nicely.

[rayner]

I've made a Python version of the script and created a pull request for it.

I didn't worry about floating-point rounding errors when porting the code, but I did notice that Matlab (or Octave, at least) seems to handle floating-point differently to Python in terms of where rounding errors popped up. Is it worth going through the code and making them as similar as possible?

[Neurophile]

Python Decimal objects can be used to control precision. It's a bit of a pain to set up, but it provides transparency and repeatability for the math. It's a built-in library so no worries about adding dependencies.

https://docs.python.org/2/library/decimal.html

[rayner]

Using Decimal objects throughout looks like a good idea - we'd have a known precision with no floating-point errors. However, it would also eliminate the floating-point errors that the original version of the sim had. Should our goal be to remove this source of error or to replicate the original simulation exactly?

(I'm leaning towards "remove this source of error", since it presumably wasn't an intended feature of the original sim in the first place, and it might cause slight differences in the output depending on what platform the code runs on)

[rayner]

I've updated the script to use Decimals. It does run 100x slower now, taking 86 seconds instead of 0.8 on my machine, but that might not be a problem if we're only expecting people to run it infrequently (e.g. to generate a graph that they can then save for later comparison).

Python 3.3 has a much faster implementation of the decimal module, so that might be a way around the slowness if someone needs to run the script frequently.

[gidili]

@rayner I think there is value in keeping the numeric "error" as in the original publication, to show exact same results are produced - maybe we can keep it in a branch?

[rayner]

@gidili That might make sense, yes. I'll take a look at the output of each version and compare it with the output of the original Matlab script, and figure out how hard it would be to make them match exactly.

For the purposes of replicating Figure 2D from the paper, the difference is small enough that it's not visible just from looking at the generated graph.

[rayner]

I've verified that the actual plotted values (i.e. the contents of I_mem) from the Python/Decimal version of the script are the same as those from the Matlab script, to at least 5 significant figures (the default number displayed by Octave).

Given that the simulation is working at about this precision anyway (the input variables from input.csv are all between 4 and 6 significant figures), I think we should be okay with just the current version of the Python script.

[pgleeson]

I'd recommend keeping both versions, maybe having vclamp_D.py etc. for the version(s) of the Matlab scripts using Decimal.

The original script will be useful for testing & rerunning with altered parameters.

Also, maybe use

f = open('../../MatlabSupport/Main_Version/data/input.csv')

as opposed to making a copy?

[pgleeson]

Also, would it be possible to extract out the

Cmem = float(data[0])
...
alphaCa = float(data[27])

block into a separate python script? This could be reused then for other scripts converted to Python, and to convert these values to the corresponding dimensional values required for the NeuroML2 files...

[rayner]

Okay, all done.

I've kept the float versions of the scripts as the defaults, for compatibility with the original Matlab results. Decimal versions of the scripts have names ending "_decimal.py".

The data import section is now in "input_vars.py", which can be imported as a module (using from input_vars import * as a drop-in replacement for the original code, or using import input_vars to avoid clashing with any other things that might share those variable names).

The Python scripts now use the same input.csv file as the Matlab scripts.

The pull request is https://github.com/openworm/muscle_model/pull/14. I'm happy to go ahead and merge it if nobody has any objections.

[pgleeson]

Thanks for that @rayner. I've merged the pull request.

[pgleeson]

It would be great to get another Python script to print out the corresponding values in the python/matlab impls vs the values used for NeuroML2. I've started a script for this here: https://github.com/openworm/muscle_model/blob/master/BoyleCohen2008/PythonSupport/Main_Version/compareToNeuroML2.py and it would be great if someone could look into this a bit more. Eventually the nml2 can be updated to make sure all the values match.

[slarson]

Closing -- moved follow on issue to #28