def genEval(self):
     text = "        bool evaluate(\n"
     args=[]
     for name,v in self.variables:
         if v.is_Matrix:
             args.append("           const Eigen::MatrixXd & %s" % name)
         else:
             args.append("           double %s" % name)
     args.append("           bool evalF=true,bool evalJ=true")
     text += ",\n".join(args) + ") {\n"
     text += "           if (evalF) {\n"
     (interm, expr) = cse(self.function,numbered_symbols("__x"));
     for dummy,exp in interm:
         text += "               double %s = %s;\n" % (str(dummy),ccode(exp))
     for i in range(self.function.rows):
         text += "               F(%d) = %s;\n" % (i,ccode(expr[0][i]))
     text += "           }\n"
     text += "           if (evalJ) {\n"
     (interm, expr) = cse(self.J,numbered_symbols("__x"));
     for dummy,exp in interm:
         text += "               double %s = %s;\n" % (str(dummy),ccode(exp))
     for i in range(self.J.rows):
         for j in range(self.J.cols):
             text += "               J(%d,%d) = %s;\n" % (i,j,ccode(expr[0][i,j]))
     text += "           }\n"
     text += "           return true;\n"
     text += "       }\n"
     return text

    def __init__(self, symbsol, V, kappa=1):

        x, y, t = smp.symbols('x[0], x[1], t')
        rhs = symbsol.diff(t) \
              + (kappa * (symbsol.diff(x))).diff(x) \
              + (kappa * (symbsol.diff(y))).diff(y) 

        from sympy.printing import ccode
        self.sol = Expression(ccode(symbsol), t=0.0)
        self.rhs = Expression(ccode(rhs), t=0.0)
        self.V = V

        v = TestFunction(V)
        u = TrialFunction(V)
        
        # Assemble system
        M = assemble(inner(u, v)*dx)
        A = assemble(kappa * inner(grad(u), grad(v)) * dx)
        # Convert DOLFIN representation to numpy arrays
        rows, cols, values = M.data()
        self.M = sps.csr_matrix((values, cols, rows))
        """csr matrix for the mass"""
        rows, cols, values = A.data()
        self.A = sps.csr_matrix((values, cols, rows))
        """csr matrix representing the weak discrete 
        :math:`-\\nabla \\cdot (\\kappa \\nabla )` operator"""
            

        # treatment of the boundary values
        nv = self.A.shape[0]
        auxu = np.zeros((nv,1))
        self.bcinds = []
        self.bc = DirichletBC(self.V, self.sol, 'on_boundary')
        self.bcdict = self.bc.get_boundary_values()
        auxu[self.bcdict.keys(),0] = self.bcdict.values()
        self.bcinds.extend(self.bcdict.keys())

        self.rhsbc = - self.A*auxu    
    
        # indices of the innernodes
        self.invinds = np.setdiff1d(range(nv),self.bcinds).astype(np.int32)
        # condense the coeff mats to the inner nodes 
        self.M = self.M[self.invinds,:][:,self.invinds]
        self.A = self.A[self.invinds,:][:,self.invinds]

        self.rhsbc = self.rhsbc[self.invinds,:]

        self.bcvals = auxu[self.bcinds]

 def rhs_xml(self):
     rhs = self.expand_integer_powers(self.rhs)
     s = ccode(rhs, user_functions=self._random_map)
     s = self.strip_L_from_rationals(s)
     s = self._ccode_print_warn_re.sub('', s)
     s = self._multiple_whitespace_re.sub(' ', s)
     return s

 def assign_str(self, lhs, rhs):
     rhs = Expression.expand_integer_powers(rhs)
     nmodl_str = ccode(rhs, user_functions=Expression._cfunc_map,
                       assign_to=lhs)
     nmodl_str = Expression.strip_L_from_rationals(nmodl_str)
     nmodl_str = nmodl_str.replace(';', '')
     return nmodl_str

    def extraToFile(self, path):

        raw = """double %sOrbitals::get_dell_alpha_phi(const Walker* walker, int qnum, int i){
    
    double dphi;
    
    __code__
    
    return dphi;
    
}""" % self.name

        shell = """if (qnum == _q_) {

    __necessities__
    
        //__simple__
        
        dphi = __expr__
        
    } else """
    
        

        Z = Symbol('Z', positive=True, real=True)

        code = "    "

        for i in range(self.maxImplemented/2):

            psi = self.orbitals[i]
            qNums = self.stateMap[i] 
            genFac = self.genericFactor(qNums, basic=False)
           
            kdiff = diff(psi, k).factor(genFac)/psi*Z
            
            simple = self.makeReadable(str(kdiff))

            nec, necList = self.getNecessities(kdiff)
            expr = printing.ccode(kdiff) + ";"
            expr = self.replaceCCode(expr, necList)

            #hack to get the right indent
            nec = "\n".join([" "*4  + nec_i for nec_i in nec.split("\n")])

            subCode = shell
            subCode = subCode.replace("\n\n    __necessities__", nec)\
                             .replace("__expr__", expr)\
                             .replace("__simple__", simple)\
                             .replace("_q_", str(i))
            
            code += subCode
   
        code = code.strip("else ")
       
        ccode = raw.replace("__code__", code)
      
        with open(pjoin(path, "%sOrbitalsAlphaDeriv.cpp" % self.name), 'w') as f:
            f.write(ccode)
            f.close()

    def getNecessities(self, expr):
        
        s = printing.ccode(expr)        
        
        nec = []
        necS = ""
        necS2 = ""

        for i, x in enumerate(self.xi):
            l = len(regxp.findall("pow\(%s\, \d+\)" % x, s))
       
            if regxp.findall("[^\w]?%s[^\w]" % x, s):
                nec.append(x)
                necS += "    %s = walker->r(i, %d);\n" % (x, i)

            if l > 0:
                x2 = x + "2"
                nec.append(x2)
                necS2 += "    %s = %s*%s;\n" % (x2, x, x)

                
        necS = ("%s\n%s" % (necS, necS2)).strip("\n")
        
        #Manually add the lineshifts so that if no nec, then no lineshift        
        if necS:
            necS = "\n\n" + necS
        
        return necS, nec

 def renderModelMethodImplementation(self):
     if self.expression is not None:
         implementation = ccode(self.expression)
     else:
         implementation = "ASSERT(False)"
     return render('model_methodImplementation.cc', dict(evalClassName=self.d['evalClassName'],
                                                         myMethod=self.d['myKeyMethod'],
                                                         myMethodDeclarationArgs=self.d['myMethodDeclarationArgs'],
                                                         myMethodImplementation=implementation))

def test_ccode_boolean():
    assert ccode(x&y) == "x&&y"
    assert ccode(x|y) == "x||y"
    assert ccode(~x) == "!x"
    assert ccode(x&y&z) == "x&&y&&z"
    assert ccode(x|y|z) == "x||y||z"
    assert ccode((x&y)|z) == "x&&y||z"
    assert ccode((x|y)&z) == "(x||y)&&z"

def GenCode(exp):
    # From http://docs.sympy.org/dev/modules/utilities/codegen.html
    [(c_name, c_code), (h_name, c_header)] = codegen(
        ("exp", exp), "C", "test", header=False, empty=False)

    print("c_code = %s" % c_code)
    print("ccode = %s" % ccode(exp))

    #!!!! TODO TODO TODO: we should check if we have symbolic denominator and generate code that takes care to check if the denominator is 0 and if so to abort the compuatation
    return None

def _str_scal_func_update(method, name, objlabel):
    expr = getattr(method, name + '_expr')
    update_h = '\nvoid {0}_update();'.format(name)
    update_cpp = '\nvoid {0}::{1}_update()'.format(objlabel, name) + '{'
    symbs = expr.free_symbols
    if any(symb in method.args for symb in symbs):
        c = ccode(expr, dereference=dereference(method))
        update_cpp += '\n_{0} = {1};'.format(name, c)
    update_cpp += '\n};'
    return update_h, update_cpp

def _str_scal_func_init_data(method, name):
    expr = getattr(method, name + '_expr')
    init_data = '{1} {0}_data = '.format(name, CONFIG_CPP['float'])
    symbs = expr.free_symbols
    if any(symb in method.args for symb in symbs):
        init_data += "0.;"
    else:
        c = ccode(expr, dereference=dereference(method))
        init_data += '{};'.format(c)
    return init_data

def _str_mat_func_update(method, name, objlabel):
    mat = types.matrix_types[0](getattr(method, name + '_expr'))
    update_h = '\nvoid {0}_update();'.format(name)
    update_cpp = '\nvoid {0}::{1}_update()'.format(objlabel, name) + '{'
    for m, n, expr in mat.row_list():
            symbs = expr.free_symbols
            if any(symb in method.args() for symb in symbs):
                c = ccode(expr, dereference=dereference(method))
                update_cpp += '\n_{0}({1}, {2}) = {3};'.format(name, m, n, c)
    update_cpp += '\n};'
    return update_h, update_cpp

def test_ccode_Piecewise():
    p = ccode(Piecewise((x,x<1),(x**2,True)))
    s = \
"""\
if (x < 1) {
x
}
else {
pow(x,2)
}\
"""
    assert p == s

    def _call_printer(self, routine):
        code_lines = []

        # Compose a list of symbols to be dereferenced in the function
        # body. These are the arguments that were passed by a reference
        # pointer, excluding arrays.
        dereference = []
        for arg in routine.arguments:
            if isinstance(arg, ResultBase) and not arg.dimensions:
                dereference.append(arg.name)

        return_val = None
        for result in routine.result_variables:
            if isinstance(result, Result):
                assign_to = routine.name + "_result"
                t = result.get_datatype('c')
                code_lines.append("{0} {1};\n".format(t, str(assign_to)))
                return_val = assign_to
            else:
                assign_to = result.result_var

            try:
                # order='none' is an optimization not in upstream
                constants, not_c, c_expr = ccode(result.expr, human=False,
                                                 assign_to=assign_to, dereference=dereference, order='none')
            except AssignmentError:
                assign_to = result.result_var
                code_lines.append(
                    "%s %s;\n" % (result.get_datatype('c'), str(assign_to)))
                constants, not_c, c_expr = ccode(result.expr, human=False,
                                                 assign_to=assign_to, dereference=dereference, order='none')

            for name, value in sorted(constants, key=str):
                code_lines.append("double const %s = %s;\n" % (name, value))
            code_lines.append("%s\n" % c_expr)

        if return_val:
            code_lines.append("   return %s;\n" % return_val)
        return code_lines

 def genEval(self):
     text = "        bool evaluate(\n"
     args=[]
     for name,v,lD in self.variables:
         if isinstance(v, Matrix):
             args.append("           const Eigen::Matrix<double, %d, 1> & %s" % (v.rows, name))
         else:
             args.append("           double %s" % name)
             
     args.append("           Eigen::Matrix<double, %d, 1> * F" % self.function.rows)
     
     for name,v,localDim in self.variables:
         if isinstance(v, Matrix):
             args.append("           Eigen::Matrix<double, %d, %d> * J%s" % (self.J.rows, localDim, name))
         else:
             args.append("           Eigen::Matrix<double, %d, 1> * J%s" % (self.J.rows, name))
     text += ",\n".join(args) + ") {\n"
     text += "           if (F) {\n"
     (interm, expr) = cse(self.function,numbered_symbols("__x"));
     for dummy,exp in interm:
         text += "               double %s = %s;\n" % (str(dummy),ccode(exp))
     for i in range(self.function.rows):
         text += "               (*F)(%d) = %s;\n" % (i,ccode(expr[0][i]))
     text += "           }\n"
     text += "           if (%s) {\n" % " && ".join([ "J" + name for name,v,lD in self.variables ])
     (interm, expr) = cse(self.J,numbered_symbols("__x"));
     for dummy,exp in interm:
         text += "               double %s = %s;\n" % (str(dummy),ccode(exp))
     colBase = 0;
     for name,v,localDim in self.variables:
         for i in range(self.J.rows):
             for j in range(0, localDim):
                 text += "               (*J%s)(%d,%d) = %s;\n" % (name, i,j,ccode(expr[0][i,colBase + j]))
         colBase+=localDim
     text += "           }\n"
     text += "           return true;\n"
     text += "       }\n"
     return text

    def __init__(self, N, omega=None, nu=None, scheme='TH'):

        self.N = N
        if scheme == 'TH':
            self.mesh = smamin_thcr_mesh.getmake_mesh(N)
            self.V = dolfin.VectorFunctionSpace(self.mesh, "CG", 2)
            self.Q = dolfin.FunctionSpace(self.mesh, "CG", 1)
        elif scheme == 'CR':
            self.mesh = dolfin.UnitSquareMesh(N, N)  # , 'crossed')
            self.V = dolfin.VectorFunctionSpace(self.mesh, "CR", 1)
            self.Q = dolfin.FunctionSpace(self.mesh, "DG", 0)
        self.velbcs = setget_velbcs_zerosq(self.mesh, self.V)
        self.Pdof = 0  # dof removed in the p approximation
        self.omega = omega
        self.nu = nu

        x, y, t, nu, om = smp.symbols('x[0], x[1], t, nu, omega')

        ft = smp.sin(om*t)

        u1 = ft*x*x*(1 - x)*(1 - x)*2*y*(1 - y)*(2*y - 1)
        u2 = ft*y*y*(1 - y)*(1 - y)*2*x*(1 - x)*(1 - 2*x)
        p = ft*x*(1 - x)*y*(1 - y)

        du1 = smp.diff(u1, t)
        du2 = smp.diff(u2, t)

        rhs1, rhs2, rhs3 = comp_symb_nserhs(u=u1, v=u2, p=p, nu=self.nu)

        from sympy.printing import ccode
        self.v = Expression((ccode(u1), ccode(u2)),
                            t=0.0, omega=self.omega)
        self.p = Expression((ccode(p)),
                            t=0.0, omega=self.omega)
        self.fv = Expression((ccode(rhs1), ccode(rhs2)),
                             t=0.0, omega=self.omega, nu=self.nu)
        self.fp = Expression((ccode(rhs3)),
                             t=0.0, omega=self.omega)
        self.vdot = Expression((ccode(du1), ccode(du2)),
                               t=0.0, omega=self.omega)

        bcinds = []
        for bc in self.velbcs:
            bcdict = bc.get_boundary_values()
            bcinds.extend(bcdict.keys())

        # indices of the inner velocity nodes
        self.invinds = np.setdiff1d(range(self.V.dim()), bcinds)

    def renderModelDerivImplementations(self):
        impls = []

        for arg,var in zip(self.args,self.vars):
            if self.expression is not None:
                print "differentiation of", self.expression, "with respect to", var
                implementation = ccode(self.expression.diff(var))
            else:
                implementation = "ASSERT(False)"
            impls.append(render('model_methodImplementation.cc',
                                dict(evalClassName=self.d['evalClassName'],
                                     myMethod="D%sD%s"%(self.d['myKeyMethod'],''.join([word[0].upper()+word[1:] for word in arg.split("_")])),
                                     myMethodDeclarationArgs=self.d['myMethodDeclarationArgs'],
                                     myMethodImplementation=implementation)))
        return '\n\n'.join(impls)

    def make_C_term(self, term, no_correct_rate, derivate=None, inverse=None, human=False, force_par=False, xify=None, set_t0=False):

        """transform a term into its ssm C expression OR the ssm C
        expression of its derivate, differentiating against the
        derivate (if derivate not None) OR compute inverse function
        """

        #prefix all the state variable and parameters by ssm___ to
        #avoid namespace collision with Sympy as QCOSINE letters are
        #used by SymPy

        myterm = self.change_user_input(term)
        safe = ''

        for r in myterm:
            if r in self.all_par:
                safe += 'ssm___' + r
            elif inverse and r == inverse:
                safe += 'ssm___' + r
            else:
                safe += r

        if derivate:
	    sy = Symbol(str('ssm___' + derivate)) if derivate != 'x' else Symbol(derivate)
	    pterm = diff(sympify(safe), sy)
        elif inverse:
            if inverse in myterm:
                sy = Symbol(str('ssm___' + inverse))
                pterm = solve(sympify(safe), sy)
                if not pterm:
                    raise SsmError("can't find a solution to " + term + "=0 solving for " + inverse)
                elif len(pterm)!=1:
                    raise SsmError("no unique solution for " + term + "=0 solving for " + inverse)
                else:
                    pterm = pterm[0]
            else:
                pterm = sympify(safe)

        else:
            pterm = sympify(safe)

        #remove the ssm___ prefix
        #term = ccode(simplify(pterm)).replace('ssm___', '') ##NOTE simplify is just too slow to be used...
        term = ccode(pterm).replace('ssm___', '')

        #make the ssm C expression
        return self.generator_C(term, no_correct_rate, force_par=force_par, xify=xify, human=human, set_t0=set_t0)

    def __init__(self, convection_field, forcing_function):
        self._forcing_function = forcing_function
        self._convection_field = convection_field

        dx0 = sg.diff(convection_field[0], x)
        dx1 = sg.diff(convection_field[1], x)
        dy0 = sg.diff(convection_field[0], y)
        dy1 = sg.diff(convection_field[1], y)

        to_ccode = lambda u: sp.ccode(
            sympy.sympify(sg.symbolic_expression(u)._sympy_()))

        self._forcing_code = _forcing_template.format(
            **{'forcing_expr': to_ccode(forcing_function)})

        self._convection_code = _convection_template.format(**
            {'value0_expr': to_ccode(convection_field[0]),
             'value1_expr': to_ccode(convection_field[1]),
             'dx0_expr': to_ccode(dx0),
             'dx1_expr': to_ccode(dx1),
             'dy0_expr': to_ccode(dy0),
             'dy1_expr': to_ccode(dy1)})

        self.so_folder = tf.mkdtemp(prefix="cfem_") + os.sep
        # even if the object is not instantiated, we can still clean it up at
        # exit.
        atexit.register(lambda folder=self.so_folder: shutil.rmtree(folder))
        shutil.copy(_module_path + "makefile", self.so_folder)
        shutil.copy(_module_path + "cfem.h", self.so_folder)
        with open(self.so_folder + "funcs.c", 'w') as fhandle:
            fhandle.write("#include <math.h>\n")
            fhandle.write("#include \"cfem.h\"\n")
            fhandle.write(self._forcing_code)
            fhandle.write(self._convection_code)

        current_directory = os.getcwd()
        try:
            os.chdir(self.so_folder)
            util.run_make(command="autogenerated_functions")
            self.so = np.ctypeslib.load_library("libfuncs.so", "./")
            self.so.cf_forcing.restype = ct.c_double
            self.so.cf_forcing.argtypes = [ct.c_double, ct.c_double,
                                           ct.c_double]
            self.so.cf_convection.restype = CConvection
            self.so.cf_convection.argtypes = [ct.c_double, ct.c_double]
        finally:
            os.chdir(current_directory)

def expr_body(expr, **kwargs):
    if not hasattr(expr, '__len__'):
        # Defined in terms of some coordinates
        xyz = set(sp.symbols('x[0], x[1], x[2]'))
        xyz_used = xyz & expr.free_symbols
        assert xyz_used <= xyz
        # Expression params which need default values
        params = (expr.free_symbols - xyz_used) & set(kwargs.keys())
        # Body
        expr = ccode(expr).replace('M_PI', 'pi')
        # Default to zero
        kwargs.update(dict((str(p), 0.) for p in params))
        # Convert
        return expr
    # Vectors, Matrices as iterables of expressions
    else:
        return [expr_body(e, **kwargs) for e in expr]