Pairwise comparisons (Chi-Square) along axis, as numpy.ndarray.

        Returns a square, symmetric matrix of test statistics for the null
        hypothesis that each vector along *axis* is equal to each other.
766021
def _pairwise_chisq(self):
        
        return self._chi_squared(self._proportions, self._margin, self._observed)
741075
def _remove_na_single(x, axis='rows'):
    
    if x.ndim == 1:
        # 1D arrays
        x_mask = ~np.isnan(x)
    else:
        # 2D arrays
        ax = 1 if axis == 'rows' else 0
        x_mask = ~np.any(np.isnan(x), axis=ax)
    # Check if missing values are present
    if ~x_mask.all():
        ax = 0 if axis == 'rows' else 1
        ax = 0 if x.ndim == 1 else ax
        x = x.compress(x_mask, axis=ax)
    return x
355573
def is_adjacency_matrix_symmetric(adjacency_matrix):
    """"""
    # Verify that the matrix is square
    num_columns = len(adjacency_matrix)
    for column in adjacency_matrix:
        # In a square matrix, every row should be the same length as the number of columns
        if len(column) != num_columns:
            return False

    # Loop through the bottom half of the matrix and compare it to the top half
    # --We do the bottom half because of how we construct adjacency matrices
    max_i = 0
    for j in range(num_columns):
        for i in range(max_i):
            # If i == j, we can skip ahead so we don't compare with ourself
            if i == j:
                continue
            # Compare the value in the bottom half with the mirrored value in the top half
            # If they aren't the same, the matrix isn't symmetric
            if adjacency_matrix[j][i] != adjacency_matrix[i][j]:
                return False
        max_i += 1

    # If we reach this far without returning false, then we know that everything matched,
    # which makes this a symmetric matrix
    return True
543712
def ucross(arr1, arr2, registry=None, axisa=-1, axisb=-1, axisc=-1, axis=None):
    
    v = np.cross(arr1, arr2, axisa=axisa, axisb=axisb, axisc=axisc, axis=axis)
    units = arr1.units * arr2.units
    arr = unyt_array(v, units, registry=registry)
    return arr
129744
def _sigmaclip_withaxis(self, data, axis=None, masked=True,
                            return_bounds=False, copy=True):
        

        # float array type is needed to insert nans into the array
        filtered_data = data.astype(float)    # also makes a copy

        # remove invalid values
        bad_mask = ~np.isfinite(filtered_data)
        if np.any(bad_mask):
            filtered_data[bad_mask] = np.nan
            warnings.warn('Input data contains invalid values (NaNs or '
                          'infs), which were automatically clipped.',
                          AstropyUserWarning)

        # remove masked values and convert to plain ndarray
        if isinstance(filtered_data, np.ma.MaskedArray):
            filtered_data = np.ma.masked_invalid(filtered_data).astype(float)
            filtered_data = filtered_data.filled(np.nan)

        # convert negative axis/axes
        if not isiterable(axis):
            axis = (axis,)
        axis = tuple(filtered_data.ndim + n if n < 0 else n for n in axis)

        # define the shape of min/max arrays so that they can be broadcast
        # with the data
        mshape = tuple(1 if dim in axis else size
                       for dim, size in enumerate(filtered_data.shape))

        nchanged = 1
        iteration = 0
        while nchanged != 0 and (iteration < self.maxiters):
            iteration += 1
            n_nan = np.count_nonzero(np.isnan(filtered_data))

            self._compute_bounds(filtered_data, axis=axis)
            if not np.isscalar(self._min_value):
                self._min_value = self._min_value.reshape(mshape)
                self._max_value = self._max_value.reshape(mshape)

            with np.errstate(invalid='ignore'):
                filtered_data[(filtered_data < self._min_value) |
                              (filtered_data > self._max_value)] = np.nan

            nchanged = n_nan - np.count_nonzero(np.isnan(filtered_data))

        self._niterations = iteration

        if masked:
            # create an output masked array
            if copy:
                filtered_data = np.ma.masked_invalid(filtered_data)
            else:
                # ignore RuntimeWarnings for comparisons with NaN data values
                with np.errstate(invalid='ignore'):
                    out = np.ma.masked_invalid(data, copy=False)

                    filtered_data = np.ma.masked_where(np.logical_or(
                        out < self._min_value, out > self._max_value),
                        out, copy=False)

        if return_bounds:
            return filtered_data, self._min_value, self._max_value
        else:
            return filtered_data
977202
def iterative_plane_errors(axes,covariance_matrix, **kwargs):
    
    sheet = kwargs.pop('sheet','upper')
    level = kwargs.pop('level',1)
    n = kwargs.pop('n',100)

    cov = N.sqrt(N.diagonal(covariance_matrix))
    u = N.linspace(0, 2*N.pi, n)

    scales = dict(upper=1,lower=-1,nominal=0)
    c1 = scales[sheet]
    c1 *= -1 # We assume upper hemisphere
    if axes[2,2] < 0:
        c1 *= -1

    def sdot(a,b):
        return sum([i*j for i,j in zip(a,b)])

    def step_func(a):
        e = [
            N.cos(a)*cov[0],
            N.sin(a)*cov[1],
            c1*cov[2]]
        d = [sdot(e,i)
            for i in axes.T]
        x,y,z = d[2],d[0],d[1]
        r = N.sqrt(x**2 + y**2 + z**2)
        lat = N.arcsin(z/r)
        lon = N.arctan2(y, x)
        return lon,lat

    # Get a bundle of vectors defining
    # a full rotation around the unit circle
    return N.array([step_func(i)
        for i in u])
1199874
def spatial_correlation_array(dmatrix, imts, correl='yes',
                              vs30clustered=True):
    
    assert correl in 'yes no full', correl
    n = len(dmatrix)
    corr = numpy.zeros((len(imts), n, n))
    for imti, im in enumerate(imts):
        if correl == 'no':
            corr[imti] = numpy.eye(n)
        if correl == 'full':
            corr[imti] = numpy.ones((n, n))
        elif correl == 'yes':
            corr[imti] = correlation.jbcorrelation(dmatrix, im, vs30clustered)
    return corr
1118567
def scoreatpercentile(data,per,axis=0):
    ''
    a = np.sort(data,axis=axis)
    idx = per/100. * (data.shape[axis]-1)

    if (idx % 1 == 0):
        return a[[slice(None) if ii != axis else idx for ii in range(a.ndim)]]
    else:
        lowerweight = 1-(idx % 1)
        upperweight = (idx % 1)
        idx = int(np.floor(idx))
        return lowerweight * a[[slice(None) if ii != axis else idx for ii in range(a.ndim)]] \
                + upperweight * a[[slice(None) if ii != axis else idx+1 for ii in range(a.ndim)]]
388968
def argsort_indices(a, axis=-1):
    
    a = np.asarray(a)
    ind = list(np.ix_(*[np.arange(d) for d in a.shape]))
    ind[axis] = a.argsort(axis)
    return tuple(ind)
57512
def K(self, X, X2=None):
                """"""        
                X,slices = X[:,:-1],index_to_slices(X[:,-1])
                if X2 is None:
                        X2,slices2 = X,slices
                        K = np.zeros((X.shape[0], X.shape[0]))
                else:
                        X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
                        K = np.zeros((X.shape[0], X2.shape[0]))

                tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
                ttdist = (X[:,0][:,None] - X2[:,0][None,:])
                
                vyt = self.variance_Yt
                
                lyt=1/(2*self.lengthscale_Yt)

                a = -self.a
                c = self.c

                kyy = lambda tdist: np.exp(-lyt*(tdist))

                k1 = lambda tdist: (2*lyt - 4*lyt**2 *(tdist) )

                k4 = lambda tdist: 2*lyt*(tdist)

                for i, s1 in enumerate(slices):
                        for j, s2 in enumerate(slices2):
                                for ss1 in s1:
                                    for ss2 in s2:
                                        if i==0 and j==0:
                                            K[ss1,ss2] = vyt*kyy(tdist[ss1,ss2])
                                        elif i==0 and j==1:
                                            K[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                            #K[ss1,ss2] = (2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                        elif i==1 and j==1:
                                            K[ss1,ss2] = ( k1(tdist[ss1,ss2]) + 1. )*vyt* kyy(tdist[ss1,ss2])+self.ubias
                                        else:
                                            K[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                                            #K[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
                #stop
                return K
1004986
def _pairwise_corr(self, columns=None, covar=None, tail='two-sided',
                   method='pearson', padjust='none', export_filename=None):
    """"""
    stats = pairwise_corr(data=self, columns=columns, covar=covar,
                          tail=tail, method=method, padjust=padjust,
                          export_filename=export_filename)
    return stats
1166196
def summary_pairwise_indices(self):
        """"""
        summary_pairwise_indices = np.empty(
            self.values[0].t_stats.shape[1], dtype=object
        )
        summary_pairwise_indices[:] = [
            sig.summary_pairwise_indices for sig in self.values
        ]
        return summary_pairwise_indices
1295143
def chisquare(n_ij, weighted):
    
    if weighted:
        m_ij = n_ij / n_ij

        nan_mask = np.isnan(m_ij)
        m_ij[nan_mask] = 0.000001  # otherwise it breaks the chi-squared test

        w_ij = m_ij
        n_ij_col_sum = n_ij.sum(axis=1)
        n_ij_row_sum = n_ij.sum(axis=0)
        alpha, beta, eps = (1, 1, 1)
        while eps > 10e-6:
            alpha = alpha * np.vstack(n_ij_col_sum / m_ij.sum(axis=1))
            beta = n_ij_row_sum / (alpha * w_ij).sum(axis=0)
            eps = np.max(np.absolute(w_ij * alpha * beta - m_ij))
            m_ij = w_ij * alpha * beta

    else:
        m_ij = (np.vstack(n_ij.sum(axis=1)) * n_ij.sum(axis=0)) / n_ij.sum().astype(float)

    dof = (n_ij.shape[0] - 1) * (n_ij.shape[1] - 1)
    chi, p_val = stats.chisquare(n_ij, f_exp=m_ij, ddof=n_ij.size - 1 - dof, axis=None)

    return (chi, p_val, dof)
469883
def leastsq_generalized(x, y, x2, y2, res_std=None, res2_std=None, fit_offset=False, fit_offset2=False, perc=None):
    
    if perc is not None:
        if not fit_offset and isinstance(perc, (list, tuple)): perc = perc[1]
        weights = csr_matrix(get_weight(x, y, perc)).astype(bool)
        x, y = weights.multiply(x).tocsr(), weights.multiply(y).tocsr()

    n_obs, n_var = x.shape
    offset, offset_ss = np.zeros(n_var, dtype="float32"), np.zeros(n_var, dtype="float32")
    gamma = np.ones(n_var, dtype="float32")

    if (res_std is None) or (res2_std is None): res_std, res2_std = np.ones(n_var), np.ones(n_var)
    ones, zeros = np.ones(n_obs), np.zeros(n_obs)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        x, y = np.vstack((make_dense(x)/res_std, x2/res2_std)), np.vstack((make_dense(y)/res_std, y2/res2_std))

    if fit_offset and fit_offset2:
        for i in range(n_var):
            A = np.c_[np.vstack((np.c_[ones/res_std[i], zeros], np.c_[zeros, ones/res2_std[i]])), x[:, i]]
            offset[i], offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset:
        for i in range(n_var):
            A = np.c_[np.hstack((ones/res_std[i], zeros)), x[:, i]]
            offset[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    elif fit_offset2:
        for i in range(n_var):
            A = np.c_[np.hstack((zeros, ones/res2_std[i])), x[:, i]]
            offset_ss[i], gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))
    else:
        for i in range(n_var):
            A = np.c_[x[:, i]]
            gamma[i] = np.linalg.pinv(A.T.dot(A)).dot(A.T.dot(y[:, i]))

    offset[np.isnan(offset)] = 0
    offset_ss[np.isnan(offset_ss)] = 0
    gamma[np.isnan(gamma)] = 0

    return offset, offset_ss, gamma
118804
def _ep_pairwise(
        n_items, comparisons, alpha, match_moments, max_iter, initial_state):
    
    # Static variable that allows to check the # of iterations after the call.
    _ep_pairwise.iterations = 0
    m = len(comparisons)
    prior_inv = alpha * np.eye(n_items)
    if initial_state is None:
        # Initially, mean and covariance come from the prior.
        mean = np.zeros(n_items)
        cov = (1 / alpha) * np.eye(n_items)
        # Initialize the natural params in the function space.
        tau = np.zeros(m)
        nu = np.zeros(m)
        # Initialize the natural params in the space of thetas.
        prec = np.zeros((n_items, n_items))
        xs = np.zeros(n_items)
    else:
        tau, nu = initial_state
        mean, cov, xs, prec = _init_ws(
                n_items, comparisons, prior_inv, tau, nu)
    for _ in range(max_iter):
        _ep_pairwise.iterations += 1
        # Keep a copy of the old parameters for convergence testing.
        tau_old = np.array(tau, copy=True)
        nu_old = np.array(nu, copy=True)
        for i in nprand.permutation(m):
            a, b = comparisons[i]
            # Update mean and variance in function space.
            f_var = cov[a,a] + cov[b,b] - 2 * cov[a,b]
            f_mean = mean[a] - mean[b]
            # Cavity distribution.
            tau_tot = 1.0 / f_var
            nu_tot = tau_tot * f_mean
            tau_cav = tau_tot - tau[i]
            nu_cav = nu_tot - nu[i]
            cov_cav = 1.0 / tau_cav
            mean_cav = cov_cav * nu_cav
            # Moment matching.
            logpart, dlogpart, d2logpart = match_moments(mean_cav, cov_cav)
            # Update factor params in the function space.
            tau[i] = -d2logpart / (1 + d2logpart / tau_cav)
            delta_tau = tau[i] - tau_old[i]
            nu[i] = ((dlogpart - (nu_cav / tau_cav) * d2logpart)
                     / (1 + d2logpart / tau_cav))
            delta_nu = nu[i] - nu_old[i]
            # Update factor params in the weight space.
            prec[(a, a, b, b), (a, b, a, b)] += delta_tau * MAT_ONE_FLAT
            xs[a] += delta_nu
            xs[b] -= delta_nu
            # Update mean and covariance.
            if abs(delta_tau) > 0:
                phi = -1.0 / ((1.0 / delta_tau) + f_var) * MAT_ONE
                upd_mat = cov.take([a, b], axis=0)
                cov = cov + upd_mat.T.dot(phi).dot(upd_mat)
            mean = cov.dot(xs)
        # Recompute the global parameters for stability.
        cov = inv_posdef(prior_inv + prec)
        mean = cov.dot(xs)
        if _converged((tau, nu), (tau_old, nu_old)):
            return mean, cov
    raise RuntimeError(
            "EP did not converge after {} iterations".format(max_iter))
279046
def X(self) -> Optional[Union[np.ndarray, sparse.spmatrix, ArrayView]]:
        """"""
        if self.isbacked:
            if not self.file.isopen: self.file.open()
            X = self.file['X']
            if self.isview:
                X = X[self._oidx, self._vidx]
            return X
        else:
            if self.n_obs == 1 and self.n_vars == 1:
                return self._X[0, 0]
            elif self.n_obs == 1 or self.n_vars == 1:
                X = self._X
                if issparse(self._X): X = self._X.toarray()
                return X.flatten()
            else:
                return self._X
103310
def chi_squared(*choices):
    """"""

    term = lambda expected, observed: float((expected - observed) ** 2) / max(expected, 1)
    mean_success_rate = float(sum([c.rewards for c in choices])) / max(sum([c.plays for c in choices]), 1)
    mean_failure_rate = 1 - mean_success_rate

    return sum([
        term(mean_success_rate * c.plays, c.rewards)
        + term(mean_failure_rate * c.plays, c.plays - c.rewards
    ) for c in choices])