dictionary support for de.jit

[sibyjackgrove]

First of all, thank you very much for 6770470 and fc32432. I am trying to solve multiple ODEs at once by concatenating them using the below code. The code without de.jit works fine.

func = de.ODEFunction(combined_model)
prob = de.ODEProblem(func, y0, (t0, dt),pvder_dict)
prob_jit = de.jit(prob ) #Error occurs if I use de.jit

def combined_model(dy,y,p,t):
	"""Concatenate ODE residuals from multiple models of same type"""

	i = 0
	pvder_dict = p
	for dss_id in pvder_dict:
		for node in pvder_dict[dss_id]:			
			for der_id in pvder_dict[dss_id][node]:
				sim = pvder_dict[dss_id][node][der_id]["sim"]
				nEqs = sim.DER_model.n_ODE
				start_index = i * nEqs
				end_index = (i + 1) * nEqs
				dy[start_index:end_index] = sim.ODE_model(y[start_index:end_index],t)
				i += 1
			
	return dy

However, if I use de.jit I am getting the following error:
@ChrisRackauckas @LilithHafner could you suggest something that I could try out? For individual ODE, I can get up to 10 times speed up with de.jit. However, my application requires solving hundreds of ODE models simultaneously.

  File "/home/tdcosim_cloud/tdcosim_cloud/model/pvder/odesolverutilities.py", line 154, in get_diffeqpy_integrator
    diffeqpy_ode = ode.jit(diffeqpy_ode)
                   ^^^^^^^^^^^^^^^^^^^^^
  File "/root/.julia/packages/PythonCall/qTEA1/src/jlwrap/any.jl", line 208, in __call__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_call)), args, kwargs)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/tdcosim_cloud/tdcosim_cloud/model/pvder/odesolverutilities.py", line 103, in funcval_diffeqpy
    for der_id in pvder_dict[dss_id][node]:
                  ~~~~~~~~~~~~~~~~~~^^^^^^
  File "/root/.julia/packages/PythonCall/qTEA1/src/jlwrap/any.jl", line 214, in __getitem__
    return self._jl_callmethod($(pyjl_methodnum(pyjlany_getitem)), k)
       ^^^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: MethodError: no method matching getindex(::Symbolics.Num, ::Symbolics.Num)

Closest candidates are:
  getindex(::Number)
   @ Base number.jl:95
  getindex(::Union{AbstractChar, Number}, !Matched::CartesianIndex{0})
   @ Base multidimensional.jl:853
  getindex(::Number, !Matched::Integer)
   @ Base number.jl:96

[LilithHafner]

I think the problem here is that de.jit (aka ModelingToolkit.modelingtoolkitize) assumes that the parameters are numeric things, and you are passing in nested dictionaries of ODE models, something that a derivative can't be computed with respect to.

Specifically, I believe that when you perform a getindex operation pvder_dict[dss_id], it thinks pvder_dict is numeric and throws.

As a workaround, perhaps don't pass the models in as parameters, instead embed them in the function. Does this work for you?

def make_combined_function(pvder_dict):
    def combined_model(dy,y,p,t):
        """Concatenate ODE residuals from multiple models of same type"""

        i = 0
        for dss_id in pvder_dict:
            for node in pvder_dict[dss_id]:			
                for der_id in pvder_dict[dss_id][node]:
                    sim = pvder_dict[dss_id][node][der_id]["sim"]
                    nEqs = sim.DER_model.n_ODE
                    start_index = i * nEqs
                    end_index = (i + 1) * nEqs
                    dy[start_index:end_index] = sim.ODE_model(y[start_index:end_index],t)
                    i += 1
                
        return dy
    return combined_model

func = de.ODEFunction(make_combined_function(pvder_dict))
prob = de.ODEProblem(func, y0, (t0, dt), None)
prob_jit = de.jit(prob)

[sibyjackgrove]

@LilithHafner Yes, the method you suggested works. Thanks a lot! However, I am getting another error from within the individual model. What could be the issue?

File "/usr/local/lib/python3.11/site-packages/pvder/DER_components.py", line 370, in Ppv_calc
    self.Ipv = (self.Np*self.Iph)-(self.Np*self.Irs*(math.exp((self.q*Vdc_actual)/(self.k*self.Tactual*self.A*self.Ns))-1))  #Getting error here
                                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/root/.julia/packages/PythonCall/qTEA1/src/jlwrap/number.jl", line 168, in __float__
    return self._jl_callmethod($(pyjl_methodnum(pyfloat)))
       ^^^^^^^^^^^^^^^^^^^^^^^^
juliacall.JuliaError: MethodError: no method matching Float64(::Symbolics.Num)

Closest candidates are:
  (::Type{T})(::Real, !Matched::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(!Matched::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50

 ...

[LilithHafner]

I think you are running into the limitations of JuliaCall/ModelingToolkit compatability. de.jit is a pretty thin wrapper around that ModelingToolkit.modelingtoolkitize.

Specifically, I'm guessing that python's Math.exp is converted into a not-completely-general Julia function that must convert its inputs to concrete Float64s. If my guess is correct, then this is probably best solved either in JuliaCall, or if that's not possible in a JuliaCall/Symbolics package extension. You could work around this particular issue by not calling math.exp, but I don't know what other Python functions fail to convert to a Symbolics.jl friendly form.

[sibyjackgrove]

The attribute Vdc_actual is obtained by multiplying one of the states with a constant. When I print it out within the method, I get this: 500.0x₁₉(t) . Could this be the problem?

[sibyjackgrove]

I think you are running into the limitations of JuliaCall/ModelingToolkit compatability. de.jit is a pretty thin wrapper around that ModelingToolkit.modelingtoolkitize.

Specifically, I'm guessing that python's Math.exp is converted into a not-completely-general Julia function that must convert its inputs to concrete Float64s. If my guess is correct, then this is probably best solved either in JuliaCall, or if that's not possible in a JuliaCall/Symbolics package extension. You could work around this particular issue by not calling math.exp, but I don't know what other Python functions fail to convert to a Symbolics.jl friendly form.

YOu are right. math.exp is the problem. Is there any other option that I can replace it with?

[LilithHafner]

I'm not sure what there is to be done on your end other than continuing to use the old version of diffeqpy if you depended on numba integration for jit-ing your ODEs. I'll look into fixing it on our end though.

[sibyjackgrove]

Ok, I understand. I guess I just have to not use de.jit for now. Appreciate your help.