def test_ccode_Pow():
    assert ccode(x**3) == "pow(x, 3)"
    assert ccode(x**(y**3)) == "pow(x, pow(y, 3))"
    assert ccode(1/(g(x)*3.5)**(x - y**x)/(x**2 + y)) == \
        "pow(3.5*g(x), -x + pow(y, x))/(pow(x, 2) + y)"
    assert ccode(x**-1.0) == '1.0/x'
    assert ccode(x**Rational(2, 3)) == 'pow(x, 2.0L/3.0L)'

def simplex_table_to_assumptions(table, symbols_used):
    common_assumps, new_assumps_by_pot = table.pots.get_common_and_not_assumptions()

    ret = ''

    common_assumpt_text = list()
    for a in common_assumps:
        process_expr_on_symbols(a.exp, symbols_used)
        common_assumpt_text.append('(' + ccode( a.exp) + ' ' + a.sign + ' 0 )')

    pots_assumpt_text = list()

    for pot in new_assumps_by_pot:

        pot_assumpt_text = list()
        for a in pot:
            process_expr_on_symbols(a.exp, symbols_used)
            pot_assumpt_text.append('(' + ccode( a.exp) + ' ' + a.sign + ' 0 )')

        pots_assumpt_text.append( '(' + " && ".join(pot_assumpt_text) + ')' )

    ret = " && ".join(common_assumpt_text)

    if len(pots_assumpt_text) > 0:
        if len(common_assumps): ret += ' && '
        ret += '(' + " || ".join(pots_assumpt_text) + ')'

    process_table_on_symbols(table, symbols_used)
    return ret

def test_ccode_reserved_words():
    x, y = symbols('x, if')
    with raises(ValueError):
        ccode(y**2, error_on_reserved=True, standard='C99')
    assert ccode(y**2) == 'pow(if_, 2)'
    assert ccode(x * y**2, dereference=[y]) == 'pow((*if_), 2)*x'
    assert ccode(y**2, reserved_word_suffix='_unreserved') == 'pow(if_unreserved, 2)'

def test_ccode_Piecewise_deep():
    p = ccode(2*Piecewise((x, x < 1), (x + 1, x < 2), (x**2, True)))
    assert p == (
            "2*((x < 1) ? (\n"
            "   x\n"
            ")\n"
            ": ((x < 2) ? (\n"
            "   x + 1\n"
            ")\n"
            ": (\n"
            "   pow(x, 2)\n"
            ")))")
    expr = x*y*z + x**2 + y**2 + Piecewise((0, x < 0.5), (1, True)) + cos(z) - 1
    assert ccode(expr) == (
            "pow(x, 2) + x*y*z + pow(y, 2) + ((x < 0.5) ? (\n"
            "   0\n"
            ")\n"
            ": (\n"
            "   1\n"
            ")) + cos(z) - 1")
    assert ccode(expr, assign_to='c') == (
            "c = pow(x, 2) + x*y*z + pow(y, 2) + ((x < 0.5) ? (\n"
            "   0\n"
            ")\n"
            ": (\n"
            "   1\n"
            ")) + cos(z) - 1;")

def test_ccode_Indexed():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    s, n, m, o = symbols('s n m o', integer=True)
    i, j, k = Idx('i', n), Idx('j', m), Idx('k', o)

    x = IndexedBase('x')[j]
    A = IndexedBase('A')[i, j]
    B = IndexedBase('B')[i, j, k]

    with warnings.catch_warnings():
        warnings.filterwarnings("ignore", category=SymPyDeprecationWarning)
        p = CCodePrinter()
        p._not_c = set()

        assert p._print_Indexed(x) == 'x[j]'
        assert p._print_Indexed(A) == 'A[%s]' % (m*i+j)
        assert p._print_Indexed(B) == 'B[%s]' % (i*o*m+j*o+k)
        assert p._not_c == set()

        A = IndexedBase('A', shape=(5,3))[i, j]
        assert p._print_Indexed(A) == 'A[%s]' % (3*i + j)

        A = IndexedBase('A', shape=(5,3), strides='F')[i, j]
        assert ccode(A) == 'A[%s]' % (i + 5*j)

        A = IndexedBase('A', shape=(29,29), strides=(1, s), offset=o)[i, j]
        assert ccode(A) == 'A[o + s*j + i]'

        Abase = IndexedBase('A', strides=(s, m, n), offset=o)
        assert ccode(Abase[i, j, k]) == 'A[m*j + n*k + o + s*i]'
        assert ccode(Abase[2, 3, k]) == 'A[3*m + n*k + o + 2*s]'

def main():
    # GENERATES WITHOUT ARGUMENTS: rev_unary_1.tex rev_unary_rate.tex rev_unary_k_b.c rev_unary_K_eq.c
    from sympy.printing.latex import latex
    open('rev_unary_1.tex', 'wt').write(as_align_env(eqs))
    open('rev_unary_rate.tex', 'wt').write(latex(rate_expr))
    open('rev_unary_k_b.c', 'wt').write('return_val = {};'.format(ccode(alt_expl_in_t)))
    open('rev_unary_K_eq.c', 'wt').write('return_val = {};'.format(
        ccode(alt_expl_in_t.subs({k_b: k_f/K_eq}))))

def test_ccode_sign():
    expr1, ref1 = sign(x) * y, 'y*(((x) > 0) - ((x) < 0))'
    expr2, ref2 = sign(cos(x)), '(((cos(x)) > 0) - ((cos(x)) < 0))'
    expr3, ref3 = sign(2 * x + x**2) * x + x**2, 'pow(x, 2) + x*(((pow(x, 2) + 2*x) > 0) - ((pow(x, 2) + 2*x) < 0))'
    assert ccode(expr1) == ref1
    assert ccode(expr1, 'z') == 'z = %s;' % ref1
    assert ccode(expr2) == ref2
    assert ccode(expr3) == ref3

    def __init__(self, model, tspan, integrator='vode', **integrator_options):

        pysb.bng.generate_equations(model)

        code_eqs = '\n'.join(['ydot[%d] = %s;' % (i, sympy.ccode(model.odes[i])) for i in range(len(model.odes))])
        code_eqs = re.sub(r's(\d+)', lambda m: 'y[%s]' % (int(m.group(1))), code_eqs)
        for e in model.expressions:
            code_eqs = re.sub(r'\b(%s)\b' % e.name,
                              sympy.ccode(e.expand_expr()), code_eqs)
        for i, p in enumerate(model.parameters):
            code_eqs = re.sub(r'\b(%s)\b' % p.name, 'p[%d]' % i, code_eqs)

        Solver._test_inline()
        # If we can't use weave.inline to run the C code, compile it as Python code instead for use with
        # exec. Note: C code with array indexing, basic math operations, and pow() just happens to also
        # be valid Python.  If the equations ever have more complex things in them, this might fail.
        if not Solver._use_inline:
            code_eqs_py = compile(code_eqs, '<%s odes>' % model.name, 'exec')
        else:
            for arr_name in ('ydot', 'y', 'p'):
                macro = arr_name.upper() + '1'
                code_eqs = re.sub(r'\b%s\[(\d+)\]' % arr_name,
                                  '%s(\\1)' % macro, code_eqs)

        def rhs(t, y, p):
            ydot = self.ydot
            # note that the evaluated code sets ydot as a side effect
            if Solver._use_inline:
                inline(code_eqs, ['ydot', 't', 'y', 'p']);
            else:
                exec code_eqs_py in locals()
            return ydot

        # build integrator options list from our defaults and any kwargs passed to this function
        options = {}
        try:
            options.update(default_integrator_options[integrator])
        except KeyError as e:
            pass
        options.update(integrator_options)

        self.model = model
        self.tspan = tspan
        self.y = numpy.ndarray((len(tspan), len(model.species)))
        self.ydot = numpy.ndarray(len(model.species))
        if len(model.observables):
            self.yobs = numpy.ndarray(len(tspan), zip(model.observables.keys(),
                                                      itertools.repeat(float)))
        else:
            self.yobs = numpy.ndarray((len(tspan), 0))
        exprs = model.expressions_dynamic()
        if len(exprs):
            self.yexpr = numpy.ndarray(len(tspan), zip(exprs.keys(),
                                                       itertools.repeat(float)))
        else:
            self.yexpr = numpy.ndarray((len(tspan), 0))
        self.yobs_view = self.yobs.view(float).reshape(len(self.yobs), -1)
        self.integrator = ode(rhs).set_integrator(integrator, **options)

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
    from sympy.printing.latex import latex

    open("irrev_binary_1.tex", "wt").write(as_align_env(eqs[:2]))
    open("irrev_binary_2.tex", "wt").write(as_align_env(eqs[2:]))
    open("irrev_binary_rate.tex", "wt").write(latex(rate_expr))
    open("irrev_binary_k_b.c", "wt").write("return_val = {};".format(ccode(alt_expl_in_t)))
    open("irrev_binary_K_eq.c", "wt").write("return_val = {};".format(ccode(alt_expl_in_t.subs({k_b: k_f / K_eq}))))

def test_ccode_user_functions():
    x = symbols('x', integer=False)
    n = symbols('n', integer=True)
    custom_functions = {
        "ceiling": "ceil",
        "Abs": [(lambda x: not x.is_integer, "fabs"), (lambda x: x.is_integer, "abs")],
    }
    assert ccode(ceiling(x), user_functions=custom_functions) == "ceil(x)"
    assert ccode(Abs(x), user_functions=custom_functions) == "fabs(x)"
    assert ccode(Abs(n), user_functions=custom_functions) == "abs(n)"

def test_ccode_sign():

    expr = sign(x) * y
    assert ccode(expr) == 'y*(((x) > 0) - ((x) < 0))'
    assert ccode(expr, 'z') == 'z = y*(((x) > 0) - ((x) < 0));'

    assert ccode(sign(2 * x + x**2) * x + x**2) == \
        'pow(x, 2) + x*(((pow(x, 2) + 2*x) > 0) - ((pow(x, 2) + 2*x) < 0))'

    expr = sign(cos(x))
    assert ccode(expr) == '(((cos(x)) > 0) - ((cos(x)) < 0))'

def test_MatrixElement_printing():
    # test cases for issue #11821
    A = MatrixSymbol("A", 1, 3)
    B = MatrixSymbol("B", 1, 3)
    C = MatrixSymbol("C", 1, 3)

    assert(ccode(A[0, 0]) == "A[0]")
    assert(ccode(3 * A[0, 0]) == "3*A[0]")

    F = C[0, 0].subs(C, A - B)
    assert(ccode(F) == "(-B + A)[0]")

def test_ccode_inline_function():
    x = symbols('x')
    g = implemented_function('g', Lambda(x, 2*x))
    assert ccode(g(x)) == "2*x"
    g = implemented_function('g', Lambda(x, 2*x/Catalan))
    assert ccode(g(x)) == "double const Catalan = %s;\n2*x/Catalan" %Catalan.n()
    A = IndexedBase('A')
    i = Idx('i', symbols('n', integer=True))
    g = implemented_function('g', Lambda(x, x*(1 + x)*(2 + x)))
    assert ccode(g(A[i]), assign_to=A[i]) == (
            "for (int i=0; i<n; i++){\n"
            "   A[i] = (1 + A[i])*(2 + A[i])*A[i];\n"
            "}"
            )

def test_ccode_loops_multiple_contractions():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m, o, p = symbols('n m o p', integer=True)
    a = IndexedBase('a')
    b = IndexedBase('b')
    y = IndexedBase('y')
    i = Idx('i', m)
    j = Idx('j', n)
    k = Idx('k', o)
    l = Idx('l', p)

    s = (
        'for (int i=0; i<m; i++){\n'
        '   y[i] = 0;\n'
        '}\n'
        'for (int i=0; i<m; i++){\n'
        '   for (int j=0; j<n; j++){\n'
        '      for (int k=0; k<o; k++){\n'
        '         for (int l=0; l<p; l++){\n'
        '            y[i] = a[%s]*b[%s] + y[i];\n' % (i*n*o*p + j*o*p + k*p + l, j*o*p + k*p + l) +\
        '         }\n'
        '      }\n'
        '   }\n'
        '}'
    )
    assert ccode(b[j, k, l]*a[i, j, k, l], assign_to=y[i]) == s

def _ccode(expr, ):
    code = sympy.ccode(expr)
    if options['unroll_square']:
        code =  re.sub(r'pow\(([^,]*), 2\)', r'((\1)*(\1))', code)
    if options['custom_sign_func']:
        code =  re.sub(r'\(\(([^,]*)\) > 0\) - \(\(([^,]*)\) < 0\)', r'Sign(\1)', code)
    return code

def run(model):
  output = StringIO()
  pysb.bng.generate_equations(model)

  output.write("% MATLAB model definition file\n")
  output.write('%% save as %s_odes.m\n' % model_name)
  output.write('function out = %s_odes(t, input, param)' % model_name)
  output.write("\n\n")

  c_code_consts = '\n'.join(['param(%d) = %s %% %s;' % (i+1, p.value, p.name) for i, p in
      enumerate(model.parameters)])
  c_code_eqs = '\n'.join(['out(%d,1) = %s;' % (i+1, sympy.ccode(model.odes[i])) for i in
      range(len(model.odes))])
  c_code_eqs = re.sub(r's(\d+)', lambda m: 'input(%s)' % (int(m.group(1))+1), c_code_eqs)
  c_code_eqs = re.sub(r'pow\(', 'power(', c_code_eqs)
  c_code = c_code_eqs

  c_code_species = '\n'.join(['%% input(%d) = %s;' % (i+1, s) for i, s in
      enumerate(model.species)])

  for i, p in enumerate(model.parameters):
    c_code = re.sub(r'\b(%s)\b' % p.name, 'param(%d)' % (i+1), c_code)

  output.write(c_code_consts + "\n\n")
  output.write(c_code_species + "\n\n")
  output.write(c_code + "\n\n")
  output.write("end\n")
  return output.getvalue()

def run(model):
    output = StringIO()
    pysb.bng.generate_equations(model)

    ic_values = [0] * len(model.odes)
    for cp, ic_param in model.initial_conditions:
        ic_values[model.get_species_index(cp)] = ic_param.value

    # list of "dynamic variables"
    pw_x = ["m = pwAddX(m, 's%d', %e);" % (i, ic_values[i]) for i in range(len(model.odes))]

    # parameters
    pw_k = ["m = pwAddK(m, '%s', %e);" % (p.name, p.value) for p in model.parameters]

    # equations (one for each dynamic variable)
    # Note that we just generate C code, which for basic math expressions
    # is identical to matlab.  We just have to change 'pow' to 'power'.
    # Ideally there would be a matlab formatter for sympy.
    pw_ode = ["m = pwAddODE(m, 's%d', '%s');" % (i, sympy.ccode(model.odes[i])) for i in range(len(model.odes))]
    pw_ode = [re.sub(r'pow(?=\()', 'power', s) for s in pw_ode]

    # observables
    pw_y = ["m = pwAddY(m, '%s', '%s');" % (obs.name, ' + '.join('%f * s%s' % t for t in zip(obs.coefficients, obs.species))) for obs in model.observables]

    output.write('% PottersWheel model definition file\n')
    output.write('%% save as %s.m\n' % model_name)
    output.write('function m = %s()\n' % model_name)
    output.write('\n')
    output.write('m = pwGetEmptyModel();\n')
    output.write('\n')
    output.write('% meta information\n')
    output.write("m.ID          = '%s';\n" % model_name)
    output.write("m.name        = '%s';\n" % model_name)
    output.write("m.description = '';\n")
    output.write("m.authors     = {''};\n")
    output.write("m.dates       = {''};\n")
    output.write("m.type        = 'PW-1-5';\n")
    output.write('\n')
    output.write('% dynamic variables\n')
    for x in pw_x:
        output.write(x)
        output.write('\n')
    output.write('\n')
    output.write('% dynamic parameters\n')
    for k in pw_k:
        output.write(k)
        output.write('\n')
    output.write('\n')
    output.write('% ODEs\n')
    for ode in pw_ode:
        output.write(ode)
        output.write('\n')
    output.write('\n')
    output.write('% observables\n')
    for y in pw_y:
        output.write(y)
        output.write('\n')
    output.write('\n')
    output.write('%% end of PottersWheel model %s\n' % model_name)
    return output.getvalue()

        def eqn_substitutions(eqns):
            """String substitutions on the sympy C code for the ODE RHS and
            Jacobian functions to use appropriate terms for variables and
            parameters."""
            # Substitute expanded parameter formulas for any named expressions
            for e in self.model.expressions:
                eqns = re.sub(r'\b(%s)\b' % e.name, '('+sympy.ccode(
                    e.expand_expr())+')', eqns)

            # Substitute sums of observable species that could've been added
            # by expressions
            for obs in self.model.observables:
                obs_string = ''
                for i in range(len(obs.coefficients)):
                    if i > 0:
                        obs_string += "+"
                    if obs.coefficients[i] > 1:
                        obs_string += str(obs.coefficients[i])+"*"
                    obs_string += "__s"+str(obs.species[i])
                if len(obs.coefficients) > 1:
                    obs_string = '(' + obs_string + ')'
                eqns = re.sub(r'\b(%s)\b' % obs.name, obs_string, eqns)

            # Substitute 'y[i]' for 'si'
            eqns = re.sub(r'_*s(\d+)', lambda m: 'y[%s]' % (int(m.group(1))),
                       eqns)

            # Substitute 'p[i]' for any named parameters
            for i, p in enumerate(self.model.parameters):
                eqns = re.sub(r'\b(%s)\b' % p.name, 'p[%d]' % i, eqns)
            return eqns

 def eval_str(self, dom, idx_str="i", qd_pt_str="QUAD_PTS"):
     qd_id_str = "(%(jac)f*%(qd_pt_str)s[%(idx_str)s]+%(shift)f)" % \
         {"qd_pt_str":qd_pt_str,
          "idx_str":idx_str,
          "jac":(dom.args[2] - dom.args[1]) / 2.0 ,
          "shift":(dom.args[2] + dom.args[1]) / 2.0 }
     return ccode(self.args[0]).replace(str(self.args[1]), qd_id_str)

def test_ccode_loops_addfactor():
    from sympy.tensor import IndexedBase, Idx
    from sympy import symbols
    n, m, o, p = symbols('n m o p', integer=True)
    a = IndexedBase('a')
    b = IndexedBase('b')
    c = IndexedBase('c')
    y = IndexedBase('y')
    i = Idx('i', m)
    j = Idx('j', n)
    k = Idx('k', o)
    l = Idx('l', p)

    s = (
        'for (int i=0; i<m; i++){\n'
        '   y[i] = 0;\n'
        '}\n'
        'for (int i=0; i<m; i++){\n'
        '   for (int j=0; j<n; j++){\n'
        '      for (int k=0; k<o; k++){\n'
        '         for (int l=0; l<p; l++){\n'
        '            y[i] = (a[%s] + b[%s])*c[%s] + y[i];\n' % (i*n*o*p + j*o*p + k*p + l, i*n*o*p + j*o*p + k*p + l, j*o*p + k*p + l) +\
        '         }\n'
        '      }\n'
        '   }\n'
        '}'
    )
    c = ccode((a[i, j, k, l] + b[i, j, k, l])*c[j, k, l], assign_to=y[i])
    assert c == s