def jointpmf2(pe, pb, n):
    # pWL(w, l) = Pr{W = w, L = l}
    #           = Pr{W = w | L = l} * Pr{L = l}
    #           = pW(w - 2; l - 2) * pL(l; n)

    T = mmbittm(pe, pb)
    Pss = mmss(T)

    pW = errpmf(pe, pb, n)

    pWL = np.zeros((n + 1, n + 1))
    pWL[0, 0] = pW[0]
    pWL[1, 1] = pW[1]

    for l in range(2, n + 1):
        Tw = mmstep(T, l - 2)
        #print('Tw =\n', Tw)
        for x in range(n - l + 1):
            y = n - l - x
            #print('0' * x + '1' + 'x' * (l - 2) + '1' + '0' * y)
            Tx = LA.matrix_power(T[0], x)
            Ty = LA.matrix_power(T[0], y)
            #Twl = Tx.dot(T[1]).dot(Tw).dot(T[1]).dot(Ty)
            #print('Tx =\n', Tx)
            #print('Ty =\n', Ty)
            #print('Twl =\n', Twl)
            #Pwl = Pss.dot(Twl).sum(axis=1)
            for w in range(2, l + 1):
                Twl = Tx.dot(T[1]).dot(Tw[w - 2]).dot(T[1]).dot(Ty)
                Pwl = Pss.dot(Twl).sum()
                #print('Twl =\n', Twl)
                #print('Pwl =\n', Pwl)
                pWL[w, l] += Pwl

    return pWL

def prob_1():
	#problem 1
	#part 1
	A = np.array([[.75, .5], [.25, .5]])
	print A.dot(A)[0,0]
	#part 2
	print la.matrix_power(A, 20)[0,0]

    def make_propagators(pb=0.0, kex=0.0, dw=0.0, r_hxy=5.0, dr_hxy=0.0,
                         r_cz=1.5, r_2hzcz=0.0, etaxy=0.0, etaz=0.0,
                         j_hc=0.0, cs_offset=0.0):

        w1 = 2.0 * pi / (4.0 * pw)

        l_free, l_w1x, l_w1y = compute_liouvillians(
            pb=pb,
            kex=kex,
            dw=dw,
            r_hxy=r_hxy,
            dr_hxy=dr_hxy,
            r_cz=r_cz,
            r_2hzcz=r_2hzcz,
            etaxy=etaxy,
            etaz=etaz,
            j_hc=j_hc,
            cs_offset=cs_offset,
            w1=w1
        )

        p_neg = expm(l_free * -2.0 * pw / pi)
        p_90px = expm((l_free + l_w1x) * pw)
        p_90py = expm((l_free + l_w1y) * pw)
        p_90mx = expm((l_free - l_w1x) * pw)
        p_180py = matrix_power(p_90py, 2)
        p_180pmx = 0.5 * (matrix_power(p_90px, 2) +
                          matrix_power(p_90mx, 2))

        ps = p_neg, p_90px, p_90py, p_90mx, p_180py, p_180pmx

        return l_free, ps

 def test_regress_mutation_probability_endpoint_conditioning(self):
     ngenerations = 10
     selection = 2.0
     mutation = 0.0001
     recombination = 0.001
     nchromosomes = 2
     npositions = 4
     K_initial = np.array([
         [1,1,1,1],
         [0,0,0,0]], dtype=np.int8)
     K_final = np.array([
         [1,1,0,0],
         [0,0,1,1]], dtype=np.int8)
     #
     ngenboundaries = ngenerations - 1
     no_mutation_prior = (1 - mutation)**(
             npositions*ngenboundaries*nchromosomes)
     #
     initial_long = popgenmarkov.bin2d_to_int(K_initial)
     final_long = popgenmarkov.bin2d_to_int(K_final)
     ci_to_short, short_to_count, sorted_chrom_lists = get_state_space_info(
             nchromosomes, npositions)
     initial_ci = chroms_to_index(
             sorted(popgenmarkov.bin_to_int(row) for row in K_initial),
             npositions)
     initial_short = ci_to_short[initial_ci]
     final_ci = chroms_to_index(
             sorted(popgenmarkov.bin_to_int(row) for row in K_final),
             npositions)
     final_short = ci_to_short[final_ci]
     # get an answer using the less efficient methods
     P_sr = popgenmarkov.get_selection_recombination_transition_matrix(
             selection, recombination, nchromosomes, npositions)
     P_m = popgenmarkov.get_mutation_transition_matrix(
             mutation, nchromosomes, npositions)
     p_b_given_a = linalg.matrix_power(np.dot(P_sr, P_m), ngenerations-1)[
             initial_long, final_long]
     p_b_given_a_no_mutation = linalg.matrix_power(P_sr, ngenerations-1)[
             initial_long, final_long]
     no_mutation_posterior = (
             no_mutation_prior * p_b_given_a_no_mutation) / p_b_given_a
     # get an answer using the more efficient methods
     P_sr_s = get_selection_recombination_transition_matrix_s(
             ci_to_short, short_to_count, sorted_chrom_lists,
             selection, recombination, nchromosomes, npositions)
     P_m_s = get_mutation_transition_matrix_s(
             ci_to_short, short_to_count, sorted_chrom_lists,
             mutation, nchromosomes, npositions)
     p_b_given_a_s = linalg.matrix_power(
             np.dot(P_sr_s, P_m_s), ngenerations-1)[
                     initial_short, final_short]
     p_b_given_a_no_mutation_s = linalg.matrix_power(
             P_sr_s, ngenerations-1)[
                     initial_short, final_short]
     no_mutation_posterior_s = (
             no_mutation_prior * p_b_given_a_no_mutation_s) / p_b_given_a_s
     #
     self.assertTrue(
             np.allclose(no_mutation_posterior, no_mutation_posterior_s))

 def test_matrix_power_operator(self):
     random.seed(1234)
     n = 5
     k = 2
     p = 3
     nsamples = 10
     for i in range(nsamples):
         A = np.random.randn(n, n)
         B = np.random.randn(n, k)
         op = MatrixPowerOperator(A, p)
         assert_allclose(op.matmat(B), matrix_power(A, p).dot(B))
         assert_allclose(op.T.matmat(B), matrix_power(A, p).T.dot(B))

    def __call__(self, options, pars):
        """Evaluate the model for a given set of parameters."""
        self.max_T = int(pars.get('max_T', 100))    # range of sample sizes
        T = np.arange(1., self.max_T + 1)
        N = map(int, np.floor(T))
        self.set_statespace(pars)                   # threshold and state space

        # which option is sampled first? [0, 1, None]
        firstoption = pars.get('firstoption', None)

        # evaluate the starting distribution
        Z = self.Z(self.m, {'theta': self.theta + 1, 'tau': pars.get('tau', .5)})

        # repeat Z for both options, and re-normalize
        if firstoption is None:
            Z = np.concatenate((Z, Z), axis=1)
        elif firstoption is 0:
            Z = np.concatenate((Z, np.matrix(np.zeros(self.m - 2))), axis=1)
        elif firstoption is 1:
            Z = np.concatenate((np.matrix(np.zeros(self.m - 2)), Z), axis=1)
        Z = Z / float(np.sum(Z))


        # transition matrix
        tm = self.transition_matrix_reflecting(options, pars)

        # min-steps
        if 'minsamplesize' in pars:
            min_ss = pars.get('minsamplesize')
            if min_ss == 2:
                Z = Z * tm
            else:
                Z = Z * matrix_power(tm, min_ss - 1)
            assert np.isclose(Z.sum(), 1.)

        # evaluate evolution in states
        M = [Z * matrix_power(tm, 0)]
        for n in N[1:]:
            M.append( np.dot(M[-1], tm) )
        p_state_t = np.array(M).reshape((len(N), self.m * 2))

        p_0 = p_state_t[:,self.theta] + p_state_t[:,self.m + self.theta]
        p_L = p_state_t[:,:self.theta].sum(axis=1) + p_state_t[:,self.m:(self.m+self.theta)].sum(axis=1) + p_0 * 0.5
        p_H = p_state_t[:,(1+self.theta):self.m].sum(axis=1) + p_state_t[:,(self.m+self.theta+1):].sum(axis=1) + p_0 * 0.5
        p_LH = np.transpose((p_L, p_H))

        assert np.isclose(p_LH.sum(), p_LH.shape[0])

        return {'T': T,
                'states_t': p_state_t,
                'p_resp_t': p_LH}

def randmatstat(t):
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i in range(t):
        a = randn(n, n)
        b = randn(n, n)
        c = randn(n, n)
        d = randn(n, n)
        P = concatenate((a, b, c, d), axis=1)
        Q = concatenate((concatenate((a, b), axis=1), concatenate((c, d), axis=1)), axis=0)
        v[i] = trace(matrix_power(dot(P.T,P), 4))
        w[i] = trace(matrix_power(dot(Q.T,Q), 4))
    return (std(v)/mean(v), std(w)/mean(w))

def DiagonalNewton(d0,d1,Q,tolerance):

    if tolerance <= 0:
            print "tolerance <= 0, algorithm won't terminate"
            return None

    #verify if Q is positive semidefinite, correct dimensions
    try:
        y = linalg.cholesky(Q)
    except linalg.LinAlgError:
        print "Matrix is either not positive semidefinite, or nonsquare"
        return None

    e = np.ones(d1.size)
    w = linalg.norm(d1 -d0)
    D = np.diag(d1)
    w = linalg.norm(np.dot(np.dot(np.dot(D,Q),D),e)-e)
    iterationcounter = 0
    while(w > tolerance):
        d0 = d1
        D = np.diag(d1)


        try:
       
            y = np.dot(linalg.matrix_power((linalg.matrix_power(D,-2)+Q),-1), np.diag( linalg.matrix_power(D,-1 ))-np.dot(Q,d1))

        except:

            return None
      
          
       # print  np.dot(np.dot(np.transpose(y) ,linalg.matrix_power(D,-2)),y)
        
        lamb = np.dot(
                    np.dot(np.transpose(y) ,linalg.matrix_power(D,-2)),
                            y) + np.dot(np.dot(np.transpose(y),Q),y) #computing lambda
        #lamb < 1, a = 1 lamb >=1 1, a = 1/2
        
        #print lamb
        if lamb > 1:
            d1 = d0 + float(1)/(1 + lamb**(0.5))*y

        else:
            d1 = d0 + y
        w = linalg.norm(np.dot(np.dot(np.dot(D,Q),D),e)-e)

        iterationcounter += 1

    return (w, d1,iterationcounter)

def directed_weighted_clustering(g,weightString):
	
	
	n = g.number_of_nodes()
	from numpy import linalg as LA
	#adjacency matrix
	A = nx.to_numpy_matrix(g,nodelist=g.nodes(),weight=None)
	A2 = LA.matrix_power(A,2)
	AT = A.T
	Asum = AT + A
	cVector = [i for i in range(n)]
	cVector = np.asmatrix(cVector)
	
	kin = {i:np.dot(AT[i],cVector.T) for i in range(n)}
	kout = {i:np.dot(A[i],cVector.T)for i in range(n)}
	kparallel = {i:np.dot(Asum[i],cVector.T)for i in range(n)}

	#print "kin"
	#print kin
	#weight matrix
	W = nx.to_numpy_matrix(g,nodelist=g.nodes(),weight=weightString)
	WT = W.T
	W2 = LA.matrix_power(W,2)
	W3 = LA.matrix_power(W,3)
			
	WWTW =  W*WT*W
	WTW2 = WT*W2
	W2WT = W2*WT
	
	ccycle = {i:0 for i in range(n)}
	cmiddle = {i:0 for i in range(n)}
	cin = {i:0 for i in range(n)}
	cout = {i:0 for i in range(n)}

	for i in range(n):
			
			if kin[i]*kout[i]  - kparallel[i] > 0:
				ccycle[i] = W3[i,i] / float((kin[i]*kout[i] - kparallel[i]))
				cmiddle[i] = WWTW[i,i] / float((kin[i]*kout[i] - kparallel[i]))
			if kin[i] > 1: 
				cin[i] = WTW2[i,i] / float((kin[i]*(kin[i]-1)))
			if kout[i] > 1: 
				cout[i] = W2WT[i,i] / float((kout[i]*(kout[i]-1))) 
	#print type((np.mean(ccycle.values()),np.mean(cmiddle.values()),np.mean(cin.values()),np.mean(cout.values())))
	#print "here"
	#input()
	#return (np.mean(ccycle.values()),np.mean(cmiddle.values()),np.mean(cin.values()),np.mean(cout.values()))
	return (ccycle,cmiddle,cin,cout)

def power_method(file_name, error_tol, u0):

    # Creates a matrix from the .dat file
    A = np.genfromtxt(file_name, delimiter=" ")

    m = int(A.shape[0]) #rows of A

    u = np.asarray(np.asmatrix(u0))
    uRows = int(u.shape[0]) #rows of u
    uCols = int(u.shape[1]) #columns of u

    # Sets the tolerance
    tol = error_tol

    # Sets the initial number of iterations
    iteration = 0

    # Sets an array with one entry 0 to use for the error in the while loop
    eigenvalues = [0]

    # While number of iterations is less than 100, the matrix
    # A is raised to a power equal to the iteration and multiplied
    # by the original u0 that was given as an input. The dominant
    # eigenvector is found and from that, the dominant eigenvalue
    # is found.
    while iteration < 100:
        copyA = LA.matrix_power(A, iteration+1)
        x = np.zeros(shape=(m, uCols))
        for i in range(m):
            for j in range(uCols):
                for k in range(uRows):
                    x[i][j] += (copyA[i][k] * u[k][j]) #Multiplies matrix A and u
        eigenvector = x / LA.norm(x) # Finds the dominant eigenvector
        eigenRows = int(eigenvector.shape[0]) #rows of dominant eigenvector
        eigenCols = int(eigenvector.shape[1]) #columns of eigenvector

        Ax = np.zeros(shape = (m, eigenCols))

        for i in range(m):
            for j in range(eigenCols):
                for k in range(eigenRows):
                    Ax[i][j] += A[i][k] * eigenvector[k][j]
        Axx = np.vdot(Ax, eigenvector)
        xx = np.vdot(eigenvector, eigenvector)
        eigenvalue = Axx / xx # Finds the dominant eigenvalue

        eigenvalues.append(eigenvalue)

        if (np.absolute(eigenvalues[iteration+1] - eigenvalues[iteration])) <= tol:
            break

        iteration += 1

    print "Dominant eigenvalue = ", eigenvalue
    print "Dominant eigenvector =\n", eigenvector

    if iteration >= 100:
        print "Did not converge after 100 iterations."
    else:
        print "Took " + str(iteration) + " iterations to converge."

def prob_3:
	#problem 3
	A = np.array([[0, 0, 1, 0, 1, 0, 1],
				  [1, 0, 0, 0, 0, 1, 0],
				  [0, 0, 0, 0, 0, 1, 0],
				  [1, 0, 0, 0, 1, 0, 0],
				  [0, 0, 0, 1, 0, 0, 0],
				  [0, 0, 1, 0, 0, 0, 1],
				  [0, 1, 0, 0, 0, 0, 0]], dtype=np.int64)
	A5 = la.matrix_power(A,5)
	coords = np.where(A5==np.max(A5))
	#note: indexing from 0
	print "maximum of 5 step connections at: ", zip(coords[0], coords[1])
	A7 = la.matrix_power(A,7)
	coords = np.where(A7==0)
	print "no 7 step connection for: ", zip(coords[0], coords[1])

def pam_matrix(d, scale=float(np.log(2)/2), as_ints=False):
	"""Creates PAM scoring matrix.

	Values calculated from (natural) log-odds ratio of PAM{d}:PAM{inf}. The
	output will be this matrix *divided* by the scale parameter (this seems
	to be the convention, but I am not sure why).

	Args:
		d: int. Mutation distance (number of time steps in the Markov chain
			model).
		scale: float. Units of returned matrix, relative to "natural" units
			of the log-odds ration. Returned matrix will be the log-odds
			values *divided* by the scale. Defaults to ln(2)/2 (half-bit
			units), because that's what I've seen everywhere.
		as_ints: bool. If true, entries will be rounded to the nearest integer.
	"""
	# Calculate matrix
	dayhoff_n = matrix_power(dayhoff_matrix, d)
	matrix = np.log(dayhoff_n / dayhoff_stationary[:, None]) / scale

	# Doesn't seem to be completely symmetrical, hopefully just due to
	# floating-point errors. Fudge it a bit.
	matrix += matrix.transpose()
	matrix *= .5

	# Round if requested
	if as_ints:
		matrix = np.round(matrix)

	return SubstitutionMatrix(aa_symbols, matrix.astype(np.float32))

def simulate_evolution(dna_sequence, alpha, time, dt, threads=1):
    tree = Tree([PhylogeneticNode(dna_sequence)], [])
    transition_matrix = build_transition_matrix(alpha)
    power_matrix = LA.matrix_power(transition_matrix, dt)

    while time > 0:
        i = 0
        fixed_length = len(tree.get_nodes())
        node_list = tree.get_nodes()
        if fixed_length > 10:
            return tree

        while i < fixed_length:
            node = node_list[i]
            origin_sequence = node.get_sequence()
            wrapper = list(origin_sequence)
            next_node = mutate_module(power_matrix, wrapper, threads=threads)

            if next_node != -1:
                tree.add_node(next_node)
                tree.add_edge(node, next_node)

            gc.collect()
            i += 1
        time -= dt
        print time

    return tree

 def spectralAnalysis(self):
     # Assumes adjacency matrix is already made
     for i in range(2,self.numNodes+1):
         m = linalg.matrix_power(self.adjMatrix,i)/i
     print "Spectral analysis, k= " + str(i) + ":\n" + str(m)
     print ""
     return m.diagonal()

    def calc_strided(self, data, power):
        """ The strided way to do the calculation """

        # Extract shapes from data.
        NE, NS, NM, NO, ND, Row, Col = data.shape

        # Make array to store results
        calc = zeros([NE, NS, NM, NO, ND, Row, Col])

        # Get the data view, from the helper function.
        data_view = self.stride_help_square_matrix(data)

        # Get the power view, from the helpwer function.
        power_view = self.stride_help_element(power)

        # The index view could be pre formed in init.
        index = self.index
        index_view = self.stride_help_array(index)

        # Zip them together and iterate over them.
        for data_i, power_i, index_i in zip(data_view, power_view, index_view):
            # Do power calculation with numpy method.
            data_power_i = matrix_power(data_i, int(power_i))

            # Extract index from index_view.
            ei, si, mi, oi, di = index_i

            # Store the result.
            calc[ei, si, mi, oi, di] = data_power_i

        return calc

def basis(H, v):

	basisset=[]
	for l in range(M):
		basisset.extend(np.dot(LA.matrix_power(H,l),v))
	basisset=np.reshape(basisset,(M,2*N))
	return (basisset)

    def calc_normal(self, data, power):
        """ The normal way to do the calculation """

        # Extract shapes from data.
        NE, NS, NM, NO, ND, Row, Col = data.shape

        # Make array to store results
        calc = zeros([NE, NS, NM, NO, ND, Row, Col])

        # Normal way, by loop of loops.
        for ei in range(NE):
            for si in range(NS):
                for mi in range(NM):
                    for oi in range(NO):
                        for di in range(ND):
                            # Get the outer square data matrix,
                            data_i = data[ei, si, mi, oi, di]

                            # Get the power.
                            power_i = power[ei, si, mi, oi, di]

                            # Do power calculation with numpy method.
                            data_power_i = matrix_power(data_i, power_i)

                            # Store result.
                            calc[ei, si, mi, oi, di] = data_power_i

        return calc

    def _calculate_unscaled_profile(self, params_local, **kwargs):
        """TODO: class docstring."""

        self.liouv.update(params_local)

        # Calculation of the propagators corresponding to all the delays
        delays = dict(zip(self.delays, self.liouv.delays(self.delays)))
        d_zeta = delays[self.t_zeta]

        # Calculation of the propagators corresponding to all the pulses
        p180_sx = self.liouv.perfect180["sx"]
        p180_ix = self.liouv.perfect180["ix"]
        p180_iy = self.liouv.perfect180["iy"]

        # Calculate starting magnetization vector
        mag0 = self.liouv.compute_mag_eq(params_local, term="2iysx")

        if self.small_protein:
            mag0 = d_zeta @ p180_sx @ p180_ix @ d_zeta @ mag0

        # Calculating the cpmg trains
        cp = {0: self.liouv.identity}

        for ncyc in set(self.data["ncycs"][~self.reference]):
            tau_cp = delays[self.tau_cps[ncyc]]
            echo = tau_cp @ p180_iy @ tau_cp
            cp_train = la.matrix_power(echo, int(ncyc))
            cp[ncyc] = cp_train @ p180_sx @ cp_train

        profile = [
            self.liouv.collapse(self.detect @ cp[ncyc] @ mag0)
            for ncyc in self.data["ncycs"]
        ]

        return np.asarray(profile)

def build_rand_sparse_diag_mat_and_multiply():
    array = diags(random.random_sample(25).tolist(), 0).toarray()
    #print array
    heat_map_data1 = [go.Heatmap(z=np.flipud(array).tolist())]
    heat_map_data2 = [go.Heatmap(z=np.flipud((LA.matrix_power(array, 2))).tolist())] #multiply with itself.
    plot_url = plotly.offline.plot(heat_map_data1, filename='basic-heatmap1.html')
    plot_url = plotly.offline.plot(heat_map_data2, filename='basic-heatmap2.html')

    def _calc_observable(pb=0.0, kex=0.0, dw=0.0, r_ixy=5.0, dr_ixy=0.0,
                         ncyc=0):
        """
        Calculate the intensity in presence of exchange during a cpmg-type pulse
        train.

        Parameters
        ----------
        i0 : float
            Initial intensity.
        pb : float
            Fractional population of state B,
            0.0 for 0%, 1.0 for 100%
        kex : float
            Exchange rate between state A and B in /s.
        dw : float
            Chemical shift difference between states A and B in rad/s.
        r_ixy : float
            Transverse relaxation rate of state a in /s.
        dr_ixy : float
            Transverse relaxation rate difference between states a and b in /s.

        Returns
        -------
        out : float
            Intensity after the CPMG block

        """

        dw *= ppm_to_rads

        mag_eq = compute_iy_eq(pb)

        if ncyc == 0:

            mag = mag_eq

        else:

            l_free = compute_liouvillians(
                pb=pb,
                kex=kex,
                dw=dw,
                r_ixy=r_ixy,
                dr_ixy=dr_ixy
            )

            t_cp = time_t2 / (4.0 * ncyc)
            p_free = expm(l_free * t_cp)

            mag = matrix_power(
                p_free
                .dot(P_180Y)
                .dot(p_free),
                2 * ncyc
            ).dot(mag_eq)

        magy_a, _ = get_iy(mag)

        return magy_a