def solve(A, b, x0, eps=1E-8, iter_max=1000):
    'Solve system Ax=b with Jacobi method.'

    # Checks of input
    m, n = A.shape
    assert m == n
    assert b.shape == (n, )

    # To prevent some integer division surprises with inverse
    A = A.astype('float', copy=False)
    b = b.astype('float', copy=False)

    # Get Jacobi matrix
    # If the diagonal is all zero, there is no way to continues
    if la.norm(A.diagonal()) < 1E-15:
        print 'Is diag(A) all zeros?'
        return x0
    else:
        diag_inv = np.array([1./d if abs(d) > 1E-14 else 0
                             for d in A.diagonal()])
        D_inv = np.diag(diag_inv)

    E = A - np.diag(A.diagonal())
    if la.norm(E) < 1E-13:
        E = np.eye(n)

    B = -D_inv.dot(E)

    # compute the rhs constant term
    z = D_inv.dot(b)

    n_iters = 0
    r = b - A.dot(x0)       # Initialize residuum
    r_norm = r.dot(r)
    r_norm0 = r_norm        # Remember size for stopping loop

    # Compare the spectra
    A_vals = la.eigvals(A)
    B_vals = la.eigvals(B)

    print 'A eigenvalues', A_vals
    print 'B eigenvalues', B_vals

    if not np.all(np.abs(B_vals) < 1):
        print '\tSome eigenvalues of B are greate than one in magnitude!'

    while n_iters < iter_max and r_norm > eps*r_norm0:
        n_iters += 1

        # Get the new solution
        x0 = B.dot(x0) + z

        # Compute norm of residuum
        r = b - A.dot(x0)
        r_norm = r.dot(r)

    return x0, n_iters

def linear_algebra():
    """ Use the `numpy.linalg` library to do Linear Algebra 
        For a reference on math, see 'Linear Algebra explained in four pages'
        http://minireference.com/static/tutorials/linear_algebra_in_4_pages.pdf
    """

    ### Setup two vectors
    x = np.array([1, 2, 3, 4])
    y = np.array([5, 6, 7, 8])

    ### Vector Operations include addition, subtraction, scaling, norm (length),
    # dot product, and cross product
    print np.vdot(x, y)  # Dot product of two vectors


    ### Setup two arrays / matrices
    a = np.array([[1, 2],
                  [3, 9]])
    b = np.array([[2, 4],
                  [5, 6]])


    ### Dot Product of two arrays
    print np.dot(a, b)


    ### Solving system of equations (i.e. 2 different equations with x and y)
    print LA.solve(a, b)


    ### Inverse of a matrix undoes the effects of the Matrix
    # The matrix multipled by the inverse matrix returns the 
    # 'identity matrix' (ones on the diagonal and zeroes everywhere else); 
    # identity matrix is useful for getting rid of the matrix in some equation
    print LA.inv(a)  # return inverse of the matrix
    print "\n"


    ### Determinant of a matrix is a special way to combine the entries of a
    # matrix that serves to check if matrix is invertible (!=0) or not (=0)
    print LA.det(a)  # returns the determinant of the array
    print "\n"  # e.g. 3, means that it is invertible


    ### Eigenvectors is a special set of input vectors for which the action of
    # the matrix is described as simple 'scaling'.  When a matrix is multiplied
    # by one of its eigenvectors, the output is the same eigenvector multipled
    # by a constant (that constant is the 'eigenvalue' of the matrix)
    print LA.eigvals(a)  # comput the eigenvalues of a general matrix
    print "\n"
    print LA.eigvalsh(a)  # Comput the eigenvalues of a Hermitian or real symmetric matrix
    print "\n"
    print LA.eig(a)  # return the eigenvalues for a square matrix
    print "\n"
    print LA.eigh(a)  # return the eigenvalues or eigenvectors of a Hermitian or symmetric matrix
    print "\n"

def ismassmatrix(M, semi=False):
    """ Check whether M is a valid mass matrix.

    :param M: the mass matrix to check
    :type  M: (6,6)-array
    :param bool semi: if set to ``True``, positive *semi*-definite matrices
                      are also considered valid.
    :return: ``True`` if M is correctly shaped and symmetric positive definite.

    """
    common = M.shape == (6, 6) and allclose(M, M.T) and allclose(M[3:6, 3:6], M[3, 3] * eye(3))
    if semi:
        return common and (eigvals(M) >= 0.0).all()
    else:
        return common and (eigvals(M) > 0.0).all()

def polyroots(cs):
    """Roots of a polynomial.

    Compute the roots of the Chebyshev series `cs`. The argument `cs` is a
    sequence of coefficients ordered from low to high. i.e., [1,2,3] is the
    polynomial ``1 + 2*x + 3*x**2``.

    Parameters
    ----------
    cs : array_like of shape(M,)
        1D array of polynomial coefficients ordered from low to high.

    Returns
    -------
    out : ndarray
        An array containing the complex roots of the polynomial series.

    Examples
    --------

    """
    # cs is a trimmed copy
    [cs] = pu.as_series([cs])
    if len(cs) <= 1 :
        return np.array([], dtype=cs.dtype)
    if len(cs) == 2 :
        return np.array([-cs[0]/cs[1]])
    n = len(cs) - 1
    cmat = np.zeros((n,n), dtype=cs.dtype)
    cmat.flat[n::n+1] = 1
    cmat[:,-1] -= cs[:-1]/cs[-1]
    roots = la.eigvals(cmat)
    roots.sort()
    return roots

    def __init__(self,ninput,nnodes,conn_input=0.4,conn_recurrent=0.2,gamma=numpy.tanh,frac_exc=0.5):
        self.ninput=ninput
        self.nnodes=nnodes
        self.gamma=gamma
        self.conn_recurrent=conn_recurrent
        self.conn_input=conn_input
        self.frac_exc=frac_exc

        w_echo = numpy.array(
            [[self.connection_weight(i,j)
              for j in range(self.nnodes)]
            for i in range(self.nnodes)])
        w_input=numpy.array(
            [[self.input_weight(i,j)
              for j in range(self.ninput)]
            for i in range(self.nnodes)])
        w_add = numpy.array(
            [self.add_weight(i)
             for i in range(self.nnodes)])
        
        eigenvalues=linalg.eigvals(w_echo)
        spectral_radius=max([abs (a) for a in eigenvalues])
        w_echo=(0.95/spectral_radius)*w_echo
        
        self.w_echo = w_echo
        self.w_input = w_input
        self.w_add = w_add

def ismassmatrix(M, semi=False):
    """Check whether M is a valid mass matrix.

    Return ``True`` if M is correctly shaped and symmetric positive
    definite.

    When ``semi`` is set to ``True``, positive *semi*-definite matrices
    are also considered valid.

    """
    common = M.shape == (6,6) and allclose(M, M.T) and \
            allclose(M[3:6, 3:6], M[3,3]*eye(3))
    if semi:
        return common and (eigvals(M) >= 0.).all()
    else:
        return common and (eigvals(M) > 0.).all()

def speigen_range(matrix, retry=True, coerce=True):
    """
    Construct the eigenrange of a potentially sparse matrix.
    """
    if spar.issparse(matrix):
        try:
            emax = spla.eigs(matrix, k=1, which='LR')[0]
        except (spla.ArpackNoConvergence, spla.ArpackError) as e:
            rowsums = np.unique(np.asarray(matrix.sum(axis=1)).flatten())
            if np.allclose(rowsums, np.ones_like(rowsums)):
                emax = np.array([1])
            else:
                Warn('Maximal eigenvalue computation failed to converge'
                     ' and matrix is not row-standardized.')
                raise e
        emin = spla.eigs(matrix, k=1, which='SR')[0]
        if coerce:
            emax = emax.real.astype(float)
            emin = emin.real.astype(float)
    else:
        try:
            eigs = nla.eigvals(matrix)
            emin, emax = eigs.min().astype(float), eigs.max().astype(float)
        except Exception as e:
            Warn('Dense eigenvector computation failed!')
            if retry:
                Warn('Retrying with sparse matrix...')
                spmatrix = spar.csc_matrix(matrix)
                speigen_range(spmatrix)
            else:
                Warn('Bailing...')
                raise e
    return emin, emax

def importance_pca(data, kpi, max_features=10):
    """
    :param data: dataframe containing training data
    :param kpi: Name of the current kpi
    :param max_features: maximum number of features to return
    :return: list of the best metrics

    The function does not use scikit-learn PCA implementation, because it is
    more difficult to associate each feature with its eigenvalue in the
    correlation matrix, also we are not interested in the projection vectors,
    but only in the eigenvalues.
    """
    columns = data[[col for col in set(data.columns) - {kpi}]].columns

    # correlation matrix and eigenvalues calculation
    corr = data[columns].corr().fillna(value=0).values
    eigenvalues = linalg.eigvals(corr)
    normalized_eigenvalues = normalize_series(eigenvalues)
    scaled_eigenvalues = normalized_eigenvalues / sum(normalized_eigenvalues)
    ranked_columns = [
        (eig, columns[n]) for n, eig in enumerate(scaled_eigenvalues)
    ]
    return [(j, i) for i, j in sorted(
        ranked_columns, reverse=True
    )][: max_features]

    def updateD_H(self, x):
        """
        Compute Hessian for update of D

        See [2] for derivation of Hessian
        """
        self.precompute(x)
        H = zeros((len(x), len(x)))
        Ai = zeros(self.A.shape[0])
        Aj = zeros(Ai.shape)
        for i in range(len(x)):
            Ai = self.A[:, i]
            ti = dot(self.AD, outer(self.R[:, i], Ai)) + dot(outer(Ai, self.R[i, :]), self.ADt)

            for j in range(i, len(x)):
                Aj = self.A[:, j]
                tj = outer(Ai, Aj)
                H[i, j] = (
                    self.E * (self.R[i, j] * tj + self.R[j, i] * tj.T) -
                    ti * (
                        dot(self.AD, outer(self.R[:, j], Aj)) +
                        dot(outer(Aj, self.R[j, :]), self.ADt)
                    )
                ).sum()
                H[j, i] = H[i, j]
        H *= -2
        e = eigvals(H).min()
        H = H + (eye(H.shape[0]) * e)
        return H

    def _hessx(xi, i):
        hess = numpy.zeros((k, k))
        for j in range(l):
            gij = snp_matrix[i, j]
            if gij == geosnp.MISSING:
                continue

            yj = Y[j]
            qj, aj, bj = yj[0], yj[1 : k + 1], yj[-1]
            qnf = (qj * sum(xi ** 2.0)) + aj.dot(xi) + bj
            fij = 1.0 / (1.0 + math.exp(qnf))
            term = 2.0 * sum(qj * xi) + aj
            hess -= fij * (1.0 - fij) * numpy.outer(term, term) + (gij - 2.0 * fij) * (2.0 * qj)

        # flip for NLL
        hess = -2.0 * hess

        # as NLL wrt X is not necessarily convex, we may need to alter the
        # hessian so that newton method still converges to local optimum
        # algorithm 6.3 from Numerical Opt Nocedal,Wright 1999
        beta = linalg.norm(hess, "fro")
        tau = 0 if min(numpy.diag(hess)) > 0 else beta
        eye = numpy.eye(k)
        while True:
            hess = hess + (tau * eye)
            # test for Positive Definiteness
            if min(linalg.eigvals(hess)) > 0:
                break
            else:
                tau = max(2 * tau, beta / 2)

        return hess

	def calculate_beta(self, column_name):
		# it doesn't make much sense to calculate beta for less than two days,
		# so return nan.
		algorithm_returns = self.prices[column_name]
		algorithm_returns = (algorithm_returns-algorithm_returns.shift(1))/algorithm_returns.shift(1)
		if (column_name==self.benchmark_code_str):
			self.prices_return[column_name]=self.prices_return[column_name]/100
		else:
			self.prices_return[column_name]=algorithm_returns
		algorithm_returns = algorithm_returns[1:]
		benchmark_returns = self.prices_return[self.benchmark_code_str]
		benchmark_returns = benchmark_returns[1:]
		if len(algorithm_returns) < 2:
			return np.nan

		returns_matrix = np.vstack([algorithm_returns,benchmark_returns])
		C = np.cov(returns_matrix, ddof=1)

		# If there are missing benchmark values, then we can't calculate the
		# beta.
		if not np.isfinite(C).all():
			return np.nan

		eigen_values = la.eigvals(C)
		condition_number = max(eigen_values) / min(eigen_values)
		algorithm_covariance = C[0][1]
		benchmark_variance = C[1][1]
		beta = algorithm_covariance / benchmark_variance
		
		#print beta
		return beta

	def A_dpf(self, theta, phi, A=None, no_psd=False, byhand=False):
		"""computes A in the dominant polarization frame. If A is supplied, it converts A to the dominant polarization frame"""
		if A==None:
			A = self.A(theta, phi, 0.0, no_psd=no_psd)

		if no_psd:
			npix, npol, npol = np.shape(A)
			a = np.empty((npix, 1, npol, npol), float)
			a[:,0,:,:] = A
			A = a

		if byhand:
			a = A[:,:,0,0]
			b = A[:,:,0,1]
			c = A[:,:,1,1]
			dpfA = np.zeros_like(A, float)
			x = ((a-c)**2 + 4*b**2)**0.5
			y = a+c
			dpfA[:,:,0,0] = 0.5*(y + x)
			dpfA[:,:,1,1] = 0.5*(y - x)
		else:
			dpfA=np.zeros_like(A, float)
			vals = linalg.eigvals(A)[:,:,::-1] ### order by decreasing eigenvalue
			for i in xrange(self.Np):
				dpfA[:,:,i,i] = vals[:,:,i]

		if no_psd:
			return dpfA[:,0,:,:]
		else:
			return dpfA

def is_chaotic(J):
    """Determines whether the Jacobian matrix is chaotic
    Arguments
    ==========
    J = string, list of lists, array of arrays, or matrix

    Returns
    =======
    bool = True if chaotic, False otherwise
    """
    # start your code here
    margin = 0.0001
    J = np.mat(J)
    eigens = linalg.eigvals(J)
    less = 0
    equal = 0
    more = 0
    for eigen in eigens:
        if abs(eigen - 1) < margin:
            equal = equal + 1
        else:
            if eigen < 1:
                less = less + 1
            if eigen > 1:
                more = more + 1
    if less == 1 and equal == 1 and more == 1:
        return True
    else:
        return False

 def _compressibility(self, A, TdAdT, dAdn, B, dBdn):
     """Calculate the compressibility from the parameters A and B"""
     c = np.array([-1, A - B * (1 + B), -(A * B)], dtype=np.float)
     self._companionmatrix[0, :] = -c
     result = la.eigvals(self._companionmatrix)
     result = result[result.imag < ERROR_TOL].real
     return np.array((result.max(), result.min()), dtype=phase_struct)

def bisect_gFB(s, tres, Q11, Q22, Q12, Q21, k1, k2):
    """
    Find number of eigenvalues of H(s) that are equal to or less than s.

    Parameters
    ----------
    s : float
        Laplace transform argument.
    tres : float
        Time resolution (dead time).
    Q11 : array_like, shape (k1, k1)
    Q22 : array_like, shape (k2, k2)
    Q21 : array_like, shape (k2, k1)
    Q12 : array_like, shape (k1, k2)
        Q11, Q12, Q22, Q21 - submatrices of Q.
    k1 : int
        A number of open/shut states in kinetic scheme.
    k2 : int
        A number of shut/open states in kinetic scheme.

    Returns
    -------
    ng : int
    """

    h = qml.H(s, tres, Q11, Q22, Q12, Q21, k2)
    eigval = nplin.eigvals(h)
    ng = (eigval <= s).sum()
    return ng

    def calculate_beta(self):
        """

        .. math::

            \\beta_a = \\frac{\mathrm{Cov}(r_a,r_p)}{\mathrm{Var}(r_p)}

        http://en.wikipedia.org/wiki/Beta_(finance)
        """
        #it doesn't make much sense to calculate beta for less than two days,
        #so return none.
        if len(self.algorithm_returns) < 2:
            return 0.0, 0.0, 0.0, 0.0, []

        returns_matrix = np.vstack([self.algorithm_returns,
                                    self.benchmark_returns])
        C = np.cov(returns_matrix)
        eigen_values = la.eigvals(C)
        condition_number = max(eigen_values) / min(eigen_values)
        algorithm_covariance = C[0][1]
        benchmark_variance = C[1][1]
        beta = C[0][1] / C[1][1]

        return (
            beta,
            algorithm_covariance,
            benchmark_variance,
            condition_number,
            eigen_values
        )

def too_small_check(ccov, eigcut = 10**(-10)):
    """Check to see if the matrix eigenvalues are too small.
    This can cause problems when computing chi^2 due to precision loss
    """
    testeig = eigvals(ccov)
    flag = 0
    for entry in testeig:
        if entry < eigcut:
            flag = 1
            print "***Warning***"
            print "Range selected has a covariance matrix with"
            print "very small eigenvalues.  This can cause problems"
            print "in computing chi^2, as well as quantities derived"
            print "from chi^2. The cuttoff is set at:", eigcut
            print "Problematic eigenvalue = ", entry
            break
    if flag == 1:
        print "List of eigenvalues of covariance matrix:"
        for entry in testeig:
            print entry
        while True:
            print "Continue? (y/n)"
            cresp = str(raw_input())
            if (cresp == "n" or cresp == "no"
                    or cresp == "No" or cresp == "N"):
                sys.exit(0)
            if (cresp == "y" or cresp == "yes"
                    or cresp == "Yes" or cresp == "Y"):
                break
            else:
                print "Sorry, I didn't understand that."
                continue
    return 0

    def point_is_stable(self, point, tol=1e-5):
        """
        returns true if a given steady state is stable
        """
        
        if not self.point_is_steady_state(point, tol):
            raise ValueError('Supplied point is not a steady state')
        
        jacobian = self.jacobian(point, tol)
        x_check, y_check = False, False #< checked direction
        
        # check special boundary cases
        if 'left' in self.region_constraint and np.abs(point[0] - self.region[0]) < 1e-8:
            x_check = True 
        elif 'right' in self.region_constraint and np.abs(point[0] - self.region[2]) < 1e-8:
            x_check = True 

        if 'bottom' in self.region_constraint and np.abs(point[1] - self.region[1]) < 1e-8:
            y_check = True
        elif 'top' in self.region_constraint and np.abs(point[1] - self.region[3]) < 1e-8:
            y_check = True
        
        # check the remaining directions
        if x_check and y_check:
            return True # both x and y direction are stable
        elif x_check:
            return jacobian[1, 1] < 0 # y direction has to be tested
        elif y_check:
            return jacobian[0, 0] < 0 # x direction has to be tested
        else:
            return all(eigvals(jacobian) < 0) # both directions have to be tested

def lanczos_eig(A):
    m,n = np.shape(A)
    Q = np.zeros((m,n))
    beta = 0
    x0 = np.random.randn(n)
    Q[:,0] = x0/(la.norm(x0))
    alpha = np.zeros(n)
    gama = []
    U = np.zeros((m,n))
    H = np.zeros((m,n))
    
    for i in range(0,n):
        U[:,i] = np.dot(A,Q[:,i])
        alpha = np.dot(Q[:,i].T,U[:,i])
        U[:,i] = U[:,i]-beta*Q[:,i-1]-alpha*Q[:,i]
        beta = la.norm(U[:,i])
        H[i,i] = alpha
        if(i+1<n):
          	H[i+1,i] = beta
          	H[i,i+1] = beta
        if beta == 0:
          	break
        if(i+1<n):
          	Q[:,i+1]=U[:,i]/beta
        gama.append(np.copy(la.eigvals(H).T))

    return gama

	def Aii_dpf(self, i, theta, phi, A=None, no_psd=False, byhand=False):
		"""computes a single component of A in the dominant polarization frame. If A is supplied, it converts to the dominant polarizatoin frame"""
		if A==None:
			A = self.A(theta, phi, 0.0, no_psd=no_psd)

		if no_psd:
			npix, npol, npol = np.shape(A)
			a = np.empty((npix,1,npol,npol),float)
			a[:,0,:,:] = A
			A = a
		if byhand:
			a = A[:,:,0,0]
                        b = A[:,:,0,1]
                        c = A[:,:,1,1]
                        dpfA = np.zeros_like(A, float)
                        x = ((a-c)**2 + 4*b**2)**0.5
                        y = a+c
			if no_psd:
				return 0.5*(y[:,0] + (-1)**i * x[:,0])
			else:
				return 0.5*(y + (-1)**i * x)
		else:
			vals = linalg.eigvals(A)[:,:,::-1] ### order by decreasing eigenvalue
			if no_psd:
				return vals[:,0,i]
			else:
				return vals[:,:,i]