Build the landmark operator

        Calculates spectral clusters on the kernel, and calculates transition
        probabilities between cluster centers by using transition probabilities
        between samples assigned to each cluster.
13690
def build_landmark_op(self):
        
        tasklogger.log_start("landmark operator")
        is_sparse = sparse.issparse(self.kernel)
        # spectral clustering
        tasklogger.log_start("SVD")
        _, _, VT = randomized_svd(self.diff_aff,
                                  n_components=self.n_svd,
                                  random_state=self.random_state)
        tasklogger.log_complete("SVD")
        tasklogger.log_start("KMeans")
        kmeans = MiniBatchKMeans(
            self.n_landmark,
            init_size=3 * self.n_landmark,
            batch_size=10000,
            random_state=self.random_state)
        self._clusters = kmeans.fit_predict(
            self.diff_op.dot(VT.T))
        # some clusters are not assigned
        tasklogger.log_complete("KMeans")

        # transition matrices
        pmn = self._landmarks_to_data()

        # row normalize
        pnm = pmn.transpose()
        pmn = normalize(pmn, norm='l1', axis=1)
        pnm = normalize(pnm, norm='l1', axis=1)
        landmark_op = pmn.dot(pnm)  # sparsity agnostic matrix multiplication
        if is_sparse:
            # no need to have a sparse landmark operator
            landmark_op = landmark_op.toarray()
        # store output
        self._landmark_op = landmark_op
        self._transitions = pnm
        tasklogger.log_complete("landmark operator")
1332648
def _nominal_kernel(x, y, out):
    """"""
    for i in range(x.shape[0]):
        for j in range(y.shape[0]):
            out[i, j] += (x[i, :] == y[j, :]).sum()

    return out
866581
def _prob_match(self, features):
        

        # compute the probabilities
        probs = self.kernel.predict_proba(features)

        # get the position of match probabilities
        classes = list(self.kernel.classes_)
        match_class_position = classes.index(1)

        return probs[:, match_class_position]
1264326
def _centroids(n_clusters: int, points: List[List[float]]) -> List[List[float]]:
    

    k_means = KMeans(n_clusters=n_clusters)
    k_means.fit(points)

    closest, _ = pairwise_distances_argmin_min(k_means.cluster_centers_, points)

    return list(map(list, np.array(points)[closest.tolist()]))
766812
def calculate_reorganize_data_output_shapes(operator):
    '''
    
    '''
    check_input_and_output_numbers(operator, input_count_range=1, output_count_range=1)
    check_input_and_output_types(operator, good_input_types=[FloatTensorType])

    output_shape = copy.deepcopy(operator.inputs[0].type.shape)

    params = operator.raw_operator.reorganizeData

    from coremltools.proto.NeuralNetwork_pb2 import ReorganizeDataLayerParams as Params

    if params.mode == Params.DEPTH_TO_SPACE:
        if output_shape[1] % (params.blockSize * params.blockSize) != 0:
            raise RuntimeError('Channel number must be divisible by the square of block size')

        output_shape = [output_shape[0], output_shape[1] / params.blockSize / params.blockSize,
                        output_shape[2] * params.blockSize, output_shape[3] * params.blockSize]
    elif params.mode == Params.SPACE_TO_DEPTH:
        if output_shape[2] % params.blockSize != 0 or output_shape[3] % params.blockSize != 0:
            raise RuntimeError('Height and weight must be divisible by block size')

        output_shape = [output_shape[0], output_shape[1] * params.blockSize * params.blockSize,
                        output_shape[2] / params.blockSize, output_shape[3] / params.blockSize]
    else:
        raise ValueError('Unsupport reorganization mode {0}'.format(params.mode))

    operator.outputs[0].type = FloatTensorType([int(i) if i != 'None' else 'None' for i in output_shape],
                                               doc_string=operator.outputs[0].type.doc_string)
1307439
def init_centers_widths(self, R):
        

        kmeans = KMeans(
            init='k-means++',
            n_clusters=self.K,
            n_init=10,
            random_state=100)
        kmeans.fit(R)
        centers = kmeans.cluster_centers_
        widths = self._get_max_sigma(R) * np.ones((self.K, 1))
        return centers, widths
50566
def spectral_entropy(data, numOfShortBlocks=10):
    """"""
    data = np.mean(data, axis=1)

    nFFT = len(data) // 2
    X = FFT(data, nFFT)
    L = len(X)                         # number of frame data
    Eol = np.sum(X ** 2)            # total spectral energy

    subWinLength = int(np.floor(L / numOfShortBlocks))   # length of sub-frame
    if L != subWinLength * numOfShortBlocks:
        X = X[0:subWinLength * numOfShortBlocks]

    # define sub-frames (using matrix reshape)
    subWindows = X.reshape(subWinLength, numOfShortBlocks, order='F').copy()

    # compute spectral sub-energies
    s = np.sum(subWindows ** 2, axis=0) / (Eol + EPSILON)

    # compute spectral entropy
    return -np.sum(s * np.log2(s + EPSILON))
14594
def k_means(
    X,
    n_clusters,
    init="k-means||",
    precompute_distances="auto",
    n_init=1,
    max_iter=300,
    verbose=False,
    tol=1e-4,
    random_state=None,
    copy_x=True,
    n_jobs=-1,
    algorithm="full",
    return_n_iter=False,
    oversampling_factor=2,
    init_max_iter=None,
):
    
    labels, inertia, centers, n_iter = _kmeans_single_lloyd(
        X,
        n_clusters,
        max_iter=max_iter,
        init=init,
        verbose=verbose,
        tol=tol,
        random_state=random_state,
        oversampling_factor=oversampling_factor,
        init_max_iter=init_max_iter,
    )
    if return_n_iter:
        return labels, centers, inertia, n_iter
    else:
        return labels, centers, inertia
619129
def inference(self, kern, X, likelihood, Y, Y_metadata=None):

        
        N = X.shape[0] #number of training points
        D = X.shape[1] #number of dimensions

        Kds = np.zeros(D, dtype=object) #vector for holding covariance per dimension
        Qs = np.zeros(D, dtype=object) #vector for holding eigenvectors of covariance per dimension
        QTs = np.zeros(D, dtype=object) #vector for holding transposed eigenvectors of covariance per dimension
        V_kron = 1 # kronecker product of eigenvalues

        # retrieve the one-dimensional variation of the designated kernel
        oneDkernel = kern.get_one_dimensional_kernel(D)

        for d in range(D):
            xg = list(set(X[:,d])) #extract unique values for a dimension
            xg = np.reshape(xg, (len(xg), 1))
            oneDkernel.lengthscale = kern.lengthscale[d]
            Kds[d] = oneDkernel.K(xg)
            [V, Q] = np.linalg.eig(Kds[d])
            V_kron = np.kron(V_kron, V)
            Qs[d] = Q
            QTs[d] = Q.T

        noise = likelihood.variance + 1e-8

        alpha_kron = self.kron_mvprod(QTs, Y)
        V_kron = V_kron.reshape(-1, 1)
        alpha_kron = alpha_kron / (V_kron + noise)
        alpha_kron = self.kron_mvprod(Qs, alpha_kron)

        log_likelihood = -0.5 * (np.dot(Y.T, alpha_kron) + np.sum((np.log(V_kron + noise))) + N*log_2_pi)

        # compute derivatives wrt parameters Thete
        derivs = np.zeros(D+2, dtype='object')
        for t in range(len(derivs)):
            dKd_dTheta = np.zeros(D, dtype='object')
            gamma = np.zeros(D, dtype='object')
            gam = 1
            for d in range(D):
                xg = list(set(X[:,d]))
                xg = np.reshape(xg, (len(xg), 1))
                oneDkernel.lengthscale = kern.lengthscale[d]
                if t < D:
                    dKd_dTheta[d] = oneDkernel.dKd_dLen(xg, (t==d), lengthscale=kern.lengthscale[t]) #derivative wrt lengthscale
                elif (t == D):
                    dKd_dTheta[d] = oneDkernel.dKd_dVar(xg) #derivative wrt variance
                else:
                    dKd_dTheta[d] = np.identity(len(xg)) #derivative wrt noise
                gamma[d] = np.diag(np.dot(np.dot(QTs[d], dKd_dTheta[d].T), Qs[d]))
                gam = np.kron(gam, gamma[d])
            
            gam = gam.reshape(-1,1)
            kappa = self.kron_mvprod(dKd_dTheta, alpha_kron)
            derivs[t] = 0.5*np.dot(alpha_kron.T,kappa) - 0.5*np.sum(gam / (V_kron + noise))

        # separate derivatives
        dL_dLen = derivs[:D]
        dL_dVar = derivs[D]
        dL_dThetaL = derivs[D+1]

        return GridPosterior(alpha_kron=alpha_kron, QTs=QTs, Qs=Qs, V_kron=V_kron), \
                log_likelihood, {'dL_dLen':dL_dLen, 'dL_dVar':dL_dVar, 'dL_dthetaL':dL_dThetaL}
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')
940414
def init_cao(X, n_clusters, dissim):
    
    n_points, n_attrs = X.shape
    centroids = np.empty((n_clusters, n_attrs), dtype='object')
    # Method is based on determining density of points.
    dens = np.zeros(n_points)
    for iattr in range(n_attrs):
        freq = defaultdict(int)
        for val in X[:, iattr]:
            freq[val] += 1
        for ipoint in range(n_points):
            dens[ipoint] += freq[X[ipoint, iattr]] / float(n_points) / float(n_attrs)

    # Choose initial centroids based on distance and density.
    centroids[0] = X[np.argmax(dens)]
    if n_clusters > 1:
        # For the remaining centroids, choose maximum dens * dissim to the
        # (already assigned) centroid with the lowest dens * dissim.
        for ik in range(1, n_clusters):
            dd = np.empty((ik, n_points))
            for ikk in range(ik):
                dd[ikk] = dissim(X, centroids[ikk]) * dens
            centroids[ik] = X[np.argmax(np.min(dd, axis=0))]

    return centroids
509936
def create_model(samples_x, samples_y_aggregation,
                 n_restarts_optimizer=250, is_white_kernel=False):
    '''
    
    '''
    kernel = gp.kernels.ConstantKernel(constant_value=1,
                                       constant_value_bounds=(1e-12, 1e12)) * \
                                                gp.kernels.Matern(nu=1.5)
    if is_white_kernel is True:
        kernel += gp.kernels.WhiteKernel(noise_level=1, noise_level_bounds=(1e-12, 1e12))
    regressor = gp.GaussianProcessRegressor(kernel=kernel,
                                            n_restarts_optimizer=n_restarts_optimizer,
                                            normalize_y=True,
                                            alpha=1e-10)
    regressor.fit(numpy.array(samples_x), numpy.array(samples_y_aggregation))

    model = {}
    model['model'] = regressor
    model['kernel_prior'] = str(kernel)
    model['kernel_posterior'] = str(regressor.kernel_)
    model['model_loglikelihood'] = regressor.log_marginal_likelihood(regressor.kernel_.theta)

    return model
641275
def _kmedoids_run(X, n_clusters, distance, max_iter, tol, rng):
    
    membs = np.empty(shape=X.shape[0], dtype=int)
    centers = kmeans._kmeans_init(X, n_clusters, method='', rng=rng)

    sse_last = 9999.9
    n_iter = 0
    for it in range(1,max_iter):
        membs = kmeans._assign_clusters(X, centers)
        centers,sse_arr = _update_centers(X, membs, n_clusters, distance)
        sse_total = np.sum(sse_arr)
        if np.abs(sse_total - sse_last) < tol:
            n_iter = it
            break
        sse_last = sse_total

    return(centers, membs, sse_total, sse_arr, n_iter)
987847
def calculate_convolution_output_shapes(operator):
    '''
    
    '''
    check_input_and_output_numbers(operator, input_count_range=1, output_count_range=1)

    params = operator.raw_operator.convolution

    input_shape = operator.inputs[0].type.shape
    operator.outputs[0].type.shape = [0, 0, 0, 0]  # Initialize output shape. It will be modified below.
    output_shape = operator.outputs[0].type.shape

    # Adjust N-axis
    output_shape[0] = input_shape[0]

    # Adjust C-axis
    output_shape[1] = params.outputChannels

    # Set up default and non-default parameters
    dilations = [1, 1]
    if len(params.dilationFactor) > 0:
        dilations = [params.dilationFactor[0], params.dilationFactor[1]]
    kernel_shape = [3, 3]
    if len(params.kernelSize) > 0:
        kernel_shape = params.kernelSize
    strides = [1, 1]
    if len(params.stride) > 0:
        strides = params.stride
    specified_output_shape = [0, 0]  # Only used with convolution transpose
    if params.isDeconvolution and len(params.outputShape) > 0:
        specified_output_shape = list(int(i) for i in params.outputShape)
    pad_mode = params.WhichOneof('ConvolutionPaddingType')
    if pad_mode == 'valid' and len(params.valid.paddingAmounts.borderAmounts) > 0:
        pad_amounts = params.valid.paddingAmounts.borderAmounts
        pad_heads = [pad_amounts[0].startEdgeSize, pad_amounts[1].startEdgeSize]
        pad_tails = [pad_amounts[0].endEdgeSize, pad_amounts[1].endEdgeSize]
    else:
        # Padding amounts are useless for same padding and valid padding uses [0, 0] by default.
        pad_heads = [0, 0]
        pad_tails = [0, 0]

    # Adjust H- and W-axes
    for i in range(2):
        if params.isDeconvolution:
            output_shape[i + 2] = calculate_convolution_transpose_1D_output_shape(
                input_shape[i + 2], kernel_shape[i], dilations[i], strides[i],
                pad_mode, pad_heads[i], pad_tails[i], specified_output_shape[i])
        else:
            output_shape[i + 2] = calculate_convolution_and_pooling_1D_output_shape(
                input_shape[i + 2], kernel_shape[i], dilations[i], strides[i],
                pad_mode, pad_heads[i], pad_tails[i])
35541
def computePCsPlink(plink_path,k,out_dir,bfile,ffile):
    
    print("Using plink to compute principal components")
    cmd = '%s --bfile %s --pca %d '%(plink_path,bfile,k)
    cmd+= '--out %s'%(os.path.join(out_dir,'plink'))
    subprocess.call(cmd,shell=True)
    plink_fn = os.path.join(out_dir, 'plink.eigenvec')
    M = sp.loadtxt(plink_fn,dtype=str)
    U = sp.array(M[:,2:],dtype=float)
    U-= U.mean(0)
    U/= U.std(0)
    sp.savetxt(ffile,U)
1231533
def _kmeans(X, n_clusters, max_iter, n_trials, tol):
    
    n_samples, n_features = X.shape[0], X.shape[1]

    centers_best = np.empty(shape=(n_clusters,n_features), dtype=float)
    labels_best  = np.empty(shape=n_samples, dtype=int)
    for i in range(n_trials):
        centers, labels, sse_tot, sse_arr, n_iter  = _kmeans_run(X, n_clusters, max_iter, tol)
        if i==0:
            sse_tot_best = sse_tot
            sse_arr_best = sse_arr
            n_iter_best = n_iter
            centers_best = centers.copy()
            labels_best  = labels.copy()
        if sse_tot < sse_tot_best:
            sse_tot_best = sse_tot
            sse_arr_best = sse_arr
            n_iter_best = n_iter
            centers_best = centers.copy()
            labels_best  = labels.copy()

    return(centers_best, labels_best, sse_arr_best, n_iter_best)
818505
def cluster_centers_(self):
        
        # TODO: Handle agglomerative clustering and LDA
        for attr in ('cluster_centers_',):
            try:
                return getattr(self.estimator, attr)
            except AttributeError:
                continue

        raise AttributeError(
            "could not find or make cluster_centers_ for {}".format(
            self.estimator.__class__.__name__
        ))