A function to order the phonon eigenvectors taken from phonopy
970897
def estimate_band_connection(prev_eigvecs, eigvecs, prev_band_order):
    
    metric = np.abs(np.dot(prev_eigvecs.conjugate().T, eigvecs))
    connection_order = []
    for overlaps in metric:
        maxval = 0
        for i in reversed(range(len(metric))):
            val = overlaps[i]
            if i in connection_order:
                continue
            if val > maxval:
                maxval = val
                maxindex = i
        connection_order.append(maxindex)

    band_order = [connection_order[x] for x in prev_band_order]

    return band_order
781914
def get_bandstructure(self):
        

        return LobsterBandStructureSymmLine(kpoints=self.kpoints_array, eigenvals=self.eigenvals, lattice=self.lattice,
                                            efermi=self.efermi, labels_dict=self.label_dict,
                                            structure=self.structure,
                                            projections=self.p_eigenvals)
489107
def blt(f: List[SYM], x: List[SYM]) -> Dict[str, Any]:
    
    J = ca.jacobian(f, x)
    nblock, rowperm, colperm, rowblock, colblock, coarserow, coarsecol = J.sparsity().btf()
    return {
        'J': J,
        'nblock': nblock,
        'rowperm': rowperm,
        'colperm': colperm,
        'rowblock': rowblock,
        'colblock': colblock,
        'coarserow': coarserow,
        'coarsecol': coarsecol
    }
360338
def energy_coefficients(m1, m2, s1z=0, s2z=0, phase_order=-1, spin_order=-1):
    
    implemented_phase_order = 7
    implemented_spin_order = 7
    if phase_order > implemented_phase_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif phase_order == -1:
        phase_order = implemented_phase_order

    if spin_order > implemented_spin_order:
        raise ValueError("pN coeffiecients of that order have not been implemented")
    elif spin_order == -1:
        spin_order = implemented_spin_order

    qmdef1 = 1.0
    qmdef2 = 1.0

    M = m1 + m2
    dm = (m1-m2)/M
    m1M = m1 / M
    m2M = m2 / M

    s1z = s1z * m1M * m1M
    s2z = s2z * m2M * m2M

    _, eta = mass1_mass2_to_mchirp_eta(m1, m2)

    ecof = numpy.zeros(phase_order+1)
    # Orbital terms
    if phase_order >= 0:
        ecof[0] = 1.0
    if phase_order >= 1:
        ecof[1] = 0
    if phase_order >= 2:
        ecof[2] = -(1.0/12.0) * (9.0 + eta)
    if phase_order >= 3:
        ecof[3] = 0
    if phase_order >= 4:
        ecof[4] = (-81.0 + 57.0*eta - eta*eta) / 24.0
    if phase_order >= 5:
        ecof[5] = 0
    if phase_order >= 6:
        ecof[6] = - 675.0/64.0 + ( 34445.0/576.0    \
              - 205.0/96.0 * lal.PI * lal.PI ) * eta  \
              - (155.0/96.0) *eta * eta - 35.0/5184.0 * eta * eta
    # Spin terms

    ESO15s1 = 8.0/3.0 + 2.0*m2/m1
    ESO15s2 = 8.0/3.0 + 2.0*m1/m2

    ESS2 = 1.0 / eta
    EQM2s1 = qmdef1/2.0/m1M/m1M
    EQM2s1L = -qmdef1*3.0/2.0/m1M/m1M
    #EQM2s2 = qmdef2/2.0/m2M/m2M
    EQM2s2L = -qmdef2*3.0/2.0/m2M/m2M

    ESO25s1 = 11.0 - 61.0*eta/9.0 + (dm/m1M) * (-3.0 + 10.*eta/3.0)
    ESO25s2 = 11.0 - 61.0*eta/9.0 + (dm/m2M) * (3.0 - 10.*eta/3.0)

    ESO35s1 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 + (dm/m1M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)
    ESO35s2 = 135.0/4.0 - 367.0*eta/4.0 + 29.0*eta*eta/12.0 - (dm/m2M) * (-27.0/4.0 + 39.0*eta - 5.0*eta*eta/4.0)

    if spin_order >=3:
        ecof[3] += ESO15s1 * s1z + ESO15s2 * s2z
    if spin_order >=4:
        ecof[4] += ESS2 * (s1z*s2z - 3.0*s1z*s2z)
        ecof[4] += EQM2s1*s1z*s1z + EQM2s1*s2z*s2z + EQM2s1L*s1z*s1z + EQM2s2L*s2z*s2z
    if spin_order >=5:
        ecof[5] = ESO25s1*s1z + ESO25s2*s2z
    if spin_order >=7:
        ecof[7] += ESO35s1*s1z + ESO35s2*s2z

    return ecof
767896
def eigenvalue_band_properties(self):
        
        vbm = -float("inf")
        vbm_kpoint = None
        cbm = float("inf")
        cbm_kpoint = None
        for spin, d in self.eigenvalues.items():
            for k, val in enumerate(d):
                for (eigenval, occu) in val:
                    if occu > self.occu_tol and eigenval > vbm:
                        vbm = eigenval
                        vbm_kpoint = k
                    elif occu <= self.occu_tol and eigenval < cbm:
                        cbm = eigenval
                        cbm_kpoint = k
        return max(cbm - vbm, 0), cbm, vbm, vbm_kpoint == cbm_kpoint
1034572
def lowpass(ts, cutoff_hz, order=3):
    """"""
    orig_ndim = ts.ndim
    if ts.ndim is 1:
        ts = ts[:, np.newaxis]
    channels = ts.shape[1]
    fs = (len(ts) - 1.0) / (ts.tspan[-1] - ts.tspan[0])
    nyq = 0.5 * fs
    cutoff = cutoff_hz/nyq
    b, a = signal.butter(order, cutoff, btype='low')
    if not np.all(np.abs(np.roots(a)) < 1.0):
        raise ValueError('Filter will not be stable with these values.')
    dtype = ts.dtype
    output = np.zeros((len(ts), channels), dtype)
    for i in range(channels):
        output[:, i] = signal.filtfilt(b, a, ts[:, i])
    if orig_ndim is 1:
        output = output[:, 0]
    return Timeseries(output, ts.tspan, labels=ts.labels)
1220853
def laplacian_reordering(G):
  ''''''
  L = G.laplacian()
  vals, vecs = np.linalg.eigh(L)
  min_positive_idx = np.argmax(vals == vals[vals>0].min())
  vec = vecs[:, min_positive_idx]
  return permute_graph(G, np.argsort(vec))
543468
def get_jac(self):
        
        if self._jac is True:
            if self.band is None:
                f = self.be.Matrix(self.nf, 1, self.exprs)
                _x = self.be.Matrix(self.nx, 1, self.x)
                return f.jacobian(_x)
            else:
                # Banded
                return self.be.Matrix(banded_jacobian(
                    self.exprs, self.x, *self.band))
        elif self._jac is False:
            return False
        else:
            return self._jac
169894
def relaxation_operators(p):
    
    k0 = np.array([[1.0, 0.0], [0.0, np.sqrt(1 - p)]])
    k1 = np.array([[0.0, np.sqrt(p)], [0.0, 0.0]])
    return k0, k1
644642
def compute_vis_noise(msname, sefd, spw_id=0):
    """"""
    from pyrap.tables import table

    tab = table(msname)
    spwtab = table(msname + "/SPECTRAL_WINDOW")

    freq0 = spwtab.getcol("CHAN_FREQ")[spw_id, 0]
    wavelength = 300e+6/freq0
    bw = spwtab.getcol("CHAN_WIDTH")[spw_id, 0]
    dt = tab.getcol("EXPOSURE", 0, 1)[0]
    dtf = (tab.getcol("TIME", tab.nrows()-1, 1)-tab.getcol("TIME", 0, 1))[0]

    # close tables properly, else the calls below will hang waiting for a lock...
    tab.close()
    spwtab.close()

    print(">>> %s freq %.2f MHz (lambda=%.2fm), bandwidth %.2g kHz, %.2fs integrations, %.2fh synthesis"%(msname, freq0*1e-6, wavelength, bw*1e-3, dt, dtf/3600))
    noise = sefd/math.sqrt(abs(2*bw*dt))
    print(">>> SEFD of %.2f Jy gives per-visibility noise of %.2f mJy"%(sefd, noise*1000))

    return noise
423540
function hm(n, v, incv, offv, tau, b, offb) {
    // Apply householder transformation b = b - tau v v' b
    // LAPACK: DLARF

    if (tau === 0.0) {
      return;
    }

    // compute v' b
    var d = b[offb + 0];
    for (var i = 1; i < n; ++i) {
      d += b[offb + i] * v[offv + incv*i];
    }

    // compute b = b - tau v v' b
    b[offb + 0] -= tau * d;
    for (var i = 1; i < n; ++i) {
      b[offb + i] -= tau * v[offv + incv*i] * d;
    }
  }
1131690
def get_principal_axes(self):
        
        # Perform eigendecomposition - returns in order P, B, T
        _ = self.eigendecompose(normalise=True)
        principal_axes = GCMTPrincipalAxes()
        # Eigenvalues
        principal_axes.p_axis = {'eigenvalue': self.eigenvalues[0]}
        principal_axes.b_axis = {'eigenvalue': self.eigenvalues[1]}
        principal_axes.t_axis = {'eigenvalue': self.eigenvalues[2]}
        # Eigen vectors
        # 1) P axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 0], True)
        principal_axes.p_axis['azimuth'] = azim
        principal_axes.p_axis['plunge'] = plun
        # 2) B axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 1], True)
        principal_axes.b_axis['azimuth'] = azim
        principal_axes.b_axis['plunge'] = plun
        # 3) T axis
        azim, plun = utils.get_azimuth_plunge(self.eigenvectors[:, 2], True)
        principal_axes.t_axis['azimuth'] = azim
        principal_axes.t_axis['plunge'] = plun
        return principal_axes
1320422
function computeHOAMatrices(matrix) {
  // We start by computing the 2nd-order matrix from the 1st-order matrix.
  for (let i = 2; i <= matrix.length; i++) {
    computeBandRotation(matrix, i);
  }
}
1238165
def trt_pmf(matrices):
    
    ntrts, nmags, ndists, nlons, nlats, neps = matrices.shape
    pmf = numpy.zeros(ntrts)
    for t in range(ntrts):
        pmf[t] = 1. - numpy.prod(
            [1. - matrices[t, i, j, k, l, m]
             for i in range(nmags)
             for j in range(ndists)
             for k in range(nlons)
             for l in range(nlats)
             for m in range(neps)])
    return pmf
1042106
def get_pmg_structure(phonopy_structure):
    

    lattice = phonopy_structure.get_cell()
    frac_coords = phonopy_structure.get_scaled_positions()
    symbols = phonopy_structure.get_chemical_symbols()
    masses = phonopy_structure.get_masses()
    mms = phonopy_structure.get_magnetic_moments()
    mms = mms or [0] * len(symbols)

    return Structure(lattice, symbols, frac_coords,
                     site_properties={"phonopy_masses": masses,
                                      "magnetic_moments": mms})
113099
def eigenvalues_samples(self):
         """
        res = np.empty((self.nsamples, self.nstates), dtype=config.dtype)
        for i in range(self.nsamples):
            res[i, :] = self._sampled_hmms[i].eigenvalues
        return res
653895
def ThermalSpectrum(self):
        
        try:
            # delegate to subclass.
            thom = _ThermalObservationMode(self._obsmode)
            self.pixscale = thom.pixscale
            return thom._getSpectrum()
        except IndexError:   # graph table is broken.
            raise IndexError("Cannot calculate thermal spectrum; graphtable may be broken")