def exact_match(test_expr, target_expr):
    """Test if the entered expression exactly matches the known expression.

       This equality checking does not expand brackets or perform much simplification,
       so is useful for checking if the submitted answer requires simplifying.
       Testing is first done using '==' which checks that the order of symbols
       matches and will not recognise 'x + 1' as equal to '1 + x'.
       The 'srepr' method outputs sympy's internal representation in a canonical form
       and thus, while performing no simplification, it allows ordering to be ignored
       in exact match checking. These two forms are treated equivalently as 'exact'
       matching.

       Returns True if the sympy expressions have the same internal structure,
       and False if not.

        - 'test_expr' should be the untrusted sympy expression to check.
        - 'target_expr' should be the trusted sympy expression to match against.
    """
    print "[EXACT TEST]"
    if test_expr == target_expr:
        print "Exact Match (with '==')"
        return True
    elif sympy.srepr(test_expr) == sympy.srepr(target_expr):
        print "Exact Match (with 'srepr')"
        return True
    else:
        return False

 def __str__(self):
     asstr = "Objective: "
     asstr += srepr(self.objective)
     asstr += "\n\n"
     asstr += "Subject to:\n"
     for c in self._constraints:
         asstr += srepr(c)
         asstr += "\n"
     return asstr

def test_complex_space():
    c1 = ComplexSpace(2)
    assert isinstance(c1, ComplexSpace)
    assert c1.dimension == 2
    assert sstr(c1) == "C(2)"
    assert srepr(c1) == "ComplexSpace(Integer(2))"

    n = Symbol("n")
    c2 = ComplexSpace(n)
    assert isinstance(c2, ComplexSpace)
    assert c2.dimension == n
    assert sstr(c2) == "C(n)"
    assert srepr(c2) == "ComplexSpace(Symbol('n'))"
    assert c2.subs(n, 2) == ComplexSpace(2)

def test_create_f():
    i, j, n, m = symbols('i,j,n,m')
    o = Fd(i)
    assert isinstance(o, CreateFermion)
    o = o.subs(i, j)
    assert o.atoms(Symbol) == {j}
    o = Fd(1)
    assert o.apply_operator(FKet([n])) == FKet([1, n])
    assert o.apply_operator(FKet([n])) == -FKet([n, 1])
    o = Fd(n)
    assert o.apply_operator(FKet([])) == FKet([n])

    vacuum = FKet([], fermi_level=4)
    assert vacuum == FKet([], fermi_level=4)

    i, j, k, l = symbols('i,j,k,l', below_fermi=True)
    a, b, c, d = symbols('a,b,c,d', above_fermi=True)
    p, q, r, s = symbols('p,q,r,s')

    assert Fd(i).apply_operator(FKet([i, j, k], 4)) == FKet([j, k], 4)
    assert Fd(a).apply_operator(FKet([i, b, k], 4)) == FKet([a, i, b, k], 4)

    assert Dagger(B(p)).apply_operator(q) == q*CreateBoson(p)
    assert repr(Fd(p)) == 'CreateFermion(p)'
    assert srepr(Fd(p)) == "CreateFermion(Symbol('p'))"
    assert latex(Fd(p)) == r'a^\dagger_{p}'

def save_sympy_matrix(matrix, filename):
    """Writes a matrix to file in the SymPy representation (srepr)."""
    num_rows, num_cols = matrix.shape
    with open(filename, 'w') as f:
        f.write(str(num_rows) + "\n")
        f.write(str(num_cols) + "\n")
        for expr in matrix:
            f.write(srepr(expr) + "\n")

def test_fixed_bosonic_basis():
    b = FixedBosonicBasis(2, 2)
    # assert b == [FockState((2, 0)), FockState((1, 1)), FockState((0, 2))]
    state = b.state(1)
    assert state == FockStateBosonKet((1, 1))
    assert b.index(state) == 1
    assert b.state(1) == b[1]
    assert len(b) == 3
    assert str(b) == '[FockState((2, 0)), FockState((1, 1)), FockState((0, 2))]'
    assert repr(b) == '[FockState((2, 0)), FockState((1, 1)), FockState((0, 2))]'
    assert srepr(b) == '[FockState((2, 0)), FockState((1, 1)), FockState((0, 2))]'

def test_annihilate_f():
    i, j, n, m = symbols('i,j,n,m')
    o = F(i)
    assert isinstance(o, AnnihilateFermion)
    o = o.subs(i, j)
    assert o.atoms(Symbol) == {j}
    o = F(1)
    assert o.apply_operator(FKet([1, n])) == FKet([n])
    assert o.apply_operator(FKet([n, 1])) == -FKet([n])
    o = F(n)
    assert o.apply_operator(FKet([n])) == FKet([])

    i, j, k, l = symbols('i,j,k,l', below_fermi=True)
    a, b, c, d = symbols('a,b,c,d', above_fermi=True)
    p, q, r, s = symbols('p,q,r,s')
    assert F(i).apply_operator(FKet([i, j, k], 4)) == 0
    assert F(a).apply_operator(FKet([i, b, k], 4)) == 0
    assert F(l).apply_operator(FKet([i, j, k], 3)) == 0
    assert F(l).apply_operator(FKet([i, j, k], 4)) == FKet([l, i, j, k], 4)
    assert str(F(p)) == 'f(p)'
    assert repr(F(p)) == 'AnnihilateFermion(p)'
    assert srepr(F(p)) == "AnnihilateFermion(Symbol('p'))"
    assert latex(F(p)) == 'a_{p}'

def findForwarBackward(formula):
    fragments = findBrackets(formula)
    
    for i in range(len(fragments)):
        #formula = re.sub("(\w)"+re.escape(fragments[0]),r"\1_fragmentNb%d_"%i,formula)
        formula = formula.replace(fragments[i],"_fragmentNb%d_"%i)
    variables = unique(re.findall("_?[A-Za-z]\w*",formula))
    command = ",".join(variables)+" = sympy.symbols(\"" + ",".join(variables) + "\")"
    exec(command)
    formula = sympy.expand(formula)
    
    backward = []
    forward = []
    additiveBlocks = sympy.srepr(formula)
    if additiveBlocks[0:3] == "Add":
        additiveBlocks = sympy.sympify(additiveBlocks[4:len(additiveBlocks)-1])
    else:
        additiveBlocks = [sympy.sympify(additiveBlocks)]
    for ee in additiveBlocks:
        ee = str(sympy.simplify(ee))
        ee = str(ee)
        if ee[0] == "-":
            ee = re.sub("^-","",ee)
            backward.append(ee)        
        else:
            forward.append(ee)
    backward = "+".join(backward)
    forward = "+".join(forward)
    for i in range(len(fragments)):
        backward = backward.replace("_fragmentNb%d_"%i,fragments[i])
        forward = forward.replace("_fragmentNb%d_"%i,fragments[i])

        
        # backward = re.sub("(?<!\w)fragmentNb%d(?!\w)"%i,fragments[i],backward)
        # forward = re.sub("(?<!\w)fragmentNb%d(?!\w)"%i,fragments[i],forward)
    return(forward,backward)

        [ omega_x_terms_expr[theta_dot_expr], omega_x_terms_expr[psi_dot_expr], omega_x_terms_expr[phi_dot_expr] ] ] )

    A_dot_expr = sympy.trigsimp(A_expr.diff(t_expr))

    syms = [ theta_expr, psi_expr, phi_expr, theta_dot_expr, psi_dot_expr, phi_dot_expr ]

    print("Dummifying sympy expressions...")

    A_expr_dummy,     A_expr_dummy_syms     = sympyutils.dummify( A_expr,     syms )
    A_dot_expr_dummy, A_dot_expr_dummy_syms = sympyutils.dummify( A_dot_expr, syms )

    print( "Saving sympy expressions...")

    current_source_file_path = pathutils.get_current_source_file_path()

    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_expr_dummy.dat",     "w") as f: f.write(sympy.srepr(A_expr_dummy))
    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_dot_expr_dummy.dat", "w") as f: f.write(sympy.srepr(A_dot_expr_dummy))

    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_expr_dummy_syms.dat",     "w") as f: f.write(sympy.srepr(sympy.Matrix(A_expr_dummy_syms)))
    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_dot_expr_dummy_syms.dat", "w") as f: f.write(sympy.srepr(sympy.Matrix(A_dot_expr_dummy_syms)))

else:

    print("Loading sympy expressions...")

    current_source_file_path = pathutils.get_current_source_file_path()

    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_expr_dummy.dat",     "r") as f: A_expr_dummy     = sympy.sympify(f.read())
    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_dot_expr_dummy.dat", "r") as f: A_dot_expr_dummy = sympy.sympify(f.read())

    with open(current_source_file_path+"/data/sympy/quadrotorcamera3d_A_expr_dummy_syms.dat",     "r") as f: A_expr_dummy_syms     = array(sympy.sympify(f.read())).squeeze()

        g = self.g
        g_inv = self.g.inv()
        scalar_curv = sympy.simplify(g_inv*Ricci)
        scalar_curv = sympy.trace(scalar_curv)
        self.scalar_curv = scalar_curv
        return self.scalar_curv
           

if __name__ == "__main__":
    import find_metric
    k = -1
    g,diff_form = find_metric.analytical(k) # k=0 gives flat space
    R = Riemann(g,dim=3,sys_title="analytical")
    print R.metric
    from sympy import srepr
    print srepr(R.system)
    RC = R.find_Christoffel_tensor()
    RR = R.find_Ricci_tensor()
    scalarRR = R.find_scalar_curvature()

    print "\nThe analytical curve element has the following metric for k=%.1f"%k
    print g
    print "\nThe Ricci tensor is given as"
    print RR
    print "\nand the scalar curvature is"
    print scalarRR 
    
    """
    from sympy.abc import r,theta, phi, u,v
    g,diff_form = find_metric.flat_sphere()
    diff = [['dv*dv','dv*dw'],['dw*dv','dw*dw']]

 def _sympyrepr(self, printer, *args):
     return "RlpSum(" + str(self.query_symbols) + " in " + str(self.query) + ", " + srepr(self.expression) + ")"

def test_hilbert_space():
    hs = HilbertSpace()
    assert isinstance(hs, HilbertSpace)
    assert sstr(hs) == "H"
    assert srepr(hs) == "HilbertSpace()"

def test_issue1689():
    assert srepr(S(1.0 + 0j)) == srepr(S(1.0)) == srepr(Real(1.0))
    assert srepr(Real(1)) != srepr(Real(1.0))

def print_conversions(*modules):
    for mod in modules:
        for x in module_conversions(mod):
            print '("' + mod.NAME + '", "%s", "%s"): "%s"' % (
                x[0], x[1], srepr(x[2]))

def test_issue1689():
    assert srepr(S(1.0 + 0J)) == srepr(S(1.0)) == srepr(Float(1.0))

def varExprNames(expression):
    variables_Str = list(set(findall(r"Symbol\(\'(\w+)\'\)", srepr(expression))))
    return [var_name for var_name in variables_Str]

 def octoutput(x, et):
     OCTCODE_INT = 1001
     OCTCODE_DOUBLE = 1002
     OCTCODE_STR = 1003
     OCTCODE_BOOL = 1005
     OCTCODE_COMPLEX = 1006
     OCTCODE_DICT = 1010
     OCTCODE_SYM = 1020
     x = objectfilter(x)
     if isinstance(x, bool):
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_BOOL)
         f = ET.SubElement(a, "f")
         f.text = str(x)
     elif x is None or isinstance(x, (sp.Basic, sp.MatrixBase)):
         # FIXME: is it weird to pretend None is a SymPy object?
         if isinstance(x, (sp.Matrix, sp.ImmutableMatrix)):
             _d = x.shape
         elif isinstance(x, sp.MatrixExpr):
             # nan for symbolic size
             _d = [float(r) if (isinstance(r, sp.Basic) and r.is_Integer)
                   else float('nan') if isinstance(r, sp.Basic)
                   else r for r in x.shape]
         elif x is None:
             _d = (1,1)
         else:
             _d = (1, 1)
         try:
             pretty_ascii = sp.pretty(x, use_unicode=False)
         except:
             # e.g., SymPy issue #10414
             pretty_ascii = str(x)
         pretty_unicode = sp.pretty(x, use_unicode=True)
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_SYM)
         f = ET.SubElement(a, "f")
         f.text = sympy.srepr(x)
         f = ET.SubElement(a, "f")
         f.text = str(_d[0])
         f = ET.SubElement(a, "f")
         f.text = str(_d[1])
         f = ET.SubElement(a, "f")
         f.text = str(x)  # esc?
         f = ET.SubElement(a, "f")
         f.text = myesc(pretty_ascii)
         f = ET.SubElement(a, "f")
         f.text = myesc(pretty_unicode)
     elif isinstance(x, (list, tuple)):
         c = ET.SubElement(et, "list")
         for y in x:
             octoutput(y, c)
     elif isinstance(x, int):
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_INT)
         f = ET.SubElement(a, "f")
         f.text = str(x)
     elif isinstance(x, float):
         # We pass IEEE doubles using the exact hex representation
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_DOUBLE)
         f = ET.SubElement(a, "f")
         f.text = d2hex(x)
     elif isinstance(x, complex):
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_COMPLEX)
         f = ET.SubElement(a, "f")
         f.text = d2hex(x.real)
         f = ET.SubElement(a, "f")
         f.text = d2hex(x.imag)
     elif isinstance(x, str) or (sys.version_info < (3, 0) and isinstance(x, unicode)):
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_STR)
         f = ET.SubElement(a, "f")
         f.text = myesc(x)
     elif isinstance(x, dict):
         # Note: the dict cannot be too complex, keys must convert to
         # strings for example.  Values can be dicts, lists.
         a = ET.SubElement(et, "item")
         f = ET.SubElement(a, "f")
         f.text = str(OCTCODE_DICT)
         # Convert keys to strings
         keystr = [str(y) for y in x.keys()]
         c = ET.SubElement(a, "list")
         octoutput(keystr, c)
         c = ET.SubElement(a, "list")
         # FIXME: bit of a kludge, use iterable instead of list, tuple above?
         if sys.version_info >= (3, 0):
             octoutput(list(x.values()), c)
         else:
             octoutput(x.values(), c)
     else:
         raise ValueError("octoutput does not know how to export type " + str(type(x)))

def test_issue_4788():
    assert srepr(S(1.0 + 0J)) == srepr(S(1.0)) == srepr(Float(1.0))

def general_symbolic(target, eqn=None, arg_map=None):
    r'''
    A general function to interpret a sympy equation and evaluate the linear
    components of the source term.

    Parameters
    ----------
    target : OpenPNM object
        The OpenPNM object where the result will be applied.

    eqn : sympy symbolic expression for the source terms
        e.g. y = a*x**b + c

    arg_map : Dict mapping the symbols in the expression to OpenPNM data
        on the target. Must contain 'x' which is the independent variable.
        e.g. arg_map={'a':'pore.a', 'b':'pore.b', 'c':'pore.c', 'x':'pore.x'}

    Example
    ----------
    >>> import openpnm as op
    >>> from openpnm.models.physics import generic_source_term as gst
    >>> import scipy as sp
    >>> import sympy as _syp
    >>> pn = op.network.Cubic(shape=[5, 5, 5], spacing=0.0001)
    >>> water = op.phases.Water(network=pn)
    >>> water['pore.a'] = 1
    >>> water['pore.b'] = 2
    >>> water['pore.c'] = 3
    >>> water['pore.x'] = sp.random.random(water.Np)
    >>> a, b, c, x = _syp.symbols('a,b,c,x')
    >>> y = a*x**b + c
    >>> arg_map = {'a':'pore.a', 'b':'pore.b', 'c':'pore.c', 'x':'pore.x'}
    >>> water.add_model(propname='pore.general',
    ...                 model=gst.general_symbolic,
    ...                 eqn=y, arg_map=arg_map,
    ...                 regen_mode='normal')
    >>> assert 'pore.general.rate' in water.props()
    >>> assert 'pore.general.S1' in water.props()
    >>> assert 'pore.general.S1' in water.props()
    '''
    # First make sure all the symbols have been allocated dict items
    for arg in _syp.postorder_traversal(eqn):
        if _syp.srepr(arg)[:6] == 'Symbol':
            key = _syp.srepr(arg)[7:].strip('(').strip(')').strip("'")
            if key not in arg_map.keys():
                raise Exception('argument mapping incomplete, missing '+key)
    if 'x' not in arg_map.keys():
        raise Exception('argument mapping must contain "x" for the ' +
                        'independent variable')
    # Get the data
    data = {}
    args = {}
    for key in arg_map.keys():
        data[key] = target[arg_map[key]]
        # Callable functions
        args[key] = _syp.symbols(key)
    r, s1, s2 = _build_func(eqn, **args)
    r_val = r(*data.values())
    s1_val = s1(*data.values())
    s2_val = s2(*data.values())
    values = {'S1': s1_val, 'S2': s2_val, 'rate': r_val}
    return values