Return the fidelity F(rho0, rho1) between two mixed quantum states.

    Note: Fidelity cannot be calculated entirely within the tensor backend.
189069
def fidelity(rho0: Density, rho1: Density) -> float:
    
    assert rho0.qubit_nb == rho1.qubit_nb   # FIXME

    rho1 = rho1.permute(rho0.qubits)
    op0 = asarray(rho0.asoperator())
    op1 = asarray(rho1.asoperator())

    fid = np.real((np.trace(sqrtm(sqrtm(op0) @ op1 @ sqrtm(op0)))) ** 2)
    fid = min(fid, 1.0)
    fid = max(fid, 0.0)     # Correct for rounding errors

    return fid
1062527
def chi_eff_from_spherical(mass1, mass2, spin1_a, spin1_polar,
                           spin2_a, spin2_polar):
    """"""
    spin1z = spin1_a * numpy.cos(spin1_polar)
    spin2z = spin2_a * numpy.cos(spin2_polar)
    return chi_eff(mass1, mass2, spin1z, spin2z)
404117
def inverse_exp_decay(max_step, min_value=0.01, step=None):
  """"""
  inv_base = tf.exp(tf.log(min_value) / float(max_step))
  if step is None:
    step = tf.train.get_global_step()
  if step is None:
    return 1.0
  step = to_float(step)
  return inv_base**tf.maximum(float(max_step) - step, 0.0)
507071
def get_inspiral_tf(tc, mass1, mass2, spin1, spin2, f_low, n_points=100,
        pn_2order=7, approximant='TaylorF2'):
    
    # handle param-dependent approximant specification
    class Params:
        pass
    params = Params()
    params.mass1 = mass1
    params.mass2 = mass2
    params.spin1z = spin1
    params.spin2z = spin2
    try:
        approximant = eval(approximant, {'__builtins__': None},
                           dict(params=params))
    except NameError:
        pass

    if approximant in ['TaylorF2', 'SPAtmplt']:
        from pycbc.waveform.spa_tmplt import findchirp_chirptime

        # FIXME spins are not taken into account
        f_high = f_SchwarzISCO(mass1 + mass2)
        track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(f_high),
                                 n_points)
        track_t = numpy.array([findchirp_chirptime(float(mass1), float(mass2),
                                        float(f), pn_2order) for f in track_f])
    elif approximant in ['SEOBNRv2', 'SEOBNRv2_ROM_DoubleSpin',
                         'SEOBNRv2_ROM_DoubleSpin_HI']:
        f_high = get_final_freq('SEOBNRv2', mass1, mass2, spin1, spin2)
        track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(f_high),
                                 n_points)
        # use HI function as it has wider freq range validity
        track_t = numpy.array([
                lalsimulation.SimIMRSEOBNRv2ROMDoubleSpinHITimeOfFrequency(f,
                    solar_mass_to_kg(mass1), solar_mass_to_kg(mass2),
                    float(spin1), float(spin2)) for f in track_f])
    elif approximant in ['SEOBNRv4', 'SEOBNRv4_ROM']:
        f_high = get_final_freq('SEOBNRv4', mass1, mass2, spin1, spin2)
        # use frequency below final freq in case of rounding error
        track_f = numpy.logspace(numpy.log10(f_low), numpy.log10(0.999*f_high),
                                 n_points)
        track_t = numpy.array([
                lalsimulation.SimIMRSEOBNRv4ROMTimeOfFrequency(
                        f, solar_mass_to_kg(mass1), solar_mass_to_kg(mass2),
                        float(spin1), float(spin2)) for f in track_f])
    else:
        raise ValueError('Approximant ' + approximant + ' not supported')
    return (tc - track_t, track_f)
511459
def _non_local_part(q0: ops.Qid,
                    q1: ops.Qid,
                    interaction_coefficients: Tuple[float, float, float],
                    atol: float = 1e-8):
    """"""

    x, y, z = interaction_coefficients

    return [
        _parity_interaction(q0, q1, x, atol),
        _parity_interaction(q0, q1, y, atol, ops.Z ** -0.5),
        _parity_interaction(q0, q1, z, atol, ops.Y ** 0.5)]
1343543
def evolve(self, rho: Density) -> Density:
        """"""
        qubits = rho.qubits
        results = [op.evolve(rho) for op in self.operators]
        tensors = [rho.tensor * w for rho, w in zip(results, self.weights)]
        tensor = reduce(add, tensors)
        return Density(tensor, qubits)
282436
def gamma(self, x, y, kwargs, diff=diff):
        
        f_xx, f_xy, f_yx, f_yy = self.hessian(x, y, kwargs, diff=diff)
        gamma1 = 1./2 * (f_xx - f_yy)
        gamma2 = f_xy
        return gamma1, gamma2
567168
def w_state(qubits: Union[int, Qubits]) -> State:
    """"""
    N, qubits = qubits_count_tuple(qubits)
    ket = np.zeros(shape=[2] * N)
    for n in range(N):
        idx = np.zeros(shape=N, dtype=int)
        idx[n] += 1
        ket[tuple(idx)] = 1 / sqrt(N)
    return State(ket, qubits)
1312500
def rho2sigma(self, rho0, Ra, Rs):
        
        return np.pi * rho0 * Ra * Rs / (Rs + Ra)
418535
def partial_trace(self, qubits: Qubits) -> 'Density':
        """"""
        vec = self.vec.partial_trace(qubits)
        return Density(vec.tensor, vec.qubits, self._memory)
863227
def getSignificance(wk1, wk2, nout, ofac):
  
  expy = exp(-wk2)
  effm = 2.0*(nout)/ofac
  sig = effm*expy
  ind = (sig > 0.01).nonzero()
  sig[ind] = 1.0-(1.0-expy[ind])**effm
  return sig
1345243
function q_f(df1, df2, f) {
        var cos2, sin2, prob, temp, i;

        if (f <= 0.0) { return 1.0; }
        if (df1 % 2 !== 0 && df2 % 2 === 0) { return 1.0 - q_f(df2, df1, 1.0 / f); }

        cos2 = 1.0 / (1.0 + df1 * f / df2);
        sin2 = 1.0 - cos2;

        if (df1 % 2 === 0) {
            prob = Math.pow(cos2, (df2 / 2.0));
            temp = prob;

            i = 2;
            while (i < df1) {
                temp = temp * ((df2 + i - 2) * sin2 / i);
                prob = prob + temp;
                i = i + 2;
            }

            return prob;
        }

        prob = Math.atan(Math.sqrt(df2 / (df1 * f)));
        temp = Math.sqrt(sin2 * cos2);
        i = 3;
        while (i <= df1) {
            prob = prob + temp;
            temp = temp * ((i - 1) * sin2 / i);
            i = i + 2.0;
        }

        temp = temp * df1;
        i = 3;
        while (i <= df2) {
            prob = prob - temp;
            temp = temp * ((df1 + i - 2) * cos2 / i);
            i = i + 2;
        }

        return prob * 2.0 / Math.PI;
    }
867246
def get_hessian(self):
        """"""
        force_const = self.fields.get("Cartesian Force Constants")
        if force_const is None:
            return None
        N = len(self.molecule.numbers)
        result = np.zeros((3*N, 3*N), float)
        counter = 0
        for row in range(3*N):
            result[row, :row+1] = force_const[counter:counter+row+1]
            result[:row+1, row] = force_const[counter:counter+row+1]
            counter += row + 1
        return result
358017
def rho_rec(self, g2):
        
        return (self.expnorm / np.sqrt(1 + self.gamma2 * g2) + (
            1 - self.expnorm) * np.exp(-0.25 * self.beta2 * g2))
1074106
def permute(self, qubits: Qubits) -> 'Density':
        """"""
        vec = self.vec.permute(qubits)
        return Density(vec.tensor, vec.qubits, self._memory)
214193
def c_rho0(self, rho0):
        
        if not hasattr(self, '_c_rho0_interp'):
            c_array = np.linspace(0.1, 10, 100)
            rho0_array = self.rho0_c(c_array)
            from scipy import interpolate
            self._c_rho0_interp = interpolate.InterpolatedUnivariateSpline(rho0_array, c_array, w=None, bbox=[None, None], k=3)
        return self._c_rho0_interp(rho0)
1118344
def _get_numdiff_hessian_element_func(objective_func, nmr_steps, step_ratio):
    """"""
    return SimpleCLFunction.from_string('''
        /**
         * Compute the Hessian using (possibly) multiple steps with various interpolations. 
         */ 
        void _numdiff_hessian_element(
                void* data, local mot_float_type* x_tmp, mot_float_type f_x_input,
                uint px, uint py, global float* initial_step, global double* derivative, 
                global double* error, local double* scratch){

            const uint nmr_steps = ''' + str(nmr_steps) + ''';
            uint nmr_steps_remaining = nmr_steps;

            local double* scratch_ind = scratch;
            local double* steps = scratch_ind;      scratch_ind += nmr_steps; 
            local double* errors = scratch_ind;     scratch_ind += nmr_steps - 1;
            local double* steps_tmp = scratch_ind;  scratch_ind += nmr_steps;
            
            if(get_local_id(0) == 0){
                for(int i = 0; i < nmr_steps - 1; i++){
                    errors[i] = 0;
                }
            }
            barrier(CLK_LOCAL_MEM_FENCE);
            
            _numdiff_hessian_steps(data, x_tmp, f_x_input, px, py, steps, initial_step);

            if(nmr_steps_remaining > 1){
                nmr_steps_remaining = _numdiff_hessian_richardson_extrapolation(steps); 
                barrier(CLK_LOCAL_MEM_FENCE);
            }
            
            if(nmr_steps_remaining >= 3){
                nmr_steps_remaining = _numdiff_wynn_extrapolation(steps, errors, nmr_steps_remaining);
                barrier(CLK_LOCAL_MEM_FENCE);                
            }
            
            if(nmr_steps_remaining > 1){
                _numdiff_find_best_step(steps, errors, steps_tmp, nmr_steps_remaining);
                barrier(CLK_LOCAL_MEM_FENCE);
            }
            
            if(get_local_id(0) == 0){
                *derivative = steps[0];
                *error = errors[0];
            }
        }
    ''', dependencies=[
        _get_numdiff_hessian_steps_func(objective_func, nmr_steps, step_ratio),
        _get_numdiff_hessian_richardson_extrapolation_func(nmr_steps, step_ratio),
        _get_numdiff_wynn_extrapolation_func(),
        _get_numdiff_find_best_step_func()
    ])