def arc_length(self, start, stop, q_type = None):
     """
     Parameters:
     ---------
     q_type:  string
         'exact', 'integral', 'numeric'
     start, stop: numeric/str
         These are the limits of integration
     """
     
     start = sym.sympify(start)
     stop = sym.sympify(stop)        
     
        
     
     if q_type is None:
         q_type = self.q_type        
     
     u = self.u
     I = sym.Integral(self.ds, (u, start, stop))
     
     if q_type == 'integral':
         return I
     elif q_type == 'exact':
         return I.doit()
     else:
         f = make_func(str(self.ds), func_params=(self.ind_var), func_type='numpy')
         value = mp.quad(f, [float(start.evalf()), float(stop.evalf())])
         return value

def test_nsimplify():
    x = Symbol("x")
    assert nsimplify(0) == 0
    assert nsimplify(-1) == -1
    assert nsimplify(1) == 1
    assert nsimplify(1+x) == 1+x
    assert nsimplify(2.7) == Rational(27, 10)
    assert nsimplify(1-GoldenRatio) == (1-sqrt(5))/2
    assert nsimplify((1+sqrt(5))/4, [GoldenRatio]) == GoldenRatio/2
    assert nsimplify(2/GoldenRatio, [GoldenRatio]) == 2*GoldenRatio - 2
    assert nsimplify(exp(5*pi*I/3, evaluate=False)) == sympify('1/2 - sqrt(3)*I/2')
    assert nsimplify(sin(3*pi/5, evaluate=False)) == sympify('sqrt(sqrt(5)/8 + 5/8)')
    assert nsimplify(sqrt(atan('1', evaluate=False))*(2+I), [pi]) == sqrt(pi) + sqrt(pi)/2*I
    assert nsimplify(2 + exp(2*atan('1/4')*I)) == sympify('49/17 + 8*I/17')
    assert nsimplify(pi, tolerance=0.01) == Rational(22, 7)
    assert nsimplify(pi, tolerance=0.001) == Rational(355, 113)
    assert nsimplify(0.33333, tolerance=1e-4) == Rational(1, 3)
    assert nsimplify(2.0**(1/3.), tolerance=0.001) == Rational(635, 504)
    assert nsimplify(2.0**(1/3.), tolerance=0.001, full=True) == 2**Rational(1, 3)
    assert nsimplify(x + .5, rational=True) == Rational(1, 2) + x
    assert nsimplify(1/.3 + x, rational=True) == Rational(10, 3) + x
    assert nsimplify(log(3).n(), rational=True) == \
           sympify('109861228866811/100000000000000')
    assert nsimplify(Float(0.272198261287950), [pi,log(2)]) == pi*log(2)/8
    assert nsimplify(Float(0.272198261287950).n(3), [pi,log(2)]) == \
        -pi/4 - log(2) + S(7)/4
    assert nsimplify(x/7.0) == x/7
    assert nsimplify(pi/1e2) == pi/100
    assert nsimplify(pi/1e2, rational=False) == pi/100.0
    assert nsimplify(pi/1e-7) == 10000000*pi

def test_product_basic():
    H, T = 'H', 'T'
    unit_line = Interval(0, 1)
    d6 = FiniteSet(1, 2, 3, 4, 5, 6)
    d4 = FiniteSet(1, 2, 3, 4)
    coin = FiniteSet(H, T)

    square = unit_line * unit_line

    assert (0, 0) in square
    assert 0 not in square
    assert (H, T) in coin ** 2
    assert (.5, .5, .5) in square * unit_line
    assert (H, 3, 3) in coin * d6* d6
    HH, TT = sympify(H), sympify(T)
    assert set(coin**2) == set(((HH, HH), (HH, TT), (TT, HH), (TT, TT)))

    assert (d4*d4).is_subset(d6*d6)

    assert square.complement(Interval(-oo, oo)*Interval(-oo, oo)) == Union(
        (Interval(-oo, 0, True, True) +
         Interval(1, oo, True, True))*Interval(-oo, oo),
         Interval(-oo, oo)*(Interval(-oo, 0, True, True) +
                  Interval(1, oo, True, True)))

    assert (Interval(-5, 5)**3).is_subset(Interval(-10, 10)**3)
    assert not (Interval(-10, 10)**3).is_subset(Interval(-5, 5)**3)
    assert not (Interval(-5, 5)**2).is_subset(Interval(-10, 10)**3)

    assert (Interval(.2, .5)*FiniteSet(.5)).is_subset(square)  # segment in square

    assert len(coin*coin*coin) == 8
    assert len(S.EmptySet*S.EmptySet) == 0
    assert len(S.EmptySet*coin) == 0
    raises(TypeError, lambda: len(coin*Interval(0, 2)))

def test_product_basic():
    H,T = 'H', 'T'
    unit_line = Interval(0,1)
    d6 = FiniteSet(1,2,3,4,5,6)
    d4 = FiniteSet(1,2,3,4)
    coin = FiniteSet(H, T)

    square = unit_line * unit_line

    assert (0,0) in square
    assert 0 not in square
    assert (H, T) in coin ** 2
    assert (.5,.5,.5) in square * unit_line
    assert (H, 3, 3) in coin * d6* d6
    HH, TT = sympify(H), sympify(T)
    assert set(coin**2) == set(((HH, HH), (HH, TT), (TT, HH), (TT, TT)))

    assert (d6*d6).subset(d4*d4)

    inf, neginf = S.Infinity, S.NegativeInfinity
    assert square.complement == Union(
       Interval(0,1) * (Interval(neginf,0,True,True)+Interval(1,inf,True,True)),
       (Interval(neginf,0,True,True)+Interval(1,inf,True,True))*Interval(0,1),
       ((Interval(neginf,0,True,True) + Interval(1,inf, True, True))
                * (Interval(neginf,0,True,True) + Interval(1,inf, True,True))))

    assert (Interval(-10,10)**3).subset(Interval(-5,5)**3)
    assert not (Interval(-5,5)**3).subset(Interval(-10,10)**3)
    assert not (Interval(-10,10)**2).subset(Interval(-5,5)**3)

    assert square.subset(Interval(.2,.5)*FiniteSet(.5)) # segment in square

def field_isomorphism(a, b, **args):
    """Construct an isomorphism between two number fields. """
    a, b = sympify(a), sympify(b)

    if not a.is_AlgebraicNumber:
        a = AlgebraicNumber(a)

    if not b.is_AlgebraicNumber:
        b = AlgebraicNumber(b)

    if a == b:
        return a.coeffs()

    n = a.minpoly.degree()
    m = b.minpoly.degree()

    if n == 1:
        return [a.root]

    if m % n != 0:
        return None

    if args.get('fast', True):
        try:
            result = field_isomorphism_pslq(a, b)

            if result is not None:
                return result
        except NotImplementedError:
            pass

    return field_isomorphism_factor(a, b)

def hyperfocal_distance(f, N, c):
    """

    Parameters
    ==========
    f: sympifiable
    Focal length of a given lens

    N: sympifiable
    F-number of a given lens

    c: sympifiable
    Circle of Confusion (CoC) of a given image format

    Example
    =======
    >>> from sympy.physics.optics import hyperfocal_distance
    >>> from sympy.abc import f, N, c
    >>> round(hyperfocal_distance(f = 0.5, N = 8, c = 0.0033), 2)
    9.47
    """

    f = sympify(f)
    N = sympify(N)
    c = sympify(c)

    return (1/(N * c))*(f**2)

 def check_output(self, want, got, optionflags):
     if IGNORE_OUTPUT & optionflags:
         return True
     # When writing tests we sometimes want to assure ourselves that the
     # results are _not_ equal
     wanted_tf = not (NOT_EQUAL & optionflags)
     # Strip dtype
     if IGNORE_DTYPE & optionflags:
         want = ignore_dtype(want)
         got = ignore_dtype(got)
     # Strip array repr from got and want if requested
     if STRIP_ARRAY_REPR & optionflags:
         # STRIP_ARRAY_REPR only matches for a line containing *only* an
         # array repr.  Use IGNORE_DTYPE to ignore a dtype specifier embedded
         # within a more complex line.
         want = strip_array_repr(want)
         got = strip_array_repr(got)
     # If testing floating point, round to required number of digits
     if optionflags & (FP_4DP | FP_6DP):
         if optionflags & FP_4DP:
             dp = 4
         elif optionflags & FP_6DP:
             dp = 6
         want = round_numbers(want, dp)
         got = round_numbers(got, dp)
     # Are the strings equal when run through sympy?
     if SYMPY_EQUAL & optionflags:
         from sympy import sympify
         res = sympify(want) == sympify(got)
         return res == wanted_tf
     # Pass tests through two-pass numpy checker
     res = NumpyOutputChecker.check_output(self, want, got, optionflags)
     # Return True if we wanted True and got True, or if we wanted False and
     # got False
     return res == wanted_tf

def sympify(arg, real=False, positive=None, cache=True):
    """Create a sympy expression."""

    if isinstance(arg, (sym.symbol.Symbol, sym.symbol.Expr)):
        return arg

    # Why doesn't sympy do this?
    if isinstance(arg, complex):
        re = sym.sympify(arg.real, rational=True)
        im = sym.sympify(arg.imag, rational=True)
        if im == 1.0:
            arg = re + sym.I
        else:
            arg = re + sym.I * im
        return arg

    if isinstance(arg, str):
        # Sympy considers E1 to be the generalized exponential integral.
        # N is for numerical evaluation.
        if symbol_pattern.match(arg) is not None:
            return symbol(arg, real=real, positive=positive, cache=cache)

        match = symbol_pattern2.match(arg)
        if match is not None:
            # Convert R_{out} to R_out for sympy to recognise.
            arg = match.groups()[0] + match.groups()[1]
            return symbol(arg, real=True)

        # Perhaps have dictionary of functions and their replacements?
        arg = arg.replace('u(t', 'Heaviside(t')
        arg = arg.replace('delta(t', 'DiracDelta(t')

    return sym.sympify(arg, rational=True, locals=symbols)

def get_sympified(instance):
    if hasattr(instance, 'iteritems'):
	return dict([(k,sympify(v)) for k,v in instance.iteritems()])
    elif hasattr(instance, '__iter__'):
	return [sympify(x) for x in instance]
    else:
	NotImplemented

    def __new__(cls, q0, q1, q2, q3):
        q0 = sympify(q0)
        q1 = sympify(q1)
        q2 = sympify(q2)
        q3 = sympify(q3)
        parent_orient = (Matrix([[q0 ** 2 + q1 ** 2 - q2 ** 2 -
                                  q3 ** 2,
                                  2 * (q1 * q2 - q0 * q3),
                                  2 * (q0 * q2 + q1 * q3)],
                                 [2 * (q1 * q2 + q0 * q3),
                                  q0 ** 2 - q1 ** 2 +
                                  q2 ** 2 - q3 ** 2,
                                  2 * (q2 * q3 - q0 * q1)],
                                 [2 * (q1 * q3 - q0 * q2),
                                  2 * (q0 * q1 + q2 * q3),
                                  q0 ** 2 - q1 ** 2 -
                                  q2 ** 2 + q3 ** 2]]))
        parent_orient = parent_orient.T

        obj = super(QuaternionOrienter, cls).__new__(cls, q0, q1, q2, q3)
        obj._q0 = q0
        obj._q1 = q1
        obj._q2 = q2
        obj._q3 = q3
        obj._parent_orient = parent_orient

        return obj

def Ruvws_from_file(filePath):
    """
    Generates nested lists where each list entry contains a again a list.
    The first entry is the Ruvw vector and the second is a list which contains all variables names that have this Ruvw vector.
    Input file must be formatted like:
    0,0,0
    variableName
    variableName
    filePath ... string location to the file which conatins the Ruvw relations
    returns ... list[np.array,list[string]] the information on all Ruvws
    """
    with open(filePath) as file:
        content = file.readlines()
    allRuvw = []
    entries = []
    r = None

    for line in content:
        if ',' in line:
            allRuvw.append([r,entries]) # if it is the first occurence, a dummy entry will be generated
            # Parsing Ruvw
            r = line.split(',')
            r = np.array([sp.sympify(r[0]).evalf(), sp.sympify(r[1]).evalf(), sp.sympify(r[2]).evalf()])
            entries = []
        else:
            entries.append(line.strip())
    # applying the last entry
    allRuvw.append([r,entries])

    del allRuvw[0] # first one is a dummy entry
    return allRuvw

 def __init__(self, name=None, mult=None, num=None, den=None, repeats=None):
     if mult == 1:
         mult = '1'
     if den == 1: 
         den = '1'
     if num == 1:
         num = '1'
     self.den = den
     self.mult = mult
     self.num = num
     self.name = name
     self.denvariables = sorted(set(re.findall(FINDVARIABLES, den)))
     self.numvariables = sorted(set(re.findall(FINDVARIABLES, num)))
     self.denominator_eq = sympy.sympify(den)
     self.numerator_eq = sympy.sympify(num)
     self.denominator = sympy.lambdify(self.denvariables, self.denominator_eq)
     self.numerator = sympy.lambdify(self.numvariables, self.numerator_eq)
     self.repeats = None
     if type(mult) == list:
         self.multvariables = sorted(set([variable for x in self.mult for variable in x.variables]))
         self.multiplier_eq = [str(x) for x in self.mult]
         self.multiplier = mult
     else:
         self.multvariables = sorted(set(re.findall(FINDVARIABLES, mult)))
         self.multiplier_eq = sympy.sympify(mult)
         self.multiplier = sympy.lambdify(self.multvariables, self.multiplier_eq)
     self.variables = self.denvariables + self.numvariables + self.multvariables

def get_gamma_keq_terms(mod, sympy_terms):
    model_map = pysces.ModelMap(mod)  # model map to get substrates, products
    # and parameters for each reaction

    messages = {}
    gamma_keq_terms = {}
    for name, terms in sympy_terms.iteritems():
        reaction_map = getattr(model_map, name)

        substrates = [sympify(substrate) for substrate in
                      reaction_map.hasSubstrates()]

        products = [sympify(product) for product in reaction_map.hasProducts()]

        if len(terms) == 2:  # condition for reversible reactions
            # make sure negative term is second in term list
            terms = sort_terms(terms)
            # divide pos term by neg term and factorise
            expressions = (-terms[0] / terms[1]).factor()
            # get substrate, product and keq terms (and strategy)
            st, pt, keq, _ = get_st_pt_keq(expressions, substrates,
                                                 products)
            if all([st, pt, keq]):
                gamma_keq_terms[name] = pt / (keq*st)
                messages[name] = 'successful generation of gamma/keq term'
            else:
                messages[name] = 'generation of gamma/keq term failed'

    return gamma_keq_terms, messages

def arctan_rule(integral):
    integrand, symbol = integral
    base, exp = integrand.as_base_exp()

    if sympy.simplify(exp + 1) == 0:
        a = sympy.Wild('a', exclude=[symbol])
        b = sympy.Wild('b', exclude=[symbol])
        match = base.match(a + b*symbol**2)
        if match:
            a, b = match[a], match[b]
            if ((isinstance(a, sympy.Number) and a < 0) or (isinstance(b, sympy.Number) and b < 0)):
                return
            if (sympy.ask(sympy.Q.negative(a) | sympy.Q.negative(b) | sympy.Q.is_true(a <= 0) | sympy.Q.is_true(b <= 0))):
                return
            #   /    dx       1  /   dx             1   /     dx                |                     |    1     1         /   du
            #  | --------- = -- | -------------- = --  | -------------------- = | sqrt(b/a)x = u      | =  -- ----------  | -------
            # /  a + bx^2    a /   1  + (b/a)x^2   a  /   1 + (sqrt(b/a)x)^2    | dx = du / sqrt(b/a) |    a   sqrt(b/a) /  1 + u^2
            if a == 1 and b == 1:
                return ArctanRule(integrand, symbol)
            if a == b:
                    constant = 1 / a
                    integrand_ = 1 / (1 + symbol**2)
                    substep = ArctanRule(integrand_, symbol)
                    return ConstantTimesRule(constant, integrand_, substep, integrand, symbol)
            u_var = new_symbol_(symbol)
            u_func = sympy.sqrt(sympy.sympify(b) / a) * symbol
            integrand_ = 1 / (1 + u_func**2)
            constant = 1 / sympy.sqrt(sympy.sympify(b) / a)
            substituted = 1 / (1 + u_var**2)
            substep = ArctanRule(substituted, u_var)
            substep = ConstantTimesRule(constant, substituted, substep, constant*substituted, u_var)
            substep = URule(u_var, u_func, constant, substep, constant*substituted, integrand_, symbol)
            return ConstantTimesRule(1/a, integrand_, substep, integrand, symbol)

def declare_functions():
    
    Expression("create_preMPF", sympify("k9/(1 + k31*OBS_p53)"))
    Expression("Hill_Mdm2", sympify("k24*OBS_Int^n / (k_m^n + OBS_Int^n)"))
    Expression("create_intermediate", sympify("k27*OBS_p53/ (1 + k26*OBS_p53*OBS_Mdm2)"))
    Expression("sig_deg", sympify("Deg_0 - k_deg*(signal-signal_damp)"))
    Expression("kdamp_DDS0", sympify("k_damp*DDS_0"))

 def test_anonymous_expression_operators(self):
     result = Expression('a + b')  # Arbitrary starting expression
     expr_iter = cycle(anonymous_expressions)
     for op in self.ops:
         if op in uniary_ops:
             ss_result = op(result.rhs)
             ee_result = op(result)
             op_str = ("{}({})".format(op.__name__, result))
         else:
             expr = next(expr_iter)
             if op in div_ops and expr.rhs == sympify(0.0):
                 expr = sympify(0.1)
             ss_result = op(result.rhs, expr.rhs)
             ee_result = op(result, expr)
             es_result = op(result, expr.rhs)
             se_result = op(result.rhs, expr)
             op_str = ("{}({}, {})".format(op.__name__, result, expr))
             self.assertEqual(
                 es_result, ss_result,
                 op_str + " not equal between Expression and sympy")
             self.assertEqual(
                 se_result, ss_result,
                 "{} not equal between Expression ({}) and sympy ({})"
                 .format(op_str, se_result, ss_result))
             self.assertIsInstance(es_result, SympyBaseClass,
                                   op_str + " did not return a Expression")
             self.assertIsInstance(se_result, SympyBaseClass,
                                   op_str + " did not return a Expression")
         self.assertEqual(
             ee_result, ss_result,
             "{} not equal between Expression ({}) and sympy ({})"
             .format(op_str, ee_result, ss_result))
         self.assertIsInstance(ee_result, SympyBaseClass,
                               op_str + " did not return a Expression")
         result = Expression(ee_result)

    def _generic_mul(q1, q2):

        q1 = sympify(q1)
        q2 = sympify(q2)

        # None is a Quaternion:
        if not isinstance(q1, Quaternion) and not isinstance(q2, Quaternion):
            return q1 * q2

        # If q1 is a number or a sympy expression instead of a quaternion
        if not isinstance(q1, Quaternion):
            if q2.real_field:
                if q1.is_complex:
                    return q2 * Quaternion(re(q1), im(q1), 0, 0)
                else:
                    return Mul(q1, q2)
            else:
                return Quaternion(q1 * q2.a, q1 * q2.b, q1 * q2.c, q1 * q2.d)


        # If q2 is a number or a sympy expression instead of a quaternion
        if not isinstance(q2, Quaternion):
            if q1.real_field:
                if q2.is_complex:
                    return q1 * Quaternion(re(q2), im(q2), 0, 0)
                else:
                    return Mul(q1, q2)
            else:
                return Quaternion(q2 * q1.a, q2 * q1.b, q2 * q1.c, q2 * q1.d)

        return Quaternion(-q1.b*q2.b - q1.c*q2.c - q1.d*q2.d + q1.a*q2.a,
                          q1.b*q2.a + q1.c*q2.d - q1.d*q2.c + q1.a*q2.b,
                          -q1.b*q2.d + q1.c*q2.a + q1.d*q2.b + q1.a*q2.c,
                          q1.b*q2.c - q1.c*q2.b + q1.d*q2.a + q1.a * q2.d)

def Conditional(cond, true_value, false_value):
    """
    Declares a conditional

    Arguments
    ---------
    cond : A conditional
        The conditional which should be evaluated
    true_value : Any model expression
        Model expression for a true evaluation of the conditional
    false_value : Any model expression
        Model expression for a false evaluation of the conditional
    """
    cond = sp.sympify(cond)

    from sympy.core.relational import Equality, Relational
    from sympy.logic.boolalg import Boolean

    # If the conditional is a bool it is already evaluated
    if isinstance(cond, bool):
        return true_value if cond else false_value

    if not isinstance(cond, (Relational, Boolean)):
        raise type_error("Cond %s is of type %s, but must be a Relational" \
                         " or Boolean." % (cond, type(cond)))

    return sp.functions.Piecewise((true_value, cond), (false_value, sp.sympify(True)),
                                  evaluate=True)

def transverse_magnification(si, so):
    """

    Calculates the transverse magnification, which is the ratio of the
    image size to the object size.

    Parameters
    ==========
    so: sympifiable
    Lens-object distance

    si: sympifiable
    Lens-image distance

    Example
    =======
    >>> from sympy.physics.optics import transverse_magnification
    >>> transverse_magnification(30, 15)
    -2

    """

    si = sympify(si)
    so = sympify(so)

    return (-(si/so))

def daubechis(N):
	# make polynomial
	q_y = [sm.binomial(N-1+k,k) for k in reversed(range(N))]

	# get polynomial roots y[k]
	y = sm.mp.polyroots(q_y, maxsteps=200, extraprec=64)

	z = []
	for yk in y:
		# subustitute y = -1/4z + 1/2 - 1/4/z to factor f(y) = y - y[k]
		f = [sm.mpf('-1/4'), sm.mpf('1/2') - yk, sm.mpf('-1/4')]

		# get polynomial roots z[k]
		z += sm.mp.polyroots(f)

	# make polynomial using the roots within unit circle
	h0z = sm.sqrt('2')
	for zk in z:
		if sm.fabs(zk) < 1:
			h0z *= sympy.sympify('(z-zk)/(1-zk)').subs('zk',zk)

	# adapt vanising moments
	hz = (sympy.sympify('(1+z)/2')**N*h0z).expand()

	# get scaling coefficients
	return [sympy.re(hz.coeff('z',k)) for k in reversed(range(N*2))]