Dear forum users,
I was trying to simulate sensitivities of my model and got the attached Error_Sens_Mem.txt. I suppose that ida_mem refers to memory and therefore I've attached the sizes of the model in sizes.txt. I have only declared one parameter as free so I guess the total DAE system to simulate has the size N_eq_F0*(p+1) = 436.
I wonder if my model is too large for sensitivity simulations or if there is a solution to this?
Best regards, Niklas
Hello Niklas,
In the call to simulate, do you specify the initialization to false?
Best
/Christian
Yes, I did. When I specify the initialization to true and simulate, the whole python session crashes and a window appears "Python has stopped working. Windows tries to find a solution".
I think I managed to find the last evaluated code and it seems to be in function make_consistent, in file implicit_ode.py.
733 if method == 'IDA_YA_YDP_INIT':
734 --> [flag, y, yd] = self.Integrator.calc_IC(method,self.tout1,self.lsoff)
735 elif method == 'IDA_Y_INIT':
736 [flag, y, yd] = self.Integrator.calc_IC(method,self.tout1,self.lsoff)
Best regards, Niklas
It seems that this is dependent of tolerances. If I have too high tolerances (atol=rtol>1e-6), python crashes. Maybe this is a known problem but I managed to reproduce it in a smaller attached SEIRS model, see below (It was not possible to upload files).
Back to the big model: If I have too low tolerances (atol=rtol=1e-9), the solver returns
"Sundials_Exception: Calculation of initial conditions failed. IDA returned flag -4
WARNING: Failure executing file: "
If I set atol=rtol=1e-8 then I get another return flag:
"Sundials_Exception: Calculation of initial conditions failed. IDA returned flag -13
WARNING: Failure executing file: "
If I set atol=2e-7, rtol=1e-7 I get a "never ending" simulation, because the inputs I'm simulating are very noisy.
So I wonder if there may be a solution to the python crashes so I can solve my problem with reasonable tolerances?
Best regards,
Niklas
============ SEIRS.py ==============
import os as O
import numpy as N
import matplotlib.pyplot as plt
import nose
from jmodelica.jmi import compile_jmu
from jmodelica.jmi import JMUModel
from jmodelica.core import TrajectoryLinearInterpolation
from scipy.optimize import *
jmu_name = compile_jmu("SEIRS", "SEIRS.mop",compiler_options={"enable_variable_scaling":True})
# Load a model instance into Python
print('loading model')
model = JMUModel(jmu_name)
# Get and set the options
opts = model.simulate_options()
opts['IDA_options']['atol'] = 1.0e-3
opts['IDA_options']['sensitivity'] = True
opts['ncp'] = 400
#Simulate
print('simulating model')
res = model.simulate(2002,2007, options=opts)
============ SEIRS.mop =================
optimization SEIRS
parameter Real p1(free=true,nominal=2.78e5) = 2.78e5; //S0
parameter Real p2(free=true) = 1.08e-1; //E0
parameter Real p3(free=true) = 1.89e-1; //I0
parameter Real p4(free=true,nominal=1.00e6) = 1.00e6; //N
parameter Real p5(free=true,nominal=5) = 5.00; //L
parameter Real p6(free=true,nominal=10e-3) = 9.59e-3; //D
parameter Real p7(free=true,nominal=5e-3) = 5.48e-3; //M
parameter Real p8(free=true,nominal=75) = 75; //P
parameter Real p9(free=true,nominal=375) = 375; //beta0
parameter Real p10(free=true,nominal=2.00e-2) = 2.00e-2;//a1
parameter Real p11(free=true,nominal=2.00e-2) = -2.00e-2;//b1
parameter Real R0 = p4-p3-p2-p1;
constant Real pi = 3.14159265;
Real S(start=p1,fixed=true,nominal=1e5);
Real E(start=p2,fixed=true,nominal=1e3);
Real I(start=p3,fixed=true,nominal=1e3);
Real R(start=R0,fixed=true,nominal=7e5);
Real Beta(start=382.5,nominal=100);
Real N(start = 1e6,nominal=1e6);
equation
Beta = p9*(1 + p10*cos(2*pi*time) + p11*sin(2*pi*time));
N = S + E + I + R;
der(S) = (1/p8)*N + (1/p5)*R - Beta*S*(I/N) - 1/p8*S;
der(E) = Beta*S*I/N - 1/p7*E - 1/p8*E;
der(I) = (1/p7)*E - (1/p6)*I - (1/p8)*I;
der(R) = (1/p6)*I - (1/p5)*R - (1/p8)*R;
end SEIRS;
Hi Niklas,
this is something that Christian will have to take a closer look at.
I just wanted to make you aware of a limitation in the sensitivity computations that makes them fail if you have free dependent parameters, i.e., parameters that depend on free parameters in their binding expressions. I have created a ticket for this: https://trac.jmodelica.org/ticket/1240.
Best
/Johan
Great Johan! That was actually the answer to the next question I was going to post! ;)
Btw, to run the above posted SEIRS model, the ticket https://trac.jmodelica.org/ticket/1239 must of course first be resolved.
Best regards, Niklas
Hi Niklas,
About the first error encountered. The initialize flag have to be set to True (when simulating sensitivities), otherwise both the writing to the file will fail and most likely also the start of the simulation. This should though be made more clearly in the documentation.
About the second problem I will have to get back to you.
Btw, ticket:1239 have been fixed.
Best
/Christian
Thanks Christian for 1239!
Then I see why the initialize flag has to be set to true. Maybe an exception should be thrown if Sensitivity=True and Initialize=False, or at least a warning.
Best regards, Niklas
Hi Niklas,
Just a status update, I have looked into the Python crashes and is seems to me that there is a bug in the Sundials suite. I have mailed the Sundials users-group for confirmation. I'll get back to you when I hear from them together with a quick fix.
Best
/Christian
Thanks for the status update. Seems to be a nasty bug.
Best
Niklas
Hi Niklas,
The bug, where Python crashes is unfortunately in Sundials. To fix it, Sundials have to be recompiled as well as Assimulo. The good thing is that only one line of code needs to be changed.
The line that needs to be added in Sundials is in /Sundials/src/idas/idas_ic.c in the method IDASensLineSrch somewhere between line 1155-1169: dtemp = phi[3];
This solves the problem with the crashes and I'm able to simulate your problem between abstol 1e-4 - 1e-8.
Best
/Christian
Great!
I've tried to recompile both sundials and assimulo with no success and when I did the nosetests in assimulo it returned the error message below.
I also tried to recompile assimulo with the working sundials-compilation in /JModelica.org-SDK/SUNDIALS resulting in the same error message.
Here is my detailed installation attempt.
1. Unzipping Assimulo-1.3b1.zip
2. Opening msys in the JModelica.org folder and changes directory to the unzipped assimulo folder.
3. Executing:
python setup_source.py install --sundials-home=/JModelica.org-SDK/SUNDIALS/
4. Checking that assimulo is created in c/Python26/Lib/site-packages/assimulo.
5. Running nosetests in assimulo/tests
Can you see anything I'm doing wrong in the steps above.
Best
Niklas
====================== Error message ======================
$ nosetests
.................E
[KINSOL ERROR] KINSol
The linear solver's setup function failed in an unrecoverable manner.
...........................F............................................................
======================================================================
ERROR: Tests if constraints are implemented correctly in kinsol.py
----------------------------------------------------------------------
Traceback (most recent call last):
File "c:\Python26\lib\site-packages\nose-0.11.3-py2.6.egg\nose\case.py", line 186, in runTest
self.test(*self.arg)
File "c:\Users\niklasa\Downloads\Assimulo-1.3b1\tests\test_KINSOL.py", line 178, in test_constraints_usage
nose.tools.assert_raises(KINSOL_Exception,KINSOL,self.pb_const)
File "c:\Python26\lib\unittest.py", line 336, in failUnlessRaises
callableObj(*args, **kwargs)
File "C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\kinsol.py", line 117, in __init__
self.problem.print_var_info(i)
File "C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\problem_algebraic.py", line 117, in print_var_info
raise ProblemAlg_Exception('The routine print_var_info is not specified.')
ProblemAlg_Exception: The routine print_var_info is not specified.
-------------------- >> begin captured stdout << ---------------------
Problem sent to KINSOL solved.
Type problem with: 1 int32
Variables breaking initial constraint:
--------------------- >> end captured stdout << ----------------------
======================================================================
FAIL: This function tests the handle result.
----------------------------------------------------------------------
Traceback (most recent call last):
File "c:\Python26\lib\site-packages\nose-0.11.3-py2.6.egg\nose\case.py", line 186, in runTest
self.test(*self.arg)
File "c:\Users\niklasa\Downloads\Assimulo-1.3b1\tests\test_implicit_ode.py", line 167, in test_handle_result
assert sim.temp == 141
AssertionError:
-------------------- >> begin captured stdout << ---------------------
Final Run Statistics: ---
Number of Steps : 67
Number of Function Evaluations : 83
Number of Jacobian Evaluations : 16
Number of F-Eval During Jac-Eval : 16
Number of Root Evaluations : 0
Number of Error Test Failures : 2
Number of Nonlinear Iterations : 83
Number of Nonlinear Convergence Failures : 0
Solver options:
Solver : IDA (BDF)
Maxord : 5
Suppress Alg : 0
Tolerances (absolute) : 1e-06
Tolerances (relative) : 1e-06
Elapsed simulation time: 0.00610756358145 seconds.
101
Final Run Statistics: ---
Number of Steps : 103
Number of Function Evaluations : 126
Number of Jacobian Evaluations : 33
Number of F-Eval During Jac-Eval : 33
Number of Root Evaluations : 0
Number of Error Test Failures : 3
Number of Nonlinear Iterations : 126
Number of Nonlinear Convergence Failures : 0
Solver options:
Solver : IDA (BDF)
Maxord : 5
Suppress Alg : 0
Tolerances (absolute) : 1e-06
Tolerances (relative) : 1e-06
Elapsed simulation time: 0.00266615640274 seconds.
138
--------------------- >> end captured stdout << ----------------------
----------------------------------------------------------------------
Ran 106 tests in 9.091s
FAILED (errors=1, failures=1)
Hi Niklas,
Assimulo is included in the source tree when you check out JModelica.org. All you should need to do is to run run
make clean
and then
make install
When you run make install, Assimulo is rebuilt. Prior to doing this, you need to recompile sundials with the change indicated by Christian.
There should be no need to unzip and install Assimulo in your Python-directory, since the Assimulo package will end up in $JMODELICA_HOME/Python when you install JModelica. Actually, the problem may be due to conflicting versions of Assimulo. Therefore, I suggest that you uninstall the Assimulo installation in your Python distro and just run make clean/make install.
Best
/Johan
Thanks,
the recompilation worked just fine now and the attached example now works with abstol: 1e-4.
Best
Niklas
Great! Sorry about the delay in finding the problem.
Best
/Christian
No worries. I think it's reasonable with delays over christmas!
Best
Niklas
Hi again!
I've tested a little more now and it seems like nothing has changed. With atol>1e-3 I still get crashes. Do you? And with atol<1e-4, the problem solves successfully. If I do not remember wrong it also solved succesfully before.
And also, the real bigger problem mentioned in the first post, still have python crashes as described above.
So if you don't get crashes with atol>1e-3, then something must be wrong with my SUNDIALS-recompilation. Any suggestions what mistake I did this time?
Best
Niklas
Here is my detailed installation attempt.
1. Unzipping sundials-2.4.0.tar
2. Steps to the folder with the configure file and runs
mkdir build
cd build
../configure
make && make install
3. copies /usr/local to /JModelica.org-SDK/SUNDIALS
4. rebuilding JModelica as usual with
cd /JModelica.org-SDK
rm -r build
rm -r install
./configure
cd build
make -j9
make install
Hi Niklas,
Did you in between step 1 and 2 modify idas_ic.c? I'm able to solve your problem upto 1e-8 where the problem seems to be very difficult to solve as I get an report that that it takes to many steps to reach the designated output.
Best
/Christian
Sorry, I forgot to mention in the list that I also modified idas_ic.c in between 1 and 2.
I managed to simulate down to 1e-7. However the results are a bit random, sometimes it solves with 1e-7 and sometimes it fails. I've seen this strange behaviour before.
But is python supposed to crash with abstol 1e-3? Does it crash for you?
Best
Niklas
Hi Niklas,
Python is not supposed to crash with abstol 1e-3. I'm able to solve the problem with 1e-3 without python crashes. I'll fire up my windows machine tonight and see if I reproduce your problem on that.
About the randomness. We have also seen some issues related to sensitivity calculations sometimes works and somethings fails. I'm still trying to investigate and see where the problem lies, but have been unable so far.
Best
/Christian
Hi Niklas,
Sorry about the delay (again). With the changes made in Sundials suggested in the previous post, my simulation does not crash python on windows. By crash I mean here that python.exe have an error. Does it always crash for you when you have abstol 1e-3?
I however, have the randomness. The randomness I have hopefully tracked down, after hours and hours, and I believe it to be in Sundials. I have posted a message on the users-group.
The easiest fix for it would be to specify the option maxord to 4 in the call to the simulate method. opts["IDA_options"]["maxord"] = 4 . I hope this will remove the randomness. Let me know if it still persists.
Best
/Christian
Hi Christian
I still and always get crashes with abstol 1e-3. And it seems like we define python crashes in the same way.
I've tried to recompile it many times now and tried a few things.
- If I add the syntax-horrible line åäö=åäö between the line dtemp = phi[3]; and compile I get some warnings in the sundials-compilation but it's still possible to simulate with JModelica with no exceptions. I expected it to totally crash because of the syntax error.
- If I delete the contents from /JModelica.org/SUNDIALS and try to build JModelica, the command ./configure returns that SUNDIALS may be corrupt. This was expected but it shows that it at least looks for SUNDIALS in the right path.
- The files that are created during the SUNDIALS-recompilation are the same as those that comes with JModelica-SDK and they get a new creation time stamp for each compilation.
I have no idea how to debug this any more. Any ideas of things to test?
I also tried to see if the randomness could be helped with changing the maximum order as proposed, but I could not reproduce any randomness at all today.
Best
Niklas
Hi Niklas,
Hmm, this is weird. I would try to run the script through a debugger, like gdb to see whats going on. Just tried gdb on my Windows machine and it seems to be working. To install gdb I just downloaded mingw-get from http://www.mingw.org/wiki/Getting_Started (there is a link under the section "Command line interface" in the second paragraph) then open a command window and browse to C:/mingw/bin (or where your mingw installation is, oh just noticed that I had a second one, the one you want is most likely under JModelica-SDK/mingw, but it should work the same) and run mingw-get install gdb.
Then from the MSYS shell you just browse to your python script and run "gdb python" and "run seris.py" then you could look at the backtrace (bt) and see what is going on and also post it here.
Best
/Christian
Hi Christian,
When I run the model with atol 1e-3 in gdb I get the output below.
Best
Niklas
============ output ======================
(gdb) run seirs.py
Starting program: c:\Python26/python.exe seirs.py
[New Thread 9564.0x19e8]
[New Thread 9564.0x9d0]
[New Thread 9564.0x21b4]
[New Thread 9564.0x265c]
[New Thread 9564.0x21dc]
[New Thread 9564.0x2158]
[New Thread 9564.0xeb8]
[New Thread 9564.0x19f8]
[New Thread 9564.0x2268]
[New Thread 9564.0x1d00]
Program received signal SIGSEGV, Segmentation fault.
0x0a6daeda in sundials_core!initsundials_core () from C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\lib\sundials_core.pyd
Hi Niklas,
What is the output if you type "bt" (backtrace) after the script have got the segmentation fault?
/Christian
Here it is.
Best
Niklas
(gdb) bt
#0 0x0a6daeda in sundials_core!initsundials_core () from C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\lib\sundials_core.pyd
#1 0x0a6f1136 in sundials_core!initsundials_core () from C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\lib\sundials_core.pyd
#2 0x0a6f29c8 in sundials_core!initsundials_core () from C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\lib\sundials_core.pyd
#3 0x0a6c15d0 in sundials_core!initsundials_core () from C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\lib\sundials_core.pyd
#4 0x1e0106fc in python26!PyObject_GetAttr () from C:\Windows\SysWOW64\python26.dll
#5 0x08f2dca0 in ?? ()
#6 0x08eec2c8 in ?? ()
Backtrace stopped: previous frame inner to this frame (corrupt stack?)
Hi Niklas,
I'm afraid that the output does not shed any light on the issue. It may be so that Python have to be built with the debugging option enable to provide more information. But I'm unsure how to do that on Windows. Maybe Jesper have an idé on providing more information from gdb? Otherwise I'm unsure on how to proceed at the moment.
Best
/Christian
Hi Christian,
I suspected that the output was of no use. I will also test this on another windows machine to see if the error exists only on this machine. Thanks for your debugging efforts so far.
Best Niklas
I've now tested it on a Vista x86 computer (the other machine was a win7 x64) and gets the same crashes and gdb outputs.
Best
Niklas
Hi Niklas,
Ok this is interesting. I have also asked Tove to try it on her machine and I might have access to a Windows 7 64bit to try it there as well. I'll get back to you with our results.
Best
/Christian
I had a J meeting today where we looked at the Sundials recompilation. He told me that it was important to always run "make clean" when rebuilding JModelica.org, which I have always skipped because I thought it was not necessary. Now I get the following error message when compiling Seirs.mop. J also said that Tove saw this bug when she tried to reproduce the python crashes.
Best
Niklas
---------------------------------------------------------------------------
JError Traceback (most recent call last)
c:\Users\niklasa\Personal\Problems\SEIRS_bugreport_pythonCrash\SEIRS.py in ()
11 from scipy.optimize import *
12
---> 13 jmu_name = compile_jmu("SEIRS", "SEIRS.mop",compiler_options={"enable_variable_scaling":True}
14 #jmu_name = 'SEIRS.jmu'
15 # Load a model instance into Python
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\jmi.py in compile_jmu(class_name, file_na
t, compiler_options, compile_to, compiler_log_level)
6124
6125 # compile jmu in Java
-> 6126 comp.compile_JMU(class_name, file_name, target, compile_to)
6127
6128 return os.path.join(compile_to, get_jmu_name(class_name))
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\compiler.py in compile_JMU(self, class_na
et, compile_to)
413 self._compiler.compileJMU(class_name, file_name, target, compile_to)
414 except jpype.JavaException, ex:
--> 415 self._handle_exception(ex)
416
417 def compile_FMU(self, class_name, file_name, target, compile_to):
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\compiler.py in _handle_exception(self, ex
659 raise JError(ex.stacktrace().encode('utf-8'))
660
--> 661 raise JError(ex.stacktrace().encode('utf-8'))
662
663 class OptimicaCompiler(ModelicaCompiler):
JError: org.jmodelica.optimica.compiler.ConstantEvaluationException: Cannot get real value of CValueU
at org.jmodelica.optimica.compiler.CValue.realValue(CValue.java:70)
at org.jmodelica.optimica.compiler.FRealVariable.startAttribute_compute(FRealVariable.java:62
at org.jmodelica.optimica.compiler.FRealVariable.startAttribute(FRealVariable.java:614)
at org.jmodelica.optimica.compiler.XMLGenerator.addRealAttributes(XMLGenerator.java:474)
at org.jmodelica.optimica.compiler.OptimicaXMLGenerator.addRealAttributes(OptimicaXMLGenerato
at org.jmodelica.optimica.compiler.XMLGenerator$DAETag_XML_modelVariables.generate(XMLGenerat
at org.jmodelica.optimica.compiler.AbstractGenerator.generate(AbstractGenerator.java:106)
at org.jmodelica.optimica.compiler.AbstractGenerator.generate(AbstractGenerator.java:50)
at org.jmodelica.optimica.compiler.ModelicaCompiler.generateCode(ModelicaCompiler.java:1153)
at org.jmodelica.optimica.compiler.ModelicaCompiler.compileModel(ModelicaCompiler.java:754)
at org.jmodelica.optimica.compiler.ModelicaCompiler.compileUnit(ModelicaCompiler.java:591)
at org.jmodelica.optimica.compiler.ModelicaCompiler.compileJMU(ModelicaCompiler.java:549)
WARNING: Failure executing file:
Hi Niklas,
I also get this error when trying to compile SEIRS, I think it was introduced in a changeset last week. Jesper is working on it and it is also because of this bug that Tove have not been able to test on her Windows computer yet. Sorry
Best
/Christian
Hi Niklas,
We have analyzed the model and came up with an explanation why there is an exception. The problem is essentially that you use a free parameter, in this case R0, as a start value for the variable R. Is this intentional, i.e., do you want to compute sensitivities for the initial condition of R? If this is the case, the problem needs to be reformulated, since this class of problems is currently not supported, I think. (Christian, is sensitivities w.r.t. initial conditions implemented in Assimulo?).
If the start value of R is indeed fixed, then I suggest you restate the problem by introducing another parameter:
parameter Real R0_start = ...;
where you inline the parameter values in the binding expression, and then use R0_start in the start attribute of R. We have discussed a fix that renders this workaround unnecessary, but until then it may enable you and Christian to debug further.
Best
/Johan
Hi again!
with the suggested correction from Johan the SEIRS model works just fine.
Back to the big model where this thread started. For some combinations of atol/rtol I get segmentation faults that crashes python. For other combinations I get MemoryError, see error message below. The memory problem happens in write_data, so I guess that the ode-solver is finished and there is some post processing that causes this. The variable simulator.y is in my case a list of length 55574 with numpy.ndarrays of size 184. I wonder if there is a workaround for this memory error?
Best regards
Niklas
c:\Users\niklasa\pcs\pydroneposts\pyobjfun.py in dynSim(self, pcsd)
--> 593 self.simres = self.sim_model.simulate(t0,tfinal,input=inputtuple,options = opts)
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\jmi.py in simulate(self, start_time, final_time, input, algo
rithm, options)
1998 input,
1999 algorithm,
-> 2000 options)
2001
2002 def optimize(self,
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\core.py in _exec_simulate_algorithm(self, start_time, final_
time, input, algorithm, options)
186 alg.solve()
187 # get and return result
--> 188 return alg.get_result()
189
190 def initialize_options(self, algorithm='IpoptInitializationAlg'):
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\algorithm_drivers.py in get_result(self)
1095 The AssimuloSimResult object.
1096 """
-> 1097 write_data(self.simulator,self.write_scaled_result,self.result_file_name)
1098 # load result file
1099 res = ResultDymolaTextual(self.result_file_name)
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\simulation\assimulo_interface.py in write_data(simulator, wr
ite_scaled_result, result_file_name)
125
126 t = N.array(simulator.t)
--> 127 y = N.array(simulator.y)
128 yd = N.array(simulator.yd)
129 if simulator._problem.input:
MemoryError:
WARNING: Failure executing file:
Hi,
When I get the MemoryError I also get two errors printed, see below, before the exception is printed out. I wonder what is causing those error messages?
Best
Niklas
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (66, 0))
ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (285, 0))
Hi Niklas,
So back to the original problem, I have got another Windows machine which I plan on testing to see if I can reproduce the crashes.
About the memoryerror, this may be solved by writing the result file incrementally http://trac.jmodelica.org/ticket/961 , I'll look at it and see if it can be implemented quickly.
Best
/Christian
Hi Niklas,
About the crashes, I tried again with a new Windows 7 64bit machine (Python 32 bit) and could replicate the crashes for various tolerances. But with the fix posted in post #10 Python stop crashes. So I'm afraid I dont know what more to do about it. I'll talk to Johan.
The memoryerror I'll try to fix in the next couple of days.
Best
/Christian
Hi Christian!
I still get compilation errors, see post #30 when running the original SEIRS model. I think these errors happens because jmodelica does not support sensitivity calculations for start values, in this case R0, according to Johan in post #32. Therefore I can't reproduce the problems anymore. When I change the model according to his suggestions, the model is simulating sensitivities with no problems (but not for start values, such as R0) and I don't get any python crashes anymore.
Thanks for looking at the memoryerror. I'm looking forward to test if it solves my problems with the big model.
Best
Niklas
Hi Niklas,
Yes, JModelica does not currently support sensitivities for start values, sorry. If we disregard that and look at the SEIRS model together with Johan's fix, can you simulate it for all tolerances without python crashes? Is it only your big model that crashes python now?
In changeset https://trac.jmodelica.org/changeset/2479 I have added the ability to write continously to file during simulation of sensitivites, the option is activated with opts["IDA_options"]["write_cont"] = True , and currently only works for simulation of sensitivites. Hopefully this solves the memory issue, if not, please report back.
Best
/Christian
Hi Christian,
After the SEIRS model is modified I have verified that tolerances within 1e-3 - 1e-6 simulates without exception or python crashes.
I've tried to get python crashes for the big model, but I was not able to reproduce it now. I've probably changed something since I first posted the problem. It could also be so that the bug fix you came up with in SUNDIALS solved it. But I'm doubtful because I know that I tried to use your debugging tool gdb also for the big model, which gave the same error message.
For the moment the sensitivity simulations with the big model works fine and I've also verified the calculated sensitivities with finite differences, which seems to match perfectly.
I will rebuild JModelica.org tomorrow and try your incremental writing to file fix. I'll let you know the result.
Best regards
Niklas
Hi Niklas,
This is good news! Lets hope it works tomorrow! Btw, I'll update your other post shortly, with the scaling problem.
Best
/Christian
Hi Christian
It's indeed very good news. Thanks for all your help! :)
I updated to r2479 and tried the write incrementally fix but there seems to be something wrong when simulating with inputs. I reproduced the bug with a modified robertson model, see error message, py- and mop-files below.
Best regards
Niklas
===== Modified robertson.py=========
import os as O
import numpy as N
import matplotlib.pyplot as plt
import nose
from jmodelica.jmi import compile_jmu
from jmodelica.jmi import JMUModel
jmu_name = compile_jmu("Robertson", "Robertson.mop")
# Load a model instance into Python
model = JMUModel(jmu_name)
# Generate input
t = N.linspace(0.,4.,20)
u1 = N.cos(t)
u_traj = N.transpose(N.vstack((t,u1)))
# Get and set the options
opts = model.simulate_options()
opts['IDA_options']['atol'] = [1.0e-8, 1.0e-14, 1.0e-6]
opts['IDA_options']['sensitivity'] = True
opts['IDA_options']['write_cont'] = True
opts['ncp'] = 400
#Simulate
res = model.simulate(final_time=4, input=('u1',u_traj), options=opts)
dy1dp1 = res['dy1/dp1']
dy2dp1 = res['dy2/dp1']
dy3dp1 = res['dy3/dp1']
time = res['time']
nose.tools.assert_almost_equal(dy1dp1[40], -0.35590, 3)
nose.tools.assert_almost_equal(dy2dp1[40], 3.9026e-04, 6)
nose.tools.assert_almost_equal(dy3dp1[40], 3.5551e-01 , 3)
==================================================
============ Robertson.mop =======================
optimization Robertson
parameter Real p1(free=true)=0.040;
parameter Real p2(free=true)=1.0e4;
parameter Real p3(free=true)=3.0e7;
Real y1(start=1.0, fixed=true);
Real y2(start=0.0, fixed=true);
Real y3(start=0.0);
input Real u1;
equation
der(y1) = -p1*y1 + p2*y2*y3;
der(y2) = p1*y1 - p2*y2*y3 - p3*(y2*y2);
0.0 = y1 + y2 + y3 - 1;
end Robertson;
===== Error message =======
In [18]: run robertson
WARNING:root:Maximum order when using IDA for simulating sensitivities is currently limited to 4.
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
C:\Users\niklasa\Personal\Problems\110413_write_increm_inputs\robertson.py in ()
44
45 #Simulate
---> 46 res = model.simulate(final_time=4, input=('u1',u_traj), options=opts)
47
48 dy1dp1 = res['dy1/dp1']
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\jmi.pyc in simulate(self, start_time, final_time, input, alg
orithm, options)
2008 input,
2009 algorithm,
-> 2010 options)
2011
2012 def optimize(self,
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\core.pyc in _exec_simulate_algorithm(self, start_time, final
_time, input, algorithm, options)
184 options)
185 # simulate
--> 186 alg.solve()
187 # get and return result
188 return alg.get_result()
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\algorithm_drivers.pyc in solve(self)
1097 if self.model.has_cppad_derivatives() and self.initialize:
1098 self.simulator.initiate()
-> 1099 self.simulator.simulate(self.final_time, self.ncp)
1100
1101 def get_result(self):
C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\ode.pyc in simulate(self, tfinal, ncp)
266
267 """
--> 268 self.__call__(tfinal,ncp)
269
270 def initiate(self):
C:\Program_Files\JModelica.org-SDK\install\Python\assimulo\implicit_ode.pyc in __call__(self, tfinal, ncp)
233 time_start = time.clock()
234
--> 235 self._problem.handle_result(self,t0,y0,yd0) #Logg the first point
236 self._flag_init = True #Reinitiate the solver
237
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\simulation\assimulo_interface.pyc in handle_result(self, sol
ver, t, y, yd)
1228 data = N.append(data, y[0:len(self._model.real_x)])
1229 if self.input!=None:
-> 1230 data = N.append(data, self.input[1].eval(t)[0,:])
1231 data = N.append(data, y[len(self._model.real_x):len(self._model.real_x)+len(self._model.real_w)])
1232
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\core.pyc in eval(self, x)
485 y = N.zeros([N.size(x),N.size(self.ordinate,1)])
486 for i in range(N.size(y,1)):
--> 487 y[:,i] = N.interp(x,self.abscissa,self.ordinate[:,i])
488 return y
489
C:\Python26\lib\site-packages\numpy\lib\function_base.pyc in interp(x, xp, fp, left, right)
1143 return compiled_interp([x], xp, fp, left, right).item()
1144 else:
-> 1145 return compiled_interp(x, xp, fp, left, right)
1146
1147
ValueError: object of too small depth for desired array
WARNING: Failure executing file:
Hi Niklas,
I'm unable to reproduce the error on my Ubuntu computer, but I believe I have seen something similar on other windows computers in the past and I will try to test later today.
Best
/Christian
Ok. It seems that x has wrong format for N.interp(), namely array(0). I get no exception if I replace t to N.array([t]): in row 1230 in assimulo_interface.py:
1230 -> data = N.append(data, self.input[1].eval(N.array([t]))[0,:]).
I run with this fix for now.
Best regards
Niklas
Hi Christian!
I updated to revision 2492 and now I get the following output when running the modified robertson example in http://jmodelica.org/358#comment-1243. I did a complete install so that my quick fix should be removed.
Best regards
Niklas
In [7]: run robertson
WARNING:root:Maximum order when using IDA for simulating sensitivities is currently limited to 4.
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
C:\Users\niklasa\Personal\Problems\110413_write_increm_inputs\robertson.py in ()
44
45 #Simulate
---> 46 res = model.simulate(final_time=4, input=('u1',u_traj), options=opts)
47
48 dy1dp1 = res['dy1/dp1']
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\jmi.py in simulate(self, start_time, final_time, input, algo
rithm, options)
2008 input,
2009 algorithm,
-> 2010 options)
2011
2012 def optimize(self,
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\core.py in _exec_simulate_algorithm(self, start_time, final_
time, input, algorithm, options)
190 # initialize algorithm
191 alg = algorithm(start_time, final_time, input, self,
--> 192 options)
193 # simulate
194 alg.solve()
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\algorithm_drivers.py in __init__(self, start_time, final_tim
e, input, model, options)
1005 TrajectoryLinearInterpolation(self.input[1][:,0], \
1006 self.input[1][:,1:])), \
-> 1007 self.result_file_name)
1008 else:
1009 if not self.input:
C:\Program_Files\JModelica.org-SDK\install\Python\jmodelica\simulation\assimulo_interface.py in __init__(self, model, in
put, result_file_name)
1086 """
1087 if input != None and not isinstance(input[0],list):
-> 1088 input[0] = [self.input[0]]
1089
1090 self._model = model
AttributeError: 'JMIDAESens' object has no attribute 'input'
WARNING: Failure executing file:
Hi Niklas,
The line should be : input = ([input[0]], input[1]) , I'll commit it in a few minutes.
Sorry
/Christian
Hi again!
With this fix the option write_cont works perfect! Thanks for fix!
Best regards
Niklas
Hi Christian!
I wonder if the bug fix in Sundials you proposed in post http://jmodelica.org/358#comment-986 is included in the latest binary installer file? If not, is there a simple way to include the sundials bug fix in an installation from the binary installer?
Best
Niklas
Hi Niklas,
We have included the bugfix in the upcoming versions of JModelica.org-SDK-1.5 and JModelica.org-1.6. We are awaiting the new Sundials release, but until then we make a local fix in our distro.
Best
/Johan