Python ODEsystem.prepDirection方法代码示例

# 需要导入模块: from PyDSTool.Generator import ODEsystem [as 别名]
# 或者: from PyDSTool.Generator.ODEsystem import prepDirection [as 别名]
    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)
#.........这里部分代码省略.........

# 需要导入模块: from PyDSTool.Generator import ODEsystem [as 别名]
# 或者: from PyDSTool.Generator.ODEsystem import prepDirection [as 别名]
    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:
#.........这里部分代码省略.........

# 需要导入模块: from PyDSTool.Generator import ODEsystem [as 别名]
# 或者: from PyDSTool.Generator.ODEsystem import prepDirection [as 别名]
    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)
#.........这里部分代码省略.........

# 需要导入模块: from PyDSTool.Generator import ODEsystem [as 别名]
# 或者: from PyDSTool.Generator.ODEsystem import prepDirection [as 别名]
    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
#.........这里部分代码省略.........