0

I have a function which should be solved using ode solver. The problems is that it has many constants (say 100) and at the same time I don't want to put them in an array because it would be hard to understand the formulas (which are chemical stuff). Some part of my code in this function looks like:

 # Component: Calcium Fluxes

    J_rel =  v1*(P_O1+P_O2)*(Ca_JSR - Ca_ss)*P_RyR # (micromolar_per_millisecond)
    J_tr = (Ca_NSR - Ca_JSR)/tau_tr #  (micromolar_per_millisecond)
    J_leak =  v2*(Ca_NSR - Ca_i) #  (micromolar_per_millisecond)
    J_up = ( v3*(Ca_i ^ 2.00000))/((K_m_up ^ 2.00000)+(Ca_i ^ 2.00000)) #  (micromolar_per_millisecond)
    J_xfer = (Ca_ss - Ca_i)/tau_xfer #  (micromolar_per_millisecond)
    P_RyR_prime =  - 0.0400000*P_RyR -  (( 0.100000*I_Ca_channels)/i_CaL_max)*(exp(( - ((V_m - 5.00000) ^ 2.00000)/648.000))) #  (micromolar_per_millisecond)

If I wanted to put them in an array it would become something like this:

J_rel =  algebraic[2]*(constants[21]+constants[3])*(constants[4] - Ca_ss)*P_RyR # (micromolar_per_millisecond)
J_tr = (Ca_NSR - Ca_JSR)/tau_tr #  (micromolar_per_millisecond)
J_leak =  constants[31]*(Ca_NSR - Ca_i) #  (micromolar_per_millisecond)
J_up = ( constants[69]*(Ca_i ^ 2.00000))/((constants[17] ^ 2.00000)+(Ca_i ^ 2.00000)) #  (micromolar_per_millisecond)
J_xfer = (Ca_ss - Ca_i)/tau_xfer #  (micromolar_per_millisecond)
P_RyR_prime =  - 0.constants[5]*P_RyR - (( 0.100000*constants[65])/constants[43])*(exp(( - ((V_m - constants[4]) ^ 2.00000)/constants[34]))) #  (micromolar_per_millisecond)

Which is extremely hard to 'see' the actual formulas. Do I have to define the constants as global variables? this functions is used by ode solver and gets iterated about 50000 times. So I need something efficient.

P.S. I'm using python(scipy library), but I think one might have the same problem in other languages too.

Thank you very much

Kamran Bigdely
  • 7,946
  • 18
  • 66
  • 86

3 Answers3

2

On second thought, why not just use global variables? It will give you a pretty, simple, and fast result. In python globals will be restricted to the module they are defined in, so you don't have to worry about them polluting the entire program as they would in say C.

goertzenator
  • 1,960
  • 18
  • 28
  • what about the efficiency? I've heard in MATLAB, each time you call the function, it copies all the global variables to the function. So if you want to run your function for 50000, it would become a problem. – Kamran Bigdely Mar 03 '11 at 22:30
  • Here, the efficiency is important to me not polluting the scope. – Kamran Bigdely Mar 03 '11 at 22:32
  • If efficiency is important, I would just profile your code using the `timeit` module using globals and some other scheme and see which gives you better performance. It's the best way to convince yourself and then you have the hard numbers to make an informed decision. – JoshAdel Mar 04 '11 at 02:08
1

If you want to avoid polluting the namespace you can use a dictionary:

c = { 'Ca_NSR' : 1.5, ... }

J_rel =  v1*(P_O1+P_O2)*(c['Ca_JSR'] - c['Ca_ss'])*P_RyR

You could also use a class:

def Constants(object):
    def __init__(self):
        self.Ca_NSR = 1.5
        # etc...

c = Constants()
J_rel =  v1*(P_O1+P_O2)*(c.Ca_JSR - c.Ca_ss)*P_RyR
Mark Byers
  • 811,555
  • 193
  • 1,581
  • 1,452
  • seems to be the answer but let me to figure out whether the scipy ode solver which is http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode can handle this or not. – Kamran Bigdely Mar 03 '11 at 21:52
1

You can use function keyword arguments to "unpack" a dict into the local namespace:

c = { 'Ca_NSR' : 1.5, ... }

def f(Ca_JSR, Ca_ss, ...):
  J_rel =  v1*(P_O1+P_O2)*(Ca_JSR - Ca_ss)*P_RyR

f(**c)

See the documentation on keyword argument for details on what f(**c) means.

goertzenator
  • 1,960
  • 18
  • 28
  • Well mybe that example helps to explain what's going on: c = { 'CA_NSR' : 1.5,... }; foo(**c) is the same as foo(CA_NSR = c['CA_NSR'],..) – Voo Mar 04 '11 at 00:44