Closed MichaelClerx closed 5 years ago
One potential issue to avoid is how we handle log(x)
, log(x, base)
, and log10
and log2
.
CellML supports log10 via log(x, 10)
, but this can cause issues with code that originally relied on the added accuracy of the specialised log10
method (esp. for small numbers).
This is from an earlier email:
We got a report that the Shannon 2004 model (I presume the second version, with the corrections) didn't run when using the exported Python code from PMR. However, if the two calls to
log(x, 10)
were replaced bylog10(x)
the code ran fine.I don't have the specific example, but I've checked and the equations:
pCa_SL = -log10(Ca_SL / 1) + 3
and
pCa_SL = -log(Ca_SL / 1) + 3
do give different results with the initial value
Ca_SL = 1.03181200000000001e-04
(I'm guessing it's not this particular initial value that causes the error though, as the difference is rather small!)
The issue is known to the Python devs: https://bugs.python.org/issue6765
Reason seems to be that the general purpose
log(x, y)
divides two logs to be generic, whilelog10(x)
has its own algorithm (it maps tolog10
in C, no idea if that's a processor instruction or it uses a custom taylor series or something).So it's really more of an `issue', as the current implementation is technically correct :-D
However, if would be good if the export could treat
log(x, 10)
as a special case and output the log10 method. (Actually, in MathMLlog(x, y)
is a special case oflog10(x)
with an added logbase element, so it'd be sticking a bit closer to the original cellml code).Does that make sense?
The referenced model is http://models.cellml.org/exposure/d72a36fe0b7e121068c96bcb1ff6044a/shannon_wang_puglisi_weber_bers_2004_b.cellml/view
At first glance, it looks like sympy handles this intelligently:
>>> from sympy import *
>>> import math
>>> log(1000, 10)
3
>>> log(1000, 10, evaluate=False)
log(1000, 10)
>>> log(1000) / log(10)
log(1000)/log(10)
>>> (log(1000) / log(10)).evalf()
3.00000000000000
>>> log(1000, 10) == (log(1000) / log(10)).evalf()
True
>>> math.log10(1000) == (log(1000, 10)).evalf()
True
For the specific example in the email:
>>> from sympy import *
>>> import math
>>> x = 1.03181200000000001e-04
>>> y1 = -math.log10(x / 1) + 3
>>> y2 = -(math.log(x / 1) / math.log(10)) + 3
>>> y1 == y2
False
>>> x = Symbol('x')
>>> y3 = -log(x / 1, 10) + 3
>>> y4 = -(log(x / 1) / log(10)) + 3
>>> y3 == y4
True
>>> y1 == (y4.subs({x:1.03181200000000001e-04})).evalf()
True
Cool!! But let's leave this open as a check to run after code generation
Adding some links here to existing Web Lab code & docs in this area:
_split_all_odes
does the ODE independent variable units conversion handlingModelModifier
class._classify_variables
, _order_variables
for analysing the equation graph prior to code generation.run()
method shows the steps in the process. CellMLToCythonTranslator
is probably the closest to what we'll want to do.
The best way to test this would probably be: