Python PyDSTool类代码示例

def event1(params):

    event_args = {'name': 'event_melt_begin',
                  'eventtol': params.RelTol,
                  'active': True,
                  'term': True}

    event_melt_begin = PyDSTool.makeZeroCrossEvent('Tmelt - Tp', 0, event_args, varnames=['Tp'], parnames=['Tmelt'])

    return event_melt_begin

def event2(params):

    event_args = {'name': 'event_melt_end',
                  'eventtol': params.RelTol,
                  'active': True,
                  'term': True}

    event_melt_end = PyDSTool.makeZeroCrossEvent('Qp / (Hf * Mp) - 1', 0, event_args, varnames=['Qp'], parnames=['Hf', 'Mp'])

    return event_melt_end

def args():
    """
    This function creates a PyDSTool 'args' object for the
    'MorrisLecar' vector field.
    """
    DSargs = PyDSTool.args()
    DSargs.name = 'MorrisLecar'
    DSargs.pars = {'gca':5.5000000000000000e+00, 'gk':8.0000000000000000e+00, 'gl':2.0000000000000000e+00, 'vca':1.1500000000000000e+02, 'vk':-8.4000000000000000e+01, 'vl':-5.5000000000000000e+01, 'c':2.0000000000000000e+01, 'phi':2.2000000000000000e-01, 'ic':9.0000000000000000e+01, 'v1':-1.2000000000000000e+00, 'v2':1.8000000000000000e+01, 'v3':2.0000000000000000e+00, 'v4':3.0000000000000000e+01}
    DSargs.varspecs = {'v':'-(1.0/2.0)*1.0/c*( 2.0*( v-vl)*gl-( vca-v)*gca*( tanh(-1.0/v2*( v1-v))+1.0)+-2.0*ic+-2.0*( vk-v)*gk*w)', 'w':'(1.0/2.0)*cosh(-(1.0/2.0)*( v3-v)/v4)*phi*( tanh(-( v3-v)/v4)+-2.0*w+1.0)'}
    DSargs.fnspecs = {'Jacobian': (['t', 'v', 'w'],
            """[[-(1.0/2.0)*1.0/c*( 2.0*gk*w+gca*( tanh(-1.0/v2*( v1-v))+1.0)+( vca-v)*gca/v2*( pow(tanh(-1.0/v2*( v1-v)),2.0)-1.0)+2.0*gl), 1.0/c*( vk-v)*gk],
                [-cosh(-(1.0/2.0)*( v3-v)/v4)*phi*( pow(tanh(-( v3-v)/v4),2.0)-1.0)/v4/2.0+sinh(-(1.0/2.0)*( v3-v)/v4)*phi*( tanh(-( v3-v)/v4)+-2.0*w+1.0)/v4/4.0, -cosh(-(1.0/2.0)*( v3-v)/v4)*phi]]""")}
    DSargs.ics = {'v':0.0, 'w':0.0}
    DSargs.tdomain = [0,10]
    return DSargs

def args():
    """
    This function creates a PyDSTool 'args' object for the
    'vanderpol' vector field.
    """
    DSargs = PyDSTool.args()
    DSargs.name = 'vanderpol'
    DSargs.pars = {'epsilon':2.0000000000000001e-01}
    DSargs.varspecs = {'x':'( x+y+-3.3333333333333331e-01*(x*x*x))/epsilon', 'y':'-x'}
    DSargs.fnspecs = {'Jacobian': (['t', 'x', 'y'],
            """[[(-(x*x)+1.0)/epsilon, 1.0/(epsilon)],
                [-1.0, 0.0]]""")}
    DSargs.ics = {'x':1.0000000000000000e-02, 'y':0.0}
    DSargs.tdomain = [0,10]
    return DSargs

def make_measure(fn_name, fn_spec, **defs):
    """Dynamically create a python function for use with calculation
    context.
    """
    all_defs = defs.copy()
    q = dst.QuantSpec('_dummy_', fn_spec, treatMultiRefs=False)
    import PyDSTool.parseUtils as pu
    mapping = pu.symbolMapClass()
    assumed_modules = []
    tokens = q.parser.tokenized
    for sym in q.freeSymbols:
        # Hack, for now: if first (therefore, assumed all)
        # occurrences of symbol are in quotes, then don't convert.
        # Better solution would be to make parser create "x" as a single
        # symbol, at least with a detect quote option
        first_ix = tokens.index(sym)
        if first_ix == 0 or (first_ix > 0 and tokens[first_ix-1] not in ['"', "'"]):
            if pu.isHierarchicalName(sym):
                parts = sym.split('.')
                if parts[0] == 'sim':
                    mapping[sym] = 'con.'+sym
##                elif parts[0] == 'bombardier':
##                    # special case as this factory function is defined in that
##                    # module so that reference will fail at runtime: remove
##                    # 'bombardier' prefix
##                    rest_sym = '.'.join(parts[1:])
##                    mapping[sym] = rest_sym
##                    scope = globals()
##                    # locals override
##                    scope.update(locals())
##                    if parts[1] in scope:
##                        all_defs[parts[1]] = scope[parts[1]]
##                    else:
##                        raise ValueError("Cannot resolve scope of symbol '%s'"%sym)
                else:
                    # assume module reference
                    assumed_modules.append(parts[0])
                    # record here to ensure inclusion in dyn_dummy
                    mapping[sym] = 'self.'+sym
            else:
                mapping[sym] = 'con.workspace.'+sym
        elif first_ix > 0 and tokens[first_ix-1] in ['"', "'"]:
            # put this symbol in the mapping values to ensure not included
            # as an argument to the function
            mapping[sym] = sym
    q.mapNames(mapping)
    import types
    for module_name in assumed_modules:
        global_scope = globals()
        # test if module name in scope
        if module_name in global_scope:
            _mod = global_scope[module_name]
            if isinstance(_mod, types.ModuleType):
                all_defs[module_name] = _mod

    # dyn_dummy contains dummy mappings but declares symbols to leave
    # evaluating until runtime
    dyn_dummy = dict(zip(mapping.values(), ['']*len(mapping)))
    funq = dst.expr2fun(q, ensure_args=['con'], ensure_dynamic=dyn_dummy,
                   for_funcspec=False, fn_name=fn_name,
                   **all_defs)

    # decorate output
    funq.attr_name = fn_name
    return funq

def model():
    K = 0.4
    E = 7.0
    M = 10.5
    N = 15
    Lambda = 0.9
    Gamma  = 12
    R  = 0.7
    PP = 20
    
    lamb_p    = (K*Gamma)/M
    eta_p     = E*K
    p_p       = (Lambda*PP)/(K*Gamma*M)
    nu_p     = N/M
    rho_p     = R
    
    # Declare names and initial values for (symbolic) parameters
    lamb   = dst.Par(lamb_p, 'lamb')
    eta    = dst.Par(eta_p, 'eta')
    p      = dst.Par(p_p, 'p')
    nu     = dst.Par(nu_p, 'nu')
    rho    = dst.Par(rho_p, 'rho')
    
    # Compute nontrivial boundary equilibrium initial condition from parameters (see reference)
    
    
    b_0 = 0.0
    w_0 = p_p/nu_p
    # Declare symbolic variables    
    b = dst.Var('b')
    w = dst.Var('w')
    t = dst.Var('t')

    # Create Symbolic Quantity objects for definitions
    brhs = dst.Fun(lamb*w*b*((1+eta*b)**2)*(1-b) - b,[b,w],'brhs')
    wrhs = dst.Fun(p - nu*w*(1-rho*b) - lamb*w*b*((1+eta*b)**2),[b,w],'wrhs')

    F   = dst.Fun([brhs(b,w),wrhs(b,w)], [b,w], 'F')
    jac = dst.Fun(dst.Diff(F,[b,w]), [t,b,w], 'Jacobian')


    # Build Generator
    DSargs = dst.args(name='fairy_circles_ode')
    DSargs.fnspecs = [jac, brhs,wrhs]
    DSargs.varspecs = {b:brhs(b,w) ,
                       w:wrhs(b,w)}
    DSargs.pars = [lamb,eta,p,nu,rho]
    
    
    # Use eval method to get a float value from the symbolic definitions given in
    # terms of parameter values
    DSargs.ics = dst.args(b=b_0, w=w_0)
    return DSargs

 def make(self, dsargs):
     # use cache if available
     for gen_ver, prev_dsargs in self.used_dsargs.items():
         if dst.filteredDict(dsargs, 'name', True) == dst.filteredDict(prev_dsargs,
                                                               'name', True):
             # compare everything but the name, but check all up to final '_ver<X>'
             parts1 = dsargs.name.split('_')
             parts2 = prev_dsargs.name.split('_') # will have one more part
             if parts1 == parts2[:-2]:
                 print("Reusing identical build")
                 return dst.loadObjects(os.path.join(self.cwd, 'models',
                                                 prev_dsargs.name+'.sav'))[0]
     # no matches
     return self.build(dsargs)

def make_shooter():
    # no friction
    # cos(atan(x)) = 1/(sqrt(1+x^2))
    Fx_str = '0' # '-speed_fn()*cos(atan2(vy,vx))'
    Fy_str = '-10'

    DSargs = dst.args()
    DSargs.varspecs = {'vx': Fx_str, 'x': 'vx',
                       'vy': Fy_str, 'y': 'vy',
                       'Fx_out': 'Fx(x,y)', 'Fy_out': 'Fy(x,y)',
                       'speed': 'speed_fn(vx, vy)',
                       'bearing': '90-180*atan2(vy,vx)/pi'}

    auxfndict = {'Fx': (['x', 'y'], Fx_str),
                 'Fy': (['x', 'y'], Fy_str),
                 'speed_fn': (['vx', 'vy'], 'sqrt(vx*vx+vy*vy)'),
                 }
    DSargs.auxvars = ['Fx_out', 'Fy_out', 'speed', 'bearing']

    DSargs.fnspecs = auxfndict
    DSargs.algparams = {'init_step':0.001,
                        'max_step': 0.1,
                        'max_pts': 20000,
                        'maxevtpts': 2,
                        'refine': 5}

    ground_event = dst.Events.makeZeroCrossEvent('y', -1,
                                                 {'name': 'ground',
                                                  'eventtol': 1e-3,
                                                  'precise': True,
                                                  'term': True},
                                                 varnames=['y'],
                                                 targetlang='python')
    peak_event = dst.Events.makeZeroCrossEvent('vy', -1,
                                                 {'name': 'peak',
                                                  'eventtol': 1e-3,
                                                  'precise': True,
                                                  'term': False},
                                                 varnames=['vy'],
                                                 targetlang='python')
    DSargs.events = [ground_event, peak_event]
    DSargs.checklevel = 2
    DSargs.ics = {'x': 0, 'y': 0,
                  'vx': 0, 'vy': 0}
    DSargs.ics.update(make_vel_ics(5,20))
    DSargs.name = 'cannon'
    DSargs.tdomain = [0, 100000]
    DSargs.tdata = [0, 10]
    return dst.embed(dst.Generator.Vode_ODEsystem(DSargs))

    def __init__(self, cwd, model_name, name_base, gen_type, gen_version=0):
        """
        Internal utility to manage versions of Generator objects within single
        session.

        cwd = current working directory (string)

        Option to set known gen version # to reuse:
        Version 0 means not yet created.
        Works across saved and restarted sessions
        """
        self.cwd = cwd
        self.model_name = model_name
        self.name_base = name_base
        self.gen_version = gen_version
        self.gen_type = gen_type
        # keyed by version
        self.used_dsargs = {}
        self.logfile = os.path.join(self.cwd, 'models', 'gen_dsargs_log.sav')
        if os.path.exists(self.logfile):
            # reload previous list of versions
            self.used_dsargs = dst.loadObjects(self.logfile)[0]
        self.classes = {'vode': dst.Generator.Vode_ODEsystem,
                        'dopri': dst.Generator.Dopri_ODEsystem,
                        'radau': dst.Generator.Radau_ODEsystem,
                        'euler': dst.Generator.Euler_ODEsystem}
        self.targetlang = self._targetlangs[gen_type]

    def build(self, dsargs, is_stiff=False):
        # re-compute in case gen type has been changed
        self.targetlang = self._targetlangs[self.gen_type]
        if is_stiff and self.targetlang == 'python' and self.gen_type == 'vode':
            dsargs.algparams['stiff'] = True

        name = dsargs.name
        if self.gen_version == 0:
            self.gen_version = 1

        # assume it's sufficient to check if .sav file there rather than .so
        found_new = False
        while not found_new:
            filename = os.path.join(self.cwd, 'models', name + '_' + \
                                    self.gen_type + \
                                    '_ver%i'%self.gen_version+'.sav')
            if not os.path.exists(filename):
                found_new = True
            else:
                print(filename + ' already exists')
                self.gen_version += 1
        dsargs.name = name+'_'+self.gen_type+'_ver%i'%self.gen_version
        gen = self.classes[self.gen_type](dsargs)
        model = dst.embed(gen, name=self.model_name, dsi_name='gen')
        self.used_dsargs[self.gen_version] = dsargs.copy()
        self.save_gen(model, name)
        return model

def simulate(args):
    modelname,ptargs,tdomain,captcnt,captincr,icdict,pardict,vardict,varspecdict,fnspecdict = args
    dsargs = pdt.args()
    dsargs.name = modelname
    dsargs.ics = icdict
    dsargs.pars = pardict
    dsargs.tdata = tdomain
    #dsargs.vars = vardict
    dsargs.varspecs = varspecdict
    #dsargs.fnspecs = fnspecdict
    dsargs.algparams = {
        'init_step':captincr/10.0,
        'atol':0.1,
            }

    dsys = pdt.Generator.Vode_ODEsystem(dsargs)
    #dsys = pdt.Generator.Radau_ODEsystem(dsargs)
    traj = dsys.compute('demo')
    pts = traj.sample()

    rshape = (len(ptargs),captcnt)
    result = numpy.zeros(shape = rshape,dtype = numpy.float)
    result[0,:] = numpy.arange(tdomain[0],tdomain[1]+0.000000001,captincr)
    for timedx in range(result.shape[1]):
        itraj = traj(result[0,timedx])
        for targdx in range(1,result.shape[0]):
            result[targdx,timedx] = itraj[ptargs[targdx]]
    return result

def test_goal(mesh_pts, goal_tol=L2_tol):
    errors_array = error_pts(mesh_pts)
    max_error = np.max(errors_array)
    result = condition(max_error, goal_tol)
    return dst.args(result=result,
                    errors=errors_array,
                    max_error=max_error)

def construct_system( I, alpha, init=(0,0), T = 10 ):
    theta0, w0 = init
    args = pd.args( name = 'Pendulum' )
    args.pars = { 'alpha' : alpha, 'I' : I }
    args.varspecs = { 'theta' : 'w', 'w' : 'I - sin(theta) - alpha * w' } 
    args.ics = { 'theta' : theta0, 'w' : w0 }
    args.tdomain = [-T, T ]
    ode = pd.Generator.Vode_ODEsystem( args )
    return ode

def build_lin():
    # make local linear system spec
    if can_cache:
        print("I'm not building this model twice!")
    DSargs = dst.args(name='lin')
    xfn_str = '(x0+yfx*y - x)/taux'
    yfn_str = '(y0+xfy*x - y)/tauy'
    DSargs.varspecs = {'x': xfn_str, 'y': yfn_str}
    DSargs.xdomain = {'x': xdom, 'y': ydom}
    DSargs.pars = {'x0': xdom_half, 'y0': ydom_half,
                   'xfy': 1, 'yfx': 1,
                   'taux': 1, 'tauy': 1}
    DSargs.algparams = {'init_step':0.001,
                        'max_step': 0.001,
                        'max_pts': 10000}
    DSargs.checklevel = 0
    DSargs.tdata = [0, 10]
    DSargs.ics = {'x': xdom_half*1.1, 'y': ydom_half*1.1}
    DSargs.fnspecs = {'Jacobian': (['t', 'x', 'y'],
                                   """[[-1/taux, yfx/taux],
                                       [xfy/tauy, -1/tauy]]""")}
    return dst.embed(dst.Generator.Vode_ODEsystem(DSargs))

def map_workspace(con, pts, *args):
    """
    Returns a list of dictionaries, each representing the state of the
    calc_context workspace for each of the points given. Optional
    positional arguments will be passed first to the calc_context when
    calling it.

    This assumes the calc_context local_init accepts `pt` as an argument.
    """
    wseq = []
    for pt in pts:
        con(*args, pt=pt)
        wseq.append(dst.filteredDict(con.workspace.__dict__, ['_name'], neg=True))
    return wseq