<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Benny,<div class=""><br class=""></div><div class="">Many thanks for your help.  It seems like perhaps my usage is a bit specialized, and maybe does not warrant a pull to SciPy.</div><div class=""><br class=""></div><div class="">In our particular case we also do some runtime Cython code generation for doing the RHS sparse matvec (if time-dependent) that I think would be a tricky thing to cast into real and imag parts without some major rewriting.  As we do quantum mechanics stuff, the vector space is naturally complex.</div><div class=""><br class=""></div><div class="">Best,</div><div class=""><br class=""></div><div class="">Paul</div><div class=""><br class=""><div><blockquote type="cite" class=""><div class="">On Mar 2, 2015, at 6:20 PM, Benny Malengier <<a href="mailto:benny.malengier@gmail.com" class="">benny.malengier@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class=""><div class="gmail_extra"><br class="Apple-interchange-newline"><br class=""><div class="gmail_quote">2015-03-02 9:59 GMT+01:00 Paul Nation<span class="Apple-converted-space"> </span><span dir="ltr" class=""><<a href="mailto:nonhermitian@gmail.com" target="_blank" class="">nonhermitian@gmail.com</a>></span>:<br class=""><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="auto" class=""><div class="">The zvode docs say that mode 5 returns the exact time answer and not an interpolated result.</div><span class=""><div class=""><br class=""></div><div class="">Paul<br class=""></div></span></div></blockquote><div class=""><br class=""></div><div class="">Yes, I did not say mode 5 was not like that. I noted that the other mode also give you a good value at the wanted output time, only interpolated. In practise this would not be very different.<br class=""><br class="">Changing RHS as you have is a valid use of stop time. Just being able to stop on stop time is quite limiting however, more complicated solver have rootfinding, and can stop on every root found (stoptime is just root of t-stoptime).<br class=""><br class=""></div><div class="">As to ZVODE, as far as I'm aware, nobody developed this further? One should take real and imaginary and convert to a CVODE problem I believe, but I never did complex problems. I do see<span class="Apple-converted-space"> </span><a href="http://netlib.sandia.gov/ode/zvode.f" class="">http://netlib.sandia.gov/ode/zvode.f</a><span class="Apple-converted-space"> </span>mentions such an approach:<span class="Apple-converted-space"> </span><br class=""><pre class="">or a complex stiff ODE system in which f is not analytic, <br class="">ZVODE is likely to have convergence
failures, and for this problem one should instead use DVODE on the
equivalent real system (in the real and imaginary parts of y).</pre></div><div class="">I don't mind somebody adding extra capabilities to VODE/ZVODE in scipy, if you do a nice PR it probably would be accepted if there is no influence on old code.<span class="Apple-converted-space"> </span><br class="">A somewhat better approach in my view would be to deprecate it and convert to code which still sees releases. The API would be more complex though ... That has been discussed before though with no movement, people seem hapy to use ode.integrate for simple things, and then use other bindings when the problem outgrows VODE.<br class=""><br class="">Note that next release of the Sundials suite should be in the coming months. Over the years, one can assume bugs in VODE/ZVODE have been fixed in CVODE only.<span class="Apple-converted-space"> </span><br class=""><br class=""></div><div class="">Benny<br class=""></div><div class=""><br class=""> <br class=""></div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div dir="auto" class=""><span class=""><div class=""><br class=""></div><div class=""><br class="">On Mar 2, 2015, at 17:43, Benny Malengier <<a href="mailto:benny.malengier@gmail.com" target="_blank" class="">benny.malengier@gmail.com</a>> wrote:<br class=""><br class=""></div></span><div class=""><div class="h5"><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><br class=""><div class="gmail_extra"><br class=""><div class="gmail_quote">2015-03-02 6:49 GMT+01:00 Paul Nation<span class="Apple-converted-space"> </span><span dir="ltr" class=""><<a href="mailto:nonhermitian@gmail.com" target="_blank" class="">nonhermitian@gmail.com</a>></span>:<br class=""><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="word-wrap: break-word;" class="">When using the single-step mode in either the VODE or ZVODE ode solver, the default mode (2) called in:<br class=""><br class="">    def step(self, *args):<br class="">       <span class="Apple-converted-space"> </span>itask = self.call_args[2]<br class="">        self.call_args[2] = 2<span class="Apple-converted-space"> </span><span style="white-space: pre-wrap;" class="">       </span># Step mode is here<br class="">       <span class="Apple-converted-space"> </span>r = self.run(*args)<br class="">        self.call_args[2] = itask<br class="">        return r<br class=""><br class=""><div class="">results in taking a single step that (typically) goes beyond the output time requested in the solver.  When doing, for example, monte carlo algorithms, this leads to a big performance hit because one must take a step back, reset the solver  and then use the normal mode to go to the requested stop time.  Instead, these solvers support a mode (5) that will never step beyond the end time.  The modified step function is in that case:</div></div></blockquote><div class=""><br class=""></div><div class="">You do obtain the output at the requested time though, it is only an interpolated value of the actually computed solutions.<span class="Apple-converted-space"> </span><br class=""></div><div class="">So, the only reason you should do above is if something is happening at a certain time and you want to change data or so. You mention monte carlo, but I don't see how that is related to such a usecase in general. I suppose in your application probably yes, but in general MC does not need this<span class="Apple-converted-space"> </span><br class=""><br class=""></div><div class="">The docs say only to use endtime if you have changing RHS or Jacobian, and to otherwise not try to outsmart the solver, as the solver needs extra work in case you set the endtime.<span class="Apple-converted-space"> </span><br class=""></div><div class=""><br class=""></div><div class="">Note that (Z)VODE was replaced by CVODE by the authors of VODE, which has many improvements and several python bindings, all of which expose setting a stop time. In my view, VODE is only present in scipy as a first attempt solver, to be replaced by more modern solvers for heavy lifting.<br class=""></div><div class=""><br class=""></div><div class="">Benny<br class=""></div><div class=""> <br class=""></div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="word-wrap: break-word;" class=""><div class=""><br class=""></div><div class=""><div class="">def step(self, *args):</div><div class="">       <span class="Apple-converted-space"> </span>itask = self.call_args[2]</div><div class="">       <span class="Apple-converted-space"> </span>self.rwork[0] = args[4]    #Set to stop time</div><div class="">       <span class="Apple-converted-space"> </span>self.call_args[2] = 5       #Set single step mode to stop at requested time.</div><div class="">       <span class="Apple-converted-space"> </span>r = self.run(*args)</div><div class="">       <span class="Apple-converted-space"> </span>self.call_args[2] = itask</div><div class="">       <span class="Apple-converted-space"> </span>return r</div></div><div class=""><br class=""></div><div class="">Currently in order to implement this, one needs to create their own ODE integrator subclass of VODE or ZVODE, overload the step function, then create an ode instance and then finally add the custom integrator using ode._integrator.  I think supporting both options natively would be a nice thing to have in SciPy.</div><div class=""><br class=""></div><div class="">In addition, often it is not necessary to do a full reset of the ode solver using ode.reset().  Often times one just needs to change the RHS vector (and possibly the time) and set the flag for the solver to start anew (ode._integrator.call_args[3] = 1).  This to results in a large performance benefit for things like monte carlo solving. Right now I need to call</div><div class=""><br class=""></div><div class=""><div class="">ode._y = new_vec</div><div class="">ode._integrator.call_args[3] = 1</div></div><div class=""><br class=""></div><div class="">when I want to accomplish this.  Adding support for a “fast reset” might also be a good thing to have in SciPy.  </div><div class=""><br class=""></div><div class="">All of the code to accomplish such things are already being used in the QuTiP monte carlo solver(<a href="https://github.com/qutip/qutip/blob/master/qutip/mcsolve.py" target="_blank" class="">https://github.com/qutip/qutip/blob/master/qutip/mcsolve.py</a>) and would therefore be fairly painless to add to SciPy.</div><div class=""><br class=""></div><div class="">Best regards,</div><div class=""><br class=""></div><div class="">Paul</div><div class=""><br class=""></div><div class=""><br class=""></div></div><br class="">_______________________________________________<br class="">SciPy-Dev mailing list<br class=""><a href="mailto:SciPy-Dev@scipy.org" target="_blank" class="">SciPy-Dev@scipy.org</a><br class=""><a href="http://mail.scipy.org/mailman/listinfo/scipy-dev" target="_blank" class="">http://mail.scipy.org/mailman/listinfo/scipy-dev</a><br class=""><br class=""></blockquote></div><br class=""></div></div></div></blockquote><blockquote type="cite" class=""><div class=""><span class="">_______________________________________________</span><br class=""><span class="">SciPy-Dev mailing list</span><br class=""><span class=""><a href="mailto:SciPy-Dev@scipy.org" target="_blank" class="">SciPy-Dev@scipy.org</a></span><br class=""><span class=""><a href="http://mail.scipy.org/mailman/listinfo/scipy-dev" target="_blank" class="">http://mail.scipy.org/mailman/listinfo/scipy-dev</a></span><br class=""></div></blockquote></div></div></div><br class="">_______________________________________________<br class="">SciPy-Dev mailing list<br class=""><a href="mailto:SciPy-Dev@scipy.org" class="">SciPy-Dev@scipy.org</a><br class=""><a href="http://mail.scipy.org/mailman/listinfo/scipy-dev" target="_blank" class="">http://mail.scipy.org/mailman/listinfo/scipy-dev</a><br class=""><br class=""></blockquote></div><br class=""></div></div><span style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px; float: none; display: inline !important;" class="">_______________________________________________</span><br style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class=""><span style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px; float: none; display: inline !important;" class="">SciPy-Dev mailing list</span><br style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class=""><a href="mailto:SciPy-Dev@scipy.org" style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class="">SciPy-Dev@scipy.org</a><br style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class=""><a href="http://mail.scipy.org/mailman/listinfo/scipy-dev" style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;" class="">http://mail.scipy.org/mailman/listinfo/scipy-dev</a></div></blockquote></div><br class=""></div></body></html>