def dampnewton(x,F,DF,q=0.5,tol=1e-10):
    cvg = []
    lup = lu_factor(DF(x))
    s = lu_solve(lup,F(x))
    xn = x-s
    lam = 1
    st = lu_solve(lup,F(xn)) # simplified Newton
    while norm(st) > tol*norm(xn):
        while norm(st) > (1-lam*0.5)*norm(s):
            lam *= 0.5
            if lam < 1e-10:
                cvg = -1
                print 'Failure of convergence'
                return x, cvg
            xn = x-lam*s
            st = lu_solve(lup,F(xn)) # simplified Newton
        cvg += [[lam, norm(xn), norm(F(xn))]]
        x = xn
        lup = lu_factor(DF(x))
        s = lu_solve(lup,F(x))
        lam = min(lam/q, 1.) # Wozu dieser Test?
        xn = x-lam*s
        st = lu_solve(lup,F(xn)) # simplified Newton
    x = xn
    return x, array(cvg)

def time_LU():
    """Print the times it takes to solve a system of equations using
    LU decomposition and (A^-1)B where A is 1000x1000 and B is 1000x500."""
    A = np.random.random((1000,1000))
    b = np.random.random((1000,500))
    start = time.time()
    L = la.lu_factor(A)
    a = time.time() - start
    start = time.time()
    A_inv = la.inv(A)
    a2 = time.time() - start
    start = time.time()
    la.lu_solve(L,b)
    a3 = time.time() - start
    start = time.time()
    np.dot(A_inv, b)
    a4 = time.time() - start

    
    time_lu_factor = a  # set this to the time it takes to perform la.lu_factor(A)
    time_inv = a2 # set this to the time it takes to take the inverse of A
    time_lu_solve = a3  # set this to the time it takes to perform la.lu_solve()
    time_inv_solve = a4  # set this to the time it take to perform (A^-1)B


    print "LU solve: " + str(time_lu_factor + time_lu_solve)
    print "Inv solve: " + str(time_inv + time_inv_solve)
    
    # What can you conclude about the more efficient way to solve linear systems?
    print "Better to use LU decomposition than inverse, cause it is NEVER a good idea to calculate an inerse"  # print your answer here.

    def solve_nonlinear(self, inputs, outputs):
        """
        Use numpy to solve Ax=b for x.

        Parameters
        ----------
        inputs : Vector
            unscaled, dimensional input variables read via inputs[key]
        outputs : Vector
            unscaled, dimensional output variables read via outputs[key]
        """
        vec_size = self.options['vec_size']
        vec_size_A = self.vec_size_A

        # lu factorization for use with solve_linear
        self._lup = []
        if vec_size > 1:
            for j in range(vec_size_A):
                lhs = inputs['A'][j] if vec_size_A > 1 else inputs['A']
                self._lup.append(linalg.lu_factor(lhs))

            for j in range(vec_size):
                idx = j if vec_size_A > 1 else 0
                outputs['x'][j] = linalg.lu_solve(self._lup[idx], inputs['b'][j])
        else:
            self._lup = linalg.lu_factor(inputs['A'])
            outputs['x'] = linalg.lu_solve(self._lup, inputs['b'])

 def updatedata(self, A):
     if self.update:
         if self.corr:
             self.data.B = self.data.w*(linalg.norm(A,1)/linalg.norm(self.data.w,1))
             self.data.C = self.data.v*(linalg.norm(A,Inf)/linalg.norm(self.data.v,1))
         else:
             # Note: Problem when singular vectors switch smallest singular value (See NewLorenz).
             #       To overcome this, I have implemented a 1e-8 random nudge.
             try:
                 ALU = linalg.lu_factor(A)
                 BC = linalg.lu_solve(ALU, c_[linalg.lu_solve(ALU, self.data.B + 1e-8*self.data.Brand), \
                                      self.data.C + 1e-8*self.data.Crand], trans=1)
                 C = linalg.lu_solve(ALU, BC[:,-1*self.data.q:])
                 B = BC[:,0:self.data.p]
             except:
                 if self.C.verbosity >= 1:
                     print 'Warning: Problem updating border vectors.  Using svd...'
                 U, S, Vh = linalg.svd(A)
                 B = U[:,-1*self.data.p:]
                 C = num_transpose(Vh)[:,-1*self.data.q:]
         
             bmult = cmult = 1
             if matrixmultiply(transpose(self.data.B), B) < 0:
                 bmult = -1
             if matrixmultiply(transpose(self.data.C), C) < 0:
                 cmult = -1
             self.data.B = bmult*B*(linalg.norm(A,1)/linalg.norm(B))
             self.data.C = cmult*C*(linalg.norm(A,Inf)/linalg.norm(C))

 def getVW(self, A):
     # V --> m, W --> n
     #print self.data
     MLU = linalg.lu_factor(c_[r_[A,transpose(self.data.C)], r_[self.data.B,self.data.D]])
     V = linalg.lu_solve(MLU,r_[zeros((self.data.n,self.data.q), Float), eye(self.data.q)])
     W = linalg.lu_solve(MLU,r_[zeros((self.data.m,self.data.p), Float), eye(self.data.p)],trans=1)
     
     return V, W

def solveLoop(A,LU,B,f):
    if f==True:
        # the fast way
        for i in xrange(B.shape[0]):
            la.lu_solve(LU,B[i])
    else:
        # the slow way
        for i in xrange(B.shape[0]):
            la.solve(A,B[i])

def Solve():
    # the students should have code to generate the random matrices, inverse, LU, and solve
    A = np.random.rand(1000,1000)
    B = np.random.rand(1000,500)
    Ai = la.inv(A)
    (lu,piv) = la.lu_factor(A)
    
    # the students should report the time for the following operations
    np.dot(Ai,B)
    la.lu_solve((lu,piv),B)

 def test(self, p):
     """ Determine which quadrilateral(s) each point is in
     Args:
         p: the points to test
     
     Returns:
         A N x len(p) boolean array, where N is the number of quadrilaterals
     """
     p1 = np.vstack((p.T, np.ones(len(p)))) # set up the points for the barycentric calculation
     bary1 = np.array(map(lambda lup: sl.lu_solve(lup, p1), self.lup1)) # calculate the barycentric coords for the first set of triangles
     bary2 = np.array(map(lambda lup: sl.lu_solve(lup, p1), self.lup2)) # ... and the second
     in1 = np.all((bary1 >=0) * (bary1 <=1), axis=1) # test that they are all in [0,1]
     in2 = np.all((bary2 >=0) * (bary2 <=1), axis=1)
     return in1+in2 # + is "or"

def do_problem_5(datafile):
    print_arr = lambda x,y: \
                print("{} =\n{}".format(y,
                                        np.array2string(x,precision = 6,
                                                        suppress_small = True,
                                                        separator=',')))
    np.set_printoptions(precision=6)
    A = loadtxt(datafile)
    (n,m) = A.shape
    (LU,p) = lup_decomp(A)
    (LU_control,p_control) = la.lu_factor(A)
    ## Check that my LU is equal to the actual LU, with a small
    ## tolerence for floating point rouding errors
    assert(np.allclose(LU,LU_control));
    
    L = np.tril(LU)
    U = np.triu(LU)
    P = np.zeros((n,n))
    for i in range(n):
        L[i,i] = 1
        P[i,p[i]] = 1
    print("Problem 5:")
    print("LUP decomposition")
    print_arr(L,"L")
    print_arr(U,"U")
    print_arr(P,"P")
    print("Solving Ax = b for various values of b")
    
    b1 = array([2,3,-1,5,7],dtype=float)
    x1 = lup_solve(LU,p,b1)
    x1_control = la.lu_solve((LU_control,p_control),b1)
    assert(np.allclose(x1,x1_control));
    print_arr(b1,"b1")
    print_arr(x1,"x1")
    
    b2 = array([15,29,8,4,-49],dtype=float)
    x2 = lup_solve(LU,p,b2)
    x2_control = la.lu_solve((LU_control,p_control),b2)
    assert(np.allclose(x2,x2_control));
    print_arr(b2,"b2")
    print_arr(x2,"x2")

    b3 = array([8,-11,3,-8,-32],dtype=float)
    x3 = lup_solve(LU,p,b3)
    x3_control = la.lu_solve((LU_control,p_control),b3 )
    assert(np.allclose(x3,x3_control));
    print_arr(b3,"b3")
    print_arr(x3,"x3")

    def solve_transpose(self, din):
        from scipy.linalg import lu_solve
        nrows = self.Alu.mshape[0]
        nblock = self.Alu.blockshape[0] #assume a square block!

        d = din.reshape((nrows,nblock))
        x = zeros(d.shape, dtype=self.Alu.dtype)
        
        #Forward substitution pass
        # b[0].T dnew[0] = d[0]
        # b[i].T dnew[i] = d[i] - c[i-1].T dnew[i-1]
        x[0] = d[0]
        for row in range(0,nrows):
            if row>0:
                x[row] = d[row] - dot(self.Alu.get_raw_block(row-1,2).T, x[row-1])
            if any(isnan(x[row])) or any(isinf(x[row])):
                print row, x[row]
            x[row] = lu_solve((self.Alu.get_raw_block(row,1),self.pivots[row]),\
                     x[row], trans=1)

        #Backward substitution
        # x[i] = d[i] - anew[i+1] x[i+1]
        for row in range(nrows-2,-1,-1):
            x[row] -= dot(self.Alu.get_raw_block(row+1,0).T, x[row+1])

        return x.reshape(din.shape)

    def SolveNextTime(self):
        r"""
        Calculate the next time (factorization)
        and update the time stack grids
        """
        
        try:
            self.tstep += 1
        except :
            self.tstep = 0
            self.LinearSystem()
            # gets the m factor from the solved system
            self.mUtfactor = ln.lu_factor(self.mUt)

        self.Source(self.tstep)
        # in time
        # As t is in [0, 1, 2] (2nd order)
        # time t in this case is Utime[2]
        v = self.Independent()
        result = ln.lu_solve(self.mUtfactor, v)
        # reshape the vector to became a matrix again
        self.Ufuture = result
        # make the update in the time stack
        # before [t-2, t-1,  t]
        # after  [t-1, t,  t+1]
        # so t-2 receive t-1 and etc.
        self.Uprevious = self.Ucurrent 
        self.Ucurrent = self.Ufuture
        
        return self.Ufuture

    def solve_linear(self, d_outputs, d_residuals, mode):
        r"""
        Back-substitution to solve the derivatives of the linear system.

        If mode is:
            'fwd': d_residuals \|-> d_outputs

            'rev': d_outputs \|-> d_residuals

        Parameters
        ----------
        d_outputs : Vector
            unscaled, dimensional quantities read via d_outputs[key]
        d_residuals : Vector
            unscaled, dimensional quantities read via d_residuals[key]
        mode : str
            either 'fwd' or 'rev'
        """
        if mode == 'fwd':
            sol_vec, forces_vec = d_outputs, d_residuals
            t = 0
        else:
            sol_vec, forces_vec = d_residuals, d_outputs
            t = 1

        sol_vec['disp_aug'] = linalg.lu_solve(self._lup, forces_vec['disp_aug'], trans=t)

def marglike(x, y, hyper, white_noise = False): # FIXME: build optional white noise into this kernel

    # Calculate covariance matrix
    K = matrifysquare(hyper, x, 0)
    K = 0.5* ( K + K.T) # Forces K to be perfectly symmetric

    # Calculate derivatives
    # dKdsigma = matrifysquare(hyper, x, 1) # Derivative w.r.t. log(sigma)
    # dKdlambda1 = matrifysquare(hyper, x, 2) # Derivative w.r.t. log(lambda1)
    # dKdh1 = matrifysquare(hyper, x, 3) # Derivative w.r.t. log(h1)

    sign, logdetK = np.linalg.slogdet(K)
    
    invKy = -0.5 * y.T *  np.mat(la.lu_solve(la.lu_factor(K),y)) \
        - 0.5 * logdetK - (y.size/2.) * np.log(2*np.pi)
    
    U = np.linalg.cholesky(K)

    n = len(x)
    L = - sum(np.log(np.diag(U))) -0.5 * y * invKy - n*0.5*np.log(2*np.pi)
    # dLdsigma = 0.5 * sum(np.diag(invKy*invKy.T*dKdsigma - (np.linalg.solve(K, dKdsigma)) ))
    # dLdlambda1 = 0.5 * sum(np.diag(invKy*invKy.T*dKdlambda1 - (np.linalg.solve(K, dKdlambda1)) ))
    # dKdh1 = 0.5 * sum(np.diag(invKy*invKy.T*dKdh1 - (np.linalg.solve(K, dKdh1)) ))

    return -L #, [-dKdsigma, -dKdlambda1, -dKdh1]

    def solve_overlap(self, b):
	"""
	x = solve_overlap(b)

	Solve for the overlap matrix: S x = b.

	Parameters
	----------
	b : 1D complex array.
	
	Returns
	-------
	x : 1D complex array.
	"""
	x = zeros(self.basis_size, dtype=complex)
	
	for i in range(self.el_basis_size):
	    
	    #Indices of a submatrix.
	    my_slice = slice(i*self.vib_basis_size, (i+1)*self.vib_basis_size)
	    
	    #Solve for the B-spline overlap matrix.
	    x[my_slice] = lu_solve(self.overlap_fact, b[my_slice])
	
	return x

def lnlike(h, X, y, covfunc):
    y = np.matrix(np.array(y).flatten()).T
    K = covfunc(X, X, h, wn = True)
    sign, logdetK = np.linalg.slogdet(K)
    alpha = np.mat(la.lu_solve(la.lu_factor(K),y))
    logL = -0.5*y.T * alpha - 0.5*logdetK - (y.size/2.)*np.log(2)
    return np.array(logL).flatten()

    def hotCmntsForTest(self, postId, nCmnts = 5):
        self.buildgraph(postId)
        
        testsizes = [shape(self.prg)[0], 800, 600, 400, 200]
        
        for size in testsizes:
            
            self.prg = self.prg[0:size,0:size]
            lil = lil_matrix(self.prg)
            
            start = clock()
            #eig  = eigs(self.prg, k=1, return_eigenvectors =False)
            eig = eigs(lil, return_eigenvectors =False, maxiter=10, tol=1E-5)
            eig = eig[0].real
            eig = 1/eig
            eigTime = clock() - start            
            print 'test_size:',size, 'eigTime:',eigTime        

            one = ones(size)
            m = eye(size) - eig*lil  
            
            start = clock()
            cmnts_ranking = lu_solve((m, one), one)
            solveTime = clock() - start
            
            print 'test_size:',size, 'solveTime:',solveTime

    def SolveNextTime(self):
        r"""
        Calculate the next time (factorization)
        and update the time stack grids
        """

        try:
            self.tstep += 1
        except :
            self.tstep = 0
            self.LinearSystem()
            # gets the m factor from the solved system
            self.mUtfactor = ln.lu_factor(self.mUt)

        # As t is in [0, 1, 2] (2nd order)
        # time t in this case is Utime[2]
        # the independent term of the matrix, due the pressure field
        v = self.Independent()

        result = ln.lu_solve(self.mUtfactor, v)
        # reshape the vector to become a matrix again
        self.Ufuture = np.reshape(result, (self.Nz, self.Nx))

        # make the update in the time stack
        # before [t-2, t-1,  t]
        # after  [t-1, t,  t+1]
        # so t-2 receive t-1 and etc.
        # make the update in the time stack
        self.Uprevious[:][:] = self.Ucurrent[:][:]
        self.Ucurrent[:][:] = self.Ufuture[:][:]        
        
        return self.Ufuture

def factiz(K):
    """
    Helper function to behave the same way scipy.sparse.factorized does, but
    for dense matrices.
    """
    luf = lu_factor(K)
    return lambda x: matrix(lu_solve(luf, x))

def invpower(e,A,B=None):
	if B is None:
		B = eye(len(A))
	K = A - B*e
	G = lu_factor(K)
	x = ones([len(A),1],complex)/len(A)
	iter = 0
	error = 10
	#for i in range(8):
	while error > 1e-8 and iter < 20:
		try:
			x = lu_solve(G,x)
		    
		except:
			print 'LU Exception'
			x = solve(K,x)
		    
		x = x/norm(x,ord=inf)
		error = norm(dot(A,x)-e*dot(B,x))
		iter = iter +1
		print 'invpower error = ',error
	x = x*conj(x[0])
	print 'Eval = ',e
	print 'Evect Real = ',norm(real(x))
	print 'Evect Imag = ',norm(imag(x))
	return x

def root_finding_newton(u, m, alpha, V, Sigma, U, G, Cov, dw):
	"""
	:param u: initial vector of u
	:param m: vector of default model for the Lehmann spectral function
	:param alpha: scalar value, controls the relative weight of maximizing entropie and minimizing kind. of least squares fit.
	:param V: part of singular SVD of K, K = V*Sigma*U.T
	:param Sigma: part of singular SVD of K, K = V*Sigma*U.T
	:param U: part of singular SVD of K, K = V*Sigma*U.T
	:param G: Vector of all average values of the Greensfunction for each time step
	:param Cov: Vector with variance of the re-binned, gaussian distributed QMC approximations for the different time-steps
	:param dw: omega step size
	:return:
	"""

	s=len(u)
	max_val = np.sum(m)
	K_s = np.dot(V,np.dot(Sigma,U.T))
	diff = 1.
	count1 = 1
	max_iter = 1000

	while diff > 1e-8 and count1 <= max_iter:
		print(count1)
		A_appr = m * np.exp(np.dot(U,u))
		A_old = A_appr
		inv_cov = (1. / np.diagonal(Cov)**2)
		inv_cov_mat = np.diag(inv_cov)
		dLdF = - inv_cov * (G - np.dot(K_s, A_appr))
		F_u = - alpha * u - np.dot(Sigma,np.dot(V.T,dLdF))
		M = np.dot(Sigma,np.dot(V.T,np.dot(inv_cov_mat,np.dot(V,Sigma))))
		T = np.dot(U.T,np.dot(np.diag(A_appr),U))
		J = alpha * np.diag(np.ones((s))) + np.dot(M,T)
		lu_and_piv = lu_factor(J)
		delta_u = lu_solve(lu_and_piv,F_u)
		A_appr = m * np.exp(np.dot(U,u + delta_u))
		count2 = 1
		while np.linalg.norm(A_appr - A_old) > max_val and count2 <= max_iter:
			J = (alpha+count2*1e10) * np.diag(np.ones((s))) + np.dot(M,T)
			lu_and_piv = lu_factor(J)
			delta_u = lu_solve(lu_and_piv,F_u)
			A_appr = m * np.exp(np.dot(U,u + delta_u))
			count2 +=1
		u_old = u 
		u = u + delta_u
		diff = np.abs(np.sum(u-u_old))
		count1 += 1
	return u