Dense [k-1]x4 transition frequency matrix
11995
def transitions(self):
        """"""
        if self._transitions is not None:
            return self._transitions
        transitions = self.array.astype(np.float)
        transitions /= transitions.sum(1)[:, np.newaxis]
        self._transitions = transitions
        return transitions
508213
def _fourier(self):
        

        freq_bin_upper = 2000
        freq_bin_lower = 40
        fs = self._metadata['fs']
        Y_transformed = {}
        for key in self.Y_dict.keys():
            # normalize by fs, bins have units strain/Hz
            fs = self._metadata["fs"]
            Y_transformed[key] = (1./fs)*np.fft.fft(self.Y_dict[key])[freq_bin_lower:freq_bin_upper]
        self.Y_transformed = Y_transformed
        # add in transform metadata
        self._metadata["dF"] = 1./self._metadata["T"]

        # because we are in the fourier domain, we will need the psd
        self.psd = load_psd()[freq_bin_lower:freq_bin_upper]
        dF = self._metadata['dF']
        self.sigma = convert_psd_to_sigma(self.psd, dF)
1070265
def _jordan_wigner_mode(n):
        
        a = []
        for i in range(n):
            xv = np.asarray([1] * i + [0] + [0] * (n - i - 1))
            xw = np.asarray([0] * i + [1] + [0] * (n - i - 1))
            yv = np.asarray([1] * i + [1] + [0] * (n - i - 1))
            yw = np.asarray([0] * i + [1] + [0] * (n - i - 1))
            a.append((Pauli(xv, xw), Pauli(yv, yw)))
        return a
389129
def stSpectralRollOff(X, c, fs):
    """"""
    totalEnergy = numpy.sum(X ** 2)
    fftLength = len(X)
    Thres = c*totalEnergy
    # Ffind the spectral rolloff as the frequency position 
    # where the respective spectral energy is equal to c*totalEnergy
    CumSum = numpy.cumsum(X ** 2) + eps
    [a, ] = numpy.nonzero(CumSum > Thres)
    if len(a) > 0:
        mC = numpy.float64(a[0]) / (float(fftLength))
    else:
        mC = 0.0
    return (mC)
884634
def s2dctmat(nfilt,ncep,freqstep):
    """"""
    melcos = numpy.empty((ncep, nfilt), 'double')
    for i in range(0,ncep):
        freq = numpy.pi * float(i) / nfilt
        melcos[i] = numpy.cos(freq * numpy.arange(0.5, float(nfilt)+0.5, 1.0, 'double'))
    melcos[:,0] = melcos[:,0] * 0.5
    return melcos
779956
def morph(clm1, clm2, t, lmax):
    """"""
    clm = (1 - t) * clm1 + t * clm2
    grid_reco = clm.expand(lmax=lmax)  # cut "high frequency" components
    agrid_reco = grid_reco.to_array()
    pts = []
    for i, longs in enumerate(agrid_reco):
        ilat = grid_reco.lats()[i]
        for j, value in enumerate(longs):
            ilong = grid_reco.lons()[j]
            th = (90 - ilat) / 57.3
            ph = ilong / 57.3
            r = value + rbias
            p = np.array([sin(th) * cos(ph), sin(th) * sin(ph), cos(th)]) * r
            pts.append(p)
    return pts
1240381
def q12d(grid, spacing=None, E=1e5, nu=0.3, dirichlet_boundary=True,
         format=None):
    
    X, Y = tuple(grid)

    if X < 1 or Y < 1:
        raise ValueError('invalid grid shape')

    if dirichlet_boundary:
        X += 1
        Y += 1

    pts = np.mgrid[0:X+1, 0:Y+1]
    pts = np.hstack((pts[0].T.reshape(-1, 1) - X / 2.0,
                     pts[1].T.reshape(-1, 1) - Y / 2.0))

    if spacing is None:
        DX, DY = 1, 1
    else:
        DX, DY = tuple(spacing)
        pts *= [DX, DY]

    # compute local stiffness matrix
    lame = E * nu / ((1 + nu) * (1 - 2*nu))  # Lame's first parameter
    mu = E / (2 + 2*nu)                   # shear modulus

    vertices = np.array([[0, 0], [DX, 0], [DX, DY], [0, DY]])
    K = q12d_local(vertices, lame, mu)

    nodes = np.arange((X+1)*(Y+1)).reshape(X+1, Y+1)
    LL = nodes[:-1, :-1]
    Id = (2*LL).repeat(K.size).reshape(-1, 8, 8)
    J = Id.copy()
    Id += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1))
    J += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1)).T
    V = np.tile(K, (X*Y, 1))

    Id = np.ravel(Id)
    J = np.ravel(J)
    V = np.ravel(V)

    # sum duplicates
    A = coo_matrix((V, (Id, J)), shape=(pts.size, pts.size)).tocsr()
    A = A.tobsr(blocksize=(2, 2))

    del Id, J, V, LL, nodes

    B = np.zeros((2 * (X+1)*(Y+1), 3))
    B[0::2, 0] = 1
    B[1::2, 1] = 1
    B[0::2, 2] = -pts[:, 1]
    B[1::2, 2] = pts[:, 0]

    if dirichlet_boundary:
        mask = np.zeros((X+1, Y+1), dtype='bool')
        mask[1:-1, 1:-1] = True
        mask = np.ravel(mask)
        data = np.zeros(((X-1)*(Y-1), 2, 2))
        data[:, 0, 0] = 1
        data[:, 1, 1] = 1
        indices = np.arange((X-1)*(Y-1))
        indptr = np.concatenate((np.array([0]), np.cumsum(mask)))
        P = bsr_matrix((data, indices, indptr),
                       shape=(2*(X+1)*(Y+1), 2*(X-1)*(Y-1)))
        Pt = P.T
        A = P.T * A * P

        B = Pt * B

    return A.asformat(format), B
1263509
def K(self, parm):
        
        return Periodic_K_matrix(self.X, parm) + np.identity(self.X.shape[0])*(10**-10)
939687
def G(self):
        
        if self.c is None:
            raise RuntimeError('Inverse cross spectral density requires '
                               'invertible noise covariance matrix c.')
        A = self.A()
        # TODO: can we do that more efficiently?
        G = np.einsum('ji..., jk... ->ik...', A.conj(), self.Cinv())
        G = np.einsum('ij..., jk... ->ik...', G, A)
        return G
32165
def xstep(self):
        """

        self.YU[:] = self.Y - self.U
        YUf = sl.rfftn(self.YU, None, self.cri.axisN)
        YUf0 = self.block_sep0(YUf)
        YUf1 = self.block_sep1(YUf)

        b = self.rho * np.sum(np.conj(self.GDf) * YUf1, axis=-1)
        if self.cri.Cd > 1:
            b = np.sum(b, axis=self.cri.axisC, keepdims=True)
        b += self.DSf + self.rho*YUf0

        # Concatenate multiple GDf components on axisC. For
        # single-channel signals, and multi-channel signals with a
        # single-channel dictionary, we end up with sl.solvemdbi_ism
        # solving a linear system of rank dimN+1 (corresponding to the
        # dictionary and a gradient operator per spatial dimension) plus
        # an identity. For multi-channel signals with a multi-channel
        # dictionary, we end up with sl.solvemdbi_ism solving a linear
        # system of rank C.d (dimN+1) (corresponding to the dictionary
        # and a gradient operator per spatial dimension for each
        # channel) plus an identity.

        # The structure of the linear system to be solved depends on the
        # number of channels in the signal and dictionary. Both branches are
        # the same in the single-channel signal case (the choice of handling
        # it via the 'else' branch is somewhat arbitrary).
        if self.cri.C > 1 and self.cri.Cd == 1:
            # Concatenate multiple GDf components on the final axis
            # of GDf (that indexes the number of gradient operators). For
            # multi-channel signals with a single-channel dictionary,
            # sl.solvemdbi_ism has to solve a linear system of rank dimN+1
            # (corresponding to the dictionary and a gradient operator per
            # spatial dimension)
            DfGDf = np.concatenate(
                [self.Df[..., np.newaxis],] +
                [np.sqrt(self.rho)*self.GDf[..., k, np.newaxis] for k
                 in range(self.GDf.shape[-1])], axis=-1)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b[..., np.newaxis],
                                          self.cri.axisM, -1)[..., 0]
        else:
            # Concatenate multiple GDf components on axisC. For multi-channel
            # signals with a multi-channel dictionary, sl.solvemdbi_ism has
            # to solve a linear system of rank C.d (dimN+1) (corresponding to
            # the dictionary and a gradient operator per spatial dimension
            # for each channel) plus an identity.
            DfGDf = np.concatenate(
                [self.Df,] + [np.sqrt(self.rho)*self.GDf[..., k] for k
                              in range(self.GDf.shape[-1])],
                axis=self.cri.axisC)
            self.Xf[:] = sl.solvemdbi_ism(DfGDf, self.rho, b, self.cri.axisM,
                                          self.cri.axisC)

        self.X = sl.irfftn(self.Xf, self.cri.Nv, self.cri.axisN)

        if self.opt['LinSolveCheck']:
            if self.cri.C > 1 and self.cri.Cd == 1:
                Dop = lambda x: sl.inner(DfGDf, x[..., np.newaxis],
                                         axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x, axis=-1)
                ax = DHop(Dop(self.Xf))[..., 0] + self.rho*self.Xf
            else:
                Dop = lambda x: sl.inner(DfGDf, x, axis=self.cri.axisM)
                DHop = lambda x: sl.inner(np.conj(DfGDf), x,
                                          axis=self.cri.axisC)
                ax = DHop(Dop(self.Xf)) + self.rho*self.Xf
            self.xrrs = sl.rrs(ax, b)
        else:
            self.xrrs = None
244228
def viterbi_tags(self,
                     logits: torch.Tensor,
                     mask: torch.Tensor) -> List[Tuple[List[int], float]]:
        
        _, max_seq_length, num_tags = logits.size()

        # Get the tensors out of the variables
        logits, mask = logits.data, mask.data

        # Augment transitions matrix with start and end transitions
        start_tag = num_tags
        end_tag = num_tags + 1
        transitions = torch.Tensor(num_tags + 2, num_tags + 2).fill_(-10000.)

        # Apply transition constraints
        constrained_transitions = (
                self.transitions * self._constraint_mask[:num_tags, :num_tags] +
                -10000.0 * (1 - self._constraint_mask[:num_tags, :num_tags])
        )
        transitions[:num_tags, :num_tags] = constrained_transitions.data

        if self.include_start_end_transitions:
            transitions[start_tag, :num_tags] = (
                    self.start_transitions.detach() * self._constraint_mask[start_tag, :num_tags].data +
                    -10000.0 * (1 - self._constraint_mask[start_tag, :num_tags].detach())
            )
            transitions[:num_tags, end_tag] = (
                    self.end_transitions.detach() * self._constraint_mask[:num_tags, end_tag].data +
                    -10000.0 * (1 - self._constraint_mask[:num_tags, end_tag].detach())
            )
        else:
            transitions[start_tag, :num_tags] = (-10000.0 *
                                                 (1 - self._constraint_mask[start_tag, :num_tags].detach()))
            transitions[:num_tags, end_tag] = -10000.0 * (1 - self._constraint_mask[:num_tags, end_tag].detach())

        best_paths = []
        # Pad the max sequence length by 2 to account for start_tag + end_tag.
        tag_sequence = torch.Tensor(max_seq_length + 2, num_tags + 2)

        for prediction, prediction_mask in zip(logits, mask):
            sequence_length = torch.sum(prediction_mask)

            # Start with everything totally unlikely
            tag_sequence.fill_(-10000.)
            # At timestep 0 we must have the START_TAG
            tag_sequence[0, start_tag] = 0.
            # At steps 1, ..., sequence_length we just use the incoming prediction
            tag_sequence[1:(sequence_length + 1), :num_tags] = prediction[:sequence_length]
            # And at the last timestep we must have the END_TAG
            tag_sequence[sequence_length + 1, end_tag] = 0.

            # We pass the tags and the transitions to ``viterbi_decode``.
            viterbi_path, viterbi_score = util.viterbi_decode(tag_sequence[:(sequence_length + 2)], transitions)
            # Get rid of START and END sentinels and append.
            viterbi_path = viterbi_path[1:-1]
            best_paths.append((viterbi_path, viterbi_score.item()))

        return best_paths
1363918
def _multitaper_spectrum(self, clm, k, convention='power', unit='per_l',
                             clat=None, clon=None, coord_degrees=True,
                             lmax=None, taper_wt=None):
        
        if lmax is None:
            lmax = clm.lmax

        if (clat is not None and clon is not None and clat == self.clat and
                clon == self.clon and coord_degrees is self.coord_degrees and
                k <= self.nwinrot):
            # use the already stored coeffs
            pass
        elif (clat is None and clon is None) and \
                (self.clat is not None and self.clon is not None and
                 k <= self.nwinrot):
            # use the already stored coeffs
            pass
        else:
            if clat is None:
                clat = self.clat
            if clon is None:
                clon = self.clon
            if (clat is None and clon is not None) or \
                    (clat is not None and clon is None):
                raise ValueError('clat and clon must both be input. ' +
                                 'clat = {:s}, clon = {:s}'
                                 .format(repr(clat), repr(clon)))
            if clat is None and clon is None:
                self.rotate(clat=90., clon=0., coord_degrees=True, nwinrot=k)
            else:
                self.rotate(clat=clat, clon=clon, coord_degrees=coord_degrees,
                            nwinrot=k)

        sh = clm.to_array(normalization='4pi', csphase=1, lmax=lmax)

        if taper_wt is None:
            mtse, sd = _shtools.SHMultiTaperMaskSE(sh, self.coeffs,
                                                   lmax=lmax, k=k)
        else:
            mtse, sd = _shtools.SHMultiTaperMaskSE(sh, self.coeffs, lmax=lmax,
                                                   k=k, taper_wt=taper_wt)

        if (unit == 'per_l'):
            pass
        elif (unit == 'per_lm'):
            degree_l = _np.arange(len(mtse))
            mtse /= (2.0 * degree_l + 1.0)
            sd /= (2.0 * degree_l + 1.0)
        else:
            raise ValueError(
                "unit must be 'per_l' or 'per_lm'." +
                "Input value was {:s}".format(repr(unit)))

        if (convention == 'power'):
            return mtse, sd
        elif (convention == 'energy'):
            return mtse * 4.0 * _np.pi, sd * 4.0 * _np.pi
        else:
            raise ValueError(
                "convention must be 'power' or 'energy'." +
                "Input value was {:s}".format(repr(convention)))
5354
def toeplitz(vect):
  
  return [[vect[abs(i-j)] for i in xrange(len(vect))]
                          for j in xrange(len(vect))]
610388
def asML(self):
        
        return newlinalg.DenseMatrix(self.numRows, self.numCols, self.values, self.isTransposed)
556852
def _opbend_transform_mean(rs, fn_low, deriv=0):
    
    v = 0.0
    d = np.zeros((4,3), float)
    dd = np.zeros((4,3,4,3), float)
    #loop over the 3 cyclic permutations
    for p in np.array([[0,1,2], [2,0,1], [1,2,0]]):
        opbend = _opbend_transform([rs[p[0]], rs[p[1]], rs[p[2]], rs[3]], fn_low, deriv)
        v += opbend[0]/3
        index0 = np.where(p==0)[0][0] #index0 is the index of the 0th atom (rs[0])
        index1 = np.where(p==1)[0][0]
        index2 = np.where(p==2)[0][0]
        index3 = 3
        if deriv>0:
            d[0] += opbend[1][index0]/3
            d[1] += opbend[1][index1]/3
            d[2] += opbend[1][index2]/3
            d[3] += opbend[1][index3]/3
        if deriv>1:
            dd[0, :, 0, :] += opbend[2][index0, :, index0, :]/3
            dd[0, :, 1, :] += opbend[2][index0, :, index1, :]/3
            dd[0, :, 2, :] += opbend[2][index0, :, index2, :]/3
            dd[0, :, 3, :] += opbend[2][index0, :, index3, :]/3

            dd[1, :, 0, :] += opbend[2][index1, :, index0, :]/3
            dd[1, :, 1, :] += opbend[2][index1, :, index1, :]/3
            dd[1, :, 2, :] += opbend[2][index1, :, index2, :]/3
            dd[1, :, 3, :] += opbend[2][index1, :, index3, :]/3

            dd[2, :, 0, :] += opbend[2][index2, :, index0, :]/3
            dd[2, :, 1, :] += opbend[2][index2, :, index1, :]/3
            dd[2, :, 2, :] += opbend[2][index2, :, index2, :]/3
            dd[2, :, 3, :] += opbend[2][index2, :, index3, :]/3

            dd[3, :, 0, :] += opbend[2][index3, :, index0, :]/3
            dd[3, :, 1, :] += opbend[2][index3, :, index1, :]/3
            dd[3, :, 2, :] += opbend[2][index3, :, index2, :]/3
            dd[3, :, 3, :] += opbend[2][index3, :, index3, :]/3
    if deriv==0:
        return v,
    elif deriv==1:
        return v, d
    elif deriv==2:
        return v, d, dd
    else:
        raise ValueError("deriv must be 0, 1 or 2.")
912943
def xmoe_tr_dense_2k():
  
  hparams = mtf_transformer2.mtf_bitransformer_base()
  hparams.encoder_layers = ["self_att", "drd"] * 4
  hparams.decoder_layers = ["self_att", "enc_att", "drd"] * 4
  hparams.batch_size = 64
  hparams.shared_embedding_and_softmax_weights = True
  hparams.mesh_shape = "batch:8"
  return hparams
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;
    }
  }