def probability(self, condition, **kwargs):
     cond_inv = False
     if isinstance(condition, Ne):
         condition = Eq(condition.args[0], condition.args[1])
         cond_inv = True
     expr = condition.lhs - condition.rhs
     rvs = random_symbols(expr)
     z = Dummy('z', real=True, Finite=True)
     dens = self.compute_density(expr)
     if any([pspace(rv).is_Continuous for rv in rvs]):
         from sympy.stats.crv import (ContinuousDistributionHandmade,
             SingleContinuousPSpace)
         if expr in self.values:
             # Marginalize all other random symbols out of the density
             randomsymbols = tuple(set(self.values) - frozenset([expr]))
             symbols = tuple(rs.symbol for rs in randomsymbols)
             pdf = self.domain.integrate(self.pdf, symbols, **kwargs)
             return Lambda(expr.symbol, pdf)
         dens = ContinuousDistributionHandmade(dens)
         space = SingleContinuousPSpace(z, dens)
         result = space.probability(condition.__class__(space.value, 0))
     else:
         from sympy.stats.drv import (DiscreteDistributionHandmade,
             SingleDiscretePSpace)
         dens = DiscreteDistributionHandmade(dens)
         space = SingleDiscretePSpace(z, dens)
         result = space.probability(condition.__class__(space.value, 0))
     return result if not cond_inv else S.One - result

    def probability(self, condition, **kwargs):
        z = Dummy('z', real=True, finite=True)
        cond_inv = False
        if isinstance(condition, Ne):
            condition = Eq(condition.args[0], condition.args[1])
            cond_inv = True
        # Univariate case can be handled by where
        try:
            domain = self.where(condition)
            rv = [rv for rv in self.values if rv.symbol == domain.symbol][0]
            # Integrate out all other random variables
            pdf = self.compute_density(rv, **kwargs)
            # return S.Zero if `domain` is empty set
            if domain.set is S.EmptySet or isinstance(domain.set, FiniteSet):
                return S.Zero if not cond_inv else S.One
            if isinstance(domain.set, Union):
                return sum(
                     Integral(pdf(z), (z, subset), **kwargs) for subset in
                     domain.set.args if isinstance(subset, Interval))
            # Integrate out the last variable over the special domain
            return Integral(pdf(z), (z, domain.set), **kwargs)

        # Other cases can be turned into univariate case
        # by computing a density handled by density computation
        except NotImplementedError:
            from sympy.stats.rv import density
            expr = condition.lhs - condition.rhs
            dens = density(expr, **kwargs)
            if not isinstance(dens, ContinuousDistribution):
                dens = ContinuousDistributionHandmade(dens)
            # Turn problem into univariate case
            space = SingleContinuousPSpace(z, dens)
            result = space.probability(condition.__class__(space.value, 0))
            return result if not cond_inv else S.One - result

def solver():
    vals = {i: base_eq(i) for i in range(1, 12)}
    coefficients = a, b, c, d, e, f, g, h, i, j, k = symbols(','.join('abcdefghijk'))
    (x, n), lhs = symbols('x n'), a
    answer = 1

    for i, element in enumerate(coefficients[1:]):
        lhs = lhs * n + element
        equation = Eq(x, lhs)
        results = solve([equation.subs(dict(n=t, x=vals[t])) for
                         t in range(1, i+3)])
        results[n] = i + 3
        if element != k:
            answer += equation.subs(results).args[1]
    return answer

def solve_it(eqn, knowns, x):
    LHS, RHS = eqn['sympy'].split('=')
    e = Eq(sympify(LHS), sympify(RHS))

    voi = Symbol(x['symbol']) # upgrade sympy to 0.7.6 if this throws a unicode error

    vals = dict()
    for k in knowns:
        vals[Symbol(k['symbol'])] = k['value']
    for c in eqn['constants']:
        vals[Symbol(c['symbol'])] = c['value']

    answers = solve(e.subs(vals),voi)
    answer = answers[0]

    return answer.evalf()

 def probability(self, condition):
     complement = isinstance(condition, Ne)
     if complement:
         condition = Eq(condition.args[0], condition.args[1])
     try:
         _domain = self.where(condition).set
         if condition == False or _domain is S.EmptySet:
             return S.Zero
         if condition == True or _domain == self.domain.set:
             return S.One
         prob = self.eval_prob(_domain)
     except NotImplementedError:
         from sympy.stats.rv import density
         expr = condition.lhs - condition.rhs
         dens = density(expr)
         if not isinstance(dens, DiscreteDistribution):
             dens = DiscreteDistributionHandmade(dens)
         z = Dummy('z', real = True)
         space = SingleDiscretePSpace(z, dens)
         prob = space.probability(condition.__class__(space.value, 0))
     if (prob == None):
         prob = Probability(condition)
     return prob if not complement else S.One - prob

def test_equality_subs1():
    f = Function("f")
    x = abc.x
    eq = Eq(f(x)**2, x)
    res = Eq(Integer(16), x)
    assert eq.subs(f(x), 4) == res

def test_equality_subs2():
    f = Function('f')
    eq = Eq(f(x)**2, 16)
    assert bool(eq.subs(f(x), 3)) is False
    assert bool(eq.subs(f(x), 4)) is True

    def set_free_surface(self, d, b, side, algo='robertsson'):
        """
        set free surface boundary condition to boundary d, at index b
        :param indices: list of indices, e.g. [t,x,y,z] for 3D
        :param d: direction of the boundary surface normal
        :param b: location of the boundary (index)
        :param algo: which algorithm to use to compute ghost cells
        algo == 'robertsson' [1]: setting all velocities at ghost cells to zero
        algo == 'levander' [2]: only valid for 4th spatial order. using 2nd order FD approximation for velocities
        side: lower boundary (0) or upper boundary (1)
        e.g. set_free_surface([t,x,y,z],1,2,0)
        set y-z plane at x=2 to be lower free surface
        ghost cells are calculated using reflection of stress fields
        store the code to populate ghost cells to self.bc
        [1] Robertsson, Johan OA. "A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography." Geophysics 61.6 (1996): 1921-1934.
        [2] Levander, Alan R. "Fourth-order finite-difference P-SV seismograms." Geophysics 53.11 (1988): 1425-1436.
        """
        idx = list(self.indices)

        if d not in self.direction:
            if (not algo == 'levander') or (not self.direction[0] == self.direction[1]):
                # shear stress, e.g. Tyz no need to recalculate at x boundary (only depends on dV/dz and dW/dy)
                self.bc[d][side] = []
                return
            else:
                # normal stress, need to recalcuate Tyy, Tzz at x boundary
                expr = self.dt
                derivatives = get_all_objects(expr, DDerivative)
                for deriv in derivatives:
                    if deriv.var == idx[d]:
                        # replacing dx at x boundary with dy, dz terms
                        expr2 = self.sfields[d].dt
                        deriv_0 = deriv
                        deriv_sub = solve(expr2, deriv)[0]
                        break
                expr = expr.subs(deriv_0, deriv_sub)
                derivatives = get_all_objects(expr, DDerivative)
                # substitution dictionary
                dict1 = {}
                for deriv in derivatives:
                    dict1[deriv] = deriv.fd[4]
                expr = expr.subs(dict1)
                eq = Eq(self.d[0][1].fd[2], expr)
                eq = eq.subs(idx[d], b)
                t = idx[0]
                idx[0] = t+hf
                idx[d] = b
                # eq = eq.subs(t, t+hf)
                # idx[0] = t+1
                # idx[d] = b
                # solve for Txx(t+1/2)
                lhs = self[idx]
                rhs = solve(eq, lhs)[0]
                rhs = self.align(rhs)
                # change t+1/2 to t+1
                lhs = lhs.subs(t, t+hf)
                rhs = rhs.subs(t, t+hf)
                eq2 = Eq(lhs, rhs)
                self.bc[d][side] = [eq2]
                return

        # use anti-symmetry to ensure stress at boundary=0
        # this applies for all algorithms

        idx = list(self.indices)  # ghost cell
        idx2 = list(self.indices)  # cell inside domain

        if not self.staggered[d]:
            # if not staggered, assign T[d]=0, assign T[d-1]=-T[d+1]
            idx[d] = b
            idx2[d] = b
            eq1 = Eq(self[idx])
        else:
            # if staggered, assign T[d-1/2]=T[d+1/2], assign T[d-3/2]=T[d+3/2]
            idx[d] = b - (1-side)
            idx2[d] = idx[d] + (-1)**side
            eq1 = Eq(self[idx], -self[idx2])
        eq1 = eq1.subs(idx[0], idx[0]+1)
        self.bc[d][side] = [eq1]

        for depth in range(self.order[d]/2-1):
            # populate ghost cells
            idx[d] -= (-1)**side
            idx2[d] += (-1)**side
            eq = Eq(self[idx], -self[idx2])
            # change t to t+1
            eq = eq.subs(idx[0], idx[0]+1)
            self.bc[d][side].append(eq)

#!/usr/bin/python

from sympy import latex
from sympy import solve
from sympy import symbols

from sympy import Eq
from sympy import Rational

a, d, t, v, v0 = symbols('a d t v v0')

kin1 = Eq(d, v0 * t + Rational(1,2) * a * t**2)
kin2 = Eq(v, v0 + a * t)

print ""
print "The basic equations for one-dimensional motion are:"
print ""
print "    {}".format(latex(kin1))
print "    {}".format(latex(kin2))
print ""
print "How long does it take an object released from rest to fall five meters?"
print ""

values = dict()
values[d] = -5
values[v0] = 0
values[a] = -9.81

print "The parameters for this problem are:"
print ""
for param, value in values.iteritems():

from sympy import symbols, Eq, Derivative, Integral, log, solve, exp, ccode

globals().update(symbs)  # see common.py: x, Y, Z, k_f, ...

eqs = []

# rate of x

rate_expr = k_f * (Y - x) * (Z - x)
rate_eq = Eq(Derivative(x, t), rate_expr)
eqs.append(rate_eq)

integrand = k_f / rate_expr
inte_eq_lhs = Integral(integrand.subs({x: chi}), (chi, 0, x))
inte_eq_rhs = Integral(k_f, (tau, 0, t))
inte_eq = Eq(inte_eq_lhs, inte_eq_rhs)
eqs.append(inte_eq)

expl_in_x_eq = inte_eq.doit(manual=True).simplify()
eqs.append(expl_in_x_eq)

expl_in_t_eq = Eq(x, solve(expl_in_x_eq, x)[0])

alt_expl_in_t = Y * (1 - exp(k_f * t * (Z - Y))) / (Y / Z - exp(k_f * t * (Z - Y)))
assert (alt_expl_in_t - expl_in_t_eq.rhs).simplify() == 0
alt_expl_in_t_eq = Eq(x, alt_expl_in_t)
eqs.append(alt_expl_in_t_eq)


def main():
    # GENERATES WITHOUT ARGUMENTS: irrev_binary_1.tex irrev_binary_2.tex irrev_binary_rate.tex irrev_binary_k_b.c irrev_binary_K_eq.c

 def solve_constants(eq, ics, d_ics):
     udiff = Eq(d_ics[0][1], eq.rhs.diff(t))
     system = [eq.subs(ics), udiff.subs(t, 0)]
     consts = solve(system, [C1, C2])
     return eq.subs(consts)

#!/usr/bin/python

from sympy import latex
from sympy import solve
from sympy import symbols

from sympy import Eq
from sympy import Rational

P, V, n, R, T = symbols('P V n R T')

igl = Eq(P * V, n * R * T)

values = {
    R : 8.3144621,
    T : 273,
    P : 101300,
    n : 1
}

answers = solve(igl.subs(values),V)
answer = answers[0]

print ""
print "The ideal gas law is ",latex(igl)
print ""
print "Where the gas constant R is {} in SI units".format(values[R])
print "Standard temperature is {} Kelvin, or 0 degrees Celsius".format(values[T])
print "Standard pressure is {} pascals, or 1 atm".format(values[P])
print ""
print "So, at standard temperature and pressure, one mole of an ideal gas"

def test_Equality_rewrite_as_Add():
    eq = Eq(x + y, y - x)
    assert eq.rewrite(Add) == 2*x
    assert eq.rewrite(Add, evaluate=None).args == (x, x, y, -y)
    assert eq.rewrite(Add, evaluate=False).args == (x, y, x, -y)

def test_equality_subs2():
    f = Function("f")
    x = abc.x
    eq = Eq(f(x)**2, 16)
    assert bool(eq.subs(f(x), 3)) == False
    assert bool(eq.subs(f(x), 4)) == True

)

globals().update(get_symbs()) # see common.py: x, Y, Z, k_f, ...

subs = {}
eqs = []

# rate of x
rate_expr = k_f*Y*(Z-x) - k_b*x
rate_eq = Eq(Derivative(x,t), rate_expr)
eqs.append(rate_eq)

integrand = 1/rate_expr
inte_eq_lhs = Integral(integrand.subs({x: chi}), (chi,0,x))
inte_eq_rhs = Integral(1, (tau,0,t))
inte_eq = Eq(inte_eq_lhs, inte_eq_rhs)
eqs.append(inte_eq)

expl_in_x_eq = inte_eq.doit().simplify()
eqs.append(expl_in_x_eq)

expl_in_t_eq = Eq(x, solve(expl_in_x_eq, x)[0])
eqs.append(expl_in_t_eq)

alt_expl_in_t = Z*k_f*Y/(k_f*Y+k_b)*(1-exp(-t*(k_f*Y+k_b)))
assert (expl_in_t_eq.rhs - alt_expl_in_t).simplify() == 0

alt_expl_in_t_eq = Eq(x, alt_expl_in_t)
eqs.append(alt_expl_in_t_eq)

###################################
##       4. Adiabatic compression
###################################

p_i   = Symbol("p_i", positive=True)
V_i   = Symbol("V_i", positive=True)
tau_i = Symbol("tau_i", positive=True)
N_i   = Symbol("N_i", positive=True)
idealgaslaw_i = Eq( p_i*V_i, N_i*tau_i)

p_f   = Symbol("p_f", positive=True)
V_f   = Symbol("V_f", positive=True)
tau_f = Symbol("tau_f", positive=True)
idealgaslaw_f = Eq( p_f*V_f, N_i*tau_f)

adia_tV = Eq( tau_i*V_i**(gamma-1) , tau_f*V_f**(gamma-1) )

##############################
#         (a)
##############################

roomtemp_K = KCconv.subs(T_C,20).rhs # room temperature in Kelvin

Prob0104ans = adia_tV.subs(gamma,1.4).subs(V_f,1).subs(V_i,15).subs(tau_i, roomtemp_K) # answer to Problem 4 of Chapter 1

Prob0104ans = N( Prob0104ans.lhs) # 866.016969686253 K 
Prob0104ansC = solve( KCconv.subs( T_K, Prob0104ans), T_C )[0] # 592.866969686253 C 
solve( FCconv.subs( T_C, Prob0104ansC ), T_F)[0] # 1099.16054543526 F

##############################
#         (b)