def main():
	A = sympy.Matrix([[3, 11], [22, 6]])
	b = sympy.Matrix([0, 0])
	print("Исходная матрица A:")
	sympy.pprint(A)
	x = gauss.gauss(A, b)
	if x:
		print("\nРезультат:", x)

def main():
	A = sympy.Matrix([[1, -3, -1], [2, -2, 1], [3, -1, 2]])
	b = sympy.Matrix([-11, -7, -4])
	print("Исходная матрица A:")
	sympy.pprint(A)
	print("\nВектор b:")
	sympy.pprint(b)
	x = gauss.gauss(A, b)
	if x:
		print("\nРезультат:", x)

def showCookToomFilter(a,n,r,fractionsIn=FractionsInG):

    AT,G,BT,f = cookToomFilter(a,n,r,fractionsIn)

    print "AT = "
    pprint(AT)
    print ""

    print "G = "
    pprint(G)
    print ""

    print "BT = "
    pprint(BT)
    print ""

    if fractionsIn != FractionsInF:
        print "AT*((G*g)(BT*d)) ="
        pprint(filterVerify(n,r,AT,G,BT))
        print ""

    if fractionsIn == FractionsInF:
        print "fractions = "
        pprint(f)
        print ""

def showCookToomConvolution(a,n,r,fractionsIn=FractionsInG):

    AT,G,BT,f = cookToomFilter(a,n,r,fractionsIn)

    B = BT.transpose()
    A = AT.transpose()
    
    print "A = "
    pprint(A)
    print ""

    print "G = "
    pprint(G)
    print ""

    print "B = "
    pprint(B)
    print ""

    if fractionsIn != FractionsInF:
        print "Linear Convolution: B*((G*g)(A*d)) ="
        pprint(convolutionVerify(n,r,B,G,A))
        print ""

    if fractionsIn == FractionsInF:
        print "fractions = "
        pprint(f)
        print ""

def tests():
    x, y, z = symbols('x,y,z')
    #print(x + x + 1)
    expr = x**2 - y**2
    factors = factor(expr)
    print(factors, " | ", expand(factors))
    pprint(expand(factors))

def quick_equation(eq = None, latex = None):
	"""
	first print the given equation. optionally, in markdown mode, generate an
	image from an equation and write a link to it.
	"""
	if eq is not None:
		sympy.pprint(eq)
	else:
		print latex
	# print an empty line
	print
	if not markdown:
		return

	global plt
	if latex is None:
		latex = sympy.latex(eq)
	img_filename = hashlib.sha1(latex).hexdigest()[: 10]

	# create the figure and hide the border. set the height and width to
	# something far smaller than the resulting image - bbox_inches will
	# expand this later
	fig = plt.figure(figsize = (0.1, 0.1), frameon = False)
	ax = fig.add_axes([0, 0, 1, 1])
	ax.axis("off")
	fig = plt.gca()
	fig.axes.get_xaxis().set_visible(False)
	fig.axes.get_yaxis().set_visible(False)
	plt.text(0, 0, r"$%s$" % latex, fontsize = 25)
	plt.savefig("img/%s.png" % img_filename, bbox_inches = "tight")
	# don't use the entire latex string for the alt text as it could be long
	quick_write("![%s](img/%s.png)" % (latex[: 20], img_filename))

def invertInteractive(f, X, Y):
  """
  Given function Y = f(X) return inverse f^{-1}(X).

  Note:
  - f could be fcn of multiple variables. This inverts in variable X
  - f must be a sympy function of sympy symbol/variable X
  - Returns function of Y and anything else f depends upon, removing
    dependence upon X
  - Accounts for non-unique inverse by prompting user to select
    from possible inverses
  """

  # Try to construct inverse
  finv = sp.solve(Y-f, X)

  # If inverse unique, return it. If not, prompt user to select one
  if len(finv) == 1:
    finv = finv[0]

  else:
    print "Please choose expression for the inverse X = f^{-1}(Y)\n"
    for opt, expr in enumerate(finv):
      print "\nOption %d:\n----------\n" % opt
      sp.pprint(expr)
    choice = int(raw_input("Your Choice: "))
    finv = finv[choice]
    os.system('clear')

  return finv

def case1():
    a, x1, x2 = symbols('a x1 x2')              # Define symbols
    b = 3 * (a - 1)                             # Relation between a and b
    y = a * x1 + b * x2 + 5                     # Given function y
    # Error of x1 and x2
    Sigx1 = 0.8
    Sigx2 = 1.5

    print "In the case that x1 and x2 are dependent:\n"
    Sigx1x2 = Sigx1 * Sigx2     # Correlation coefficient equals to 1
    SigmaXX = np.matrix([[Sigx1**2, Sigx1x2],   # Define covariance
                         [Sigx1x2, Sigx2**2]])  # matrix of x1 and x2
    print "SigmaXX =\n"
    pprint(Matrix(SigmaXX))
    print

    # Jacobian matrix of y function with respect to x1 and x2
    Jyx = np.matrix([diff(y, x1), diff(y, x2)])

    SigmaYY = Jyx * SigmaXX * Jyx.T   # Compute covariance matrix of y function

    # Compute solution in witch the SigmaYY is minimum
    a = round(solve(diff(expand(SigmaYY[0, 0]), a))[0], 8)
    b = round(b.evalf(subs={'a': a}), 8)

    # Compute sigma y
    SigmaY = sqrt(SigmaYY[0, 0].evalf(subs={'a': a, 'b': b}))

    print "a = %.4f\nb = %.4f" % (a, b)
    print "SigmaY = %.4f" % float(SigmaY)

def test1():
    """
    Wiedemann's Formel hat einen Fehler!
    """
    print("\n\n...................",    """Wiedemann's Formel hat einen Fehler!""")
    fs,ld = S.symbols('fs ld')
    QF  = S.Matrix([[1,0],[-1/fs,1]])
    DR  = S.Matrix([[1,ld],[0,1]])
    QD  = S.Matrix([[1,0],[1/fs,1]])
    FODO1 = QD*DR*QF
    FODO2 = FODO1.subs(fs, -fs)
    
    FODO = FODO1*FODO2
    FODO = S.expand(FODO)
    print('... this is the correct one:')
    S.pprint(FODO)

    T = 25.           # kin. energy [MeV]
    T=T*1.e-3         # [GeV] kin.energy 
    betakin=M.sqrt(1.-(1+T/E0)**-2)   # beta kinetic
    E=E0+T            # [GeV] energy
    k=0.2998*Bg/(betakin*E)    # [1/m^2]
    L=0.596                    # distance between quads [m]
    Ql=0.04                    # full quad length [m]
    f = 1./(k*Ql)
    fodo1 = FODO.subs(fs,2*f).subs(ld,L/2.)
    fodo1 = np.array(fodo1)
    fodo2 = Mfodo(2*f,L/2.)         # Wiedemann  (the book is wrong!!)
    my_debug('matrix probe: fodo1-fodo2 must be zero matrix')
    zero = fodo1-fodo2
    my_debug(f'{abs(zero[0][0])}  {abs(zero[0][1])}')
    my_debug(f'{abs(zero[1][0])}  {abs(zero[1][1])}')
    return

    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 main():

    print 'Initial metric:'
    pprint(gdd)
    print '-'*40
    print 'Christoffel symbols:'
    for i in [0,1,2,3]:
    	for k in [0,1,2,3]:
    		for l in [0,1,2,3]:
    			if Gamma.udd(i,k,l) != 0 :
		    		pprint_Gamma_udd(i,k,l)
    print'-'*40
    print'Ricci tensor:'
    for i in [0,1,2,3]:
    	for j in [0,1,2,3]:
    		if Rmn.dd(i,j) !=0:
    			pprint_Rmn_dd(i,j)
    
    print '-'*40
    #Solving EFE for A and B
    s = ( Rmn.dd(1,1)/ A(r) ) + ( Rmn.dd(0,0)/ B(r) )
    pprint (s)
    t = dsolve(s, A(r))
    pprint(t)
    metric = gdd.subs(A(r), t)
    print "metric:"
    pprint(metric)
    r22 = Rmn.dd(3,3).subs( A(r), 1/B(r))
    h = dsolve( r22, B(r) )
    pprint(h)

def p1_5(part=None):
    """
    Complete solution to problem 1.5

    Parameters
    ----------
    part: str, optional(default=None)
        The part number you would like evaluated. If this is left blank
        the default value of None is used and the entire problem will be
        solved.

    Returns
    -------
    i-dont-know-yet
    """
    f, fp, fp2, fp3, fp4, h, fplus, fminus = symbols('f, fp, fp2, fp3, fp4, \
                                                     h, fplus, fminus')
    eqns = (Eq(fplus, f + fp * h + fp2 * h ** 2 / 2 + \
                      fp3 * h ** 3 / 6 + fp4 * h ** 4 / 24),
            Eq(fminus, f - fp * h + fp2 * h ** 2 / 2 - \
                       fp3 * h ** 3 / 6 + fp4 * h ** 4 / 24))
    soln = solve(eqns, fp, fp2)
    pprint(soln)

    pass  # return nothing

 def display(self,opt='repr'):
     """
     Procedure Name: display
     Purpose: Displays the random variable in an interactive environment
     Arugments:  1. self: the random variable
     Output:     1. A print statement for each piece of the distribution
                     indicating the function and the relevant support
     """
     if self.ftype[0] in ['continuous','Discrete']:
         print ('%s %s'%(self.ftype[0],self.ftype[1]))
         for i in range(len(self.func)):
             cons_list=['0<'+str(cons) for cons in self.constraints[i]]
             cons_string=', '.join(cons_list)
             print('for x,y enclosed in the region:')
             print(cons_string)
             print('---------------------------')
             pprint(self.func[i])
             print('---------------------------')
             if i<len(self.func)-1:
                 print(' ');print(' ')
         
     if self.ftype[0]=='discrete':
         print '%s %s where {x->f(x)}:'%(self.ftype[0],
                                         self.ftype[1])
         for i in range(len(self.support)):
             if i!=(len(self.support)-1):
                 print '{%s -> %s}, '%(self.support[i],
                                       self.func[i]),
             else:
                 print '{%s -> %s}'%(self.support[i],
                                     self.func[i])

def getRelError(variables,func,showFunc = False):
	"""getError generates a function to calculate the error of a function. I.E. a function that
	gives you the error in its result given the error in its input. Output function is numpy ready. 
	Output function will take twice as many args as variables, one for the var and one for the error in that var.

	arguments
	variables      : a list of sympy symbols in func
	errorVariables : list of sympy ymbols representing the error in each value of variables. must be the same length as variables
	func           : a function containing your variables that you want the error of"""
	ErrorFunc = 0 # need to set function to start value
	erVars    = {}
	for i in range(len(variables)): #run through all variables in the function
		v                  = variables[i]
		dv                 =  Symbol('d'+str(v), positive = True)
		erVars['d'+str(v)] = dv
		D                  = (diff(func,v)*dv)**2
		ErrorFunc          += D

	ErrorFunc = sqrt(ErrorFunc)/func
	if showFunc:
		pprint(ErrorFunc)
		
	variables.extend(erVars.values()) #create a list of all sympy symbols involved
	func      = lambdify(tuple(variables), ErrorFunc ,"numpy") #convert ErrorFunc to a numpy rteady python lambda
	return(func)

def print_formulas(G):
    for node in G.nodes():
        print("Node {}:".format(node))
        pprint(G.formula_conj(node))
        #pprint(simplify(G.formula_conj(node)))
        #pprint(to_cnf(G.formula_conj(node), simplify=True))
    print("\n")

def setup_loss(loss_type, native_lang = "C", verbose=False):
    """
    set up symbolic derivations
    """

    print "setting up autowrapped loss", loss_type

    x = x = sympy.Symbol("x")
    y = y = sympy.Symbol("y")
    z = z = sympy.Symbol("z")
    cx = cx = sympy.Symbol("cx")
    cy = cy = sympy.Symbol("cy")
    wx = wx = sympy.Symbol("wx")
    wy = wy = sympy.Symbol("wy")
    rx = rx = sympy.Symbol("rx")
    ry = ry = sympy.Symbol("ry")

    if loss_type == "algebraic_squared":
        fun = ((x-(cx + z*wx))**2/(rx*rx) + (y-(cy + z*wy))**2/(ry*ry) - 1)**2

    if loss_type == "algebraic_abs":
        fun = sympy.sqrt(((x-(cx + z*wx))**2/(rx*rx) + (y-(cy + z*wy))**2/(ry*ry) - 1)**2 + 0.01)

    #TODO replace x**2 with x*x
    fun = fun.expand(deep=True)

    if verbose:
        sympy.pprint(fun)

    d_cx = fun.diff(cx).expand(deep=True)
    d_cy = fun.diff(cy).expand(deep=True)
    d_wx = fun.diff(wx).expand(deep=True)
    d_wy = fun.diff(wy).expand(deep=True)
    d_rx = fun.diff(rx).expand(deep=True)
    d_ry = fun.diff(ry).expand(deep=True)

    if native_lang == "fortran":
        c_fun = autowrap(fun, language="F95", backend="f2py")
        c_d_cx = autowrap(d_cx)
        c_d_cy = autowrap(d_cy)
        w_d_wx = autowrap(d_wx)
        w_d_wy = autowrap(d_wy)
        c_d_rx = autowrap(d_rx)
        c_d_ry = autowrap(d_ry)
    else:
        c_fun = autowrap(fun, language="C", backend="Cython", tempdir=".")
        c_d_cx = autowrap(d_cx, language="C", backend="Cython", tempdir=".")
        c_d_cy = autowrap(d_cy, language="C", backend="Cython", tempdir=".")
        c_d_wx = autowrap(d_wx, language="C", backend="Cython", tempdir=".")
        c_d_wy = autowrap(d_wy, language="C", backend="Cython", tempdir=".")
        c_d_rx = autowrap(d_rx, language="C", backend="Cython", tempdir=".")
        c_d_ry = autowrap(d_ry, language="C", backend="Cython", tempdir=".")

    grads = {"cx": d_cx, "cy": d_cy, "wx": d_wx, "wy": d_wy, "rx": d_rx, "ry": d_ry}
    c_grads = {"cx": c_d_cx, "cy": c_d_cy, "wx": c_d_wx, "wy": c_d_wy, "rx": c_d_rx, "ry": c_d_ry}

    print "done setting up autowrapped loss"

    return c_fun, c_grads

def print_formulas(G, label=None):
    if label:
        print("{}".format(label))
        print("----------------------")
    for node in G.nodes():
        print("Node {}:".format(node))
        pprint(G.formula_conj(node))
    print("\n")

def main():
    A = sympy.Matrix([[3, -1],
                      [-6, 2]])
    b = sympy.Matrix([0, 0])
    print("\nInitial matrix:")
    sympy.pprint(A)
    x = gauss.sym_gauss(A, b)
    sympy.pprint(x)

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)