def pprint(self, **settings):
        """
        Print the PHS structure :math:`\\mathbf{b} = \\mathbf{M} \\cdot \\mathbf{a}`
        using sympy's pretty printings.

        Parameters
        ----------

        settings : dic
            Parameters for sympy.pprint function.

        See Also
        --------

        sympy.pprint

        """

        sympy.init_printing()

        b = types.matrix_types[0](self.dtx() +
                                  self.w +
                                  self.y +
                                  self.cy)

        a = types.matrix_types[0](self.dxH() +
                                  self.z +
                                  self.u +
                                  self.cu)

        sympy.pprint([b, self.M, a], **settings)

def hst(n=16):
    # %pylab inline # command in ipython that imports sci modulues
    import sympy as sym
    sym.init_printing()
    for k in range(n+1):
        n = float(n)
        print('HST63%02d = %s' %  ( k, sym.Rational(k,n*2) ) )  

def in_mm(n=16):
    # %pylab inline # command in ipython that imports sci modulues
    import sympy as sym
    sym.init_printing()
    for k in range(n+1):
        n = float(n)
        print(' %5s in - %.6f in - %.6f mm ' %  (sym.Rational(k,n) , k/n, 25.4*k/n ) )  

def setup_pprint():
    from sympy import pprint_use_unicode, init_printing

    # force pprint to be in ascii mode in doctests
    pprint_use_unicode(False)

    # hook our nice, hash-stable strprinter
    init_printing(pretty_print=False)

def triangle_length(l):
    sympy.init_printing()
    for n in range(1, 17):
        print("Triangle: ", n)
        print (l)
        sympy.pprint (l * sympy.sqrt(n), use_unicode=True) 
        sympy.pprint (l * sympy.sqrt(n+1), use_unicode=True) 
    print ("-------------------------")

def print_series(n):
    # initialize printing system with
    # reverse order
    init_printing(order='rev-lex')
    x = Symbol('x')
    series = x
    for i in range(2, n+1):
        series = series + (x**i)/i
    pprint(series)

def main():
    if sympy.__version__ != "1.0":
        sys.exit("The doctests must be run against SymPy 1.0. Please install SymPy 1.0 and run them again.")
    full_text = ""

    for file in files:
        with open(file, 'r') as f:
            s = f.read()
            st = s.find(begin)
            while st != -1:
                if not (st >= len(skip)) or s[st - len(skip) : st] != skip:
                    full_text += s[st + len(begin) : s.find(end, st)]
                st = s.find(begin, st+ len(begin))

    full_text = full_text.replace(r'\end{verbatim}', '')

    with open(output_file, "w") as f:
        f.write("'''\n %s \n'''" % full_text)

    # force pprint to be in ascii mode in doctests
    pprint_use_unicode(False)

    # hook our nice, hash-stable strprinter
    init_printing(pretty_print=False)

    import test_full_paper

    # find the doctest
    module = pdoctest._normalize_module(test_full_paper)
    tests = SymPyDocTestFinder().find(module)
    test = tests[0]

    runner = SymPyDocTestRunner(optionflags=pdoctest.ELLIPSIS |
            pdoctest.NORMALIZE_WHITESPACE |
            pdoctest.IGNORE_EXCEPTION_DETAIL)
    runner._checker = SymPyOutputChecker()

    old = sys.stdout
    new = StringIO()
    sys.stdout = new

    future_flags = __future__.division.compiler_flag | __future__.print_function.compiler_flag

    try:
        f, t = runner.run(test, compileflags=future_flags,
                          out=new.write, clear_globs=False)
    except KeyboardInterrupt:
        raise
    finally:
        sys.stdout = old

    if f > 0:
        print(new.getvalue())
        return 1
    else:
        return 0

def print_series(n, x_value):
	init_printing(order = 'rev-lex')
	x = Symbol('x')
	series = x
	for i in range(2, n+1):
		series = series + (x**i)/i
	pprint(series)

	series_value = series.subs({x:x_value})
	print('Value of the series at {0}: {1}'.format(x_value, series_value))

def printSeries(n, val):
    init_printing(order='rev-lex')
    x = Symbol('x')
    expr = x
    for i in range(2, n + 1):
        expr += (x**i / i)
    pprint(expr)

    res = expr.subs({x: val})
    print(res)

	def new_MM(self):
		init_printing()
		
		m, n, k, self.tag = self.__choose_problem_type__()
		
		self.A = self.random_integer_matrix( m, k )
		self.x = self.random_integer_matrix( k, n )
		self.answer = self.A * self.x
		
		return Math( "$$" + latex( MatMul( self.A, self.x ), mat_str = "matrix" ) + "=" + "?" + "$$" )

	def new_problem(self):
		init_printing()

		m = self.random_integer( 3, 4 )
		n = self.random_integer( 2, 3 )

		self.A = self.random_integer_matrix( m, n )
		self.x = self.random_integer_matrix( n, 1 )
		self.answer = self.A * self.x

		return Math( "$$" + latex( MatMul( self.A, self.x ), mat_str = "matrix" ) + "=" + "?" + "$$" )

def print_series(n, x_value):
    # initialize printing system with
    # reverse order
    init_printing(order='rev-lex')
    x = Symbol('x')
    series = x
    for i in range(2, n+1):
        series = series + (x**i)/i
    pprint(series)
    # evaluate the series at x_value
    series_value = series.subs({x:x_value})
    print('Value of the series at {0}: {1}'.format(x_value, series_value))

 def init_sympy_printing(self, line):
     from sympy import init_printing
     init_printing()
     ip = get_ipython()
     latex = ip.display_formatter.formatters["text/latex"]
     for key in latex.type_printers.keys():
         if key.__module__ == "__builtin__":
             del latex.type_printers[key]
     png = ip.display_formatter.formatters["image/png"]
     for key in png.type_printers.keys():
         if key.__module__ == "__builtin__":
             del png.type_printers[key]

def exer_mult_step(expression):
    import re
    import sympy as sym
    from sympy import init_printing; init_printing() 
    from IPython.display import display, Math, Latex
    a, b, c, x, y, z = sym.symbols("a b c x y z")
    expression=str(expression)
    evexpr=sym.expand(eval(expression))
    print 'Expand the following expression step by step:'
    display(Math(str(expression).replace('**','^').replace('*','')))
    print 'Multiply ALL the monomials:'
    factor=re.split('\)\s*\*\s*\(',str(expression))
    refactor=[ re.split('\(|\)',factor[i])[1-i] for i in range(len(factor))]
    terms=[ re.split(':',re.sub('\s*(\+|\-)\s*',r':\1',refactor[i])) for i in range(len(refactor)) ]
    expandednr=''
    for i in range(len(terms[0])):
    #nci=''
        if terms[0][i]:
            nci=terms[0][i]
            terms[0][i]=r'{\color{red}{%s}}' %terms[0][i]
            for j in range(len(terms[1])):
            #ncj=''
                if terms[1][j]:
                    ncj=terms[1][j]
                    terms[1][j]=r'\color{red}{%s}' %terms[1][j]
        
                    l=jointerms(terms)
                    display(Math(l.replace('**','^').replace('*','')))
                    terms[1][j]=ncj
                    newterm=raw_input()
                    if re.search('^-',newterm):
                        newterm=newterm
                    else:
                        newterm='+'+newterm
                    expandednr=expandednr+newterm
                    expandednr=re.sub('^\+','',expandednr)
                
            terms[0][i]=nci  
        
    print 'Reduce the expression:' 
    display(Math(expandednr.replace('**','^').replace('*','')))
    result=raw_input()
    if evexpr== eval(result): 
        print "Good job!, final result:"
        return eval(result)
        
    else:
        print('WRONG!!!\nRight result:')
        return sym.expand(expression)

def ItosLemma(f, drift, diffusion, finv=None, replaceX=True, pretty_print=False):
  """
    This function applies Ito's lemma given a twice-differentiable,
    invertible function (f) and drift and diffusion coefficients for a
    stochastic process (X_t) satisfying SDE

      dX = drift(X,t) dt + diffusion(X,t) dW

    This routine then prints out a new the SDE for Y=f(X)

      dY = drift_new(X,t) dt + diffusion_new(X,t) dW_t

    If replaceX=True, this will replace instances of X with f^{-1}(Y)

    NOTE: Inputs f, drift, diffusion _all_ expected to be sympy fcns.
    Derivatives are taken with respect to sympy symbols/variables X and
    t in those functions.
  """

  # Define symbols for display later
  t, Y, dt, dW, dY, dX, Y_t = sp.symbols('t Y dt dW dY dX Y_t')

  # Define differentiation functions
  DX, Dt = map(lambda var: (lambda fcn: sp.diff(fcn, var)), [X, t])

  # Compute drift_new and diffusion_new according to Ito's Lemma
  drift_new     = Dt(f) + DX(f)*b + 0.5*DX(DX(f))*(s**2)
  diffusion_new = DX(f)*s

  # If finv not given, invert Y=f(X) to get Y = f^{-1}(X)
  if finv == None:
    finv = invertInteractive(f, X, Y)

  # Substitute out X with f^{-1}(Y)
  if replaceX:
    drift_new, diffusion_new = map(lambda fcn: fcn.subs(X,finv).simplify(),
                                   [drift_new, diffusion_new])

  # Print the results to screen
  if pretty_print:
    sp.init_printing()
  print "\nOriginal SDE\n------------\n"
  sp.pprint(sp.Eq(dX, (drift)*dt + (diffusion)*dW))
  print "\n\nNew SDE for Y = f(X)\n--------------------\n"
  sp.pprint(sp.Eq(dY,(drift_new)*dt + (diffusion_new)*dW))
  print "\nwhere"
  sp.pprint(sp.Eq(Y,f))
  print("\n\n")

def exer_mult(expression):
    import sympy as sym
    from sympy import init_printing; init_printing() 
    from IPython.display import display, Math, Latex
    a, b, c, x, y, z = sym.symbols("a b c x y z")
    expression=str(expression)
    evexpr=sym.expand(eval(expression))
    print 'Expand the expression:'
    display(Math(str(expression).replace('**','^').replace('*','')))
    result=raw_input()
    if evexpr== eval(result): 
        print "Good job! Final result:"
        return eval(result)
    else:
        print('WRONG!!!\nRight result:')
        return sym.expand(eval(expression))

def init_from_matmodlab_magic(p):
    if p == BOKEH:
        from bokeh.plotting import output_notebook
        output_notebook()
        environ.plotter = BOKEH
        i = 2
    elif p == MATPLOTLIB:
        environ.plotter = MATPLOTLIB
        i = 1

    environ.notebook = i
    environ.verbosity = 0
    try:
        from sympy import init_printing
        init_printing()
    except ImportError:
        pass

def Format(Fmode=True, Dmode=True, dop=1, inverse='full'):
    """
    Set modes for latex printer -

        Fmode:  Suppress function arguments (True)          Use sympy latex for functions (False)
        Dmode:  Use compact form of derivatives (True)      Use sympy latex for derivatives (False)

    and redirects printer output so that latex compiler can capture it.
    """
    global Format_cnt

    GaLatexPrinter.Dmode = Dmode
    GaLatexPrinter.Fmode = Fmode
    GaLatexPrinter.inv_trig_style = inverse

    if Format_cnt == 0:
        Format_cnt += 1

        """
        if metric.in_ipynb():
            GaLatexPrinter.ipy = True
        else:
            GaLatexPrinter.ipy = False

        GaLatexPrinter.dop = dop
        GaLatexPrinter.latex_flg = True

        if not GaLatexPrinter.ipy:
            GaLatexPrinter.redirect()
        """
        GaLatexPrinter.dop = dop
        GaLatexPrinter.latex_flg = True
        GaLatexPrinter.redirect()

        Basic.__str__ = lambda self: GaLatexPrinter().doprint(self)
        Matrix.__str__ = lambda self: GaLatexPrinter().doprint(self)
        Basic.__repr_ = lambda self: GaLatexPrinter().doprint(self)
        Matrix.__repr__ = lambda self: GaLatexPrinter().doprint(self)

        if isinteractive():
            init_printing(use_latex= 'mathjax')

    return

def exer_mult_step_old(expression):
    import re
    import sympy as sym
    from sympy import init_printing; init_printing() 
    from IPython.display import display, Math, Latex
    a, b, c, x, y, z = sym.symbols("a b c x y z")
    expression=str(expression)
    evexpr=sym.expand(eval(expression))
    print 'Expand the following expression step by step:'
    display(Math(str(expression).replace('**','^').replace('*','')))
    print 'Multiply ALL the monomials:'
    factor=re.split('\)\s*\*\s*\(',str(expression))
    refactor=[ re.split('\(|\)',factor[i])[1-i] for i in range(len(factor))]
    terms=[ re.split(':',re.sub('\s*(\+|\-)\s*',r':\1',refactor[i])) for i in range(len(refactor)) ]
    expandednr=''
    for i in terms[0]:
        for j in terms[1]:
            if i and j: #mo empty strings
                print '(%s)*(%s)' %(i,j)
                newterm=raw_input()
                if re.search('^-',newterm):
                    newterm=newterm
                else:
                    newterm='+'+newterm
                expandednr=expandednr+newterm
                expandednr=re.sub('^\+','',expandednr)
    print 'Reduce the expression:' 
    display(Math(expandednr.replace('**','^').replace('*','')))
    result=raw_input()
    if evexpr== eval(result): 
        print "Good job!, final result:"
        return eval(result)
        
    else:
        print('WRONG!!!\nRight result:')
        return sym.expand(expression)
           

 def __init__(self, coordinates, metric):
     sp.init_printing()
     if len(coordinates)==metric.rows and len(coordinates)==metric.cols:
         self.gdd = metric #stands for down-down components of the metric g
         self.x = coordinates
         self.dim = len(self.x)
         self.index = range(self.dim)
         self.guu = metric.inv()
         
         self.connection = [[[0 for i in self.index] for i in self.index] for i in self.index]
         for i in self.index:
             for j in self.index:
                 for k in self.index:
                     for m in self.index:
                         self.connection[i][j][k] += self.guu[i,m]*(self.gdd[m,j].diff(self.x[k]) + self.gdd[m,k].diff(self.x[j]) - self.gdd[k,j].diff(self.x[m]))/2
         self.Riemanntensor = [[[[0for i in self.index]for i in self.index]for i in self.index]for i in self.index]
         for i in self.index:
             for j in self.index:
                 for k in self.index:
                     for m in self.index:
                         self.Riemanntensor[i][j][k][m] +=  self.connection[i][j][k].diff(self.x[m]) - self.connection[i][j][m].diff(self.x[k])
                         for n in self.index:
                             self.Riemanntensor[i][j][k][m] += self.connection[n][j][k]*self.connection[i][m][n] - self.connection[n][j][m]*self.connection[i][k][n]
         
         self.Riccitensor = [[0 for i in self.index]for i in self.index]
         for i in self.index:
             for j in self.index:
                 for k in self.index:
                     self.Riccitensor[i][j] += self.Riemanntensor[k][i][k][j]
                     
         self.Riemannscalar = 0 
         for i in self.index:
             for j in self.index:
                 self.Riemannscalar += self.guu[i,j]*self.Riccitensor[i][j]           
     else:
         print "Creation failed. Please use an n-tuple of coordinate names with a square, n times n matrix"