Python scipy.spatial.distance.cdist() Examples

def query(self, X, k, return_dist=False):
        if self.compact:
            X = self._validate_input(X)
            cdist = hamming_cdist(X, self._fit_X)
        else:
            X = self._validate_input(X, False)
            cdist = distance.cdist(X, self._fit_X, 'hamming')
        ind = np.argsort(cdist, 1)[:, :k]
        if return_dist:
            rows = np.arange(ind.shape[0])[:, np.newaxis]
            dist = cdist[rows, ind]
            if not self.compact:
                dist = (dist * X.shape[1]).astype(int)
            return ind, dist
        else:
            return ind 

def __init__(self, cities, **kwargs):
        """
        A two-dimensional traveling salesman problem (TSP)

        Parameters
        ----------
        cities : numpy.array
            The cities with 2-dimensional coordinates provided by a matrix where where city is represented by a row.

        """
        n_cities, _ = cities.shape

        self.cities = cities
        self.D = cdist(cities, cities)

        super(TravelingSalesman, self).__init__(
            n_var=n_cities,
            n_obj=1,
            xl=0,
            xu=n_cities,
            type_var=np.int,
            elementwise_evaluation=True,
            **kwargs
        ) 

def vna(self, coords, **kw):
      """
      Compute the neutral-atom potential V_NA(coords) for a set of Cartesian
      coordinates coords.
      The subroutine could be also used for computing the non-linear core
      corrections or some other atom-centered fields.
      """

      sp2v = kw['sp2v'] if 'sp2v' in kw else self.ao_log.sp2vna
      sp2rcut = kw['sp2rcut'] if 'sp2rcut' in kw else self.ao_log.sp2rcut_vna
      atom2coord = kw['atom2coord'] if 'atom2coord' in kw else self.atom2coord

      nc = coords.shape[0]
      vna = np.zeros(nc)
      for ia,(R,sp) in enumerate(zip(atom2coord, self.atom2sp)):
          if sp2v[sp] is None: # This can be done better via preparation of a special atom2sp excluding ghost atoms
              continue
          #print(__name__, ia, sp, sp2rcut[sp])
          dd = cdist(R.reshape((1,3)), coords).reshape(nc)
          vnaa = self.ao_log.interp_rr(sp2v[sp], dd, rcut=sp2rcut[sp])
          vna = vna + vnaa
      return vna 

def kernel_matrix_xX(svm_model, original_x, original_X):

        if (svm_model.svm_kernel == 'polynomial_kernel' or svm_model.svm_kernel == 'soft_polynomial_kernel'):
            K = (svm_model.zeta + svm_model.gamma * np.dot(original_x, original_X.T)) ** svm_model.Q
        elif (svm_model.svm_kernel == 'gaussian_kernel' or svm_model.svm_kernel == 'soft_gaussian_kernel'):
            K = np.exp(-svm_model.gamma * (cdist(original_X, np.atleast_2d(original_x), 'euclidean').T ** 2)).ravel()

        '''
        K = np.zeros((svm_model.data_num, svm_model.data_num))

        for i in range(svm_model.data_num):
            for j in range(svm_model.data_num):
                if (svm_model.svm_kernel == 'polynomial_kernel' or svm_model.svm_kernel == 'soft_polynomial_kernel'):
                    K[i, j] = Kernel.polynomial_kernel(svm_model, original_x, original_X[j])
                elif (svm_model.svm_kernel == 'gaussian_kernel' or svm_model.svm_kernel == 'soft_gaussian_kernel'):
                    K[i, j] = Kernel.gaussian_kernel(svm_model, original_x, original_X[j])
        '''

        return K 

def _similarity_smooth(presented, recalled, features, distance):
    lists = presented.index.get_values()
    res = np.empty((len(lists), len(features), recalled.iloc[0].shape[0], presented.iloc[0].shape[0]))*np.nan
    for li, l in enumerate(lists):
        p_list = presented.loc[l]
        r_list = recalled.loc[l]
        for i, feature in enumerate(features):
            get_feature = lambda x: np.array(x[feature]) if np.array(pd.notna(x['item'])).any() else np.nan
            p = np.vstack(p_list.apply(get_feature).get_values())
            r = r_list.dropna().apply(get_feature).get_values()
            r = np.vstack(list(filter(lambda x: x is not np.nan, r)))
            tmp = 1 - cdist(r, p, distance)
            res[li, i, :tmp.shape[0], :] =  tmp
    if distance == 'correlation':
        return np.nanmean(res, 1)
    else:
        return np.mean(res, 1) 

def _cold_start_initialize(self):

        # Initialize the models for the untrained arms with the data from the most similar arm
        for arm in self.untrained_arms:

            distances = {}
            for n in self.arms:
                if n not in self.untrained_arms:
                    distances[n] = cdist(np.asarray([self.features[n]]), np.asarray([self.features[arm]]),
                                         metric='cosine')

            # Identify the closest arm
            closest_arm = min(distances, key=distances.get)
            print('Cold Start Arm:', arm, 'Closest arm:', closest_arm)

            # Set the arm to use the values of the closets arm
            self.arm_to_model[arm].beta = self.arm_to_model[closest_arm].beta.copy()
            self.arm_to_model[arm].A = self.arm_to_model[closest_arm].A.copy()
            self.arm_to_model[arm].A_inv = self.arm_to_model[closest_arm].A_inv.copy()
            self.arm_to_model[arm].Xty = self.arm_to_model[closest_arm].Xty.copy()

            self.untrained_arms.remove(arm) 

def compute_pairwise_distances(first_xyz, second_xyz):
  """Computes pairwise distances between two molecules.

  Takes an input (m, 3) and (n, 3) numpy arrays of 3D coords of
  two molecules respectively, and outputs an m x n numpy
  array of pairwise distances in Angstroms between the first and
  second molecule. entry (i,j) is dist between the i"th 
  atom of first molecule and the j"th atom of second molecule.

  Parameters
  ----------
  first_xyz: np.ndarray
    Of shape (m, 3)
  seocnd_xyz: np.ndarray
    Of shape (n, 3)

  Returns
  -------
  np.ndarray of shape (m, n)
  """

  pairwise_distances = cdist(first_xyz, second_xyz, metric='euclidean')
  return pairwise_distances 

def execute(cls, ctx, op):
        from scipy.spatial.distance import cdist

        inputs, device_id, xp = as_same_device(
            [ctx[inp.key] for inp in op.inputs], device=op.device, ret_extra=True)

        if xp is cp:  # pragma: no cover
            raise NotImplementedError('`cdist` does not support running on GPU yet')

        with device(device_id):
            inputs_iter = iter(inputs)
            xa = next(inputs_iter)
            xb = next(inputs_iter)
            kw = dict()
            if op.p is not None:
                kw['p'] = op.p
            if op.w is not None:
                kw['w'] = next(inputs_iter)
            if op.v is not None:
                kw['V'] = next(inputs_iter)
            if op.vi is not None:
                kw['VI'] = next(inputs_iter)

        ctx[op.outputs[0].key] = cdist(xa, xb, metric=op.metric, **kw) 

def tracklet_classify(A, pca, D, knn, clf_coding):
    encode_fea = np.zeros((len(A),len(D)))
    for n in range(len(A)):
        pca_fea = pca.transform(A[n])
        dist = distance.cdist(pca_fea, D, 'euclidean')
        x = np.zeros((len(pca_fea),len(D)))
        for k in range(len(dist)):
            sort_idx = np.argsort(dist[k,:])
            temp_D = D[sort_idx[0:knn],:]
            temp_coder = SparseCoder(dictionary=temp_D, transform_n_nonzero_coefs=10, 
                                 transform_alpha=0.05, transform_algorithm='lasso_lars')
            #import pdb; pdb.set_trace()
            xx = np.zeros((1,D.shape[1]))
            xx[:,:] = pca_fea[k,:]
            temp_x = temp_coder.transform(xx)
            x[k,sort_idx[0:knn]] = temp_x

        encode_fea[n,:] = np.max(x, axis=0)
    pred_set_label = clf_coding.predict(encode_fea)
    return pred_set_label 

def gen_gausprocess_se(ntrain, ntest, noise=1., lenscale=1., scale=1.,
                       xmin=-10, xmax=10):
    """
    Generate a random (noisy) draw from a Gaussian Process with a RBF kernel.
    """

    # Xtrain = np.linspace(xmin, xmax, ntrain)[:, np.newaxis]
    Xtrain = np.random.rand(ntrain)[:, np.newaxis] * (xmin - xmax) - xmin
    Xtest = np.linspace(xmin, xmax, ntest)[:, np.newaxis]
    Xcat = np.vstack((Xtrain, Xtest))

    K = scale * np.exp(-cdist(Xcat, Xcat, metric='sqeuclidean') /
                       (2 * lenscale**2))
    U, S, V = np.linalg.svd(K)
    L = U.dot(np.diag(np.sqrt(S))).dot(V)
    f = np.random.randn(ntrain + ntest).dot(L)

    ytrain = f[0:ntrain] + np.random.randn(ntrain) * noise
    ftest = f[ntrain:]

    return Xtrain, ytrain, Xtest, ftest 

def transform(self, X, lenscale=None):
        r"""
        Apply the sigmoid basis function to X.

        Parameters
        ----------
        X: ndarray
            (N, d) array of observations where N is the number of samples, and
            d is the dimensionality of X.
        lenscale: float
            the length scale (scalar) of the RBFs to apply to X. If not input,
            this uses the value of the initial length scale.

        Returns
        -------
        ndarray:
            of shape (N, D) where D is number of centres.
        """
        N, d = X.shape
        lenscale = self._check_dim(d, lenscale)

        return expit(cdist(X / lenscale, self.C / lenscale, 'euclidean')) 

def _spatial_sort(glyph):
  from scipy.spatial.distance import cdist
  from numpy import argsort
  from numpy import argmin

  curr = argmin(glyph[:,0])
  visited = set([curr])
  order = [curr]

  dd = cdist(glyph, glyph)

  while len(visited)<len(glyph):
    row = dd[curr,:]

    for i in argsort(row):
      if row[i]<=0.0 or i==curr or i in visited:
        continue
      order.append(i)
      visited.add(i)
      break
  glyph[:,:] = glyph[order,:] 

def analysis_KMeans():
	mean_distortions = []
	K = len(labels_idx)
	K_range = range(320, 1000)
	for k in K_range:
		print("Cluster k is {}".format(k))
		kmeans_model = KMeans(n_clusters=k, init="k-means++", n_jobs=-1)
		kmeans_model.fit(np_features)
		t_distortions = sum(
			np.min(cdist(np_features, kmeans_model.cluster_centers_, 'euclidean'), axis=1)) / np_features.shape[0]
		mean_distortions.append(t_distortions)

	with open("./kmeans_cluster.csv", "a+") as wh:
		for idx in range(len(K_range)):
			wh.write("{},{}\n".format(K_range[idx], mean_distortions[idx]))

	# plt.plot(K_range, mean_distortions, 'bx-')
	# plt.xlabel('k')
	# plt.ylabel(u'Avgerage distortion degree')
	# plt.title(u'Elbows rule to select the best K value')
	# plt.savefig("kmeans_cluster.png") 

def _find_site_neighbors(site_locs, n_neighbors, shank_map):
    from scipy.spatial.distance import cdist

    if np.unique(shank_map).size <= 1:
        pass

    n_sites = site_locs.shape[0]
    n_neighbors = int(min(n_neighbors, n_sites))

    neighbors = np.zeros((n_sites, n_neighbors), dtype=np.int)
    for i in range(n_sites):
        i_loc = site_locs[i, :][np.newaxis, :]
        dists = cdist(i_loc, site_locs).ravel()
        neighbors[i, :] = dists.argsort()[:n_neighbors]

    return neighbors 

def tsne_to_grid(X_2d):
    from lapjv import lapjv
    from scipy.spatial.distance import cdist

    out_dim = np.sqrt(len(X_2d))
    out_dim = int(out_dim)
    to_plot = np.square(out_dim)
    grid = np.dstack(np.meshgrid(np.linspace(0, 1, out_dim), np.linspace(0, 1, out_dim))).reshape(-1, 2)
    cost_matrix = cdist(grid, X_2d[:to_plot], "sqeuclidean").astype(np.float32)
    cost_matrix = cost_matrix * (100000 / cost_matrix.max())
    row_asses, col_asses, _ = lapjv(cost_matrix)
    grid_jv = grid[col_asses]

    return grid_jv, to_plot
    # out = np.ones((out_dim*out_res, out_dim*out_res, 3))

    # to_plot = np.square(out_dim)
    # for pos, img in zip(grid_jv, img_collection[0:to_plot]):
    #     h_range = int(np.floor(pos[0]* (out_dim - 1) * out_res))
    #     w_range = int(np.floor(pos[1]* (out_dim - 1) * out_res))
    #     out[h_range:h_range + out_res, w_range:w_range + out_res]  = image.img_to_array(img)

    # im = image.array_to_img(out)
    # im.save(out_dir + out_name, quality=100) 

def test_euclidean_distances(dtype, x_array_constr, y_array_constr):
    # check that euclidean distances gives same result as scipy cdist
    # when X and Y != X are provided
    rng = np.random.RandomState(0)
    X = rng.random_sample((100, 10)).astype(dtype, copy=False)
    X[X < 0.8] = 0
    Y = rng.random_sample((10, 10)).astype(dtype, copy=False)
    Y[Y < 0.8] = 0

    expected = cdist(X, Y)

    X = x_array_constr(X)
    Y = y_array_constr(Y)
    distances = euclidean_distances(X, Y)

    # the default rtol=1e-7 is too close to the float32 precision
    # and fails due too rounding errors.
    assert_allclose(distances, expected, rtol=1e-6)
    assert distances.dtype == dtype 

def transform(self, X, lenscale=None):
        """
        Apply the RBF to X.

        Parameters
        ----------
        X: ndarray
            (N, d) array of observations where N is the number of samples, and
            d is the dimensionality of X.
        lenscale: scalar or ndarray, optional
            scalar or array of shape (d,) length scales (one for each dimension
            of X). If not input, this uses the value of the initial length
            scale.

        Returns
        -------
        ndarray:
            of shape (N, D) where D is number of RBF centres.
        """
        N, d = X.shape
        lenscale = self._check_dim(d, lenscale)

        den = (2 * lenscale**2)
        return np.exp(- cdist(X / den, self.C / den, 'sqeuclidean')) 

def test_euclidean_distances_upcast(batch_size, x_array_constr,
                                    y_array_constr):
    # check batches handling when Y != X (#13910)
    rng = np.random.RandomState(0)
    X = rng.random_sample((100, 10)).astype(np.float32)
    X[X < 0.8] = 0
    Y = rng.random_sample((10, 10)).astype(np.float32)
    Y[Y < 0.8] = 0

    expected = cdist(X, Y)

    X = x_array_constr(X)
    Y = y_array_constr(Y)
    distances = _euclidean_distances_upcast(X, Y=Y, batch_size=batch_size)
    distances = np.sqrt(np.maximum(distances, 0))

    # the default rtol=1e-7 is too close to the float32 precision
    # and fails due too rounding errors.
    assert_allclose(distances, expected, rtol=1e-6) 

def _shrink(X, npoints):
        """
        When designs are generated that are larger than the requested number of points (N* > N), resize them.
        If the size was correct all along, the LHD is returned unchanged.

        :param X: Generated LHD, size N* x D, with N* >= N
        :param npoints: What size to resize to (N)
        :return: LHD data matrix, size N x D
        """
        npStar, nv = X.shape

        # Pick N samples nearest to centre of X
        centre = npStar * np.ones((1, nv)) / 2.
        distances = cdist(X, centre).ravel()
        idx = np.argsort(distances)
        X = X[idx[:npoints], :]

        # Translate to origin
        X -= np.min(X, axis=0) - 1

        # Collapse gaps in the design to assure all cell projections onto axes have 1 sample
        Xs = np.argsort(X, axis=0)
        X[Xs, np.arange(nv)] = np.tile(np.arange(1, npoints + 1), (nv, 1)).T
        assert (X.shape[0] == npoints)
        return X 

def tracklet_classify(A, pca, D, knn, clf_coding):
    encode_fea = np.zeros((len(A),len(D)))
    for n in range(len(A)):
        pca_fea = pca.transform(A[n])
        dist = distance.cdist(pca_fea, D, 'euclidean')
        x = np.zeros((len(pca_fea),len(D)))
        for k in range(len(dist)):
            sort_idx = np.argsort(dist[k,:])
            temp_D = D[sort_idx[0:knn],:]
            temp_coder = SparseCoder(dictionary=temp_D, transform_n_nonzero_coefs=10, 
                                 transform_alpha=0.05, transform_algorithm='lasso_lars')
            #import pdb; pdb.set_trace()
            xx = np.zeros((1,D.shape[1]))
            xx[:,:] = pca_fea[k,:]
            temp_x = temp_coder.transform(xx)
            x[k,sort_idx[0:knn]] = temp_x

        encode_fea[n,:] = np.max(x, axis=0)
    pred_set_label = clf_coding.predict(encode_fea)
    return pred_set_label 

def test_pdist_bool_metrics(metric):
    D_true = cdist(X1_bool, X1_bool, metric)
    check_pdist_bool(metric, D_true) 

def test_pdist(metric):
    argdict = METRICS_DEFAULT_PARAMS[metric]
    keys = argdict.keys()
    for vals in itertools.product(*argdict.values()):
        kwargs = dict(zip(keys, vals))
        D_true = cdist(X1, X1, metric, **kwargs)
        check_pdist(metric, kwargs, D_true) 

def __init__(self,
                 ref_dirs,
                 n_neighbors=20,
                 decomposition='auto',
                 prob_neighbor_mating=0.9,
                 display=MultiObjectiveDisplay(),
                 **kwargs):
        """

        Parameters
        ----------
        ref_dirs
        n_neighbors
        decomposition
        prob_neighbor_mating
        display
        kwargs
        """

        self.n_neighbors = n_neighbors
        self.prob_neighbor_mating = prob_neighbor_mating
        self.decomposition = decomposition

        set_if_none(kwargs, 'pop_size', len(ref_dirs))
        set_if_none(kwargs, 'sampling', FloatRandomSampling())
        set_if_none(kwargs, 'crossover', SimulatedBinaryCrossover(prob=1.0, eta=20))
        set_if_none(kwargs, 'mutation', PolynomialMutation(prob=None, eta=20))
        set_if_none(kwargs, 'survival', None)
        set_if_none(kwargs, 'selection', None)

        super().__init__(display=display, **kwargs)

        # initialized when problem is known
        self.ref_dirs = ref_dirs

        if self.ref_dirs.shape[0] < self.n_neighbors:
            print("Setting number of neighbours to population size: %s" % self.ref_dirs.shape[0])
            self.n_neighbors = self.ref_dirs.shape[0]

        # neighbours includes the entry by itself intentionally for the survival method
        self.neighbors = np.argsort(cdist(self.ref_dirs, self.ref_dirs), axis=1, kind='quicksort')[:, :self.n_neighbors] 

def grad(self, X, lenscale=None):
        r"""
        Get the gradients of this basis w.r.t.\ the length scale.

        Parameters
        ----------
        X: ndarray
            (N, d) array of observations where N is the number of samples, and
            d is the dimensionality of X.
        lenscale: scalar or ndarray, optional
            scalar or array of shape (d,) length scales (one for each dimension
            of X). If not input, this uses the value of the initial length
            scale.

        Returns
        -------
        ndarray:
            of shape (N, D) where D is number of RBF centres. This is
            :math:`\partial \Phi(\mathbf{X}) / \partial l`
        """
        N, d = X.shape
        lenscale = self._check_dim(d, lenscale)

        Phi = self.transform(X, lenscale)
        dPhi = []
        for i, l in enumerate(lenscale):
            ldist = cdist(X[:, [i]] / l**3, self.C[:, [i]] / l**3,
                          'sqeuclidean')
            dPhi.append(Phi * ldist)

        return np.dstack(dPhi) if len(lenscale) != 1 else dPhi[0] 

def grad(self, X, lenscale=None):
        r"""
        Get the gradients of this basis w.r.t.\  the length scale.

        Parameters
        ----------
        X: ndarray
            (N, d) array of observations where N is the number of samples, and
            d is the dimensionality of X.
        lenscale: float, optional
            the length scale (scalar) of the RBFs to apply to X. If not input,
            this uses the value of the initial length scale.

        Returns
        -------
        ndarray:
            of shape (N, D) where D is number of centres. This is
            :math:`\partial \Phi(\mathbf{X}) / \partial l`
        """
        N, d = X.shape
        lenscale = self._check_dim(d, lenscale)

        Phi = self.transform(X, lenscale)
        dPhi = []
        for i, l in enumerate(lenscale):
            ldist = cdist(X[:, [i]] / l**2, self.C[:, [i]] / l**2, 'euclidean')
            dPhi.append(- ldist * Phi * (1 - Phi))

        return np.dstack(dPhi) if len(lenscale) != 1 else dPhi[0] 

def compute_distanceMatrix(ft1,ft2):
    distMat = cdist(ft1,ft2,"cosine")
    return distMat


# In[ ]: 

def test_cdist_bool_metric(metric):
    D_true = cdist(X1_bool, X2_bool, metric)
    check_cdist_bool(metric, D_true) 

def calc_pts_diameter2(pts):
    """
    Calculates diameter of a set of points (i.e. the maximum distance between
    any two points in the set). Faster but requires more memory than
    calc_pts_diameter.

    :param pts: nx3 ndarray with 3D points.
    :return: Diameter.
    """
    dists = distance.cdist(pts, pts, 'euclidean')
    diameter = np.max(dists)
    return diameter 

def mine_triplets(anchor, targets, negative_samples):
    distances = cdist(anchor, targets, 'cosine')
    distances = cdist(anchor, targets, 'cosine').tolist()
    QnQ_duplicated = [
        [target_index for target_index, dist in enumerate(QnQ_dist) if dist == QnQ_dist[query_index]]
        for query_index, QnQ_dist in enumerate(distances)]
    for i, QnT_dist in enumerate(QnT_dists):
        for j in QnQ_duplicated[i]:
            QnT_dist.itemset(j, np.inf)

    QnT_dists_topk = QnT_dists.argsort(axis=1)[:, :negative_samples]
    top_k_index = np.array([np.insert(QnT_dist, 0, i) for i, QnT_dist in enumerate(QnT_dists_topk)])
    return top_k_index 

def calc_pts_diameter2(pts):
    """
    Calculates diameter of a set of points (i.e. the maximum distance between
    any two points in the set). Faster but requires more memory than
    calc_pts_diameter.

    :param pts: nx3 ndarray with 3D points.
    :return: Diameter.
    """
    dists = distance.cdist(pts, pts, 'euclidean')
    diameter = np.max(dists)
    return diameter