Python Generator.ODEsystem类代码示例

 def addMethods(self):
     # override to add _solver function
     ODEsystem.addMethods(self)
     # Jacobian ignored
     self._solver = euler_solver(getattr(self,self.funcspec.spec[1]))
     self._funcreg['_solver'] = ('self', 'euler_solver(getattr(self,' \
                                               + 'self.funcspec.spec[1]))')

    def __init__(self, kw):
        ODEsystem.__init__(self, kw)

        self.diagnostics._errorcodes = \
                {0: 'Unrecognized error code returned (see stderr output)',
                -1: 'Excess work done on this call. (Perhaps wrong method MF.)',
                -2: 'Excess accuracy requested. (Tolerances too small.)',
                -3: 'Illegal input detected. (See printed message.)',
                -4: 'Repeated error test failures. (Check all input.)',
                -5: 'Repeated convergence failures. (Perhaps bad'
                    ' Jacobian supplied or wrong choice of method MF or tolerances.)',
                -6: 'Error weight became zero during problem. (Solution'
                    ' component i vanished, and ATOL or ATOL(i) = 0.)'
                }
        self.diagnostics.outputStatsInfo = {'errorStatus': 'Error status on completion.'}
        # note: VODE only supports array atol, not rtol.
        algparams_def = {'poly_interp': False,
                         'rtol': 1e-9,
                         'atol': [1e-12 for dimix in range(self.dimension)],
                         'stiff': False,
                         'max_step': 0.0,
                         'min_step': 0.0,
                         'init_step': 0.01,
                         'max_pts': 1000000,
                         'strictdt': False,
                         'use_special': False,
                         'specialtimes': []
                         }
        for k, v in algparams_def.items():
            if k not in self.algparams:
                self.algparams[k] = v

 def __init__(self, kw):
     ODEsystem.__init__(self, kw)
     assert self.funcspec.targetlang == 'python', \
            ('Wrong target language for functional specification. '
             'Python needed for this class')
     assert isinstance(self.funcspec, RHSfuncSpec), ('Vode '
                                 'requires RHSfuncSpec type to proceed')
     self._paraminfo = {'init_step': 'Fixed step size for time mesh.',
                        'strictdt': 'Boolean determining whether to evenly space time mesh (default=False), or to use exactly dt spacing.',
                        'stiff': 'Boolean to activate the BDF method, otherwise Adams method used. Default False.'}
     self._errorcodes = {0: 'Unrecognized error code returned (see stderr output)',
             -1: 'Excess work done on this call. (Perhaps wrong method MF.)',
             -2: 'Excess accuracy requested. (Tolerances too small.)',
             -3: 'Illegal input detected. (See printed message.)',
             -4: 'Repeated error test failures. (Check all input.)',
             -5: 'Repeated convergence failures. (Perhaps bad'
                 ' Jacobian supplied or wrong choice of method MF or tolerances.)',
             -6: 'Error weight became zero during problem. (Solution'
                 ' component i vanished, and ATOL or ATOL(i) = 0.)'
             }
     self._outputinfo = {'errorStatus': 'Error status on completion.'}
     # note: VODE only supports array atol, not rtol.
     algparams_def = {'rtol': 1e-9,
                      'atol': [1e-12 for dimix in xrange(self.dimension)],
                      'stiff': False,
                      'max_step': 0.0,
                      'min_step': 0.0,
                      'init_step': 0.01,
                      'max_pts': 1000000,
                      'strictdt': False
                      }
     for k, v in algparams_def.iteritems():
         if k not in self._algparams:
             self._algparams[k] = v
     self.outputstats = {}

 def addMethods(self):
     # override to add _solver function
     ODEsystem.addMethods(self, usePsyco=False)
     if self.haveJacobian():
         self._solver = ode(getattr(self,self.funcspec.spec[1]),
                         getattr(self,self.funcspec.auxfns["Jacobian"][1]))
         self._funcreg['_solver'] = ('self',
              'ode(getattr(self,self.funcspec.spec[1]),' \
                + 'getattr(self,self.funcspec.auxfns["Jacobian"][1]))')
     else:
         self._solver = ode(getattr(self,self.funcspec.spec[1]))
         self._funcreg['_solver'] = ('self', 'ode(getattr(self,' \
                                               + 'self.funcspec.spec[1]))')

 def __init__(self, kw):
     if 'user_func_beforestep' in kw:
         self.ufunc_before = kw['user_func_beforestep']
         # delete because not covered in ODEsystem
         del kw['user_func_beforestep']
     else:
         self.ufunc_before = _dummy_userfunc
     if 'user_func_afterstep' in kw:
         self.ufunc_after = kw['user_func_afterstep']
         # delete because not covered in ODEsystem
         del kw['user_func_afterstep']
     else:
         self.ufunc_after = _dummy_userfunc
     ODEsystem.__init__(self, kw)
     self._paraminfo = {'init_step': 'Fixed step size for time mesh.'}
     self.diagnostics._errorcodes = {1: 'Step OK'}
     self.diagnostics.outputStatsInfo = {'errorStatus': 'Error status on completion.'}
     algparams_def = {'poly_interp': False,
                      'init_step': 0.01,
                      'max_pts': 100000
                      }
     for k, v in algparams_def.items():
         if k not in self.algparams:
             self.algparams[k] = v

 def __del__(self):
     ODEsystem.__del__(self)

    def compute(self, trajname, dirn='f', ics=None):
        continue_integ = ODEsystem.prepDirection(self, dirn)
        if self._dircode == -1:
            raise NotImplementedError('Backwards integration is not implemented')
        if ics is not None:
            self.set(ics=ics)
        self.validateICs()
        self.diagnostics.clearWarnings()
        self.diagnostics.clearErrors()
        pnames = sortedDictKeys(self.pars)
        xnames = self._var_ixmap  # ensures correct order
        # Check i.c.'s are well defined (finite)
        self.checkInitialConditions()
        if self.algparams['stiff']:
            methstr = 'bdf'
            methcode = 2
        else:
            methstr = 'adams'
            methcode = 1
        haveJac = int(self.haveJacobian())
        if isinstance(self.algparams['atol'], list):
            if len(self.algparams['atol']) != self.dimension:
                raise ValueError('atol list must have same length as phase '
                                 'dimension')
        else:
            atol = self.algparams['atol']
            self.algparams['atol'] = [atol for dimix in range(self.dimension)]
        indepdom0, indepdom1 = self.indepvariable.depdomain.get()
        if continue_integ:
            if indepdom0 > self._solver.t:
                print("Previous end time is %f"%self._solver.t)
                raise ValueError("Start time not correctly updated for "
                                 "continuing orbit")
            x0 = self._solver.y
            indepdom0 = self._solver.t
            self.indepvariable.depdomain.set((indepdom0,indepdom1))
        else:
            x0 = sortedDictValues(self.initialconditions,
                                        self.funcspec.vars)
        t0 = indepdom0
        if self.algparams['use_special']:
            tmesh_special = list(self.algparams['specialtimes'])
            if continue_integ:
                if self._solver.t not in tmesh_special:
                    raise ValueError("Invalid time to continue integration:"
                                     "it is not in 'special times'")
                tmesh_special = tmesh_special[tmesh_special.index(self._solver.t):]
            try:
                dt = min([tmesh_special[1]-t0, self.algparams['init_step']])
            except:
                raise ValueError("Invalid list of special times")
            if not isincreasing(tmesh_special):
                raise ValueError("special times must be given in increasing "
                                 "order")
            if self.indepvariable.depdomain.contains(t0) is notcontained or \
              self.indepvariable.depdomain.contains(tmesh_special[-1]) is notcontained:
                raise PyDSTool_BoundsError("special times were outside of independent "
                                 "variable bounds")
        else:
            tmesh_special = []
            dt = self.algparams['init_step']
        # speed up repeated access to solver by making a temp name for it
        solver = self._solver
        if solver._integrator is None:
            # Banded Jacobians not yet supported
            #
            # start a new integrator, because method may have been
            # switched
            solver.set_integrator('vode', method=methstr,
                                 rtol=self.algparams['rtol'],
                                 atol=self.algparams['atol'],
                                 nsteps=self.algparams['max_pts'],
                                 max_step=self.algparams['max_step'],
                                 min_step=self.algparams['min_step'],
                                 first_step=dt/2.,
                                 with_jacobian=haveJac)
        else:
            solver.with_jacobian = haveJac
#            self.mu = lband
#            self.ml = uband
            solver.rtol = self.algparams['rtol']
            solver.atol = self.algparams['atol']
            solver.method = methcode
#            self.order = order
            solver.nsteps = self.algparams['max_pts']
            solver.max_step = self.algparams['max_step']
            solver.min_step = self.algparams['min_step']
            solver.first_step = dt/2.
        solver.set_initial_value(x0, t0)
        ## Broken code for going backwards (doesn't respect 'continue' option
        ##  either)
##            if self._dircode == 1:
##                solver.set_initial_value(x0, t0)
##            else:
##                solver.set_initial_value(x0, indepdom1)
        # wrap up each dictionary initial value as a singleton list
        alltData = [t0]
        allxDataDict = dict(zip(xnames, map(listid, x0)))
        plist = sortedDictValues(self.pars)
        extralist = copy(plist)
#.........这里部分代码省略.........

    def __init__(self, kw):
        """Use the nobuild key to postpone building of the library, e.g. in
        order to provide additional build options to makeLibSource and
        compileLib methods or to make changes to the C code by hand.
        No build options can be specified otherwise."""

        # Building is just doing make
        if 'nobuild' in kw:
            nobuild = kw['nobuild']
            del kw['nobuild']
        else:
            nobuild = False
        ODEsystem.__init__(self, kw)
        self._solver = None
        assert self.funcspec.targetlang == 'matlab', \
               ('Wrong target language for functional specification. '
                'matlab needed for this class')
        assert isinstance(self.funcspec, RHSfuncSpec), ('ADMC++ '
                                    'requires RHSfuncSpec type to proceed')
        assert not self.inputs, \
                        'ADMC++ does not support external inputs feature'
        self._errorcodes = {}
        self._paraminfo = {}

        self.vftype = 'vfieldts'

        # currently the final four of these params are for event handling
        # NEED TO CHECK WHICH ONES ARE SUPPORTED BY ADMC -- LOOKS LIKE EVTOLS ONLY FOR NOW
        # HACK: vftype is alg param for now, tells us whether parent class is hybridvf, vfieldts, etc.
        algparams_def = {'evtols' : 0.0001, 'vftype' : 'vfieldts'}

        # Remove this later
        for k, v in algparams_def.iteritems():
            if k not in self.algparams:
                self.algparams[k] = v

        # verify that no additional keys are present in algparams, after
        # defaults are added above
        if len(self.algparams) != len(algparams_def):
            raise ValueError("Invalid keys present in algparams argument: " \
                     + str(remain(self.algparams.keys(),algparams_def.keys())))

        thisplatform = platform.system()

        self._compilation_tempdir = os.path.join(os.getcwd(),
                                                      "admcpp_temp")
        if not os.path.isdir(self._compilation_tempdir):
            try:
                assert not os.path.isfile(self._compilation_tempdir), \
                     "A file already exists with the same name"
                os.mkdir(self._compilation_tempdir)
            except:
                print "Could not create compilation temp directory " + \
                      self._compilation_tempdir
                raise

        # ADMC targets must go in their own directories with appropriate names
        self._model_dir = "@"+self.name
        self._target_dir = os.path.join(self._compilation_tempdir,self._model_dir)
        # Make the target directory
        if not os.path.isdir(self._target_dir):
            try:
                assert not os.path.isfile(self._target_dir), \
                       "A file already exists with the same name"
                os.mkdir(self._target_dir)
            except:
                print "Could not creat target ADMC model directory " + \
                      self._target_dir
                raise


        """ An ADMC model has the following files:
        vfield.m -- contains code for the RHS of the vector field
        set.m -- a generic method that overload matlab's set method; only need to insert vfield name
        get.m -- a generic method that overloads matlab's get method; only need to insert appropriate parent name
        """

        # model.m, get.m, set.m, vfield.m are minimal files required. TO DO: EVENTS
        self._model_file = self.name+".m"
        self._ic_file = self.name+"_ics.m"
        self._param_file = self.name+"_params.m"
        self._set_file = "set.m"
        self._get_file = "get.m"
        self._vfield_file = "vfield.m"
        self._events_file = "events.m"

        self._vf_filename_ext = "_"+self._model_file[:-2]

        if not nobuild:
            self.makeLibSource()
        else:
            print "Build the library using the makeLib method, or in "
            print "stages using the makeLibSource and compileLib methods."

 def __del__(self):
     genDB.unregister(self)
     ODEsystem.__del__(self)

    def compute(self, trajname, dirn='f', ics=None):
        continue_integ = ODEsystem.prepDirection(self, dirn)
        if ics is not None:
            self.set(ics=ics)
        self.validateICs()
        self.diagnostics.clearWarnings()
        self.diagnostics.clearErrors()
        if isinstance(self.algparams['rtol'], list):
            if len(self.algparams['rtol']) != self.dimension:
                raise ValueError('rtol list must have same length as phase dimension')
        else:
            rtol = self.algparams['rtol']
            self.algparams['rtol'] = [rtol for dimix in range(self.dimension)]
        if isinstance(self.algparams['atol'], list):
            if len(self.algparams['atol']) != self.dimension:
                raise ValueError('atol list must have same length as phase dimension')
        else:
            atol = self.algparams['atol']
            self.algparams['atol'] = [atol for dimix in range(self.dimension)]
        anames = self.funcspec.auxvars
        # Check i.c.'s are well defined (finite)
        self.checkInitialConditions()
        self.setEventICs(self.initialconditions, self.globalt0)
        # update event params in case changed since last run
        self._prepareEventSpecs()
        # Main integration
        t0 = self.indepvariable.depdomain[0]
        t1 = self.indepvariable.depdomain[1]
        plist = sortedDictValues(self.pars)
        self.algparams['hasJac'] = self.haveJacobian()
        self.algparams['hasJacP'] = self.haveJacobian_pars()
        if self._solver is None:
            self._solver = dopri(self.modname,
                                 rhs=self.name, phaseDim=self.dimension,
                                 paramDim=len(plist), nAux=len(anames),
                                 nEvents=len(self._eventNames),
                                 nExtInputs=len(self.inputs),
                                 hasJac=self.algparams['hasJac'],
                                 hasJacP=self.algparams['hasJacP'],
                                 hasMass=self.haveMass(),
                                 extraSpace=self.algparams['extraspace'],
                                 )
            try:
                genDB.register(self)
            except PyDSTool_KeyError:
                errstr = "Generator " + self.name + ": this vector field's " +\
                         "DLL is already in use"
                raise RuntimeError(errstr)
        if self._dircode == 1:
            tbegin = t0
            tend = t1
        elif self._dircode == -1:
            # dopri does reverse time integration simply by switching t0 and t1
            # and using negative steps
            tbegin = t1
            tend = t0
        if len(self.algparams['specialtimes'])>0:
            use_special = self.algparams['use_special']
        else:
            use_special = 0
        bounds = [[],[]]  # lower, then upper
        for v in self.funcspec.vars:
            bds = self.xdomain[v]
            bounds[0].append(bds[0])
            bounds[1].append(bds[1])
        for p in self.funcspec.pars:
            bds = self.pdomain[p]
            bounds[0].append(bds[0])
            bounds[1].append(bds[1])
        if continue_integ:
            x0 = self._solver.lastPoint
            # overwrite t0 from self.indepvariable.domain, but use its t1
            tbegin = self._solver.lastTime
            if abs(self._solver.lastStep) < abs(self.algparams['init_step']):
                self.algparams['init_step'] = self._solver.lastStep
            if abs(t1-tbegin) < abs(self.algparams['init_step']):
                raise ValueError("Integration end point too close to initial "
                                 "point")
#            if self.inputs and self._extInputsChanged:
#                self._extInputsChanged = False
#                self._solver.setContParams(tend, plist,
#                                           use_special,
#                                           self.algparams['verbose'],
#                                           True, deepcopy(self._inputVarList),
#                                           deepcopy(self._inputTimeList))
        else:
            if self._solver.numRuns > 0:
                self._solver.clearAll()
            x0 = sortedDictValues(self.initialconditions, self.funcspec.vars)
            self._solver.setInteg(maxpts=self.algparams['max_pts'],
                rtol=self.algparams['rtol'], atol=self.algparams['atol'])
            self._solver.setRunParams(ic=x0, params=plist,
                                  t0=tbegin, tend=tend, gt0=self.globalt0,
                                  refine=self.algparams['refine'],
                                  specTimes=self.algparams['specialtimes'],
                                  bounds=bounds)
        if self.inputs:
            # self._extInputsChanged if global t0 changed so that can
            # adjust times given to the integrator (it is blind to global t0
            # when accesses input variable times)
#.........这里部分代码省略.........

    def __init__(self, kw):
        """Use the nobuild key to postpone building of the library, e.g. in
        order to provide additional build options to makeLibSource and
        compileLib methods or to make changes to the C code by hand.
        No build options can be specified otherwise."""

        # delete because not covered in ODEsystem
        nobuild = kw.pop('nobuild', False)
        ODEsystem.__init__(self, kw)
        self.diagnostics.outputStatsInfo = {
            'last_step': 'Predicted step size of the last accepted step (useful for a subsequent call to dop853).',
            'num_steps': 'Number of used steps.',
            'num_accept': 'Number of accepted steps.',
            'num_reject': 'Number of rejected steps.',
            'num_fcns': 'Number of function calls.',
            'errorStatus': 'Error status on completion.'
                        }
        self.diagnostics._errorcodes = {
             0 : 'Unrecognized error code returned (see stderr output)',
            -1 : 'input is not consistent',
            -2 : 'larger nmax is needed',
            2 : 'larger nmax or maxevtpts is probably needed (error raised by solout)',
            -3 : 'step size becomes too small',
            -4 : 'the problem is probably stiff (interrupted)',
            -8 : 'The solution exceeded a magbound (poor choice of initial step)'}
        self._solver = None
        algparams_def = {'poly_interp': False,
                        'init_step': 0,
                        'max_step': 0,
                        'rtol': [1e-9 for i in range(self.dimension)],
                        'atol': [1e-12 for i in range(self.dimension)],
                        'fac1': 0.2,
                        'fac2': 10.0,
                        'safety': 0.9,
                        'beta': 0.04,
                        'max_pts': 10000,
                        'refine': 0,
                        'maxbisect': [], # for events
                        'maxevtpts': 1000, # for events
                        'eventInt': [],  # set using setEventInterval only
                        'eventDelay': [], # set using setEventDelay only
                        'eventTol': [], # set using setEventTol only
                        'use_special': 0,
                        'specialtimes': [],
                        'check_aux': 1,
                        'extraspace': 100,
                        'verbose': 0,
                        'hasJac': 0,
                        'hasJacP': 0,
                        'magBound': 1e7,
                        'boundsCheckMaxSteps': 1000,
                        'checkBounds': self.checklevel
                        }
        for k, v in algparams_def.items():
            if k not in self.algparams:
                self.algparams[k] = v
        # verify that no additional keys are present in algparams, after
        # defaults are added above
        if len(self.algparams) != len(algparams_def):
            raise ValueError("Invalid keys present in algparams argument: " \
                     + str(remain(self.algparams.keys(),algparams_def.keys())))

        if self.haveMass():
            raise ValueError("Mass matrix declaration is incompatible "
                             "with Dopri853 integrator system specification")

        self._prepareEventSpecs()
        self._inputVarList = []
        self._inputTimeList = []

        if nobuild:
            print("Build the library using the makeLib method, or in ")
            print("stages using the makeLibSource and compileLib methods.")
        else:
            self.makeLib()

    def compute(self, trajname, dirn='f', ics=None):
        continue_integ = ODEsystem.prepDirection(self, dirn)
        if self._dircode == -1:
            raise NotImplementedError('Backwards integration is not implemented')
        if ics is not None:
            self.set(ics=ics)
        self.validateICs()
        self.diagnostics.clearWarnings()
        self.diagnostics.clearErrors()
        pnames = sortedDictKeys(self.pars)
        xnames = self._var_ixmap  # ensures correct order
        # Check i.c.'s are well defined (finite)
        self.checkInitialConditions()
        haveJac = int(self.haveJacobian())
        indepdom0, indepdom1 = self.indepvariable.depdomain.get()
        if continue_integ:
            if indepdom0 > self._solver.t:
                print("Previous end time is %f"%self._solver.t)
                raise ValueError("Start time not correctly updated for "
                                 "continuing orbit")
            x0 = self._solver.y
            indepdom0 = self._solver.t
            self.indepvariable.depdomain.set((indepdom0,indepdom1))
        else:
            x0 = sortedDictValues(self.initialconditions,
                                        self.funcspec.vars)
        t0 = indepdom0
        dt = self.algparams['init_step']
        # speed up repeated access to solver by making a temp name for it
        solver = self._solver
        solver.set_initial_value(x0, t0)
        solver.dt = dt
        # wrap up each dictionary initial value as a singleton list
        alltData = [t0]
        allxDataDict = dict(zip(xnames, map(listid, x0)))
        plist = sortedDictValues(self.pars)
        extralist = copy(plist)
        if self.inputs:
            # inputVarList is a list of Variables
            inames = sortedDictKeys(self.inputs)
            listend = self.numpars + len(self.inputs)
            inputVarList = sortedDictValues(self.inputs)
            ilist = _pollInputs(inputVarList, alltData[0]+self.globalt0,
                               self.checklevel)
        else:
            ilist = []
            inames = []
            listend = self.numpars
            inputVarList = []
        extralist.extend(ilist)
        solver.set_f_params(extralist)
        if haveJac:
            solver.set_jac_params(extralist)
        do_poly = self.algparams['poly_interp']
        if do_poly:
            rhsfn = getattr(self, self.funcspec.spec[1])
            dx0 = rhsfn(t0, x0, extralist)
            alldxDataDict = dict(zip(xnames, map(listid, dx0)))
        auxvarsfn = getattr(self,self.funcspec.auxspec[1])
        # Make t mesh if it wasn't given as 'specialtimes'
        if not all(isfinite(self.indepvariable.depdomain.get())):
            print("Time domain was: %r" % self.indepvariable.depdomain.get())
            raise ValueError("Ensure time domain is finite")
        if dt == indepdom1 - indepdom0:
            # single-step integration required
            # special times will not have been set (unless trivially
            # they are [indepdom0, indepdom1])
            tmesh = [indepdom0, indepdom1]
        else:
            notDone = True
            while notDone:
                tmesh = self.indepvariable.depdomain.sample(dt,
                                          strict=True,
                                          avoidendpoints=True)
                notDone = False
        first_found_t = None
        eventslist = self.eventstruct.query(['active', 'notvarlinked'])
        termevents = self.eventstruct.query(['term'], eventslist)
        tmesh.pop(0)  # get rid of first entry for initial condition
        xnames = self.funcspec.vars
        # storage of all auxiliary variable data
        allaDataDict = {}
        anames = self.funcspec.auxvars
        avals = auxvarsfn(t0, x0, extralist)
        for aix in range(len(anames)):
            aname = anames[aix]
            try:
                allaDataDict[aname] = [avals[aix]]
            except IndexError:
                print("\nEuler generator: There was a problem evaluating " \
                      + "an auxiliary variable")
                print("Debug info: avals (length %d) was %r" % (len(avals), avals))
                print("Index out of range was %d" % aix)
                print(self.funcspec.auxspec[1])
                print(hasattr(self, self.funcspec.auxspec[1]))
                print("Args were: %r" % [t0, x0, extralist])
                raise
        # Initialize signs of event detection objects at IC
        self.setEventICs(self.initialconditions, self.globalt0)
        if self.inputs:
#.........这里部分代码省略.........

    def compute(self, trajname, dirn='f'):
        continue_integ = ODEsystem.prepDirection(self, dirn)
        if self._dircode == -1:
            raise NotImplementedError, ('Backwards integration is not implemented')
        # validate spec if there exists a prior trajectory computation
        if self.defined:
            self.validateSpec()
            self.validateICs()
            self.clearWarnings()
            self.clearErrors()
        pnames = sortedDictKeys(self.pars)
        xnames = self._var_ixmap  # ensures correct order
        # Check i.c.'s are well defined (finite)
        self.checkInitialConditions()
        if self._algparams['stiff']:
            methstr = 'bdf'
            methcode = 2
        else:
            methstr = 'adams'
            methcode = 1
        if self.haveJacobian():
            haveJac = 1
        else:
            haveJac = 0
        if isinstance(self._algparams['atol'], list):
            if len(self._algparams['atol']) != self.dimension:
                raise ValueError, 'atol list must have same length as phase dimension'
        else:
            atol = self._algparams['atol']
            self._algparams['atol'] = [atol for dimix in xrange(self.dimension)]
        indepdom0 = self.indepvariable.depdomain.get(0)
        indepdom1 = self.indepvariable.depdomain.get(1)
        if continue_integ:
            if self._tdata[0] != self._solver.t:
                print "Previous end time is %f"%self._solver.t
                raise ValueError, \
                      "Start time not correctly updated for continuing orbit"
            x0 = self._solver.y
            indepdom0 = self._solver.t
        else:
            x0 = sortedDictValues(self.initialconditions,
                                        self.funcspec.vars)
        if self._solver._integrator is None:
            # Banded Jacobians not yet supported
            #
            # start a new integrator, because method may have been
            # switched
            self._solver.set_integrator('vode', method=methstr,
                                 rtol=self._algparams['rtol'],
                                 atol=self._algparams['atol'],
                                 nsteps=self._algparams['max_pts'],
                                 max_step=self._algparams['max_step'],
                                 min_step=self._algparams['min_step'],
                                 first_step=self._algparams['init_step'],
                                 with_jacobian=haveJac)
            # speed up repeated access to solver by making a temp name for it
            solver = self._solver
        else:
            # speed up repeated access to solver by making a temp name for it
            solver = self._solver
            solver.with_jacobian = haveJac
#            self.mu = lband
#            self.ml = uband
            solver.rtol = self._algparams['rtol']
            solver.atol = self._algparams['atol']
            solver.method = methcode
#            self.order = order
            solver.nsteps = self._algparams['max_pts']
            solver.max_step = self._algparams['max_step']
            solver.min_step = self._algparams['min_step']
            solver.first_step = self._algparams['init_step']
        solver.set_initial_value(x0, indepdom0)
##            if self._dircode == 1:
##                solver.set_initial_value(x0, indepdom0)
##            else:
##                solver.set_initial_value(x0, indepdom1)
        # wrap up each dictionary initial value as a singleton list
        alltData = [indepdom0]
        allxDataDict = dict(zip(xnames, map(listid, x0)))
        plist = sortedDictValues(self.pars)
        extralist = copy(plist)
        if self.inputs:
            # inputVarList is a list of Variables
            inames = sortedDictKeys(self.inputs)
            listend = self.numpars + len(self.inputs)
            inputVarList = sortedDictValues(self.inputs)
            ilist = _pollInputs(inputVarList, alltData[0]+self.globalt0,
                               self.checklevel)
        else:
            ilist = []
            inames = []
            listend = self.numpars
            inputVarList = []
        extralist.extend(ilist)
        solver.set_f_params(extralist)
        if haveJac:
            solver.set_jac_params(extralist)
        dt = self._algparams['init_step']
        strict = self._algparams['strictdt']
        # Make t mesh
#.........这里部分代码省略.........