[Scipy-tickets] [SciPy] #478: A new infodict() method for integrate.ode

SciPy scipy-tickets@scipy....
Sat Aug 11 16:55:41 CDT 2007


#478: A new infodict() method for integrate.ode
-----------------------------+----------------------------------------------
 Reporter:  friisj           |       Owner:  somebody
     Type:  enhancement      |      Status:  new     
 Priority:  normal           |   Milestone:          
Component:  scipy.integrate  |     Version:  devel   
 Severity:  normal           |    Keywords:  ode     
-----------------------------+----------------------------------------------
 This patch adds a new infodict() method to trunk/Lib/integrate/ode.py,
 making it possigle to access information about the solution returned by
 the vode solver.

 A important note about how the matrix returned by the Jacobian function
 should be constructed for banded systems is also added to the module
 documentation.


 ----

 {{{
 --- trunk/Lib/integrate/ode.py.orig     2007-08-11 22:36:18.000000000
 +0200
 +++ trunk/Lib/integrate/ode.py  2007-08-11 22:54:07.000000000 +0200
 @@ -26,6 +26,7 @@
             integrator = integrator.set_jac_params(*args)
             y1 = integrator.integrate(t1,step=0,relax=0)
             flag = integrator.successful()
 +           infodict = integrator.infodict()

  Supported integrators:
    vode - Variable-coefficient Ordinary Differential Equation solver,
 @@ -39,12 +40,36 @@
             atol=float|seq
             rtol=float|seq
             lband=None|int
 -           rband=None|int
 +           uband=None|int
             method='adams'|'bdf'
             with_jacobian=0|1
             nsteps = int
             (first|min|max)_step = float
             order = int        # <=12 for adams, <=5 for bdf
 +         When used (with_jacobian=1), the Jacobian should return a
 +         matrix pd of size nrowpd by len(y). For a full system
 +         (lband=uband=None) nrowpd=len(y) and pd[i,j]=df[i]/dy[j]. For
 +         a banded system nrowpd=lband*2+uband+1 and
 +         pd[uband+i-j][i]=df[i]/dy[j], i.e. the diagonals are returned
 +         along the rows in pd from top down (fortran indexing).
 +
 +         After the first call to integrate(), the infodict() method
 +         returns a dictionary with the following elements:
 +           hu     Last (successfully) used step size.
 +           hcur   The step size to be attempted on the next step.
 +           tcur   Current value of t which the solver has actually
 +                  reached, i.e. the C current internal mesh point in t.
 +           tolsf  Tolerance scale factor.
 +           nst    Number of steps taken so far.
 +           nfs    Number of f evaluations so far.
 +           nje    Number of Jacobian evaluations so far.
 +           nqu    Last (successfully) used order.
 +           nqcur  Order to be attempted on the next step.
 +           imxer  Index of the largest component of the weighted local
 error.
 +           nlu    Number of matrix LU decompositions so far.
 +           nni    Number of nonlinear (Newton) iterations so far.
 +           ncfn   Number of convergence failures so far.
 +           netf   Number of error test failures so far.
  """
  """
  XXX: Integrators must have:
 @@ -105,6 +130,7 @@
      integrator = integrator.set_jac_params(*args)
      y1 = integrator.integrate(t1,step=0,relax=0)
      flag = integrator.successful()
 +    infodict = integrator.infodict()

    Typical usage:
      r = ode(f,jac).set_integrator('vode').set_initial_value(y0,t0)
 @@ -197,6 +223,13 @@
          self.jac_params = args
          return self

 +    def infodict(self):
 +        """Return a dictionary of optional output."""
 +        try: self._integrator
 +        except AttributeError: self.set_integrator('')
 +        return self._integrator.infodict()
 +
 +
  #############################################################
  #### Nothing interesting for an end-user in what follows ####
  #############################################################
 @@ -240,6 +273,12 @@
          raise NotImplementedError,'%s does not support run_relax()
 method' %\
                (self.__class__.__name__)

 +    def infodict(self):
 +        """Return a dictionary of optional output."""
 +        raise NotImplementedError,'%s does not support infodict() method'
 %\
 +              (self.__class__.__name__)
 +
 +
      #XXX: __str__ method for getting visual state of the integrator

  class vode(IntegratorBase):
 @@ -250,7 +289,9 @@
          _vode = None
      runner = getattr(_vode,'dvode',None)

 -    messages = {-1:'Excess work done on this call. (Perhaps wrong MF.)',
 +    messages = {1: 'Nothing was done as TOUT was equal to T.',
 +                2: 'The integration was performed successfully.',
 +                -1:'Excess work done on this call. (Perhaps wrong MF.)',
                  -2:'Excess accuracy requested. (Tolerances too small.)',
                  -3:'Illegal input detected. (See printed message.)',
                  -4:'Repeated error test failures. (Check all input.)',
 @@ -354,6 +395,7 @@

      def run(self,*args):
          y1,t,istate =
 self.runner(*(args[:5]+tuple(self.call_args)+args[5:]))
 +        self.istate = istate
          if istate <0:
              print 'vode:',self.messages.get(istate,'Unexpected
 istate=%s'%istate)
              self.success = 0
 @@ -375,6 +417,29 @@
          self.call_args[2] = itask
          return r

 +    def infodict(self):
 +        if hasattr(self, 'istate'):
 +            return {'message': self.messages.get(self.istate,'Unexpected
 istate=%s'%self.istate),
 +                    'hu':    self.rwork[10],
 +                    'hcur':  self.rwork[11],
 +                    'tcur':  self.rwork[12],
 +                    'tolsf': self.rwork[13],
 +                    'nst':   self.iwork[10],
 +                    'nfe':   self.iwork[11],
 +                    'nje':   self.iwork[12],
 +                    'nqu':   self.iwork[13],
 +                    'nqcur': self.iwork[14],
 +                    'imxer': self.iwork[15],
 +                    'lenrw': self.iwork[16],
 +                    'lwniw': self.iwork[17],
 +                    'nlu':   self.iwork[18],
 +                    'nni':   self.iwork[19],
 +                    'ncfn':  self.iwork[20],
 +                    'netf':  self.iwork[21],
 +                    }
 +        else:
 +            return {'message': 'the integrate() method has not yet been
 called'}
 +
  if vode.runner:
      IntegratorBase.integrator_classes.append(vode)

 @@ -391,6 +456,8 @@
      while ode_runner.successful() and ode_runner.t < 50:
          y1 = ode_runner.integrate(ode_runner.t+2)
          print ode_runner.t,y1[:3]
 +    for k,v in ode_runner.infodict().iteritems():
 +        print "    %-7s = %s"%(k,v)

  def test2():
      # Stiff problem. Requires analytic Jacobian.
 @@ -416,10 +483,14 @@
          r.integrate(tout)
          print 'At t=%s  y=%s'%(r.t,r.y)
          tout *= 10
 +    for k,v in r.infodict().iteritems():
 +        print "    %-7s = %s"%(k,v)

  if __name__ == "__main__":
      print 'Integrators available:',\
            ', '.join(map(lambda c:c.__name__,
                          IntegratorBase.integrator_classes))
 +
 +    from scipy.integrate import ode
      test1()
      test2()
 }}}

-- 
Ticket URL: <http://projects.scipy.org/scipy/scipy/ticket/478>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list