Creates a plot of a catalogue showing where particular clusters exist
329240
def add_catalogue_cluster(self, catalogue, vcl, flagvector,
                              cluster_id=None, overlay=True):
        
        # Create simple magnitude scaled point basemap
        self.add_size_scaled_points(catalogue.data['longitude'],
                                    catalogue.data['latitude'],
                                    catalogue.data['magnitude'],
                                    shape="o",
                                    alpha=0.8,
                                    colour=(0.5, 0.5, 0.5),
                                    smin=1.0,
                                    sscale=1.5,
                                    overlay=True)
        # If cluster ID is not specified just show mainshocks
        if cluster_id is None:
            idx = flagvector == 0
            self.add_size_scaled_points(catalogue.data['longitude'][idx],
                                        catalogue.data['latitude'][idx],
                                        catalogue.data['magnitude'][idx],
                                        shape="o",
                                        colour="r",
                                        smin=1.0,
                                        sscale=1.5,
                                        overlay=overlay)
            return
        if not isinstance(cluster_id, collections.Iterable):
            cluster_id = [cluster_id]
        for iloc, clid in enumerate(cluster_id):
            if iloc == (len(cluster_id) - 1):
                # On last iteration set overlay to function overlay
                temp_overlay = overlay
            else:
                temp_overlay = True
            idx = vcl == clid
            self.add_size_scaled_points(
                catalogue.data["longitude"][idx],
                catalogue.data["latitude"][idx],
                catalogue.data["magnitude"][idx],
                shape="o",
                colour=DISSIMILAR_COLOURLIST[(iloc + 1) % NCOLS],
                smin=1.0,
                sscale=1.5,
                overlay=temp_overlay)
416481
def output_clusters(labels, cluster_centers_indices):
    
            
    here = os.getcwd()
    try:
        output_directory = os.path.join(here, 'concurrent_AP_output')
        os.makedirs(output_directory)
    except OSError:
        if not os.path.isdir(output_directory):
            print("ERROR: concurrent_AP: output_clusters: cannot create a directory "
                  "for storage of the results of Affinity Propagation clustering "
                  "in your current working directory")
            sys.exit(1)
            
    if any(np.isnan(labels)):
        fmt = '%.1f'
    else:
        fmt = '%d'
                
    with open(os.path.join(output_directory, 'labels.tsv'), 'w') as fh:
        np.savetxt(fh, labels, fmt = fmt, delimiter = '\t')
        
    if cluster_centers_indices is not None:    
        with open(os.path.join(output_directory, 'cluster_centers_indices.tsv'), 'w') as fh:
            np.savetxt(fh, cluster_centers_indices, fmt = '%.1f', 
                       delimiter = '\t')
575256
def _get_gosubdagplotnts(self, hdrgo, usrgos, pltargs, go2parentids):
        """"""
        dotgraphs = []
        go2color = pltargs.get_go2color_inst(hdrgo)
        # namedtuple fields: hdrgo gosubdag tot_usrgos parentcnt desc
        ntpltgo0 = self._get_pltdag_ancesters(hdrgo, usrgos, desc="")
        ntpltgo1 = self._get_pltdag_path_hdr(hdrgo, usrgos, desc="pruned")
        num_go0 = len(ntpltgo0.gosubdag.go2obj)
        num_go1 = len(ntpltgo1.gosubdag.go2obj)
        title = "{GO} {NAME}; {N} GO sources".format(
            GO=hdrgo,
            NAME=self.gosubdag.go2obj[hdrgo].name,
            N=len(ntpltgo0.gosubdag.go_sources))
        # print("PltGroupedGos::_get_gosubdagplotnts TITLE", title)
        if num_go0 < pltargs.max_gos:
            # print("PltGroupedGos::_get_gosubdagplotnts ntpltgo0 ALWAYS IF NOT TOO BIG")
            dotgraphs.append(self._get_gosubdagplotnt(ntpltgo0, title, go2color, pltargs))
        # PLOT A: Plot the entire GO ID group under GO header, hdrgo, if not too big
        if num_go0 < pltargs.max_gos and \
            ntpltgo0.tot_usrgos == ntpltgo1.tot_usrgos:
            # print("PltGroupedGos::_get_gosubdagplotnts ntpltgo0 AAAAAAAAAAAAAAAAA")
            dotgraphs.append(self._get_gosubdagplotnt(ntpltgo0, title, go2color, pltargs))
        # PLOT B: Plot only the GO ID group passing thru the GO header, hdrgo, if not too big
        elif num_go1 < pltargs.max_gos and ntpltgo0.tot_usrgos != ntpltgo1.tot_usrgos:
            # print("PltGroupedGos::_get_gosubdagplotnts ntpltgo1(pruned) BBBBBBBBBBBBBBBBB")
            dotgraphs.append(self._get_gosubdagplotnt(ntpltgo1, title, go2color, pltargs))
        # PLOT C: If DAG is large, just print upper portion
        # PLOT D: If DAG is large, just print upper portion passing through the GO header
        elif num_go1 >= pltargs.upper_trigger:
            # print("PltGroupedGos::_get_gosubdagplotnts upper_pruned CCCCCCCCCCCCCCCCC")
            gos_upper = self._get_gos_upper(ntpltgo1, pltargs.max_upper, go2parentids)
            #ntpltgo2 = self._get_pltdag_ancesters(hdrgo, gos_upper, "{BASE}_upper.png")
            ntpltgo3 = self._get_pltdag_path_hdr(hdrgo, gos_upper, "upper_pruned")
            # Middle GO terms chosen to be start points will be green unless reset back
            for goid in gos_upper:
                if goid not in go2color:
                    go2color[goid] = 'white'
            dotgraphs.append(self._get_gosubdagplotnt(ntpltgo3, title, go2color, pltargs))
        else:
            # print("PltGroupedGos::_get_gosubdagplotnts EEEEEEEEEEEEEEEEE")
            self._no_ntplt(ntpltgo0)
        return dotgraphs
254207
def _bisect_kmeans(X, n_clusters, n_trials, max_iter, tol):
    
    membs = np.empty(shape=X.shape[0], dtype=int)
    centers = dict() #np.empty(shape=(n_clusters,X.shape[1]), dtype=float)
    sse_arr = dict() #-1.0*np.ones(shape=n_clusters, dtype=float)

    ## data structure to store cluster hierarchies
    tree = treelib.Tree()
    tree = _add_tree_node(tree, 0, ilev=0, X=X) 

    km = _kmeans.KMeans(n_clusters=2, n_trials=n_trials, max_iter=max_iter, tol=tol)
    for i in range(1,n_clusters):
        sel_clust_id,sel_memb_ids = _select_cluster_2_split(membs, tree)
        X_sub = X[sel_memb_ids,:]
        km.fit(X_sub)

        #print("Bisecting Step %d    :"%i, sel_clust_id, km.sse_arr_, km.centers_)
        ## Updating the clusters & properties
        #sse_arr[[sel_clust_id,i]] = km.sse_arr_
        #centers[[sel_clust_id,i]] = km.centers_
        tree = _add_tree_node(tree, 2*i-1, i, \
                              size=np.sum(km.labels_ == 0), center=km.centers_[0], \
                              sse=km.sse_arr_[0], parent= sel_clust_id)
        tree = _add_tree_node(tree, 2*i,   i, \
                             size=np.sum(km.labels_ == 1), center=km.centers_[1], \
                             sse=km.sse_arr_[1], parent= sel_clust_id)

        pred_labels = km.labels_
        pred_labels[np.where(pred_labels == 1)[0]] = 2*i
        pred_labels[np.where(pred_labels == 0)[0]] = 2*i - 1
        #if sel_clust_id == 1:
        #    pred_labels[np.where(pred_labels == 0)[0]] = sel_clust_id
        #    pred_labels[np.where(pred_labels == 1)[0]] = i
        #else:
        #    pred_labels[np.where(pred_labels == 1)[0]] = i
        #    pred_labels[np.where(pred_labels == 0)[0]] = sel_clust_id

        membs[sel_memb_ids] = pred_labels


    for n in tree.leaves():
        label = n.data['label']
        centers[label] = n.data['center']
        sse_arr[label] = n.data['sse']

    return(centers, membs, sse_arr, tree)
1262481
def cluster(self, dist_type='cosine', run_clustering=True,
                 dendro=True, views=['N_row_sum', 'N_row_var'],
                 linkage_type='average', sim_mat=False, filter_sim=0.1,
                 calc_cat_pval=False, run_enrichr=None, enrichrgram=None):
    '''
    
    '''
    initialize_net.viz(self)

    make_clust_fun.make_clust(self, dist_type=dist_type, run_clustering=run_clustering,
                                   dendro=dendro,
                                   requested_views=views,
                                   linkage_type=linkage_type,
                                   sim_mat=sim_mat,
                                   filter_sim=filter_sim,
                                   calc_cat_pval=calc_cat_pval,
                                   run_enrichr=run_enrichr,
                                   enrichrgram=enrichrgram)
1160792
def read_clusters(self, min_cluster_size):
        """"""
        num_found = 0
        clusters = []
        with open(self.get_path(OslomRunner.OUTPUT_FILE), "r") as reader:
            # Read the output file every two lines
            for line1, line2 in itertools.izip_longest(*[reader] * 2):
                info = OslomRunner.RE_INFOLINE.match(line1.strip()).groups()
                nodes = line2.strip().split(" ")
                if len(nodes) >= min_cluster_size: # Apply min_cluster_size
                    clusters.append({
                        "id": int(info[0]),
                        "bs": float(info[2]),
                        "nodes": [{"id": self.id_remapper.get_str_id(int(n))} for n in nodes],
                    })
                num_found += 1
        return {"num_found": num_found, "clusters": clusters}
1140107
def _group_cluster(self, clusters, adj_list, counts):
        ''' '''

        groups = []
        for cluster in clusters:
            groups.append(sorted(cluster, key=lambda x: counts[x],
                                 reverse=True))

        return groups
509832
def create_gp(self):
        """"""

        categories_order = [
            ("{U}", "#ee82ee", 'Unmapped correctly'),
            ("{u}", "#ff0000", 'Unmapped incorrectly'),
            ("{T}", "#00ff00", 'Thresholded correctly'),
            ("{t}", "#008800", 'Thresholded incorrectly'),
            ("{P}", "#ffff00", 'Multimapped'),
            ("{w}+{x}", "#7f7f7f", 'Mapped, should be unmapped'),
            ("{m}", "#000000", 'Mapped to wrong position'),
            ("{M}", "#0000ff", 'Mapped correctly'),
        ]

        plot_lines = [
            '"{roc_fn}" using (( ({x}) )):({y}) lt rgb "{color}" with filledcurve x1 title "{title}", \\'.format(
                roc_fn=self._roc_fn,
                x=rnftools.lavender._format_xxx(self.default_x_axis),
                y=rnftools.lavender._format_xxx(
                    '({sum})*100/{{all}}'.format(sum="+".join([c[0] for c in categories_order[i:]]))
                ),
                color=categories_order[i][1],
                title=categories_order[i][2],
            ) for i in range(len(categories_order))
        ]

        plot = os.linesep.join((["plot \\"] + plot_lines + [""]))

        with open(self._gp_fn, "w+") as gp:
            gp_content = """
					set title "{{/:Bold=16 {title}}}"
	
					set x2lab "{x_lab}"
					set log x
					set log x2
	
					set format x "10^{{%L}}"
					set format x2 "10^{{%L}}"
					set xran  [{xran}]
					set x2ran [{xran}]
					set x2tics
					unset xtics
	
					set ylab "Part of all reads (%)"
	
					set format y "%g %%"
					set yran [{yran}]
					set y2ran [{yran}]
	
					set pointsize 1.5
	
					set grid ytics lc rgb "#777777" lw 1 lt 0 front
					set grid x2tics lc rgb "#777777" lw 1 lt 0 front
	
					set datafile separator "\\t"
					set palette negative
	
					set termin svg size {svg_size} enhanced
					set out "{svg_fn}"
					set key spacing 0.8 opaque width -5
				""".format(
                svg_fn=self._svg_fn,
                xran="{:.10f}:{:.10f}".format(self.report.default_x_run[0], self.report.default_x_run[1]),
                yran="{:.10f}:{:.10f}".format(self.report.default_y_run[0], self.report.default_y_run[1]),
                svg_size="{},{}".format(self.report.default_svg_size_px[0], self.report.default_svg_size_px[1]),
                title=os.path.basename(self._bam_fn)[:-4],
                x_lab=self.default_x_label,
            )
            gp_content = textwrap.dedent(gp_content) + "\n" + plot
            # gp_lines=gp_content.split("\n")
            # gp_lines=[x.strip() for x in gp_lines]
            # gp_content=gp_lines.join("\n")
            gp.write(gp_content)
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
38686
def _assign_clusters(self):
        
        self.clusters = np.array([self._closest_centroid(x) for x in self._X])
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)
1078687
def _single_cluster(c, data, out_file, args):
    
    valid, ann = 0, 0
    raw_file = None
    freq = defaultdict()
    [freq.update({s.keys()[0]: s.values()[0]}) for s in data[0][c]['freq']]
    names = [s.keys()[0] for s in data[0][c]['seqs']]
    seqs = [s.values()[0] for s in data[0][c]['seqs']]
    loci = data[0][c]['loci']

    if loci[0][3] - loci[0][2] > 500:
        logger.info("locus bigger > 500 nt, skipping: %s" % loci)
        return valid, ann, {}
    if not file_exists(out_file):
        if args.razer:
            logger.debug("map with razer all sequences to all loci %s " % loci)
            map_to_precursors(seqs, names, {loci[0][0]: [loci[0][0:5]]}, out_file, args)
        else:
            logger.debug("map with biopython fn all sequences to all loci %s " % loci)
            if args.debug:
                raw_file = out_file
            out_file = map_to_precursor_biopython(seqs, names, loci[0][0:5], args)

    logger.debug("plot sequences on loci")
    df = _convert_to_df(out_file, freq, raw_file)
    if df:
        valid, ann = _make(data[0][c])

    return valid, ann, df
796480
def run_tutorial(plot=False, multiplex=True, return_streams=False, cores=4,
                 verbose=False):
    
    client = Client("GEONET", debug=verbose)
    cat = client.get_events(
        minlatitude=-40.98, maxlatitude=-40.85, minlongitude=175.4,
        maxlongitude=175.5, starttime=UTCDateTime(2016, 5, 1),
        endtime=UTCDateTime(2016, 5, 20))
    print("Downloaded a catalog of %i events" % len(cat))
    # This gives us a catalog of events - it takes a while to download all
    # the information, so give it a bit!
    # We will generate a five station, multi-channel detector.
    cat = filter_picks(catalog=cat, top_n_picks=5)
    stachans = list(set(
        [(pick.waveform_id.station_code, pick.waveform_id.channel_code)
         for event in cat for pick in event.picks]))
    # In this tutorial we will only work on one cluster, defined spatially.
    # You can work on multiple clusters, or try to whole set.
    clusters = space_cluster(catalog=cat, d_thresh=2, show=False)
    # We will work on the largest cluster
    cluster = sorted(clusters, key=lambda c: len(c))[-1]
    # This cluster contains 32 events, we will now download and trim the
    # waveforms.  Note that each chanel must start at the same time and be the
    # same length for multiplexing.  If not multiplexing EQcorrscan will
    # maintain the individual differences in time between channels and delay
    # the detection statistics by that amount before stacking and detection.
    client = Client('GEONET')
    design_set = []
    st = Stream()
    for event in cluster:
        print("Downloading for event {0}".format(event.resource_id.id))
        bulk_info = []
        t1 = event.origins[0].time
        t2 = t1 + 25.1  # Have to download extra data, otherwise GeoNet will
        # trim wherever suits.
        t1 -= 0.1
        for station, channel in stachans:
            bulk_info.append(('NZ', station, '*', channel[0:2] + '?', t1, t2))
        st += client.get_waveforms_bulk(bulk=bulk_info)
    print("Downloaded %i channels" % len(st))
    for event in cluster:
        t1 = event.origins[0].time
        t2 = t1 + 25
        design_set.append(st.copy().trim(t1, t2))
    # Construction of the detector will process the traces, then align them,
    # before multiplexing.
    print("Making detector")
    detector = subspace.Detector()
    detector.construct(
        streams=design_set, lowcut=2.0, highcut=9.0, filt_order=4,
        sampling_rate=20, multiplex=multiplex, name='Wairarapa1', align=True,
        reject=0.2, shift_len=6, plot=plot).partition(9)
    print("Constructed Detector")
    if plot:
        detector.plot()
    # We also want the continuous stream to detect in.
    t1 = UTCDateTime(2016, 5, 11, 19)
    t2 = UTCDateTime(2016, 5, 11, 20)
    # We are going to look in a single hour just to minimize cost, but you can
    # run for much longer.
    bulk_info = [('NZ', stachan[0], '*',
                  stachan[1][0] + '?' + stachan[1][-1],
                  t1, t2) for stachan in detector.stachans]
    print("Downloading continuous data")
    st = client.get_waveforms_bulk(bulk_info)
    st.merge().detrend('simple').trim(starttime=t1, endtime=t2)
    # We set a very low threshold because the detector is not that great, we
    # haven't aligned it particularly well - however, at this threshold we make
    # two real detections.
    print("Computing detections")
    detections, det_streams = detector.detect(
        st=st, threshold=0.3, trig_int=2, extract_detections=True,
        cores=cores)
    if return_streams:
        return detections, det_streams
    else:
        return detections
471086
def cluster(list_of_texts, num_clusters=3):
    
    pipeline = Pipeline([
        ("vect", CountVectorizer()),
        ("tfidf", TfidfTransformer()),
        ("clust", KMeans(n_clusters=num_clusters))
    ])

    try:
        clusters = pipeline.fit_predict(list_of_texts)
    except ValueError:
        clusters = list(range(len(list_of_texts)))

    return clusters
422536
def _plot_stars(catalog, observer, project, ax, mag1, mag2, margin=1.25):
    """"""

    art = []

    # from astropy import wcs
    # w = wcs.WCS(naxis=2)
    # w.wcs.crpix = [-234.75, 8.3393]
    # w.wcs.cdelt = np.array([-0.066667, 0.066667])
    # w.wcs.crval = [0, -90]
    # w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
    # w.wcs.set_pv([(2, 1, 45.0)])

    # import matplotlib.pyplot as plt

    # plt.subplot(projection=wcs)
    # #plt.imshow(hdu.data, vmin=-2.e-5, vmax=2.e-4, origin='lower')
    # plt.grid(color='white', ls='solid')
    # plt.xlabel('Galactic Longitude')
    # plt.ylabel('Galactic Latitude')

    xmin, xmax = ax.get_xlim()
    ymin, ymax = ax.get_xlim()
    lim = max(abs(xmin), abs(xmax), abs(ymin), abs(ymax)) * margin
    lims = (-lim, lim)
    ax.set_xlim(lims)
    ax.set_ylim(lims)
    ax.set_aspect('equal')

    o = observer[0]

    # Dim stars: points of with varying gray levels.

    c = catalog
    c = c[c['magnitude'] > mag1]
    c = c[c['magnitude'] <= mag2]
    #print('Second star group:', len(c))
    c = c.sort_values('magnitude', ascending=False)
    s = Star(ra_hours=c.ra_hours, dec_degrees=c.dec_degrees)
    spos = o.observe(s)
    x, y = project(spos)
    m = (mag2 - c['magnitude']) / (mag2 - mag1)
    # Note that "gray_r" is white for 0.0 and black for 1.0
    art.append(ax.scatter(
        x, y, s=1.0,
        c=1 - 0.8 * m, cmap='gray_r', vmin=0.0, vmax=1.0,
    ))

    # Bright stars: black circles of varying radius, surrounded by a
    # white gap in case stars are touching.  Draw the brightest stars
    # first to stop them from completely occluding smaller companions.

    def mag_to_radius(m):
        return (mag1 - m) * scale + 1.0

    c = catalog
    c = c[c['magnitude'] <= mag1]
    c = c.sort_values('magnitude', ascending=True)
    #print('First star group:', len(c))
    s = Star(ra_hours=c.ra_hours, dec_degrees=c.dec_degrees)
    spos = o.observe(s)
    x, y = project(spos)
    scale = 1.5
    radius = mag_to_radius(c['magnitude'])

    x2 = np.repeat(x, 2)
    y2 = np.repeat(y, 2)
    radius2 = (radius[:,None] + (3.0,0.0)).flatten()
    c2 = ('w', 'k')
    c2 = ('k', 'w')

    art.append(ax.scatter(x2, y2, s=radius2 ** 2.0, c=c2))

    return art, mag_to_radius
831646
def get_cluster_labels(hdf5_file, N_processes):
    

    with Worker.hdf5_lock:
        with tables.open_file(hdf5_file, 'r+') as fileh:
            A = fileh.root.aff_prop_group.availabilities
            R = fileh.root.aff_prop_group.responsibilities
    
            N = A.nrows
            diag_ind = np.diag_indices(N)
            
            a = A[diag_ind]
            r = R[diag_ind]
    
    I = np.where(np.add(a[:], r[:]) > 0)[0]
    K = I.size
    
    if K == 0:
        labels = np.empty((N, 1))
        labels.fill(np.nan)
        cluster_centers_indices = None
        
    else:
        c_array = multiprocessing.Array(c_int, N, lock = True)
    
        chunk_size = get_chunk_size(N, 3 * N_processes)
        numpy_args = c_array, N, np.int32
    
        with closing(multiprocessing.Pool(N_processes, 
                     initializer = cluster_labels_init, 
                     initargs = (hdf5_file, I, c_array.get_lock()) 
                               + numpy_args)) as pool:
            pool.map_async(multiprocessing_cluster_labels_A, 
                  chunk_generator(N, 3 * N_processes), chunk_size)
        
        pool.close()
        pool.join()
        
        gc.collect()
        
        c = to_numpy_array(*numpy_args)
        c[I] = np.arange(K)
        # determine the exemplars of clusters, applying some 
        # cosmetics to our results before returning them
        for k in range(K):
            ii = np.where(c == k)[0]
            
            iix = np.full(N, False, dtype = bool)
            iix[ii] = True
            
            s_reduced_array = multiprocessing.Array(c_double, ii.size, 
                                                    lock = True)
        
            chunk_size = get_chunk_size(N, 3 * N_processes)
            numpy_args = s_reduced_array, ii.size, np.float64
            
            with closing(multiprocessing.Pool(N_processes,
                         initializer = cluster_labels_init_B,
                         initargs = (hdf5_file, I, ii, iix, 
                         s_reduced_array.get_lock())
                         + numpy_args)) as pool:
                pool.map_async(multiprocessing_cluster_labels_B,
                      chunk_generator(N, 3 * N_processes), chunk_size)
                               
            pool.close()
            pool.join()
            
            s_reduced = to_numpy_array(*numpy_args)
            
            j = np.argmax(s_reduced)
            I[k] = ii[j]
            
            gc.collect()
            
        c_array = multiprocessing.Array(c_int, N, lock = True)
        
        chunk_size = get_chunk_size(N, 3 * N_processes)
        numpy_args = c_array, N, np.int32
    
        with closing(multiprocessing.Pool(N_processes, 
                     initializer = cluster_labels_init, 
                     initargs = (hdf5_file, I, c_array.get_lock()) 
                               + numpy_args)) as pool:
            pool.map_async(multiprocessing_cluster_labels_C, 
                  chunk_generator(N, 3 * N_processes), chunk_size)
        
        pool.close()
        pool.join()
        
        c = to_numpy_array(*numpy_args)
        c[I] = np.arange(K)
        labels = I[c]
        
        gc.collect()

        cluster_centers_indices = np.unique(labels)
        labels = np.searchsorted(cluster_centers_indices, labels)
        
    return cluster_centers_indices, labels
57761
def plot_results(fout_png, goea_results, *args, **kws):
    """"""
    if "{NS}" not in fout_png:
        plt_goea_results(fout_png, goea_results, *args, **kws)
    else:
        # Plot separately by NS: BP, MF, CC
        ns2goea_results = cx.defaultdict(list)
        for rec in goea_results:
            ns2goea_results[rec.NS].append(rec)
        for ns_name, ns_res in ns2goea_results.items():
            png = fout_png.format(NS=ns_name)
            plt_goea_results(png, ns_res, *args, **kws)