def align_magnetism(m, vectors):
    """ Rotates a matrix, to align its components with the direction
  of the magnetism """
    if not len(m) == 2 * len(vectors):  # stop if they don't have
        # compatible dimensions
        raise
    # pauli matrices
    from scipy.sparse import csc_matrix, bmat

    sx = csc_matrix([[0.0, 1.0], [1.0, 0.0]])
    sy = csc_matrix([[0.0, -1j], [1j, 0.0]])
    sz = csc_matrix([[1.0, 0.0], [0.0, -1.0]])
    n = len(m) / 2  # number of sites
    R = [[None for i in range(n)] for j in range(n)]  # rotation matrix
    from scipy.linalg import expm  # exponenciate matrix

    for (i, v) in zip(range(n), vectors):  # loop over sites
        vv = np.sqrt(v.dot(v))  # norm of v
        if vv > 0.000001:  # if nonzero scale
            u = v / vv
        else:  # if zero put to zero
            u = np.array([0.0, 0.0, 0.0])
        #    rot = u[0]*sx + u[1]*sy + u[2]*sz
        uxy = np.sqrt(u[0] ** 2 + u[1] ** 2)  # component in xy plane
        phi = np.arctan2(u[1], u[0])
        theta = np.arctan2(uxy, u[2])
        r1 = phi * sz / 2.0  # rotate along z
        r2 = theta * sy / 2.0  # rotate along y
        # a factor 2 is taken out due to 1/2 of S
        rot = expm(1j * r2) * expm(1j * r1)
        R[i][i] = rot  # save term
    R = bmat(R)  # convert to full sparse matrix
    mout = R * csc_matrix(m) * R.H  # rotate matrix
    return mout.todense()  # return dense matrix

def make_random_su(N):
    
    A = spla.expm(2 * np.pi * 1j * npr.random((N, N)))
    B = (A - np.trace(A) * np.identity(N) / N)
    C = 0.5 * (B - np.conj(B.T))
    
    return spla.expm(C)

    def _solve_numericaly(self, u, x0, t, t_list, t0, dps=2):
        """ returns the numeric evaluation of the system for input u, know state x0 at time t0 and times t_list
        """
        result = []
        for t_i in t_list:
            # we use the arbitrary precision module mpmath for numercial evaluation of the matrix exponentials
            first = (
                np.array(np.array(self.represent[2]), np.float)
                .dot(expm(np.array(np.array((self.represent[0] * (t_i - t0)).evalf()), np.float)))
                .dot(np.array(np.array(x0), np.float))
            )

            second = np.array(np.array((self.represent[3] * u.subs(t, t_i)).evalf()), np.float)

            integrand = (
                lambda tau: np.array(np.array(self.represent[2]), np.float)
                .dot(expm(np.array(np.array((self.represent[0] * (t_i - tau)).evalf()), np.float)))
                .dot(np.array(np.array(self.represent[1]), np.float))
                .dot(np.array(np.array(u.subs(t, tau).evalf()), np.float))
            )

            # the result must have the same shape as D:
            integral = zeros(self.represent[2].rows, 1)

            # Loop through every entry and evaluate the integral using mpmath.quad()
            for row_idx in xrange(self.represent[2].rows):

                integral[row_idx, 0] = quad(lambda x: integrand(x)[row_idx, 0], t0, t_i)[0]

            result.append(Matrix(first) + Matrix(second) + integral)

        # return sum of results
        return result

    def make_propagators(pb=0.0, kex=0.0, dw=0.0, r_coxy=5.0, dr_coxy=0.0,
                         r_nz=1.5, r_2coznz=0.0, etaxy=0.0, etaz=0.0,
                         j_nco=0.0, dj_nco=0.0, cs_offset=0.0):

        w1 = 2.0 * pi / (4.0 * pwco90)
        l_free, l_w1x, l_w1y = compute_liouvillians(pb=pb, kex=kex, dw=dw,
                                                    r_coxy=r_coxy, dr_coxy=dr_coxy,
                                                    r_nz=r_nz, r_2coznz=r_2coznz,
                                                    etaxy=etaxy, etaz=etaz,
                                                    j_nco=j_nco, dj_nco=dj_nco,
                                                    cs_offset=cs_offset, w1=w1)

        p_equil = expm(l_free * time_equil)
        p_taucc = expm(l_free * taucc)
        p_neg = expm(l_free * -2.0 * pwco90 / pi)
        p_90py = expm((l_free + l_w1y) * pwco90)
        p_90my = expm((l_free - l_w1y) * pwco90)
        p_180px = P180X  # Perfect 180 for CPMG blocks
        p_180py = matrix_power(p_90py, 2)
        p_180my = matrix_power(p_90py, 2)

        ps = (p_equil, p_taucc, p_neg, p_90py, p_90my,
              p_180px, p_180py, p_180my)

        return l_free, ps

    def _learnStep(self):
        """ Main part of the algorithm. """
        I = eye(self.numParameters)
        self._produceSamples()
        utilities = self.shapingFunction(self._currentEvaluations)
        utilities /= sum(utilities)  # make the utilities sum to 1
        if self.uniformBaseline:
            utilities -= 1./self.batchSize
        samples = array(map(self._base2sample, self._population))

        dCenter = dot(samples.T, utilities)
        covGradient = dot(array([outer(s,s) - I for s in samples]).T, utilities)
        covTrace = trace(covGradient)
        covGradient -= covTrace/self.numParameters * I
        dA = 0.5 * (self.scaleLearningRate * covTrace/self.numParameters * I
                    +self.covLearningRate * covGradient)

        self._lastLogDetA = self._logDetA
        self._lastInvA = self._invA

        self._center += self.centerLearningRate * dot(self._A, dCenter)
        self._A = dot(self._A, expm(dA))
        self._invA = dot(expm(-dA), self._invA)
        self._logDetA += 0.5 * self.scaleLearningRate * covTrace
        if self.storeAllDistributions:
            self._allDistributions.append((self._center.copy(), self._A.copy()))

def transfer_f(dw,aas,aai,eps,deltaw,f):
    """
    Args:
    dw: size of the grid spacing
    aas=relative slowness of the signal mode
    aai=relative slowness of the idler mode
    lnl=inverse of the strength of the nonlinearity
    deltaw:  specifies the size of the frequency grid going from
    -deltaw to deltaw for each frequency
    f: shape of the pump function
    """
    ddws=np.arange(-deltaw-dw/2,deltaw+dw/2,dw)
    deltaks=aas*ddws
    ddwi=np.arange(-deltaw-dw/2,deltaw+dw/2,dw)
    deltaki=aai*ddwi
    ds=np.diag(deltaks)
    di=np.diag(deltaki)


    def ff(x,y):
        return f(x+y)
    
    v=eps*(dw)*ff(ddwi[:,None],ddws[None,:])
    G=1j*np.concatenate((np.concatenate((ds,v),axis=1),np.concatenate((-v,-di),axis=1)),axis=0)
    z=1;
    dsi=np.concatenate((deltaks,-deltaki),axis=0)
    U0=linalg.expm(-1j*np.diag(dsi)*z/2)
    GG=np.dot(np.dot(U0,linalg.expm(G)),U0)
    n=len(ddws)
    return (GG[0:n,0:n],GG[n:2*n,0:n],GG[0:n,n:2*n],GG[n:2*n,n:2*n])

def ghz_simult_trajectory(stages):
    rhos = []
    for stage in stages:
        U = np.eye(8, dtype=complex)
        I = pyle.tomo.sigmaI
        X = pyle.tomo.sigmaX
        Y = pyle.tomo.sigmaY
        couple = (pyle.tensor((X,X,I)) + pyle.tensor((Y,Y,I)) +
                  pyle.tensor((X,I,X)) + pyle.tensor((Y,I,Y)) +
                  pyle.tensor((I,X,X)) + pyle.tensor((I,Y,Y))) / 2.0
        
        if stage > 0:
            fraction = np.clip(stage-0, 0, 1)
            H = -1j * (np.pi/2)/2 * (pyle.tensor((Y,I,I)) + pyle.tensor((I,Y,I)) + pyle.tensor((I,I,Y)))
            U = np.dot(expm(fraction * H), U)
    
        if stage > 1:
            fraction = np.clip(stage-1, 0, 1)
            H = -1j * np.pi/2 * couple
            U = np.dot(expm(fraction * H), U)
        
        if stage > 2:
            fraction = np.clip(stage-2, 0, 1)
            H = -1j * (np.pi/2)/2 * (pyle.tensor((X,I,I)) + pyle.tensor((I,X,I)) + pyle.tensor((I,I,X)))
            U = np.dot(expm(fraction * H), U)
    
        psi0 = np.array([1,0,0,0,0,0,0,0], dtype=complex)
        psi_th = np.dot(U, psi0)
        rho_th = pyle.ket2rho(psi_th)
        rhos.append(rho_th)
    return np.array(rhos)

def linear_ode_discretation(F, L=None, Q=None, dt=1):
    n = F.shape[0]

    if L is None:
        L = eye(n)

    if Q is None:
        Q = zeros((n,n))

    A = expm(F*dt)

    phi = zeros((2*n, 2*n))

    phi[0:n,     0:n] = F
    phi[0:n,   n:2*n] = L.dot(Q).dot(L.T)
    phi[n:2*n, n:2*n] = -F.T

    zo = vstack((zeros((n,n)), eye(n)))

    CD = expm(phi*dt).dot(zo)

    C = CD[0:n,:]
    D = CD[n:2*n,:]
    q = C.dot(inv(D))

    return (A, q)

def Iterate_LS_Strain(xys0_pix, xys1_pix, num = 10):
    solve = Get_LS_DGradient(xys0_pix, xys1_pix)
    for j in range(0,num):
        xy_pix = Simulate_Shifted_Spots(xys0_pix, linalg.expm(solve), det2lab_mat)
        dsolve = Get_LS_DGradient(xy_pix, xys1_pix)
        solve = solve+dsolve
    return linalg.expm(solve)

def Char_Gate(NV,res ,B_field=400):
    """
    Characterize the gate, take the NV centre, the resonance paramters and the Bfield as input
    returns the fidelity with which an x-gate can be implemented.
    """


    #data = np.loadtxt("NV_Sim_8.dat") #Placeholder data to test the script
    #NV = np.vstack((data[:,3],data[:,4]))
    #physical constants
    gamma_c = 1.071e3 #g-factor for C13 in Hz/G
    #Model parameters
    omega_larmor = 2*np.pi*gamma_c*B_field
    tau_larmor = 2*np.pi/omega_larmor
    tau = res[0]
    n_pulses = int(res[1]*2) #So that we do a pi -pulse

    Ix = 0.5 * np.array([[0,1],[1,0]])
    Iz = 0.5* np.array([[1,0],[0,-1]])
    H0 = (omega_larmor)*Iz
    exH0 =linalg.expm(-1j*H0*tau)

    S_final =1
    for idC in range(np.shape(NV)[1]):
        A= 2*np.pi*NV[0,idC]
        B= 2*np.pi*NV[1,idC]  #Converts to radial frequency in Hz/G
        H1 = (A+omega_larmor) *Iz +B*Ix
        exH1 = linalg.expm(-1j*H1*tau)
        V0 = exH0.dot(exH1.dot(exH1.dot(exH0)))
        V1 = exH1.dot(exH0.dot(exH0.dot(exH1)))
        S = np.real(np.trace(np.dot(np.linalg.matrix_power(V0,n_pulses),np.linalg.matrix_power(V1,n_pulses)))/2)
        S_final = S_final *S
    F = (1-(S_final+1)/2) #Converting from probability of measuring +X to fidelity of -X (x-rotation)
    return F

def discretize(F,G,Q,Ts):
    Phi = sp_linalg.expm(F*Ts)

    # Control matrix Lambda, not to be used
    L = np.zeros((F.shape[0],1))

    A_z = np.zeros((L.shape[1], F.shape[1]+L.shape[1]))
    A1 = np.vstack((np.hstack((F,L)),A_z))

    Loan1 = sp_linalg.expm(A1*Ts)
    Lambda = Loan1[0:L.shape[0], F.shape[1]:F.shape[1]+L.shape[1]]

    # Covariance
    Qc = symmetrize(np.dot(G, np.dot(Q, G.T)))
    dim = F.shape[0]
    A2 = np.vstack((np.hstack((-F, Qc)), np.hstack((np.zeros((dim,dim)), F.T))))
    Loan2 = sp_linalg.expm(A2*Ts)
    G2 = Loan2[0:dim, dim:2*dim]
    F3 = Loan2[dim:2*dim, dim:2*dim]

    # Calculate Gamma*Gamma.T
    Qd = symmetrize(np.dot(F3.T, G2))
    L = np.linalg.cholesky(Qd)
    
    return Phi, Qd, L

    def pdf(self, x):
        """ probability density function """
        if not np.isscalar(x):
            x = np.asarray(x)
            res = np.zeros_like(x)
            nz = (x > 0)
            if np.any(nz):
                if self.method == 'sum':
                    factor = np.exp(-x[nz, None] * self.rates[..., :]) \
                                / self.rates[..., :]
                    res[nz] = np.sum(self._terms[..., :] * factor, axis=1)
                else:
                    Theta = (np.diag(-self.rates, 0) + 
                             np.diag(self.rates[:-1], 1))
                    for i in np.flatnonzero(nz):
                        res.flat[i] = \
                                1 - linalg.expm(x.flat[i]*Theta)[0, :].sum()
 
        elif x == 0:
            res = 0
        else:
            if self.method == 'sum':
                factor = np.exp(-x*self.rates)/self.ratesx
                res[nz] = np.sum(self._terms * factor)
            else:
                Theta = np.diag(-self.rates, 0) + np.diag(self.rates[:-1], 1)
                res = 1 - linalg.expm(x*Theta)[0, :].sum()
        return res

 def test_padecases_float(self):
     # test single-precision cases
     a1 = eye(3, dtype=float)*1e-1; e1 = exp(1e-1)*eye(3)
     a2 = eye(3, dtype=float);      e2 = exp(1.0)*eye(3)
     a3 = eye(3, dtype=float)*10;   e3 = exp(10.)*eye(3)
     assert_array_almost_equal(expm(a1),e1)
     assert_array_almost_equal(expm(a2),e2)
     assert_array_almost_equal(expm(a3),e3)

    def test_consistency(self):
        a = array([[0.,1],[-1,0]])
        assert_array_almost_equal(expm(a), expm2(a))
        assert_array_almost_equal(expm(a), expm3(a))

        a = array([[1j,1],[-1,-2j]])
        assert_array_almost_equal(expm(a), expm2(a))
        assert_array_almost_equal(expm(a), expm3(a))

def test_QobjExpmExplicitDense():
    "Qobj expm (explicit dense)"
    data = np.random.random(
        (15, 15)) + 1j * np.random.random((15, 15)) - (0.5 + 0.5j)
    A = Qobj(data)
    B = A.expm(method='dense')
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())
    B = A.expm(method='scipy-delse')
    assert_((B.data.todense() - np.matrix(la.expm(data)) < 1e-10).all())

def KPMF(input_matrix, approx=50, iterations=30, learning_rate=.001, adjacency_width=5, adjacency_strength=.5):
    A = input_matrix
    Z = np.asarray(A > 0,dtype=np.int)
    A1d = np.ravel(A)
    mean = np.mean(A1d)
    A = A-mean
    K = approx
    R = itr = iterations
    l = learning_rate
    N = A.shape[0]
    M = A.shape[1]
    U = np.random.randn(N,K)
    V = np.random.randn(K,M)
    #KPMF using gradient descent as per paper
    #Kernelized Probabilistic Matrix Factorization: Exploiting Graphs and Side Information
    #T. Zhou, H. Shan, A. Banerjee, G. Sapiro
    #Using diffusion kernel
    #U are the rows, we use an adjacency matrix CU to reprent connectivity
    #This matrix connects rows +-adjacency_width
    #V are the columns, connected columns are CV
    #Operate on graph laplacian L, which is the degree matrix D - C
    #Applying the diffusion kernel to L, this forms a spatial smoothness graph
    bw = adjacency_width
    #Use scipy.sparse.diags to generate band matrix with bandwidth = 2*adjacency_width+1
    #Example of adjacency_width = 1, N = 4
    #[1 1 0 0]
    #[1 1 1 0]
    #[0 1 1 1]
    #[0 0 1 1]
    print "Running KPMF with:"
    print "learning rate=" + `l`
    print "bandwidth=" + `bw`
    print "beta=" + `b`
    print "approximation rank=" + `K`
    print "iterations=" + `R`
    print ""
    CU = sp.diags([1]*(2*bw+1),range(-bw,bw+1),shape=(N,N)).todense()
    DU = np.diagflat(np.sum(CU,1))
    CV = sp.diags([1]*(2*bw+1),range(-bw,bw+1),shape=(M,M)).todense()
    DV = np.diagflat(np.sum(CV,1))
    LU = DU - CU
    LV = DV - CV
    beta = adjacency_strength
    KU = sl.expm(beta*LU)
    KV = sl.expm(beta*LV)
    SU = np.linalg.pinv(KU)
    SV = np.linalg.pinv(KV)
    for r in range(R):
        for i in range(N):
            for j in range(M):
                if Z[i,j] > 0:
                    e = A[i,j] - np.dot(U[i,:],V[:,j])
                    U[i,:] = U[i,:] + l*(e*V[:,j] - np.dot(SU[i,:],U))
                    V[:,j] = V[:,j] + l*(e*U[i,:] - np.dot(V,SV[:,j]))
    A_ = np.dot(U,V)
    return A_+mean

 def get_symmetry(self,x):
     """Return the symmetry operator"""
     if self.real: 
         M = v2O(x,exp=self.discrete)
         if self.discrete: return M # return the operator
         else: return lg.expm(M) # return the operator
     else:
         M = v2U(x,exp=self.discrete)
         if self.discrete: return M # return the operator
         else: return lg.expm(1j*M) # return the operator

    def test_consistency(self):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore", DeprecationWarning)
            a = array([[0.,1],[-1,0]])
            assert_array_almost_equal(expm(a), expm2(a))
            assert_array_almost_equal(expm(a), expm3(a))

            a = array([[1j,1],[-1,-2j]])
            assert_array_almost_equal(expm(a), expm2(a))
            assert_array_almost_equal(expm(a), expm3(a))

def get_plot_array(N_diploid, Nr, theta_values, Ns_values):
    """
    Compute expected hitting times.
    Theta is 4*N*mu, and the units of time are 4*N*mu generations.
    @param N_diploid: diploid population size
    @param Nr: recombination rate
    @param theta_values: mutation rates
    @param Ns_values: selection values
    @return: arr[i][j] gives time for Ns_values[i] and theta_values[j]
    """
    # set up the state space
    k = 4
    M = multinomstate.get_sorted_states(2 * N_diploid, k)
    T = multinomstate.get_inverse_map(M)
    nstates = M.shape[0]
    lmcs = wfengine.get_lmcs(M)
    # precompute rate matrices
    R_rate = wfcompens.create_recomb(M, T)
    M_rate = wfcompens.create_mutation(M, T)
    # precompute a recombination probability matrix
    R_prob = linalg.expm(Nr * R_rate / float((2 * N_diploid) ** 2))
    #
    arr = []
    for theta in theta_values:
        # Compute the expected number of mutation events per generation.
        mu = theta / 2
        # Precompute the mutation matrix
        # and the product of mutation and recombination.
        M_prob = linalg.expm(mu * M_rate / float(2 * 2 * N_diploid))
        MR_prob = np.dot(M_prob, R_prob)
        #
        row = []
        for Ns in Ns_values:
            s = Ns / float(N_diploid)
            lps = wfcompens.create_selection(s, M)
            S_prob = np.exp(wfengine.create_genic(lmcs, lps, M))
            P = np.dot(MR_prob, S_prob)
            # compute the stationary distribution
            v = MatrixUtil.get_stationary_distribution(P)
            # compute the transition matrix limit at time infinity
            # P_inf = np.outer(np.ones_like(v), v)
            # compute the fundamental matrix Z
            # Z = linalg.inv(np.eye(nstates) - (P - P_inf)) - P_inf
            #
            # Use broadcasting instead of constructing P_inf.
            Z = linalg.inv(np.eye(nstates) - (P - v)) - v
            # compute the hitting time from state AB to state ab.
            i = 0
            j = 3
            hitting_time_generations = (Z[j, j] - Z[i, j]) / v[j]
            hitting_time = hitting_time_generations * theta
            row.append(hitting_time)
        arr.append(row)
    return arr

 def gibbsstate(self):
     nb=self.p.nbar
     hbar=self.p.hbar
     kb=self.p.kb
     nu=self.p.nu
     T=hbar*nu/(kb*np.log((nb + 1.)/nb))
     Z=np.trace(expm(-hbar*nu*self.ada/(kb*T)))
     result = np.kron(expm(-hbar*nu*self.ada/(kb*T)),[[0,0],[0,1]])/Z
     if result[2*(self.p.nmax+1)-1,2*(self.p.nmax+1)-1]>0.0001:
         print 'Warning: nmax may not be high enough for chosen value of nbar'
     return result