Python numpy.diag_indices() Examples

def compute_outflow_pdf(self):
        """
        Lifetime model. The method compute outflow_pdf returns an array year-by-cohort of the probability of a item added to stock in year m (aka cohort m) leaves in in year n. This value equals pdf(n,m).
        This is the only method for the inflow-driven model where the lifetime distribution directly enters the computation. All other stock variables are determined by mass balance.
        The shape of the output pdf array is NoofYears * NoofYears, but the meaning is years by age-cohorts.
        The method does nothing if the pdf alreay exists.
        """
        if self.pdf is None:
            self.compute_sf() # computation of pdfs moved to this method: compute survival functions sf first, then calculate pdfs from sf.
            self.pdf   = np.zeros((len(self.t), len(self.t)))
            self.pdf[np.diag_indices(len(self.t))] = np.ones(len(self.t)) - self.sf.diagonal(0)
            for m in range(0,len(self.t)):
                self.pdf[np.arange(m+1,len(self.t)),m] = -1 * np.diff(self.sf[np.arange(m,len(self.t)),m])            
            return self.pdf
        else:
            # pdf already exists
            return self.pdf 

def diag_indices_from(arr):
    """
    Return the indices to access the main diagonal of an n-dimensional array.
    See `diag_indices` for full details.

    Args:
        arr (cupy.ndarray): At least 2-D.

     .. seealso:: :func:`numpy.diag_indices_from`

    """
    if not isinstance(arr, cupy.ndarray):
        raise TypeError("Argument must be cupy.ndarray")

    if not arr.ndim >= 2:
        raise ValueError("input array must be at least 2-d")
    # For more than d=2, the strided formula is only valid for arrays with
    # all dimensions equal, so we check first.
    if not cupy.all(cupy.diff(arr.shape) == 0):
        raise ValueError("All dimensions of input must be of equal length")

    return diag_indices(arr.shape[0], arr.ndim) 

def test_field_gradients(self):
        mol = gto.M(atom='H1 0.5 -0.6 0.4; H2 -0.5, 0.4, -0.3; H -0.4 -0.3 0.5; H 0.3 0.5 -0.6',
                    unit='B',
                    basis={'H': [[0,[2., 1]]], 'H1':[[1,[.5, 1]]], 'H2':[[1,[1,1]]]})
        grids = dft.gen_grid.Grids(mol)
        grids.build()
        ao = mol.eval_gto('GTOval', grids.coords)
        r0 = mol.atom_coord(0)
        dr = grids.coords - r0
        dd = numpy.linalg.norm(dr, axis=1)
        rr = 3 * numpy.einsum('ix,iy->ixy', dr, dr)
        for i in range(3):
            rr[:,i,i] -= dd**2
        h1ref = lib.einsum('i,ixy,ip,iq->xypq', grids.weights/dd**5, rr, ao, ao)

        h1ao = rhf_ssc._get_integrals_fcsd(mol, 0)
        h1ao[numpy.diag_indices(3)] += rhf_ssc._get_integrals_fc(mol, 0)
        self.assertAlmostEqual(abs(h1ref - h1ao).max(), 0, 5) 

def energy_nuc(mol):
    atom_charges = mol.atom_charges()
    atom_coords = mol.atom_coords()

    distances = numpy.linalg.norm(atom_coords[:,None,:] - atom_coords, axis=2)
    distances_in_AA = distances * lib.param.BOHR
    # numerically exclude atomic self-interaction terms
    distances_in_AA[numpy.diag_indices(mol.natm)] = 1e60

    # one atom is H, another atom is N or O
    where_NO = (atom_charges == 7) | (atom_charges == 8)
    mask = (atom_charges[:,None] == 1) & where_NO
    mask = mask | mask.T
    scale = alpha = _get_alpha(atom_charges[:,None], atom_charges)
    scale[mask] *= numpy.exp(-distances_in_AA[mask])
    scale[~mask] = numpy.exp(-alpha[~mask] * distances_in_AA[~mask])

    gamma = _get_gamma(mol)
    z_eff = mopac_param.CORE[atom_charges]
    e_nuc = .5 * numpy.einsum('i,ij,j->', z_eff, gamma, z_eff)
    e_nuc += .5 * numpy.einsum('i,j,ij,ij->', z_eff, z_eff, scale,
                               mopac_param.E2/distances_in_AA - gamma)
    return e_nuc 

def test_square(self):
        m, n = (3, 3)

        # Obtain a random matrix of orthonormal rows
        Q = random_unitary_matrix(n)
        Q = Q[:m, :]
        Q = Q[:m, :]

        # Get Givens decomposition of Q
        givens_rotations, V, diagonal = givens_decomposition(Q)

        # There should be no Givens rotations
        self.assertEqual(givens_rotations, list())

        # Compute V * Q * U^\dagger
        W = V.dot(Q)

        # Construct the diagonal matrix
        D = numpy.zeros((m, n), dtype=complex)
        D[numpy.diag_indices(m)] = diagonal

        # Assert that W and D are the same
        for i in range(m):
            for j in range(n):
                self.assertAlmostEqual(D[i, j], W[i, j]) 

def test_antidiagonal(self):
        m, n = (3, 3)
        Q = numpy.zeros((m, n), dtype=complex)
        Q[0, 2] = 1.
        Q[1, 1] = 1.
        Q[2, 0] = 1.
        givens_rotations, V, diagonal = givens_decomposition(Q)

        # There should be no Givens rotations
        self.assertEqual(givens_rotations, list())

        # VQ should equal the diagonal
        VQ = V.dot(Q)
        D = numpy.zeros((m, n), dtype=complex)
        D[numpy.diag_indices(m)] = diagonal
        for i in range(n):
            for j in range(n):
                self.assertAlmostEqual(VQ[i, j], D[i, j]) 

def dist_diff(self, other=None):
        r"""
        Calculate distances and vectors between atoms.

        Args:
            other (:class:`~berny.Geometry`): calculate distances between two
                geometries if given or within a geometry if not

        Returns:
            :math:`R_{ij}:=|\mathbf R_i-\mathbf R_j|` and
            :math:`R_{ij\alpha}:=(\mathbf R_i)_\alpha-(\mathbf R_j)_\alpha`.
        """
        if other is None:
            other = self
        diff = self.coords[:, None, :] - other.coords[None, :, :]
        dist = np.sqrt(np.sum(diff ** 2, 2))
        dist[np.diag_indices(len(self))] = np.inf
        return dist, diff 

def __init__(self, one_body, two_body, constant=0.):
        if two_body.dtype != numpy.float:
            raise ValueError('Two-body tensor has invalid dtype. Expected {} '
                             'but was {}'.format(numpy.float, two_body.dtype))
        if not numpy.allclose(two_body, two_body.T):
            raise ValueError('Two-body tensor must be symmetric.')
        if not numpy.allclose(one_body, one_body.T.conj()):
            raise ValueError('One-body tensor must be Hermitian.')

        # Move the diagonal of two_body to one_body
        diag_indices = numpy.diag_indices(one_body.shape[0])
        one_body[diag_indices] += two_body[diag_indices]
        numpy.fill_diagonal(two_body, 0.)

        self.one_body = one_body
        self.two_body = two_body
        self.constant = constant 

def __init__(self, features, using_cache=False, identity_init=True, eps=1e-3):
        super().__init__(features, using_cache)

        self.eps = eps

        self.lower_indices = np.tril_indices(features, k=-1)
        self.upper_indices = np.triu_indices(features, k=1)
        self.diag_indices = np.diag_indices(features)

        n_triangular_entries = ((features - 1) * features) // 2

        self.lower_entries = nn.Parameter(torch.zeros(n_triangular_entries))
        self.upper_entries = nn.Parameter(torch.zeros(n_triangular_entries))
        self.unconstrained_upper_diag = nn.Parameter(torch.zeros(features))

        self._initialize(identity_init) 

def _get_knmat(exog, xcov, sl):
    # Utility function, see equation 2.2 of Barber & Candes.

    nobs, nvar = exog.shape

    ash = np.linalg.inv(xcov)
    ash *= -np.outer(sl, sl)
    i, j = np.diag_indices(nvar)
    ash[i, j] += 2 * sl

    umat = np.random.normal(size=(nobs, nvar))
    u, _ = np.linalg.qr(exog)
    umat -= np.dot(u, np.dot(u.T, umat))
    umat, _ = np.linalg.qr(umat)

    ashr, xc, _ = np.linalg.svd(ash, 0)
    ashr *= np.sqrt(xc)
    ashr = ashr.T

    ex = (sl[:, None] * np.linalg.solve(xcov, exog.T)).T
    exogn = exog - ex + np.dot(umat, ashr)

    return exogn 

def _get_knmat(exog, xcov, sl):
    # Utility function, see equation 2.2 of Barber & Candes.

    nobs, nvar = exog.shape

    ash = np.linalg.inv(xcov)
    ash *= -np.outer(sl, sl)
    i, j = np.diag_indices(nvar)
    ash[i, j] += 2 * sl

    umat = np.random.normal(size=(nobs, nvar))
    u, _ = np.linalg.qr(exog)
    umat -= np.dot(u, np.dot(u.T, umat))
    umat, _ = np.linalg.qr(umat)

    ashr, xc, _ = np.linalg.svd(ash, 0)
    ashr *= np.sqrt(xc)
    ashr = ashr.T

    ex = (sl[:, None] * np.linalg.solve(xcov, exog.T)).T
    exogn = exog - ex + np.dot(umat, ashr)

    return exogn 

def _initialize_error_cov_diagonal(self, scalar=False):
        # Initialize the parameters
        self.parameters['error_cov'] = 1 if scalar else self.k_endog

        # Setup fixed components of state space matrices

        # Setup indices of state space matrices
        k_endog = self.k_endog
        k_factors = self.k_factors
        idx = np.diag_indices(k_endog)
        if self.error_order > 0:
            matrix = 'state_cov'
            idx = (idx[0] + k_factors, idx[1] + k_factors)
        else:
            matrix = 'obs_cov'
        self._idx_error_cov = (matrix,) + idx 

def test_get_data_dst_nodata(rast, dst_nodata, dst_arr):
    fp = rast.fp.dilate(1)
    inner_slice = rast.fp.slice_in(fp)
    rast.set_data(dst_arr, channels=None)

    arr = rast.get_data(band=[-1], dst_nodata=dst_nodata, fp=fp)

    arr2 = rast.get_data(band=[-1], dst_nodata=dst_nodata, fp=fp.erode(1))
    assert np.all(arr2 == arr[inner_slice])

    if rast.nodata is not None:
        assert np.all(arr[inner_slice][np.diag_indices(dst_arr.shape[0])] == dst_nodata)
        arr[inner_slice][np.diag_indices(dst_arr.shape[0])] = rast.nodata
    assert np.all(arr[inner_slice] == dst_arr)

    outer_mask = np.ones(fp.shape, bool)
    outer_mask[inner_slice] = False
    assert np.all(arr[outer_mask] == dst_nodata) 

def _posDef_adjustment(mat, scaling_factor=0.99,max_it=1000):
    '''
    Checks whether the provided is pos semidefinite. If it is not, then it performs the the adjustment procedure descried in 1.2.2 of the Supplementary Note

    scaling_factor: the multiplicative factor that all off-diagonal elements of the matrix are scaled by in the second step of the procedure.
    max_it: max number of iterations set so that
    '''
    logging.info('Checking for positive definiteness ..')
    assert mat.ndim == 2
    assert mat.shape[0] == mat.shape[1]
    is_pos_semidef = lambda m: np.all(np.linalg.eigvals(m) >= 0)
    if is_pos_semidef(mat):
        return mat
    else:
        logging.info('matrix is not positive definite, performing adjustment..')
        P = mat.shape[0]
        for i in range(P):
            for j in range(i,P):
                if np.abs(mat[i,j]) > np.sqrt(mat[i,i] * mat[j,j]):
                    mat[i,j] = scaling_factor*np.sign(mat[i,j])*np.sqrt(mat[i,i] * mat[j,j])
                    mat[j,i] = mat[i,j]
        n=0
        while not is_pos_semidef(mat) and n < max_it:
            dg = np.diag(mat)
            mat = scaling_factor * mat
            mat[np.diag_indices(P)] = dg
            n += 1
        if n == max_it:
            logging.info('Warning: max number of iterations reached in adjustment procedure. Sigma matrix used is still non-positive-definite.')
        else:
            logging.info('Completed in {} iterations'.format(n))
        return mat 

def flatten_out_omega(omega_est):
    # stacks the lower part of the cholesky decomposition ROW_WISE [(0,0) (1,0) (1,1) (2,0) (2,1) (2,2) ...]
    P_c = len(omega_est)
    x_chol = np.linalg.cholesky(omega_est)

    # transform components of cholesky decomposition for better optimization
    lowTr_ind = np.tril_indices(P_c)
    x_chol_trf = np.zeros((P_c,P_c))
    for i in range(P_c):
        for j in range(i): # fill in lower triangular components not on diagonal
            x_chol_trf[i,j] = x_chol[i,j]/np.sqrt(x_chol[i,i]*x_chol[j,j])
    x_chol_trf[np.diag_indices(P_c)] = np.log(np.diag(x_chol))  # replace with log transformation on the diagonal
    return tuple(x_chol_trf[lowTr_ind]) 

def test_real_numbers(self):
        for m, n in self.test_dimensions:
            # Obtain a random real matrix of orthonormal rows
            Q = random_unitary_matrix(n, real=True)
            Q = Q[:m, :]

            # Get Givens decomposition of Q
            givens_rotations, V, diagonal = givens_decomposition(Q)

            # Compute U
            U = numpy.eye(n, dtype=complex)
            for parallel_set in givens_rotations:
                combined_givens = numpy.eye(n, dtype=complex)
                for i, j, theta, phi in reversed(parallel_set):
                    c = numpy.cos(theta)
                    s = numpy.sin(theta)
                    phase = numpy.exp(1.j * phi)
                    G = numpy.array([[c, -phase * s],
                                    [s, phase * c]], dtype=complex)
                    givens_rotate(combined_givens, G, i, j)
                U = combined_givens.dot(U)

            # Compute V * Q * U^\dagger
            W = V.dot(Q.dot(U.T.conj()))

            # Construct the diagonal matrix
            D = numpy.zeros((m, n), dtype=complex)
            D[numpy.diag_indices(m)] = diagonal

            # Assert that W and D are the same
            for i in range(m):
                for j in range(n):
                    self.assertAlmostEqual(D[i, j], W[i, j]) 

def make_corr(dimension, offdiag=0.0):

    corr=np.array([[offdiag]*dimension]*dimension)
    corr[np.diag_indices(dimension)]=1.0
    return corr 

def test_main_procedure(self):
        for m, n in self.test_dimensions:
            # Obtain a random matrix of orthonormal rows
            Q = random_unitary_matrix(n)
            Q = Q[:m, :]

            # Get Givens decomposition of Q
            givens_rotations, V, diagonal = givens_decomposition(Q)

            # Compute U
            U = numpy.eye(n, dtype=complex)
            for parallel_set in givens_rotations:
                combined_givens = numpy.eye(n, dtype=complex)
                for i, j, theta, phi in reversed(parallel_set):
                    c = numpy.cos(theta)
                    s = numpy.sin(theta)
                    phase = numpy.exp(1.j * phi)
                    G = numpy.array([[c, -phase * s],
                                     [s, phase * c]], dtype=complex)
                    givens_rotate(combined_givens, G, i, j)
                U = combined_givens.dot(U)

            # Compute V * Q * U^\dagger
            W = V.dot(Q.dot(U.T.conj()))

            # Construct the diagonal matrix
            D = numpy.zeros((m, n), dtype=complex)
            D[numpy.diag_indices(m)] = diagonal

            # Assert that W and D are the same
            for i in range(m):
                for j in range(n):
                    self.assertAlmostEqual(D[i, j], W[i, j]) 

def test_sds_and_corr_to_cov():
    sds = [1, 2, 3]
    corr = np.ones((3, 3)) * 0.2
    corr[np.diag_indices(3)] = 1
    calculated = sds_and_corr_to_cov(sds, corr)
    expected = np.array([[1.0, 0.4, 0.6], [0.4, 4.0, 1.2], [0.6, 1.2, 9.0]])
    aaae(calculated, expected) 

def get_modified_adj(self, ori_adj):
        adj_changes_square = self.adj_changes - torch.diag(torch.diag(self.adj_changes, 0))
        ind = np.diag_indices(self.adj_changes.shape[0])
        adj_changes_symm = torch.clamp(adj_changes_square + torch.transpose(adj_changes_square, 1, 0), -1, 1)
        modified_adj = adj_changes_symm + ori_adj
        return modified_adj 

def test_cov_to_sds_and_corr():
    cov = np.array([[1.0, 0.4, 0.6], [0.4, 4.0, 1.2], [0.6, 1.2, 9.0]])
    calc_sds, calc_corr = cov_to_sds_and_corr(cov)
    exp_sds = [1, 2, 3]
    exp_corr = np.ones((3, 3)) * 0.2
    exp_corr[np.diag_indices(3)] = 1
    aaae(calc_sds, exp_sds)
    aaae(calc_corr, exp_corr) 

def unvech(v):
    # quadratic formula, correct fp error
    rows = .5 * (-1 + np.sqrt(1 + 8 * len(v)))
    rows = int(np.round(rows))

    result = np.zeros((rows, rows))
    result[np.triu_indices(rows)] = v
    result = result + result.T

    # divide diagonal elements by 2
    result[np.diag_indices(rows)] /= 2

    return result 

def _diag_indices(n):
    rows, cols = np.diag_indices(n)
    return rows * n + cols 

def get_mvn_data(total_rvs, dimensionality=2, scale_sigma_offdiagonal_by=1., total_samples=1000):
    data_space_size = total_rvs * dimensionality

    # initialise distribution
    mu = np.random.randn(data_space_size)
    sigma = np.random.rand(data_space_size, data_space_size)
    # sigma = 1. + 0.5*np.random.randn(data_space_size, data_space_size)

    # ensures that sigma is positive semi-definite
    sigma = np.dot(sigma.transpose(), sigma)

    # scale off-diagonal entries -- might want to change that to block diagonal entries
    # diag = np.diag(sigma).copy()
    # sigma *= scale_sigma_offdiagonal_by
    # sigma[np.diag_indices(len(diag))] = diag

    # scale off-block diagonal entries
    d = dimensionality
    for ii, jj in itertools.product(list(range(total_rvs)), repeat=2):
        if ii != jj:
            sigma[d*ii:d*(ii+1), d*jj:d*(jj+1)] *= scale_sigma_offdiagonal_by

    # get samples
    samples = multivariate_normal(mu, sigma).rvs(total_samples)

    return [samples[:,ii*d:(ii+1)*d] for ii in range(total_rvs)] 

def _initialize_error_transition_individual(self):
        k_endog = self.k_endog
        _error_order = self._error_order

        # Initialize the parameters
        self.parameters['error_transition'] = _error_order

        # Fixed components already setup above

        # Setup indices of state space matrices
        # Here we want to set only the diagonal elements of the coefficient
        # matrices, and we want to set them in order by equation, not by
        # matrix (i.e. set the first element of the first matrix's diagonal,
        # then set the first element of the second matrix's diagonal, then...)

        # The basic setup is a tiled list of diagonal indices, one for each
        # coefficient matrix
        idx = np.tile(np.diag_indices(k_endog), self.error_order)
        # Now we need to shift the rows down to the correct location
        row_shift = self._factor_order
        # And we need to shift the columns in an increasing way
        col_inc = self._factor_order + np.repeat(
            [i * k_endog for i in range(self.error_order)], k_endog)
        idx[0] += row_shift
        idx[1] += col_inc

        # Make a copy (without the row shift) so that we can easily get the
        # diagonal parameters back out of a generic coefficients matrix array
        idx_diag = idx.copy()
        idx_diag[0] -= row_shift
        idx_diag[1] -= self._factor_order
        idx_diag = idx_diag[:, np.lexsort((idx_diag[1], idx_diag[0]))]
        self._idx_error_diag = (idx_diag[0], idx_diag[1])

        # Finally, we want to fill the entries in in the correct order, which
        # is to say we want to fill in lexicographically, first by row then by
        # column
        idx = idx[:, np.lexsort((idx[1], idx[0]))]
        self._idx_error_transition = np.s_['transition', idx[0], idx[1]] 

def loo_negloglikelihood(self):
		# get all pdf values of the training data (including 'self interaction')
		pdfs = self._individual_pdfs(self.data)
		
		# get indices to remove diagonal values for LOO part :)
		indices = np.diag_indices(pdfs.shape[0])

		# combine values based on fully_dimensional!
		if self.fully_dimensional:
			pdfs[indices] = 0 # remove self interaction		
			
			pdfs2 = np.sum(np.prod(pdfs, axis=-1), axis=-1)
			sm_return_value = -np.log(pdfs2).sum()
			pdfs = np.prod(pdfs, axis=-1)
			
			# take weighted average (accounts for LOO!)
			lhs = np.sum(pdfs*self.weights, axis=-1)/(1-self.weights)
		else:
			#import pdb; pdb.set_trace()
			pdfs[indices] = 0 # we sum first so 0 is the appropriate value
			pdfs *= self.weights[:,None,None]
			
			pdfs = pdfs.sum(axis=-2)/(1-self.weights[:,None])
			lhs = np.prod(pdfs, axis=-1)
			
		return(-np.sum(self.weights*np.log(lhs))) 

def test_naf_layer_diag():
    batch_size = 2
    for nb_actions in (1, 3):
        # Construct single model with NAF as the only layer, hence it is fully deterministic
        # since no weights are used, which would be randomly initialized.
        L_flat_input = Input(shape=(nb_actions,))
        mu_input = Input(shape=(nb_actions,))
        action_input = Input(shape=(nb_actions,))
        x = NAFLayer(nb_actions, mode='diag')([L_flat_input, mu_input, action_input])
        model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
        model.compile(loss='mse', optimizer='sgd')
        
        # Create random test data.
        L_flat = np.random.random((batch_size, nb_actions)).astype('float32')
        mu = np.random.random((batch_size, nb_actions)).astype('float32')
        action = np.random.random((batch_size, nb_actions)).astype('float32')

        # Perform reference computations in numpy since these are much easier to verify.
        P = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
        for p, l_flat in zip(P, L_flat):
            p[np.diag_indices(nb_actions)] = l_flat
        print(P, L_flat)
        A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
        A_ref *= -.5

        # Finally, compute the output of the net, which should be identical to the previously
        # computed reference.
        A_net = model.predict([L_flat, mu, action]).flatten()
        assert_allclose(A_net, A_ref, rtol=1e-5) 

def test_naf_layer_full():
    batch_size = 2
    for nb_actions in (1, 3):
        # Construct single model with NAF as the only layer, hence it is fully deterministic
        # since no weights are used, which would be randomly initialized.
        L_flat_input = Input(shape=((nb_actions * nb_actions + nb_actions) // 2,))
        mu_input = Input(shape=(nb_actions,))
        action_input = Input(shape=(nb_actions,))
        x = NAFLayer(nb_actions, mode='full')([L_flat_input, mu_input, action_input])
        model = Model(inputs=[L_flat_input, mu_input, action_input], outputs=x)
        model.compile(loss='mse', optimizer='sgd')
        
        # Create random test data.
        L_flat = np.random.random((batch_size, (nb_actions * nb_actions + nb_actions) // 2)).astype('float32')
        mu = np.random.random((batch_size, nb_actions)).astype('float32')
        action = np.random.random((batch_size, nb_actions)).astype('float32')

        # Perform reference computations in numpy since these are much easier to verify.
        L = np.zeros((batch_size, nb_actions, nb_actions)).astype('float32')
        LT = np.copy(L)
        for l, l_T, l_flat in zip(L, LT, L_flat):
            l[np.tril_indices(nb_actions)] = l_flat
            l[np.diag_indices(nb_actions)] = np.exp(l[np.diag_indices(nb_actions)])
            l_T[:, :] = l.T
        P = np.array([np.dot(l, l_T) for l, l_T in zip(L, LT)]).astype('float32')
        A_ref = np.array([np.dot(np.dot(a - m, p), a - m) for a, m, p in zip(action, mu, P)]).astype('float32')
        A_ref *= -.5

        # Finally, compute the output of the net, which should be identical to the previously
        # computed reference.
        A_net = model.predict([L_flat, mu, action]).flatten()
        assert_allclose(A_net, A_ref, rtol=1e-5) 

def __init__(self, mu_tra=np.zeros(3), mu_rot=np.eye(3), 
                 sigma_tra=np.eye(3), sigma_rot=np.eye(3),
                 from_frame='world', to_frame='world',
                 *args, **kwargs):
        if np.abs(np.linalg.det(mu_rot) - 1.0) > 1e-3:
            raise ValueError('Illegal rotation. Must have determinant == 1.0')
        if not is_positive_semi_definite(sigma_tra):
            raise ValueError('Translation covariance is not positive semi-definite!')
        if not is_positive_semi_definite(sigma_rot):
            raise ValueError('Rotation covariance is not positive semi-definite!')

        # read params
        self._mu_tra = mu_tra.copy()
        self._mu_rot = mu_rot.copy()
        self._sigma_tra = sigma_tra.copy()
        self._sigma_rot = sigma_rot.copy()

        diag_idx = np.diag_indices(3)
        self._sigma_tra[diag_idx] = np.clip(np.diag(self._sigma_tra), 1e-10, np.inf)
        self._sigma_rot[diag_idx] = np.clip(np.diag(self._sigma_rot), 1e-10, np.inf)
        
        self._from_frame = from_frame
        self._to_frame = to_frame

        # setup random variables
        self._t_rv = scipy.stats.multivariate_normal(self._mu_tra, self._sigma_tra)
        self._r_xi_rv = scipy.stats.multivariate_normal(np.zeros(3), self._sigma_rot)
        super(GaussianRigidTransformRandomVariable, self).__init__(*args, **kwargs) 

def make_corr(dimension, offdiag=0.0):

    corr=np.array([[offdiag]*dimension]*dimension)
    corr[np.diag_indices(dimension)]=1.0
    return corr