Python sympy.solve() Examples

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def solve(jarvis, s):
    """
    Prints where expression equals zero
    -- Example:
        solve x**2 + 5*x + 3
        solve x + 3 = 5
    """
    x = sympy.Symbol('x')

    def _format(solutions):
        if solutions == 0:
            return "No solution!"
        ret = ''
        for count, point in enumerate(solutions):
            if x not in point:
                return "Please use 'x' in expression."
            x_value = point[x]
            ret += "{}. x: {}\n".format(count, x_value)
        return ret

    def _calc(expr):
        return sympy.solve(expr, x, dict=True)

    s = remove_equals(jarvis, s)
    calc(jarvis, s, calculator=_calc, formatter=_format, do_evalf=False) 

def calc(jarvis, s, calculator=sympy.sympify, formatter=None, do_evalf=True):
    s = format_expression(s)
    try:
        result = calculator(s)
    except sympy.SympifyError:
        jarvis.say("Error: Something is wrong with your expression", Fore.RED)
        return
    except NotImplementedError:
        jarvis.say("Sorry, cannot solve", Fore.RED)
        return

    if formatter is not None:
        result = formatter(result)

    if do_evalf:
        result = result.evalf()

    jarvis.say(str(result), Fore.BLUE) 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def testAlgebraInverse(self):
    dataset_objects = algorithmic_math.math_dataset_init(26)
    counter = 0
    for d in algorithmic_math.algebra_inverse(26, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      solve_var, expression = decoded_input.split(":")
      lhs, rhs = expression.split("=")

      # Solve for the solve-var.
      result = sympy.solve("%s-(%s)" % (lhs, rhs), solve_var)
      target_expression = dataset_objects.int_decoder(d["targets"])

      # Check that the target and sympy's solutions are equivalent.
      self.assertEqual(
          0, sympy.simplify(str(result[0]) + "-(%s)" % target_expression))
    self.assertEqual(counter, 10) 

def findomega(stab_fh):
    assert np.array_equal(np.shape(stab_fh), [2, 2]), 'Not 2x2 matrix...'
    omega = sympy.Symbol('omega')
    func = (sympy.exp(-1j * omega) - stab_fh[0, 0]) * (sympy.exp(-1j * omega) - stab_fh[1, 1]) - \
        stab_fh[0, 1] * stab_fh[1, 0]
    solsym = sympy.solve(func, omega)
    sol0 = complex(solsym[0])
    sol1 = complex(solsym[1])
    if sol0.real >= 0:
        sol = sol0
    elif sol1.real >= 0:
        sol = sol1
    else:
        print("Two roots with real part of same sign...")
        sol = sol0
    return sol 

def DFE(ode, disease_state):
    '''
    Returns the disease free equilibrium from an ode object

    Parameters
    ----------
    ode: :class:`.BaseOdeModel`
        a class object from pygom
    diseaseState: array like
        name of the disease states

    Returns
    -------
    e: array like
        disease free equilibrium

    '''

    eqn = ode.get_ode_eqn()
    index = ode.get_state_index(disease_state)
    states = [s for s in ode._iterStateList()]
    states_subs = {states[i]: 0 for i in index}
    eqn = eqn.subs(states_subs)

    DFE_solution = sympy.solve(eqn, states)
    if len(DFE_solution) == 0: DFE_solution = {}

    for s in states:
        if s not in states_subs.keys() and s not in DFE_solution.keys():
            DFE_solution.setdefault(s, 0)
    return DFE_solution 

def solve_movie_equations():
    from sympy import Eq, solve, symbols

    hash, title, year, release_group = symbols('hash title year release_group')
    format, audio_codec, resolution, video_codec = symbols('format audio_codec resolution video_codec')
    hearing_impaired = symbols('hearing_impaired')

    equations = [
        # hash is best
        Eq(hash, title + year + release_group + format + audio_codec + resolution + video_codec),

        # title counts for the most part in the total score
        Eq(title, year + release_group + format + audio_codec + resolution + video_codec + 1),

        # year is the second most important part
        Eq(year, release_group + format + audio_codec + resolution + video_codec + 1),

        # release group is the next most wanted match
        Eq(release_group, format + audio_codec + resolution + video_codec + 1),

        # format counts as much as audio_codec, resolution and video_codec
        Eq(format, audio_codec + resolution + video_codec),

        # audio_codec is more valuable than video_codec
        Eq(audio_codec, video_codec + 1),

        # resolution counts as much as video_codec
        Eq(resolution, video_codec),

        # video_codec is the least valuable match but counts more than the sum of all scoring increasing matches
        Eq(video_codec, hearing_impaired + 1),

        # hearing impaired is only used for score increasing, so put it to 1
        Eq(hearing_impaired, 1),
    ]

    return solve(equations, [hash, title, year, release_group, format, audio_codec, resolution, hearing_impaired,
                             video_codec]) 

def solve(eq, target, **kwargs):
    """
    Algebraically rearrange an Eq w.r.t. a given symbol.

    This is a wrapper around ``sympy.solve``.

    Parameters
    ----------
    eq : expr-like
        The equation to be rearranged.
    target : symbol
        The symbol w.r.t. which the equation is rearranged. May be a `Function`
        or any other symbolic object.
    **kwargs
        Symbolic optimizations applied while rearranging the equation. For more
        information. refer to ``sympy.solve.__doc__``.
    """
    if isinstance(eq, Eq):
        eq = eq.lhs - eq.rhs if eq.rhs != 0 else eq.lhs
    sols = []
    for e, t in zip(as_tuple(eq), as_tuple(target)):
        # Try first linear solver
        try:
            cc = linear_coeffs(e.evaluate, t)
            sols.append(-cc[1]/cc[0])
        except ValueError:
            warning("Equation is not affine w.r.t the target, falling back to standard"
                    "sympy.solve that may be slow")
            kwargs['rational'] = False  # Avoid float indices
            kwargs['simplify'] = False  # Do not attempt premature optimisation
            sols.append(sympy.solve(e.evaluate, t, **kwargs)[0])
    # We need to rebuild the vector/tensor as sympy.solve outputs a tuple of solutions
    if len(sols) > 1:
        return target.new_from_mat(sols)
    else:
        return sols[0] 

def test_solve(so):
    """
    Test that our solve produces the correct output and faster than sympy's
    default behavior for an affine equation (i.e. PDE time steppers).
    """
    grid = Grid((10, 10, 10))
    u = TimeFunction(name="u", grid=grid, time_order=2, space_order=so)
    v = Function(name="v", grid=grid, space_order=so)
    eq = u.dt2 - div(v * grad(u))

    # Standard sympy solve
    t0 = time.time()
    sol1 = sympy.solve(eq.evaluate, u.forward, rational=False, simplify=False)[0]
    t1 = time.time() - t0

    # Devito custom solve for linear equation in the target ax + b (most PDE tie steppers)
    t0 = time.time()
    sol2 = solve(eq.evaluate, u.forward)
    t12 = time.time() - t0

    diff = sympy.simplify(sol1 - sol2)
    # Difference can end up with super small coeffs with different evaluation
    # so zero out everything very small
    assert diff.xreplace({k: 0 if abs(k) < 1e-10 else k
                          for k in diff.atoms(sympy.Float)}) == 0
    # Make sure faster (actually much more than 10 for very complex cases)
    assert t12 < t1/10 

def test_events():
    # use bouncing ball to test events work

    # simulate in block diagram
    int_opts = block_diagram.DEFAULT_INTEGRATOR_OPTIONS.copy()
    int_opts['rtol'] = 1E-12
    int_opts['atol'] = 1E-15
    int_opts['nsteps'] = 1000
    int_opts['max_step'] = 2**-3
    x = x1, x2 = Array(dynamicsymbols('x_1:3'))
    mu, g = sp.symbols('mu g')
    constants = {mu: 0.8, g: 9.81}
    ic = np.r_[10, 15]
    sys = SwitchedSystem(
        x1, Array([0]),
        state_equations=r_[x2, -g],
        state_update_equation=r_[sp.Abs(x1), -mu*x2],
        state=x,
        constants_values=constants,
        initial_condition=ic
    )
    bd = BlockDiagram(sys)
    res = bd.simulate(5, integrator_options=int_opts)

    # compute actual impact time
    tvar = dynamicsymbols._t
    impact_eq = (x2*tvar - g*tvar**2/2 + x1).subs(
        {x1: ic[0], x2: ic[1], g: 9.81}
    )
    t_impact = sp.solve(impact_eq, tvar)[-1]

    # make sure simulation actually changes velocity sign around impact
    abs_diff_impact = np.abs(res.t - t_impact)
    impact_idx = np.where(abs_diff_impact == np.min(abs_diff_impact))[0]
    assert np.sign(res.x[impact_idx-1, 1]) != np.sign(res.x[impact_idx+1, 1]) 

def equilibrium_points(self, input_=None):
        return sp.solve(self.state_equation, self.state, dict=True) 

def compute_fg(self):
        n_states = len(self.x)
        n_vars = len(self.v)
        n_eqs = len(self.eqs)
        if n_states + n_vars != n_eqs:
            raise RuntimeError('# states: {:d} + # variables: {:d} != # equations {:d}'.format(
                n_states, n_vars, n_eqs))
        lhs = list(self.x.diff(self.t)) + list(self.v) + list(self.y)
        fg_sol = sympy.solve(self.eqs, lhs, dict=True)[0]
        self.f = self.x.diff(self.t).subs(fg_sol)
        assert len(self.x) == len(self.f)
        self.g = self.y.subs(fg_sol)
        assert len(self.y) == len(self.g) 

def test_solveset():
	x = Symbol('x')
	A = Matrix([[x,2,x*x],[4,5,x],[x,8,9]])
	
	solns = solve(det(A), x)
	solns_set = list(solveset(det(A), x))
	
	print(solns)
	print('\n')
	print(solns_set)
	
	print('\n\n\n')
	print((solns[0]))
	print('\n')
	print((solns_set[0]))
	
	soln_sub = solns[0].subs(x, 1)
	solnset_sub = solns_set[0].subs(x, 1)
	
	s1 = soln_sub.evalf()
	s1set = solnset_sub.evalf()
	
	s2set = solns_set[1].subs(x, 1).evalf()
	
	print(s1)
	print(s1set)
	print(s2set) 

def solow_residual(self):
        """
        Symbolic expression for the Solow residual which is used as a
        measure of technology.

        :getter: Return the symbolic expression.
        :type: sym.Basic

        """
        return sym.solve(Y - self.output, A)[0] 

def solve_y(s):
    if 'y' in s:
        y = sympy.Symbol('y')
        try:
            results = sympy.solve(s, y)
        except NotImplementedError:
            return 'unknown'
        if len(results) == 0:
            return '0'
        else:
            return results[0]
    else:
        return solve_y("({}) -y".format(s)) 

def equations(jarvis, term):
    """
    Solves linear equations system

    Use variables: a, b, c, ..., x, y,z

    Example:

    ~> Hi, what can I do for you?
    equations
    1. Equation: x**2 + 2y - z = 6
    2. Equation: (x-1)(y-1) = 0
    3. Equation: y**2 - x -10 = y**2 -y
    4. Equation:
    [{x: -9, y: 1, z: 77}, {x: 1, y: 11, z: 17}]

    """
    a, b, c, d, e, f, g, h, i, j, k, l, m = sympy.symbols(
        'a,b,c,d,e,f,g,h,i,j,k,l,m')
    n, o, p, q, r, s, t, u, v, w, x, y, z = sympy.symbols(
        'n,o,p,q,r,s,t,u,v,w,x,y,z')

    equations = []
    count = 1
    user_input = jarvis.input('{}. Equation: '.format(count))
    while user_input != '':
        count += 1
        user_input = format_expression(user_input)
        user_input = remove_equals(jarvis, user_input)
        equations.append(user_input)
        user_input = jarvis.input('{}. Equation: '.format(count))

    calc(
        jarvis,
        term,
        calculator=lambda expr: sympy.solve(
            equations,
            dict=True)) 

def solow_residual(self):
        """
        Symbolic expression for the Solow residual which is used as a measure
        of technology.

        :getter: Return the symbolic expression.
        :type: sympy.Basic

        """
        return sym.solve(Y - self.output, A)[0] 

def testLinearSystem(self, degree):
    for _ in range(100):  # test a few times
      target = [random.randint(-100, 100) for _ in range(degree)]
      variables = [sympy.Symbol(chr(ord('a') + i)) for i in range(degree)]
      system = linear_system.linear_system(
          variables=variables,
          solutions=target,
          entropy=10.0)
      solved = sympy.solve(system, variables)
      solved = [solved[symbol] for symbol in variables]
      self.assertEqual(target, solved) 

def enforce_volumetric_constraint(left_volume, inshape):
        left_symbols = set()
        if is_symbolic(left_volume):
            left_symbols = left_volume.free_symbols
        # Generally, we want to solve for right in terms of left. But this
        # is kinda annoying actually.
        shape = list(inshape)

        # Handling when reshape is given 0 instead of actual input
        # input tensor shape: [4, 3, 2], reshape:[0, -1], output tensor shape: [4, 6]
        if shape.count(-1) > 1:
            raise ValueError(
                "Reshape op supports only one dimension to be -1. Given {}".format(
                    shape.count(-1)
                )
            )

        infer_dim_index = shape.index(-1) if -1 in shape else None
        right_volume = 1
        for i in shape:
            if i != -1:
                right_volume = right_volume * i

        if infer_dim_index:
            shape[infer_dim_index] = left_volume // right_volume

        if not is_symbolic(right_volume):
            return shape

        constraints = [left_volume - right_volume]
        solve_for = [s for s in shape if is_symbolic(s)]

        for rightsym in solve_for:
            sol = sm.solve(constraints, [rightsym], dict=True)
            if not isinstance(sol, list):
                sol = [sol]
            # look for an acceptable solution
            for s in sol:
                if 0 in s.values():
                    continue
                for i in range(len(shape)):
                    if shape[i] in s:
                        v = s[shape[i]]
                        if len(v.free_symbols - left_symbols) > 0:
                            continue
                        try:
                            shape[i] = int(v)
                        except:
                            shape[i] = v
        return shape 

def control_systems(request):
    ct_sys, ref = request.param
    Ac, Bc, Cc = ct_sys.data
    Dc = np.zeros((Cc.shape[0], 1))

    Q = np.eye(Ac.shape[0])
    R = np.eye(Bc.shape[1] if len(Bc.shape) > 1 else 1)

    Sc = linalg.solve_continuous_are(Ac, Bc.reshape(-1, 1), Q, R,)
    Kc = linalg.solve(R, Bc.T @ Sc).reshape(1, -1)
    ct_ctr = LTISystem(Kc)

    evals = np.sort(np.abs(
        linalg.eig(Ac, left=False, right=False, check_finite=False)
    ))
    dT = 1/(2*evals[-1])

    Tsim = (8/np.min(evals[~np.isclose(evals, 0)])
            if np.sum(np.isclose(evals[np.nonzero(evals)], 0)) > 0
            else 8
            )

    dt_data = signal.cont2discrete((Ac, Bc.reshape(-1, 1), Cc, Dc), dT)
    Ad, Bd, Cd, Dd = dt_data[:-1]
    Sd = linalg.solve_discrete_are(Ad, Bd.reshape(-1, 1), Q, R,)
    Kd = linalg.solve(Bd.T @ Sd @ Bd + R, Bd.T @ Sd @ Ad)

    dt_sys = LTISystem(Ad, Bd, dt=dT)
    dt_sys.initial_condition = ct_sys.initial_condition
    dt_ctr = LTISystem(Kd, dt=dT)

    yield ct_sys, ct_ctr, dt_sys, dt_ctr, ref, Tsim 

def limit(jarvis, s):
    """
    Prints limit to +/- infinity or to number +-. Use 'x' as variable.
    -- Examples:
        limit 1/x
        limit @1 1/(1-x)
        limit @1 @2 1/((1-x)(2-x))
    """
    def try_limit(term, x, to, directory=''):
        try:
            return sympy.Limit(term, x, to, directory).doit()
        except sympy.SympifyError:
            return 'Error'
        except NotImplementedError:
            return "Sorry, cannot solve..."

    if s == '':
        jarvis.say("Usage: limit TERM")
        return

    s_split = s.split()
    limit_to = []
    term = ""
    for token in s_split:
        if token[0] == '@':
            if token[1:].isnumeric():
                limit_to.append(int(token[1:]))
            else:
                jarvis.say("Error: {} Not a number".format(
                    token[1:]), Fore.RED)
        else:
            term += token

    term = remove_equals(jarvis, term)
    term = format_expression(term)

    try:
        term = solve_y(term)
    except (sympy.SympifyError, TypeError):
        jarvis.say('Error, not a valid term')
        return

    x = sympy.Symbol('x')

    # infinity:
    jarvis.say("lim ->  ∞\t= {}".format(try_limit(term,
                                                  x, +sympy.S.Infinity)), Fore.BLUE)
    jarvis.say("lim -> -∞\t= {}".format(try_limit(term,
                                                  x, -sympy.S.Infinity)), Fore.BLUE)

    for limit in limit_to:
        limit_plus = try_limit(term, x, limit, directory="+")
        limit_minus = try_limit(term, x, limit, directory="-")

        jarvis.say("lim -> {}(+)\t= {}".format(limit, limit_plus), Fore.BLUE)
        jarvis.say("lim -> {}(-)\t= {}".format(limit, limit_minus), Fore.BLUE) 

def solve_episode_equations():
    from sympy import Eq, solve, symbols

    hash, series, year, season, episode, release_group = symbols('hash series year season episode release_group')
    format, audio_codec, resolution, video_codec = symbols('format audio_codec resolution video_codec')
    hearing_impaired = symbols('hearing_impaired')

    equations = [
        # hash is best
        Eq(hash, series + year + season + episode + release_group + format + audio_codec + resolution + video_codec),

        # series counts for the most part in the total score
        Eq(series, year + season + episode + release_group + format + audio_codec + resolution + video_codec + 1),

        # year is the second most important part
        Eq(year, season + episode + release_group + format + audio_codec + resolution + video_codec + 1),

        # season is important too
        Eq(season, release_group + format + audio_codec + resolution + video_codec + 1),

        # episode is equally important to season
        Eq(episode, season),

        # release group is the next most wanted match
        Eq(release_group, format + audio_codec + resolution + video_codec + 1),

        # format counts as much as audio_codec, resolution and video_codec
        Eq(format, audio_codec + resolution + video_codec),

        # audio_codec is more valuable than video_codec
        Eq(audio_codec, video_codec + 1),

        # resolution counts as much as video_codec
        Eq(resolution, video_codec),

        # video_codec is the least valuable match but counts more than the sum of all scoring increasing matches
        Eq(video_codec, hearing_impaired + 1),

        # hearing impaired is only used for score increasing, so put it to 1
        Eq(hearing_impaired, 1),
    ]

    return solve(equations, [hash, series, year, season, episode, release_group, format, audio_codec, resolution,
                             hearing_impaired, video_codec])