From Bern to 4-quark Fierz basis for Classes III, IV and V.
    `qqqq` should be of the form 'sbuc', 'sdcc', 'ucuu' etc.
188801
def _Bern_to_Fierz_III_IV_V(C, qqqq):
    """"""
    # 2nd != 4th, color-octet redundant
    if qqqq in ['sbss', 'dbdd', 'dbds', 'sbsd', 'bsbd', 'dsdd']:
        return {
                'F' + qqqq + '1': C['1' + qqqq] + 16 * C['3' + qqqq],
                'F' + qqqq + '1p': C['1p' + qqqq] + 16 * C['3p' + qqqq],
                'F' + qqqq + '3': C['1' + qqqq] + 4 * C['3' + qqqq],
                'F' + qqqq + '3p': C['1p' + qqqq] + 4 * C['3p' + qqqq],
                'F' + qqqq + '5': C['5p' + qqqq] + 64 * C['9p' + qqqq],
                'F' + qqqq + '5p': C['5' + qqqq] + 64 * C['9' + qqqq],
                'F' + qqqq + '7': C['5p' + qqqq] + 16 * C['9p' + qqqq],
                'F' + qqqq + '7p': C['5' + qqqq] + 16 * C['9' + qqqq],
                'F' + qqqq + '9': C['7p' + qqqq] - 16 * C['9p' + qqqq],
                'F' + qqqq + '9p': C['7' + qqqq] - 16 * C['9' + qqqq],
                }
    if qqqq in ['dbbb', 'sbbb', 'dsss']:  # 2nd = 4th, color-octet redundant
        return {
                'F' + qqqq + '1': C['1' + qqqq] + 16 * C['3' + qqqq],
                'F' + qqqq + '1p': C['1p' + qqqq] + 16 * C['3p' + qqqq],
                'F' + qqqq + '3': C['1' + qqqq] + 4 * C['3' + qqqq],
                'F' + qqqq + '3p': C['1p' + qqqq] + 4 * C['3p' + qqqq],
                'F' + qqqq + '5': C['5' + qqqq] + 64 * C['9' + qqqq],
                'F' + qqqq + '5p': C['5p' + qqqq] + 64 * C['9p' + qqqq],
                'F' + qqqq + '7': C['5' + qqqq] + 16 * C['9' + qqqq],
                'F' + qqqq + '7p': C['5p' + qqqq] + 16 * C['9p' + qqqq],
                'F' + qqqq + '9': C['7' + qqqq] - 16 * C['9' + qqqq],
                'F' + qqqq + '9p': C['7p' + qqqq] - 16 * C['9p' + qqqq],
                }
    # generic case
    if qqqq in ['sbuu', 'sbdd', 'sbuu', 'sbuc', 'sbcu', 'sbcc',
                'dbuu', 'dbss', 'dbuu', 'dbuc', 'dbcu', 'dbcc',
                'dsuu', 'dsbb', 'dsuu', 'dsuc', 'dscu', 'dscc',]:
        return {
                'F' + qqqq + '1': C['1' + qqqq] - C['2' + qqqq] / 6 + 16 * C['3' + qqqq] - (8 * C['4' + qqqq]) / 3,
                'F' + qqqq + '10': -8 * C['10' + qqqq] + C['8' + qqqq] / 2,
                'F' + qqqq + '10p': -8 * C['10p' + qqqq] + C['8p' + qqqq] / 2,
                'F' + qqqq + '1p': C['1p' + qqqq] - C['2p' + qqqq] / 6 + 16 * C['3p' + qqqq] - (8 * C['4p' + qqqq]) / 3,
                'F' + qqqq + '2': C['2' + qqqq] / 2 + 8 * C['4' + qqqq],
                'F' + qqqq + '2p': C['2p' + qqqq] / 2 + 8 * C['4p' + qqqq],
                'F' + qqqq + '3': C['1' + qqqq] - C['2' + qqqq] / 6 + 4 * C['3' + qqqq] - (2 * C['4' + qqqq]) / 3,
                'F' + qqqq + '3p': C['1p' + qqqq] - C['2p' + qqqq] / 6 + 4 * C['3p' + qqqq] - (2 * C['4p' + qqqq]) / 3,
                'F' + qqqq + '4': C['2' + qqqq] / 2 + 2 * C['4' + qqqq],
                'F' + qqqq + '4p': C['2p' + qqqq] / 2 + 2 * C['4p' + qqqq],
                'F' + qqqq + '5': -((32 * C['10' + qqqq]) / 3) + C['5' + qqqq] - C['6' + qqqq] / 6 + 64 * C['9' + qqqq],
                'F' + qqqq + '5p': -((32 * C['10p' + qqqq]) / 3) + C['5p' + qqqq] - C['6p' + qqqq] / 6 + 64 * C['9p' + qqqq],
                'F' + qqqq + '6': 32 * C['10' + qqqq] + C['6' + qqqq] / 2,
                'F' + qqqq + '6p': 32 * C['10p' + qqqq] + C['6p' + qqqq] / 2,
                'F' + qqqq + '7': -((8 * C['10' + qqqq]) / 3) + C['5' + qqqq] - C['6' + qqqq] / 6 + 16 * C['9' + qqqq],
                'F' + qqqq + '7p': -((8 * C['10p' + qqqq]) / 3) + C['5p' + qqqq] - C['6p' + qqqq] / 6 + 16 * C['9p' + qqqq],
                'F' + qqqq + '8': 8 * C['10' + qqqq] + C['6' + qqqq] / 2,
                'F' + qqqq + '8p': 8 * C['10p' + qqqq] + C['6p' + qqqq] / 2,
                'F' + qqqq + '9': (8 * C['10' + qqqq]) / 3 + C['7' + qqqq] - C['8' + qqqq] / 6 - 16 * C['9' + qqqq],
                'F' + qqqq + '9p': (8 * C['10p' + qqqq]) / 3 + C['7p' + qqqq] - C['8p' + qqqq] / 6 - 16 * C['9p' + qqqq],
                }
    raise ValueError("Case not implemented: {}".format(qqqq))
668608
def _run_sm_scale_in(self, C_out, scale_sm=91.1876):
        """"""
        # initialize an empty SMEFT instance
        smeft_sm = SMEFT()
        C_in_sm = beta.C_array2dict(np.zeros(9999))
        # set the SM parameters to the values obtained from smpar.smeftpar
        C_SM = smpar.smeftpar(scale_sm, self.scale_high, C_out, basis='Warsaw')
        C_SM = {k: v for k, v in C_SM.items() if k in definitions.SM_keys}
        # set the Wilson coefficients at the EW scale to C_out
        C_in_sm.update(C_out)
        C_in_sm.update(C_SM)
        smeft_sm.set_initial(C_in_sm, scale_sm, scale_high=self.scale_high)
        # run up (with 1% relative precision, ignore running of Wilson coefficients)
        C_SM_high = smeft_sm.rgevolve(self.scale_in, newphys=False, rtol=0.01, atol=1)
        return {k: v for k, v in C_SM_high.items() if k in definitions.SM_keys}
280319
def quatInv(q):
    """"""
    qConj = -q
    qConj[0] = -qConj[0]
    normSqr = multiplyQuats(q, qConj)[0]
    return qConj/normSqr
1208923
def warsaw_to_warsawmass(C, parameters=None, sectors=None):
    
    p = default_parameters.copy()
    if parameters is not None:
        # if parameters are passed in, overwrite the default values
        p.update(parameters)
    # start out with a 1:1 copy
    C_out = C.copy()
    # rotate left-handed up-type quark fields in uL-uR operator WCs
    C_rotate_u = ['uphi', 'uG', 'uW', 'uB']
    for name in C_rotate_u:
        _array = smeft_toarray(name, C)
        V = ckmutil.ckm.ckm_tree(p["Vus"], p["Vub"], p["Vcb"], p["delta"])
        UuL = V.conj().T
        _array = UuL.conj().T @ _array
        _dict = smeft_fromarray(name, _array)
        C_out.update(_dict)
    # diagonalize dimension-5 Weinberg operator
    _array = smeft_toarray('llphiphi', C)
    _array = np.diag(ckmutil.diag.msvd(_array)[1])
    _dict = smeft_fromarray('llphiphi', _array)
    C_out.update(_dict)
    return C_out
1290013
def as_qubit_order(val: 'qubit_order_or_list.QubitOrderOrList'
                       ) -> 'QubitOrder':
        
        if isinstance(val, collections.Iterable):
            return QubitOrder.explicit(val)
        if isinstance(val, QubitOrder):
            return val
        raise ValueError(
            "Don't know how to interpret <{}> as a Basis.".format(val))
473464
def canonical(self):
        
        Ac1 = np.hstack((self.deltah, np.zeros((self.nh, self.nz))))
        Ac2 = np.hstack((np.zeros((self.nz, self.nh)), self.a22))
        Ac = np.vstack((Ac1, Ac2))
        Bc = np.vstack((self.thetah, np.zeros((self.nz, self.nc))))
        Cc = np.vstack((np.zeros((self.nh, self.nw)), self.c2))
        Rc1 = np.hstack((self.llambda.T.dot(self.llambda), -
                         self.llambda.T.dot(self.ub)))
        Rc2 = np.hstack((-self.ub.T.dot(self.llambda), self.ub.T.dot(self.ub)))
        Rc = np.vstack((Rc1, Rc2))
        Qc = self.pih.T.dot(self.pih)
        Nc = np.hstack(
            (self.pih.T.dot(self.llambda), -self.pih.T.dot(self.ub)))

        lq_aux = LQ(Qc, Rc, Ac, Bc, N=Nc, beta=self.beta)

        P1, F1, d1 = lq_aux.stationary_values()

        self.F_b = F1[:, 0:self.nh]
        self.F_f = F1[:, self.nh:]

        self.pihat = np.linalg.cholesky(self.pih.T.dot(
            self.pih) + self.beta.dot(self.thetah.T).dot(P1[0:self.nh, 0:self.nh]).dot(self.thetah)).T
        self.llambdahat = self.pihat.dot(self.F_b)
        self.ubhat = - self.pihat.dot(self.F_f)

        return
1352597
def symbolic_rotation_matrix(phi, theta, symbolic_psi):
    """"""
    return sympy.Matrix(Rz_matrix(phi)) * sympy.Matrix(Rx_matrix(theta)) * symbolic_Rz_matrix(symbolic_psi)
1302442
def xi2_from_mass1_mass2_spin2x_spin2y(mass1, mass2, spin2x, spin2y):
    
    q = q_from_mass1_mass2(mass1, mass2)
    a1 = 2 + 3 * q / 2
    a2 = 2 + 3 / (2 * q)
    return a1 / (q**2 * a2) * chi_perp_from_spinx_spiny(spin2x, spin2y)
893692
def Fierz_to_Bern_lep(C, ddll):
    """"""
    ind = ddll.replace('l_','').replace('nu_','')
    dic = {
        '1' + ind : 5 * C['F'+ ind + '10'] / 3 + C['F'+ ind + '9'],
        '3' + ind : -C['F' + ind + '10'] / 6,
        '5' + ind : C['F' + ind + 'S'] - 5 * C['F' + ind + 'P'] / 3,
        '7' + ind : 2 * C['F' + ind + 'P'] / 3 + C['F' + ind + 'T']
                    + C['F' + ind + 'T5'],
        '9' + ind : C['F' + ind + 'P'] / 24,
        '1p' + ind : C['F' + ind + '9p'] - 5 * C['F' + ind + '10p'] / 3,
        '3p' + ind : C['F' + ind + '10p'] / 6,
        '5p' + ind : 5 * C['F' + ind + 'Pp'] / 3 + C['F' + ind + 'Sp'],
        '7p' + ind : -2 * C['F' + ind + 'Pp'] / 3 + C['F' + ind + 'T']
                    - C['F' + ind + 'T5'],
        '9p' + ind : -C['F' + ind + 'Pp'] / 24,
    }
    return dic
882987
def _Bern_to_JMS_I(C, qq):
    """'"""
    if qq in ['sb', 'db', 'ds']:
        dd = 'dd'
        ij = '{}{}'.format(dflav[qq[0]] + 1, dflav[qq[1]] + 1)
    elif qq == 'cu':
        dd = 'uu'
        ij = '{}{}'.format(uflav[qq[0]] + 1, uflav[qq[1]] + 1)
    else:
        raise ValueError("not in Bern_I: ".format(qq))
    ji = ij[1] + ij[0]
    d = {"V{}LL_{}{}".format(dd, ij, ij): C['1' + 2 * qq],
         "S1{}RR_{}{}".format(dd, ji, ji): C['2' + 2 * qq].conjugate() + C['3' + 2 * qq].conjugate() / 3,
         "S8{}RR_{}{}".format(dd, ji, ji): 2 * C['3' + 2 * qq].conjugate(),
         "V1{}LR_{}{}".format(dd, ij, ij): -C['4' + 2 * qq] / 6 - C['5' + 2 * qq] / 2,
         "V8{}LR_{}{}".format(dd, ij, ij): -C['4' + 2 * qq],
         "V{}RR_{}{}".format(dd, ij, ij): C['1p' + 2 * qq],
         "S1{}RR_{}{}".format(dd, ij, ij): C['2p' + 2 * qq] + C['3p' + 2 * qq] / 3,
         "S8{}RR_{}{}".format(dd, ij, ij): 2 * C['3p' + 2 * qq],
         }
    if qq == 'cu':
        # here we need to convert some operators that are not in the basis
        for VXY in ['VuuRR', 'V1uuLR', 'V8uuLR', 'VuuLL']:
            d[VXY + '_1212'] = d.pop(VXY + '_2121').conjugate()
    return d
1164261
def mixing_phases(U):
    """"""
    f = {}
    # angles
    f['t13'] = asin(abs(U[0,2]))
    if U[0,0] == 0:
        f['t12'] = pi/2
    else:
        f['t12'] = atan(abs(U[0,1])/abs(U[0,0]))
    if U[2,2] == 0:
        f['t23'] = pi/2
    else:
        f['t23'] = atan(abs(U[1,2])/abs(U[2,2]))
    s12 = sin(f['t12'])
    c12 = cos(f['t12'])
    s13 = sin(f['t13'])
    c13 = cos(f['t13'])
    s23 = sin(f['t23'])
    c23 = cos(f['t23'])
    # standard phase
    if (s12*s23) == 0 or (c12*c13**2*c23*s13) == 0:
        f['delta'] = 0
    else:
        f['delta'] = -phase((U[0,0].conj()*U[0,2]*U[2,0]*U[2,2].conj()/(c12*c13**2*c23*s13) + c12*c23*s13)/(s12*s23))
    # Majorana phases
    f['delta1']  = phase(exp(1j*f['delta']) * U[0, 2])
    f['delta2']  = phase(U[1, 2])
    f['delta3']  = phase(U[2, 2])
    f['phi1'] = 2*phase(exp(1j*f['delta1']) * U[0, 0].conj())
    f['phi2'] = 2*phase(exp(1j*f['delta1']) * U[0, 1].conj())
    return f
876339
def numpartition_qaoa(n_step, nums, minimizer=None, sampler=None):
    
    hamiltonian = pauli.Expr.zero()
    for i, x in enumerate(nums):
        hamiltonian += pauli.Z[i] * x
    hamiltonian = (hamiltonian ** 2).simplify()

    return vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, n_step), minimizer, sampler)
465470
def spin2x_from_mass1_mass2_xi2_phi_a_phi_s(mass1, mass2, xi2, phi_a, phi_s):
    
    chi_perp = chi_perp_from_mass1_mass2_xi2(mass1, mass2, xi2)
    phi2 = phi2_from_phi_a_phi_s(phi_a, phi_s)
    return chi_perp * numpy.cos(phi2)
929826
def phi_a(mass1, mass2, spin1x, spin1y, spin2x, spin2y):
    
    phi1 = phi_from_spinx_spiny(primary_spin(mass1, mass2, spin1x, spin2x),
                                primary_spin(mass1, mass2, spin1y, spin2y))
    phi2 = phi_from_spinx_spiny(secondary_spin(mass1, mass2, spin1x, spin2x),
                                secondary_spin(mass1, mass2, spin1y, spin2y))
    return (phi1 - phi2) % (2 * numpy.pi)
988117
def mf_aBl(self):
        
        mf_aBl = self._mf_aBl = np.zeros((self.T, self.num_states))
        ids, dds, eds = self.init_dynamics_distns, self.dynamics_distns, \
            self.emission_distns

        for idx, (d1, d2, d3) in enumerate(zip(ids, dds, eds)):
            mf_aBl[0,idx] = d1.expected_log_likelihood(
                stats=self.E_init_stats)
            mf_aBl[:-1,idx] += d2.expected_log_likelihood(
                stats=self.E_dynamics_stats)
            mf_aBl[:,idx] += d3.expected_log_likelihood(
                stats=self.E_emission_stats)

        mf_aBl[np.isnan(mf_aBl).any(1)] = 0.
        return mf_aBl
1230665
def Bern_to_Fierz_nunu(C,ddll):
    """"""
    ind = ddll.replace('l_','').replace('nu_','')
    return {
            'F' + ind + 'nu': C['nu1' + ind],
            'F' + ind + 'nup': C['nu1p' + ind],
            }
1175004
def _BernI_to_Flavio_I(C, qq):
    """'"""
    qqf = qq[::-1]  # flavio uses "bs" instead of "sb" etc.
    if qq in ['sb', 'db', 'ds', 'cu']:
        return {
            'CVLL_' + 2*qqf: C["1" + 2*qq],
            'CSLL_' + 2*qqf: C["2" + 2*qq] - 1 / 2 * C["3" + 2*qq],
            'CTLL_' + 2*qqf: -1 / 8 * C["3" + 2*qq],
            'CVLR_' + 2*qqf: -1 / 2 * C["5" + 2*qq],
            'CVRR_' + 2*qqf: C["1p" + 2*qq],
            'CSRR_' + 2*qqf: C["2p" + 2*qq] - 1 / 2 * C["3p" + 2*qq],
            'CTRR_' + 2*qqf: -1 / 8 * C["3p" + 2*qq],
            'CSLR_' + 2*qqf: C["4" + 2*qq]
            }
    else:
        raise ValueError("not in Flavio_I: ".format(qq))