r"""Calculate the matrices omega_ij, gamma_ij, r_pij.

    This function calculates the matrices omega_ij, gamma_ij and r_pij given a
    list of atomic states. The states can be arbitrarily in their fine,
    hyperfine or magnetic detail.
333608
def calculate_matrices(states, Omega=1):
    
    """
    # We check that all states belong to the same element and the same isotope.
    iso = states[0].isotope
    element = states[0].element
    for state in states[1:]:
        if state.element != element:
            raise ValueError('All states must belong to the same element.')
        if state.isotope != iso:
            raise ValueError('All states must belong to the same isotope.')

    # We find the fine states involved in the problem.
    fine_states = find_fine_states(states)

    # We find the full magnetic states. The matrices will be first calculated
    # for the complete problem and later reduced to include only the states of
    # interest.
    full_magnetic_states = make_list_of_states(fine_states, 'magnetic',
                                               verbose=0)

    # We calculate the indices corresponding to each sub matrix of fine and
    # hyperfine levels.

    # We calculate the frequency differences between states.
    omega_full = calculate_omega_matrix(full_magnetic_states, Omega)
    # We calculate the matrix gamma.

    gamma_full = calculate_gamma_matrix(full_magnetic_states, Omega)
    # We calculate the reduced matrix elements
    reduced_matrix_elements = calculate_reduced_matrix_elements(fine_states)

    # We calculate the matrices r_-1, r_0, r_1
    r_full = calculate_r_matrices(fine_states, reduced_matrix_elements)

    # Reduction to be implemented
    omega = omega_full
    r = r_full
    gamma = gamma_full
    return omega, gamma, r
879773
def getMatrix(self):
        
        if self.mflag:
            m0 = 9.10938215e-31
            e0 = 1.602176487e-19
            c0 = 299792458.0
            rho = np.sqrt(self.gamma ** 2 - 1) * m0 * c0 / self.bend_field / e0
            theta = np.arcsin(self.bend_length / rho)
            self.bangle = theta

            m_rb_1 = transRbend(theta, rho, self.gamma, -1)
            m_rb_2 = transRbend(-theta, -rho, self.gamma, 1)
            m_rb_3 = transRbend(-theta, -rho, self.gamma, -1)
            m_rb_4 = transRbend(theta, rho, self.gamma, 1)
            m_df_12 = transDrift(self.dflist[0], self.gamma)
            m_df_23 = transDrift(self.dflist[1], self.gamma)
            m_df_34 = transDrift(self.dflist[2], self.gamma)

            self.transM = reduce(np.dot, [m_rb_1, m_df_12, m_rb_2, m_df_23, m_rb_3, m_df_34, m_rb_4])

            return self.transM
444924
def nekmin(omega_in,q,x0=0.5,z0=0.5):

    ''''''

    def Omega_xz(q,x,z):
        return 1./np.sqrt(x**2+z**2)+q/np.sqrt((1-x)**2+z**2)+(q+1)*x**2/2.-q*x

    def Omega_xy(q,x,y):
        return 1./np.sqrt(x**2+y**2)+q/np.sqrt((1-x)**2+y**2)+(q+1)*(x**2+y**2)/2.-q*x

    def dOmegadx_z(q,x,z):
        return -x/(x**2+z**2)**(3./2)+q*(1-x)/((1-x)**2+z**2)**(3./2.)+(q+1)*x-q

    def dOmegadx_y(q,x,y):
        return -x/(x**2+y**2)**(3./2)+q*(1-x)/((1-x)**2+y**2)**(3./2.)+(q+1)*x-q

    def dOmegadz(q,x,z):
        return -z/(x**2+z**2)**(3./2)-q*z/((1-x)**2+z**2)**(3./2.)

    def dOmegady(q,x,y):
        return -y/(x**2+y**2)**(3./2)-q*y/((1-x)**2+y**2)**(3./2.)+(q+1)*y

    def d2Omegadx2_z(q,x,z):
        return (2*x**2-z**2)/(x**2+z**2)**(5./2)+q*(2*(1-x)**2-z**2)/((1-x)**2+z**2)**(5./2)+(q+1)

    def d2Omegadx2_y(q,x,y):
        return (2*x**2-y**2)/(x**2+y**2)**(5./2)+q*(2*(1-x)**2-y**2)/((1-x)**2+y**2)**(5./2)+(q+1)

    def d2Omegadxdz(q,x,z):
        return 3*x*z/(x**2+z**2)**(5./2)-3*q*x*(1-x)/((1-x)**2+z**2)**(5./2)

    def d2Omegadxdy(q,x,y):
        return 3*x*y/(x**2+y**2)**(5./2)-3*q*x*(1-x)/((1-x)**2+y**2)**(5./2)

    xz,z = x0,z0
    dxz, dz = 1.,1.

    # find solution in xz plane
    while abs(dxz)>1e-8 and abs(dz)>1e-8:

        delz = 1.
        z=0.05
        while abs(delz) > 0.000001:
            delom = omega_in - Omega_xz(q,xz,z)
            delz = delom/dOmegadz(q,xz,z)
            z = abs(z+delz)

        DN = np.array([[dOmegadx_z(q,xz,z),dOmegadz(q,xz,z)],[d2Omegadx2_z(q,xz,z),d2Omegadxdz(q,xz,z)]])
        EN = np.array([omega_in-Omega_xz(q,xz,z),(-1)*dOmegadx_z(q,xz,z)])

        a,b,c,d = DN[0][0],DN[0][1],DN[1][0],DN[1][1]

        if (a*d-b*c)!=0.:
            DNINV = 1./(a*d-b*c)*np.array([[d,(-1)*b],[(-1)*c,d]])
            #DNINV = inv(DN)

            dd = np.dot(DNINV,EN)
            dxz,dz = dd[0],dd[1]
            xz=xz+dxz
            z=z+dz
        else:
            xz = xz+0.5
            z = z+0.5
            dxz = 1.
            dz = 1.

    return xz,z
42529
def calc_state_matrix(self):
        
        system = self.system

        Gyx = matrix(system.dae.Gx)
        self.solver.linsolve(system.dae.Gy, Gyx)

        self.As = matrix(system.dae.Fx - system.dae.Fy * Gyx)

        # ------------------------------------------------------
        # TODO: use scipy eigs
        # self.As = sparse(self.As)
        # I = np.array(self.As.I).reshape((-1,))
        # J = np.array(self.As.J).reshape((-1,))
        # V = np.array(self.As.V).reshape((-1,))
        # self.As = csr_matrix((V, (I, J)), shape=self.As.size)
        # ------------------------------------------------------
        return self.As
1320778
def calculate_iI_correspondence(omega):
    """
    Ne = len(omega[0])
    om = omega[0][0]
    correspondence = []
    I = 0
    for i in range(Ne):
        if omega[i][0] != om:
            om = omega[i][0]
            I += 1
        correspondence += [(i+1, I+1)]
    Nnd = I+1

    def I_nd(i):
        return correspondence[i-1][1]

    def i_d(I):
        for i in range(Ne):
            if correspondence[i][1] == I:
                return correspondence[i][0]

    return i_d, I_nd, Nnd
326048
def convertToRateMatrix(self, Q):
        
        rowSums             = Q.sum(axis=1).getA1()
        idxRange            = np.arange(Q.shape[0])
        Qdiag               = coo_matrix((rowSums,(idxRange,idxRange)),shape=Q.shape).tocsr()
        return Q-Qdiag
70614
def invariants(mol):
    """"""
    atom_ids = {}
    for a in mol.atoms:
        components = []
        components.append(a.number)
        components.append(len(a.oatoms))
        components.append(a.hcount)
        components.append(a.charge)
        components.append(a.mass)
        if len(a.rings) > 0:
            components.append(1)

        atom_ids[a.index] = gen_hash(components)

    return atom_ids
924100
function qrJs() {
                var m = this.rows();
                var n = this.cols();
                var Q = Matrix$b.I(m);
                var A = this;
                
                for(var k = 1; k < Math.min(m, n); k++) {
            	var ak = A.slice(k, 0, k, k).col(1);
            	var oneZero = [1];
            	
            	while(oneZero.length <=  m - k)
            	    oneZero.push(0);
            	
            	oneZero = $V(oneZero);
            	var vk = ak.add(oneZero.x(ak.norm() * Math.sign(ak.e(1))));
            	var Vk = $M(vk);
            	var Hk = Matrix$b.I(m - k + 1).subtract(Vk.x(2).x(Vk.transpose()).div(Vk.transpose().x(Vk).e(1, 1)));
            	var Qk = identSize(Hk, m, n, k);
            	A = Qk.x(A);
            	// slow way to compute Q
            	Q = Q.x(Qk);
                }
                
                return {Q: Q, R: A};
            }
212048
function(a, b, c, d) {
            	var elements = [], ni = c, i, nj, j;
            	var rows = this.elements.length, cols = this.elements[0].length;
            	while (ni--) {
            	    i = c - ni - 1;
            	    elements[i] = [];
            	    nj = d;
            	    while (nj--) {
            		j = d - nj - 1;
            		elements[i][j] = this.elements[(a + i - 1) % rows][(b + j - 1) % cols];
            	    }
            	}
            	return Matrix$b.create(elements);
                }
470058
def get_matrix(theta, phi, psi, dx, dy, dz,
                   matrix=np.zeros((4, 4), dtype=DTYPE),
                   angles=np.zeros(3, dtype=DTYPE)):
        
        # NB!: matrix and angles by default are being overwritten on each call
        # thus, only created once at compile time.

        angles[0] = theta
        angles[1] = phi
        angles[2] = psi

        cx, cy, cz = np.cos(angles)
        sx, sy, sz = np.sin(angles)

        rotation = matrix[:3, :3]
        rotation.flat = (cx * cz - sx * cy * sz,
                         cx * sz + sx * cy * cz, sx * sy,
                         -sx * cz - cx * cy * sz,
                         -sx * sz + cx * cy * cz, cx * sy,
                         sy * sz,
                         -sy * cz, cy)

        # Translation component
        matrix[:3, 3] = dx, dy, dz
        matrix[3, 3] = 1.
        return matrix
756248
def calc_z(bands, filling, interaction, hund_cu, name):
    """"""

    while True:
        try:
            data = np.load(name+'.npz')
            break
        except IOError:
            dopout = []
            for dop in filling:
                slsp = Spinon(slaves=2*bands, orbitals=bands, \
                             hopping=[0.5]*2*bands, populations=[dop]*2*bands)
                dopout.append(solve_loop(slsp, interaction, hund_cu)[0][0])
            np.savez(name, zeta=dopout, u_int=interaction, doping=filling, hund=hund_cu)

    return data
432953
def Theta(i,j,theta,omega_rescaled,omega_min,
			detunings,detuningsij,combinations,detuning_indices,
			Lij,i_d,I_nd,Nnd,
			states=None,verbose=1,other_the=None):
	""""""

	if i==j:
		return ''
	elif j>i:
		raise ValueError,'i should never be less than j.'

	Ne=len(omega_rescaled[0]); Nl=len(theta[0])
	if other_the!=None:
		the=other_the[:]
	else:
		the=[theta[j-1][l]-theta[i-1][l] for l in range(Nl)]

	###########################################################################################
	#This part is about finding detunings that reduce to -omega_ij
	if the==[0 for l in range(Nl)]:
		if omega_rescaled[i-1][j-1]==0:
			return ''
		#We need a combination of detunings that yields -omega_i,j.
		#According to equation labeled "detuning-exception1"
		#-omega_ij= delta^l_ik - delta^l_jk
		#so we seek a k that is lower than i and j
		#and an l such that l is in L_ik and in L_jk

		ii=min(i,j)
		band1=False
		band2=False
		for k in range(1,ii):
			for l in range(1,Nl+1):
				if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
					s   ='I found that -omega_'+str(i)+','+str(j)
					s+='= delta^'+str(l)+'_'+str(i)+','+str(k)
					s+= '-delta^'+str(l)+'_'+str(j)+','+str(k)
					s+='\n and also l='+str(l)
					s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
					s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'

					#print s
					band1=True
					break
			if band1:
				break

		if not band1:
			#If no such a k exists then we follow te equation labeled "detuning-exception2"
			#-omega_ij= delta^l_kj - delta^l_ki
			#to look for a k greater than i and j
			#and an l such that l is in L_ik and in L_jk
			for k in range(ii+1,Ne+1):
				for l in range(1,Nl+1):
					if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
						s   ='I found that -omega_'+str(i)+','+str(j)
						s+='= delta^'+str(l)+'_'+str(k)+','+str(j)
						s+='-delta^'+str(l)+'_'+str(k)+','+str(i)
						s+='\n and also l='+str(l)
						s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
						s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'
						#print s
						band2=True
						break

				if band2:
					break

		if band1:
			The=''
			#We need to find which detunings are delta^l_i,k and delta^l_j,k.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(i) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(j) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		elif band2:
			The=''
			#We need to find which detunings are delta^l_k,j and delta^l_k,i.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(j):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(i):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		else:
			if verbose>1: print 'WARNING: Optical frequencies will be used instead for -omega_',i,j,'=',omega_rescaled[i-1][j-1],'\n'
			return format_double(-omega_rescaled[i-1][j-1])
	###########################################################################################
	#This part is about finding detunings that reduce to (theta_j -theta_i -omega_ij)

	#We establish to which omegaij the detunings should reduce to
	omegaij=[0 for k in range(Nnd)]
	omegaij[I_nd(i)-1]=1;omegaij[I_nd(j)-1]=-1

	#We search for a combination of detunings that reduces to omegaij
	#That is a linear combination of detunings with coefficients the
	#that reduces to omegaij. Here
	detuning_failure=True
	for comb in combinations:
		suma=[sum([the[l]*detunings[l][comb[l]][ii] for l in range(Nl)])  for ii in range(Nnd)]
		if suma==omegaij:
			detuning_failure=False
			break

	#We stop if no combination reduces to the needed expression.
	if detuning_failure:

		if verbose>1: print 'We will see if it is possible to express Theta_'+str(i)+','+str(j)+'=theta_'+str(j)+'-theta_'+str(i)+'-omega_'+str(i)+','+str(j)
		if verbose>1: print 'in terms of other indices a,b such that omega_ab=omega_ij and transition i -> j is allowed by Lij.'

		band3=False
		for a in range(Ne):
			for b in range(a):
				if omega_rescaled[a][b]==omega_rescaled[i-1][j-1] and Lij[a][b]!=[]:
					band3=True
					break
			if band3: break
		if verbose>1: print band3,a,b,i,j

		a=a+1; b=b+1
		if band3 and (i!=a or j!=b):
			if verbose>1: print omega_rescaled[i-1][j-1],omega_rescaled[a-1][b-1]
			if verbose>1: print the
			if verbose>1: print 'This was possible for omega_'+str(a)+','+str(b)
			return Theta(a,b,theta,omega_rescaled,omega_min,
					detunings,detuningsij,combinations,detuning_indices
					,Lij,other_the=the,verbose=verbose,states=states)
		else:
			#verbose=2
			#print 111
			#qi=states[i-1].quantum_numbers
			#qj=states[j-1].quantum_numbers
			#print verbose
			if verbose>0: print 'WARNING: It was impossible to express laser frequencies',the
			if verbose>0: print 'and atomic transition -omega_',i,j,'=',-omega_rescaled[i-1][j-1]
			if verbose>0: print 'in terms of detunings from the transitions given by Lij'
			#for i in range(Ne):
			#	print i+1,Lij[i]
			#print
			if verbose>0: print 'I Will use optical frequencies instead.'
			The=''
			#We give the optical frequencies
			for l in range(Nl):
				a=the[l]
				if a==1:
					The+='+'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				elif a==-1:
					The+='-('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'
				elif a==0:
					The+=''
				elif a>0:
					The+='+'+str(a)+'*'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				else:
					The+=    str(a)+'*('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'

			#We substract omega_ij
			The+=format_double(-omega_rescaled[i-1][j-1])
			if verbose>1: print The
			if verbose>1: print
			return The

	#For each optical frequency in the, we write the corresponding detuning.
	#This way of assigining a global index ll to the detunings ammounts to
	#   ll=   number_of_previous_detunings
	#       + number_of_detuning_ordered_by_row_and_from_left_to_right_column
	The=''
	acum=0
	####################################################################
	#print 'comb',comb,'the',the,'i,j',i,j
	for l in range(Nl):
		if the[l]==1:
			The+='+detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==-1:
			The+='-detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==0:
			pass
		elif the[l]>0:
			The+='+'+str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		else:
			The+= str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		acum+=detuning_indices[l]

	if The[:1]=='+':
		The=The[1:]

	####################################################################
	#print 'The',The
	return The
676261
def inter_spin_hamiltonian(self, u_int, J_coup):
        """"""
        J_coup *= u_int
        h_int  = (u_int - 2*J_coup)/2.*self.oper['sumSz2']
        h_int += J_coup*self.oper['sumSz-sp2']
        h_int -= J_coup/2.*self.oper['sumSz-or2']
        h_int -= J_coup*self.oper['Sfliphop']

        return h_int
1241877
def update(self, Pi, Tij):
         """
        from msmtools import analysis as msmana

        # update transition matrix by copy
        self._Tij = np.array(Tij)
        assert msmana.is_transition_matrix(self._Tij), 'Given transition matrix is not a stochastic matrix'
        assert self._Tij.shape[0] == self._nstates, 'Given transition matrix has unexpected number of states '
        # reset spectral decomposition
        self._spectral_decomp_available = False

        # check initial distribution
        assert np.all(Pi >= 0), 'Given initial distribution contains negative elements.'
        assert np.any(Pi > 0), 'Given initial distribution is zero'
        self._Pi = np.array(Pi) / np.sum(Pi)
1283103
def get_pij_matrix(t, diag, A, A_inv):
    
    return A.dot(np.diag(np.exp(diag * t))).dot(A_inv)
439727
def spin_z_op(param, oper):
    """"""
    slaves = param['slaves']
    oper['Sz'] = np.array([spin_z(slaves, spin) for spin in range(slaves)])
    oper['Sz+1/2'] = oper['Sz'] + 0.5*np.eye(2**slaves)
    oper['sumSz2'] = oper['Sz'].sum(axis=0)**2  # because Sz is diagonal
    Sz_mat_shape = oper['Sz'].reshape(param['orbitals'], 2, 2**slaves, 2**slaves)
    oper['sumSz-sp2'] = (Sz_mat_shape.sum(axis=1)**2).sum(axis=0)
    oper['sumSz-or2'] = (Sz_mat_shape.sum(axis=0)**2).sum(axis=0)
443585
def U_ij(q_vars: List[fl.Var], mass: np.ndarray, i: int, j: int):
    """"""
    # Check that the lengths are consistent
    assert len(q_vars) == 3 * len(mass)
    
    # Masses of the bodies i and j
    mi = mass[i]
    mj = mass[j]
    
    # Gravitational potential is -G * m1 * m2 / r
    U = -(G * mi * mj) / flux_r(q_vars, i, j)
    return U