%prog ld map

    Calculate pairwise linkage disequilibrium given MSTmap.
322514
def ld(args):
    
    import numpy as np
    from random import sample

    from jcvi.algorithms.matrix import symmetrize

    p = OptionParser(ld.__doc__)
    p.add_option("--subsample", default=1000, type="int",
                 help="Subsample markers to speed up [default: %default]")
    opts, args, iopts = p.set_image_options(args, figsize="8x8")

    if len(args) != 1:
        sys.exit(not p.print_help())

    mstmap, = args
    subsample = opts.subsample
    data = MSTMap(mstmap)

    markerbedfile = mstmap + ".subsample.bed"
    ldmatrix = mstmap + ".subsample.matrix"
    # Take random subsample while keeping marker order
    if subsample < data.nmarkers:
        data = [data[x] for x in \
                sorted(sample(xrange(len(data)), subsample))]
    else:
        logging.debug("Use all markers, --subsample ignored")

    nmarkers = len(data)
    if need_update(mstmap, (ldmatrix, markerbedfile)):
        fw = open(markerbedfile, "w")
        print("\n".join(x.bedline for x in data), file=fw)
        logging.debug("Write marker set of size {0} to file `{1}`."\
                        .format(nmarkers, markerbedfile))
        fw.close()

        M = np.zeros((nmarkers, nmarkers), dtype=float)
        for i, j in combinations(range(nmarkers), 2):
            a = data[i]
            b = data[j]
            M[i, j] = calc_ldscore(a.genotype, b.genotype)

        M = symmetrize(M)

        logging.debug("Write LD matrix to file `{0}`.".format(ldmatrix))
        M.tofile(ldmatrix)
    else:
        nmarkers = len(Bed(markerbedfile))
        M = np.fromfile(ldmatrix, dtype="float").reshape(nmarkers, nmarkers)
        logging.debug("LD matrix `{0}` exists ({1}x{1})."\
                        .format(ldmatrix, nmarkers))

    from jcvi.graphics.base import plt, savefig, Rectangle, draw_cmap

    plt.rcParams["axes.linewidth"] = 0

    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])
    ax = fig.add_axes([.1, .1, .8, .8])  # the heatmap

    ax.matshow(M, cmap=iopts.cmap)

    # Plot chromosomes breaks
    bed = Bed(markerbedfile)
    xsize = len(bed)
    extent = (0, nmarkers)
    chr_labels = []
    ignore_size = 20

    for (seqid, beg, end) in bed.get_breaks():
        ignore = abs(end - beg) < ignore_size
        pos = (beg + end) / 2
        chr_labels.append((seqid, pos, ignore))
        if ignore:
            continue
        ax.plot((end, end), extent, "w-", lw=1)
        ax.plot(extent, (end, end), "w-", lw=1)

    # Plot chromosome labels
    for label, pos, ignore in chr_labels:
        pos = .1 + pos * .8 / xsize
        if not ignore:
            root.text(pos, .91, label,
                ha="center", va="bottom", rotation=45, color="grey")
            root.text(.09, pos, label,
                ha="right", va="center", color="grey")

    ax.set_xlim(extent)
    ax.set_ylim(extent)
    ax.set_axis_off()

    draw_cmap(root, "Pairwise LD (r2)", 0, 1, cmap=iopts.cmap)

    root.add_patch(Rectangle((.1, .1), .8, .8, fill=False, ec="k", lw=2))
    m = mstmap.split(".")[0]
    root.text(.5, .06, "Linkage Disequilibrium between {0} markers".format(m), ha="center")

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    image_name = m + ".subsample" + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
394987
def _M(self, k, t, tips=None, gaps=None):
        """"""
        if self._distributionmodel:
            return self.model.M(k, t, tips, gaps)
        else:
            return self.model.M(t, tips, gaps)
227310
def estimate_map(interface, state, label, inp):
    """"""
    out = interface.output(0)
    centers = {}
    for row in inp:
        row = row.strip().split(state["delimiter"])
        if len(row) > 1:
            x = [(0 if row[i] in state["missing_vals"] else float(row[i])) for i in state["X_indices"]]
            cluster = min((state['dist'](c, x), i) for i, c in state['centers'])[1]
            vertex = state['create'](x, 1.0)
            centers[cluster] = vertex if cluster not in centers else state["update"](centers[cluster], vertex)

    for cluster, values in centers.iteritems():
        out.add(cluster, values)
324484
def get_linkage(self, line):
        """"""
        conf1, id1, chain1, pos1 = line[16].strip(), line[17:20].strip(), line[21].strip(), int(line[22:26])
        conf2, id2, chain2, pos2 = line[46].strip(), line[47:50].strip(), line[51].strip(), int(line[52:56])
        return self.covlinkage(id1=id1, chain1=chain1, pos1=pos1, conf1=conf1,
                               id2=id2, chain2=chain2, pos2=pos2, conf2=conf2)
1363111
def _compute_diff(map1, map2):
        
        data = map1.data - map2.data
        return HpxMap(data, map1.hpx)
1180438
def run_gap_k_selection(data, k_min=1, k_max=50, B=5,
        skip=5, **kwargs):
    
    if k_min == k_max:
        return k_min
    gap_vals = []
    sk_vals = []
    k_range = list(range(k_min, k_max, skip))
    min_k = 0
    min_i = 0
    for i, k in enumerate(k_range):
        km = KMeans(k)
        clusters = km.fit_predict(data)
        gap, sk = calculate_gap(data, clusters, km, B=B)
        if len(gap_vals) > 1:
            if gap_vals[-1] >= gap - (skip+1)*sk:
                min_i = i
                min_k = k_range[i-1]
                break
                #return k_range[-1], gap_vals, sk_vals
        gap_vals.append(gap)
        sk_vals.append(sk)
    if min_k == 0:
        min_k = k_max
    if skip == 1:
        return min_k, gap_vals, sk_vals
    gap_vals = []
    sk_vals = []
    for k in range(min_k - skip, min_k + skip):
        km = KMeans(k)
        clusters = km.fit_predict(data)
        gap, sk = calculate_gap(data, clusters, km, B=B)
        if len(gap_vals) > 1:
            if gap_vals[-1] >= gap - sk:
                min_k = k-1
                return min_k, gap_vals, sk_vals
        gap_vals.append(gap)
        sk_vals.append(sk)
    return k, gap_vals, sk_vals
764780
def assign_order(self):
        
        linkage_groups = self.linkage_groups
        for mlg in linkage_groups:
            mapname = mlg.mapname
            if mapname == self.pivot:
                pivot_position = mlg.position

        for mlg in linkage_groups:
            position = mlg.position
            # Flip order if path goes in the opposite direction to the pivot
            common = []
            for a, ap in position.items():
                if a not in pivot_position:
                    continue
                pp = pivot_position[a]
                common.append((ap, pp))

            mlg.rho = get_rho(common)
            if mlg.rho < 0:
                mlg.path = mlg.path[::-1]

            mlg.populate_pairwise_distance()

        # Preparation of TSP
        distances = defaultdict(list)
        for mlg in linkage_groups:
            mapname = mlg.mapname
            position = mlg.position
            length = mlg.length
            path = mlg.path
            rho = mlg.rho
            dd = mlg.distances
            for a, b in combinations(path, 2):
                d = dd[a, b]
                distances[a, b].append((d, mapname))
            for p in path:
                adist, bdist = position[p], length - position[p]
                if rho < 0:
                    adist, bdist = bdist, adist
                distances[START, p].append((adist, mapname))
                distances[p, END].append((bdist, mapname))

        self.distances = distances
        tour = self.distances_to_tour()
        return tour
84532
def lazy_binmap(f, xs):
    
    return (f(x, y) for x, y in zip(xs, xs[1:]))
1362232
def find_map( x1, x2, eps = 1e-8, allow_double = False, join = True ):
    
    off, dim = x1.shape
    ir = nm.zeros( (off + x2.shape[0],), dtype = nm.int32 )
    ir[off:] = off

    x1 = nm.round( x1.T / eps ) * eps
    x2 = nm.round( x2.T / eps ) * eps
    xx = nm.c_[x1, x2]

    keys = [xx[ii] for ii in range( dim )]
    iis = nm.lexsort( keys = keys )

    xs = xx.T[iis]
##     import scipy.io as io
##     io.write_array( 'sss1', x1.T )
##     io.write_array( 'sss2', x2.T )
##     io.write_array( 'sss', xs, precision = 16 )
##     pause()
    xd = nm.sqrt( nm.sum( nm.diff( xs, axis = 0 )**2.0, axis = 1 ) )

    ii = nm.where( xd < eps )[0]
    off1, off2 = ir[iis][ii], ir[iis][ii+1]
    i1, i2 = iis[ii] - off1, iis[ii+1] - off2
    dns = nm.where( off1 == off2 )[0]
    if dns.size:
        print 'double node(s) in:'
        for dn in dns:
            if off1[dn] == 0:
                print 'x1: %d %d -> %s %s' % (i1[dn], i2[dn],
                                              x1[:,i1[dn]], x1[:,i2[dn]])
            else:
                print 'x2: %d %d -> %s %s' % (i1[dn], i2[dn],
                                              x2[:,i1[dn]], x2[:,i2[dn]])
        if not allow_double:
            raise ValueError

    if join:
        cmap = nm.c_[i1, i2]
        return cmap
    else:
        return i1, i2
292947
def _map_parallel(function, args, n_jobs):
    """"""
    # Following the error checking found in joblib
    multiprocessing = int(os.environ.get('JOBLIB_MULTIPROCESSING', 1)) or None
    if multiprocessing:
        try:
            import multiprocessing
            import multiprocessing.pool
        except ImportError:
            multiprocessing = None
        if sys.platform.startswith("win") and PY2:
            msg = "Multiprocessing is not supported on Windows with Python 2.X. Setting n_jobs=1"
            logger.warning(msg)
            n_jobs = 1
    # 2nd stage: validate that locking is available on the system and
    #            issue a warning if not
    if multiprocessing:
        try:
            _sem = multiprocessing.Semaphore()
            del _sem  # cleanup
        except (ImportError, OSError) as e:
            multiprocessing = None
            logger.warning('{}. _map_parallel will operate in serial mode'.format(e))
    if multiprocessing and int(n_jobs) not in (0, 1):
        if n_jobs == -1:
            n_jobs = None
        try:
            pool = multiprocessing.Pool(processes=n_jobs)
            map_result = pool.map(function, args)
        finally:
            pool.close()
            pool.join()
    else:
        map_result = list(map(function, args))
    return map_result
534870
def connectivity_map(dset,prefix,x,y,z,radius=2):
    ''''''
    seed_series = nl.sphere_average(dset,x,y,z,radius)
    with tempfile.NamedTemporaryFile(delete=False) as temp:
        temp.write('\n'.join([str(x) for x in seed_series]))
    decon = nl.Decon()
    decon.input_dsets = dset
    decon.stim_files = {'seed':temp.name}
    decon.prefix = prefix
    decon.run()
    try:
        os.remove(temp.name)
    except:
        pass
634149
def cplus(*args):
    '''
    
    '''
    n = len(args)
    if   n == 0: return np.asarray(0)
    elif n == 1: return np.asarray(args[0])
    elif n >  2: return reduce(plus, args)
    (a,b) = args
    if sps.issparse(a):
        if not sps.issparse(b):
            b = np.asarray(b)
            if len(b.shape) == 0: b = np.reshape(b, (1,1))
    elif sps.issparse(b):
        a = np.asarray(a)
        if len(a.shape) == 0: a = np.reshape(a, (1,1))
    else:
        a = np.asarray(a)
        b = np.asarray(b)
    return a + b
1050941
def run_diff_umap(self,use_rep='X_pca', metric='euclidean', n_comps=15,
                      method='gauss', **kwargs):
                     
        import scanpy.api as sc
        sc.pp.neighbors(self.adata,use_rep=use_rep,n_neighbors=self.k,
                                       metric=self.distance,method=method)
        sc.tl.diffmap(self.adata, n_comps=n_comps)
        
        sc.pp.neighbors(self.adata,use_rep='X_diffmap',n_neighbors=self.k,
                                       metric='euclidean',method=method)
                
        if 'X_umap' in self.adata.obsm.keys():
                self.adata.obsm['X_umap_sam'] = self.adata.obsm['X_umap']
                
        sc.tl.umap(self.adata,min_dist=0.1,copy=False)
15715
def gen_ldm(src, dst):
        
        return ReilBuilder.build(ReilMnemonic.LDM, src, ReilEmptyOperand(), dst)
242229
def populate_iteration(self, iteration):
        """"""

        cur_idx = iteration.cur_idx
        genotypes = self.genotype_file.next().split()
        iteration.chr, iteration.rsid, junk, iteration.pos = genotypes[0:4]
        iteration.chr = int(iteration.chr)
        iteration.pos = int(iteration.pos)

        if DataParser.boundary.TestBoundary(iteration.chr, iteration.pos, iteration.rsid):
            try:
                [iteration.genotype_data,
                    iteration.major_allele,
                    iteration.minor_allele,
                    iteration.hetero_count,
                    iteration.maj_allele_count,
                    iteration.min_allele_count,
                    iteration.missing_allele_count,
                    iteration.allele_count2] = self.process_genotypes(genotypes)
                return iteration.maf >= DataParser.min_maf and iteration.maf <= DataParser.max_maf
            except TooFewAlleles:
                print "\n\n\nSkipping %s:%s %s %s" % (iteration.chr, iteration.pos, iteration.rsid, cur_idx)

        return False
710602
def mape(y, p):
    

    filt = np.abs(y) > EPS
    return np.mean(np.abs(1 - p[filt] / y[filt]))
485365
def mapper_graph(df, lens_data=None, lens='pca', resolution=10, gain=0.5, equalize=True, clust='kmeans', stat='db',
                 max_K=5):
    
    if lens_data is None:
        lens_data = apply_lens(df, lens=lens)

    patch_clusterings = {}
    counter = 0
    patches = covering_patches(lens_data, resolution=resolution, gain=gain, equalize=equalize)
    for key, patch in patches.items():
        if len(patch) > 0:
            patch_clusterings[key] = optimal_clustering(df, patch, method=clust, statistic=stat, max_K=max_K)
            counter += 1
    print 'total of {} patches required clustering'.format(counter)

    all_clusters = []
    for key in patch_clusterings:
        all_clusters += patch_clusterings[key]
    num_nodes = len(all_clusters)
    print 'this implies {} nodes in the mapper graph'.format(num_nodes)

    A = np.zeros((num_nodes, num_nodes))
    for i in range(num_nodes):
        for j in range(i):
            overlap = set(all_clusters[i]).intersection(set(all_clusters[j]))
            if len(overlap) > 0:
                A[i, j] = 1
                A[j, i] = 1

    G = nx.from_numpy_matrix(A)
    total = []
    all_clusters_new = []
    mapping = {}
    cont = 0
    for m in all_clusters:
        total += m
    for n, m in enumerate(all_clusters):
        if len(m) == 1 and total.count(m) > 1:
            G.remove_node(n)
        else:
            all_clusters_new.append(m)
            mapping[n] = cont
            cont += 1
    H = nx.relabel_nodes(G, mapping)
    return H, all_clusters_new, patches