text:
from __future__ import print_function
import argparse
import numpy as np
import datetime
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import glob
import sys
import os
import re
import concurrent.futures
import warnings
import subprocess
from collections import defaultdict
from operator import itemgetter
from argparse import RawTextHelpFormatter
from io import StringIO

# Plotting imports
import bokeh.io as bk_io
from bokeh.plotting import figure, output_file, show, save
from bokeh.layouts import gridplot
from bokeh.models import ColumnDataSource, CDSView, BooleanFilter
from bokeh.transform import factor_cmap, factor_mark
from bokeh.palettes import Spectral6, Dark2, cividis, Viridis, viridis

def get_bedgraph_data(params):

    (bedgraph, genome) = params
    subsample_blacklist_vals =  ['_001', '_002', '_003', '_004', '_005', '_006', '_007', '_008', '_009', \
                                 '_01', '_02', '_03', '_04', '_05', '_06', '_07', '_08', '_09', \
                                 '_1', '_2', '_3', '_4', '_5', '_6', '_7', '_8', '_9']
    sample_name = os.path.splitext(os.path.basename(bedgraph))[0]
    splt = sample_name.split('_')
    sample_id = '_'.join(splt[:-1]) if '_{}'.format(splt[-1:][0]) in subsample_blacklist_vals else '_'.join(splt)

    coverage_file = pd.read_csv(bedgraph, comment='#', names=['chr','start','stop','coverage'], sep='\t', header=None \
                                 ).assign(sample=sample_name)
    coverage_variance = coverage_file['coverage'].var(ddof=1)

# Made a mistake previously where I did not take the abs -- resulted in linear relationship -- should explore this
    log_coverage_variance = (np.log2(abs(coverage_file['coverage']))).var(ddof=1)

    complement = subprocess.Popen(['complementBed', '-i', bedgraph, '-g', genome],
               stdout=subprocess.PIPE,
               stderr=subprocess.STDOUT)
    compout,comperr = complement.communicate()
    compdata = StringIO(str(compout,'utf-8'))
    gaps=pd.read_csv(compdata, sep='\t', header=None, \
                  names = ['chr','start','end'])

    gaps['gap_size'] = gaps.end - gaps.start

    gap_std = gaps['gap_size'].std(ddof=1)
    log_gap_std = np.log2(gaps['gap_size']).std(ddof=1)
    gap_var = gaps['gap_size'].var(ddof=1)
    log_gap_var = np.log2(gaps['gap_size']).var(ddof=1)
    gap_median_size = gaps['gap_size'].median()
    log_gap_median_size = np.log2((gaps['gap_size']).median())
    number_of_1bp_gaps = (gaps['gap_size'] == 0).sum()
    log_number_of_1bp_gaps = np.log2((gaps['gap_size'] == 0).sum())

    coverage_stats = { 'sample': sample_name, \
                       'sample_id': sample_id, \
                       'gap_std': float(gap_std), \
                       'log_gap_std': float(log_gap_std), \
                       'gap_var': float(gap_var), \
                       'log_gap_var': float(log_gap_var), \
                       'gap_median_size': int(gap_median_size), \
                       'log_gap_median_size': float(log_gap_median_size), \
                       'coverage_variance': float(coverage_variance), \
                       'log_coverage_variance': float(log_coverage_variance), \
                       'number_of_1bp_gaps': int(number_of_1bp_gaps), \
                       'log_number_of_1bp_gaps': float(log_number_of_1bp_gaps) }

    print("Done with sample %s bedGraph." % sample_name)

    return coverage_stats

def get_mapped_complexity_data(dup_file):

    sample_name = os.path.splitext(os.path.basename(dup_file))[0]
    if sample_name.endswith('.marked_dup_metrics'):
        sample_name = sample_name[:-19]
    dup_file_df = pd.read_csv(dup_file, comment='#', sep='\t', usecols = range(1,9), skiprows=[6], header=None, \
                                names=['unpaired_reads','paired_reads','secondary_reads', \
                                'unmapped_reads', 'unpaired_read_dups', 'paired_read_dups', \
                                'read_pair_optical_dups', 'percent_dup'] \
                                 ).assign(sample=sample_name)

    print("Done with sample %s mapped duplication file." % sample_name)

    dup_stats = { 'sample': sample_name, \
                         'unpaired_reads': int(dup_file_df.iloc[0]['unpaired_reads']), \
                         'paired_reads': int(dup_file_df.iloc[0]['paired_reads']), \
                         'secondary_reads': int(dup_file_df.iloc[0]['secondary_reads']), \
                         'unmapped_reads': int(dup_file_df.iloc[0]['unmapped_reads']), \
                         'unpaired_read_dups': int(dup_file_df.iloc[0]['unpaired_read_dups']), \
                         'paired_read_dups': int(dup_file_df.iloc[0]['paired_read_dups']), \
                         'read_pair_optical_dups': int(dup_file_df.iloc[0]['read_pair_optical_dups']), \
                         'percent_dup': float(dup_file_df.iloc[0]['percent_dup']) }
    return dup_stats

def get_read_complexity_data(dup_file):

    sample_name = os.path.splitext(os.path.basename(dup_file))[0]
    if sample_name.endswith('.fastqc_stats'):
        sample_name = sample_name[:-13]
    read_dup_file_df = pd.read_csv(dup_file, comment='#', sep='\t')

    print("Done with sample %s read duplication file." % sample_name)

    read_dup_stats = { 'sample': read_dup_file_df.iloc[0]['SRR'], \
                         'percent_gc': int(read_dup_file_df.iloc[0]['%GC']), \
                         'total_sequences': int(read_dup_file_df.iloc[0]['Total_Sequences']), \
                         'percent_deduplicated': float(read_dup_file_df.iloc[0]['%Total_Deduplicated']) }
    return read_dup_stats

def main():
    parser = argparse.ArgumentParser(description='Nascent QC\n\n===============================================================================================================\n\nAn algorithm which provides a metric for nascent data quality.', epilog='@Dowell Lab, Margaret Gruca, margaret.gruca@colorado.edu\nFor questions and issues, see https://github.com/Dowell-Lab/nqc', usage='%(prog)s --bedgraph /my/sample/dir/*.bedGraph --output /my/out/dir/', formatter_class=RawTextHelpFormatter)

    required = parser.add_argument_group('Required Arguments')
    optional = parser.add_argument_group('Optional Arguments')

    required.add_argument('-b', '--bedgraph', dest='bedgraphdir', \
                        help='Path to bedgraph file(s) of interest. Wildcard may be used to process multiple bedGraphs simultaneously with one output report.', required=True)

    required.add_argument('-d', '--mapped-duplicates', dest='mappeddupdir', \
                        help='Path to picard duplication summary file(s). Wildcard may be used to process multiple duplication files simultaneously.', required=True)

    required.add_argument('-r', '--read-duplicates', dest='readdupdir', \
                        help='Path to fastqc duplication summary file(s). Wildcard may be used to process multiple duplication files simultaneously.', required=True)

    required.add_argument('-o', '--output', dest='output', \
                        help='Path to where output stats file and plots will be saved.', required=False, default='./')

    required.add_argument('-g', '--genome', dest='genome', \
                        help='Chromosome sizes file -- bedGraph file sort must match genome file. See README for details on generating this file for your specific genome.', \
                        default='', required=False, type=str)

    optional.add_argument('-png', '--save-png', dest='save_png', action='store_true', \
                        help='Save png in addition to the default html output.', default=False, required=False)

    optional.add_argument('-t', '--threads', type=int, dest='threads', metavar='<THREADS>', \
                        help='Number of threads for multi-processing. Default=1', default=1, required=False)

    args = parser.parse_args()

    genome = args.genome
    warnings.simplefilter(action='ignore', category=(FutureWarning,RuntimeWarning))

    threads = args.threads
    bedgraph_list = glob.glob(args.bedgraphdir + '*.bedGraph')
    mapped_duplication_file_list = glob.glob(args.mappeddupdir + '*.marked_dup_metrics.txt')
    read_duplication_file_list = glob.glob(args.readdupdir + '*.fastqc_stats.txt')

    #subsample_blacklist_vals =  ['_001', '_002', '_003', '_004', '_005', '_006', '_007', '_008', '_009', \
    #                             '_01', '_02', '_03', '_04', '_05', '_06', '_07', '_08', '_09', \
    #                             '_1', '_2', '_3', '_4', '_5', '_6', '_7', '_8', '_9']

    print("Starting to process bedGraphs...\n" + str(datetime.datetime.now()))

    with concurrent.futures.ProcessPoolExecutor(threads) as executor:
        jobs = [executor.submit(get_bedgraph_data, [file,genome])
                for file in bedgraph_list]
        coverage_stats = [r.result() for r in jobs]

    with concurrent.futures.ProcessPoolExecutor(threads) as executor:
        jobs = [executor.submit(get_mapped_complexity_data, file)
                for file in mapped_duplication_file_list]
        mapped_duplication_stats = [r.result() for r in jobs]

    with concurrent.futures.ProcessPoolExecutor(threads) as executor:
        jobs = [executor.submit(get_read_complexity_data, file)
                for file in read_duplication_file_list]
        read_duplication_stats = [r.result() for r in jobs]

    d = defaultdict(dict)
    for l in (coverage_stats, mapped_duplication_stats):
        for elem in l:
            d[elem['sample']].update(elem)
    mapdup_bg_stats = d.values()

    e = defaultdict(dict)
    for l in (mapdup_bg_stats, read_duplication_stats):
        for elem in l:
            e[elem['sample']].update(elem)
    nqc_stats = e.values()

    sorted_nqc_stats = sorted(nqc_stats, key=itemgetter('sample'))

    nqc_stats_file = open(args.output + '/nqc_stats.txt', 'w')
    for stat in sorted_nqc_stats:
        nqc_stats_file.write("%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % \
                        (stat['sample'], \
                         stat['sample_id'], \
                         stat['gap_std'], stat['log_gap_std'], \
                         stat['gap_var'], stat['log_gap_var'], \
                         stat['gap_median_size'], stat['log_gap_median_size'], \
                         stat['coverage_variance'], stat['log_coverage_variance'], \
                         stat['number_of_1bp_gaps'], stat['log_number_of_1bp_gaps'], \
                         stat['unpaired_reads'], stat['paired_reads'], stat['secondary_reads'], \
                         stat['unmapped_reads'], stat['unpaired_read_dups'], \
                         stat['paired_read_dups'], stat['read_pair_optical_dups'], \
                         stat['percent_dup'], stat['percent_gc'], stat['total_sequences'], stat['percent_deduplicated']))
    nqc_stats_file.close()

    #sample = [d['sample'] for d in sorted_nqc_stats]

    #identifier = []
    #for i in sample:
    #    splt = i.split('_')
    #    value = '_'.join(splt[:-1]) if '_{}'.format(splt[-1:][0]) in subsample_blacklist_vals else '_'.join(splt)
    #    identifier.append(value)
    #
    #df["id"] = identifier

    df = pd.DataFrame(sorted_nqc_stats).replace([np.inf, -np.inf], np.nan).fillna(0)
    print(list(df.columns.values))
    df['log2_unique_mapped_reads'] = np.log2(df['unpaired_reads'] - df['unpaired_read_dups'])
    df['percent_raw_read_dup'] = 1 - df['percent_deduplicated']
    df['log2_unique_raw_reads'] = np.log2(df['total_sequences'] * df['percent_deduplicated'])
    df['log2_unmapped_reads'] = np.log2(df['unmapped_reads'])

    output_file(args.output + '/nqc_plot.html')
    source = ColumnDataSource(df)
    SAMPLE = df['sample_id'].unique()
    pal=viridis(len(SAMPLE))

    # create a view of the source for one plot to use
    #view = CDSView(source=source, filters=[BooleanFilter([True if y > 250 or y < 100 else False for y in y1])])

    TOOLS = "box_select,lasso_select,hover,help,save"
    #colors = factor_cmap('sample_id', palette=Category20b_20, factors=sample_id.unique())

    # create a new plot and add a renderer
    s1 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s1.xaxis.axis_label = 'Log\u2082 Number of Unique Reads'
    s1.yaxis.axis_label = 'Log\u2082 Gap Standard Deviation'
    s1.circle('log2_unique_mapped_reads', 'log_gap_std', size=6, hover_color="firebrick", legend='sample_id',
                color=factor_cmap('sample_id', pal, SAMPLE),
                source=source)

    s2 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s2.xaxis.axis_label = 'Log\u2082 Number of Unique Reads'
    s2.yaxis.axis_label = "Log\u2082 Gap Variance"
    s2.circle('log2_unique_mapped_reads', 'log_gap_var', size=6, hover_color="firebrick", legend='sample_id',
                  color=factor_cmap('sample_id', pal, SAMPLE),
                  source=source)

    s3 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s3.xaxis.axis_label = 'Log\u2082 Number of Unique Reads'
    s3.yaxis.axis_label = "Log\u2082 Gap Median Size"
    s3.circle('log2_unique_mapped_reads', 'log_gap_median_size', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)

    s4 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s4.xaxis.axis_label = 'Log\u2082 Unmapped Reads'
    s4.yaxis.axis_label = "Log\u2082 Coverage Variance"
    s4.circle('log2_unmapped_reads', 'log_coverage_variance', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)

    s5 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s5.xaxis.axis_label = 'Log\u2082 Number of Unique Reads'
    s5.yaxis.axis_label = "Log\u2082 Coverage Variance"
    s5.circle('log2_unique_mapped_reads', 'log_coverage_variance', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)

    s6 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s6.xaxis.axis_label = 'Log\u2082 Number of Unique Reads'
    s6.yaxis.axis_label = "Log\u2082 Number of 1bp Gaps"
    s6.circle('log2_unique_raw_reads', 'log_number_of_1bp_gaps', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)
    s6.legend.location = "bottom_right"

    s7 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s7.xaxis.axis_label = 'Log\u2082 Number of Unique Raw Reads'
    s7.yaxis.axis_label = "Log\u2082 Gap Median Size"
    s7.circle('log2_unique_raw_reads', 'log_gap_median_size', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)

    s8 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s8.xaxis.axis_label = 'Log\u2082 Number of Unique Raw Reads'
    s8.yaxis.axis_label = "Log\u2082 Coverage Variance"
    s8.circle('log2_unique_raw_reads', 'log_coverage_variance', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)

    s9 = figure(tools=TOOLS, plot_width=600, plot_height=600, title=None)
    s9.xaxis.axis_label = 'Log\u2082 Number of Unique Raw Reads'
    s9.yaxis.axis_label = "Log\u2082 Number of Unique Mapped Reads"
    s9.circle('log2_unique_raw_reads', 'log2_unique_mapped_reads', size=6, hover_color="firebrick", legend='sample_id',
                    color=factor_cmap('sample_id', pal, SAMPLE),
                    source=source)
    s9.legend.location = "bottom_right"

    # make a grid
    grid = gridplot([s1, s2, s3, s4, s5, s6, s7, s8, s9], ncols=3, plot_width=600, plot_height=600)
    save(grid, browser=args.output + 'nqc_plot.html')
    if args.save_png:
        bk_io.export_png(grid, filename='./nqc_plot.png')

    #fig = plt.figure()
    #ax = fig.add_subplot(111)
    ##ax2 = fig.add_subplot(3, 3, 2)
    ## We use ax parameter to tell seaborn which subplot to use for this plot
    #sns.scatterplot(y='log_gap_std', x='log2_unique_reads', hue='id', data=df)
    #sns.set_style("white")
    #sns.set_style("ticks")
    #plt.tight_layout()
    ##ax.legend(loc='center left', bbox_to_anchor=(1, 0.5));
    #plt.savefig('./plot_test' + '.png')

    print("Done.\n" + str(datetime.datetime.now()))

    sys.exit(0)

if __name__ == '__main__':
    main()

max_stars_repo_name: Dowell-Lab/Nascent-Flow

max_issues_repo_name: Dowell-Lab/Nascent-Flow

max_forks_repo_name: Dowell-Lab/Nascent-Flow


text:
import os
import time
import torch
import torch.nn as nn
import utils
from torch.autograd import Variable
import datetime
import numpy as np
import sys

def LangMCriterion(input, target):
    target = target.view(-1, 1)
    logprob_select = torch.gather(input, 1, target)

    mask = target.data.gt(0)  # generate the mask
    if isinstance(input, Variable):
        mask = Variable(mask, volatile=input.volatile)
    out = torch.masked_select(logprob_select, mask)
    loss = -torch.sum(out)  # get the average loss.
    return loss

def train(model, train_loader, eval_loader, args):
    t = datetime.datetime.now()
    cur_time = '%s-%s-%s-%s-%s' % (t.year, t.month, t.day, t.hour, t.minute)
    save_path = os.path.join(args.output, cur_time)
    args.save_path = save_path
    utils.create_dir(save_path)
    optim = torch.optim.Adam(model.parameters(), lr=args.lr, betas=(args.beta1, 0.999))

    lr_default = args.lr if eval_loader is not None else 7e-4
    lr_decay_step = 2
    lr_decay_rate = .25
    lr_decay_epochs = range(10, 30, lr_decay_step) if eval_loader is not None else range(10, 20, lr_decay_step)
    gradual_warmup_steps = [0.5 * lr_default, 1.0 * lr_default, 1.5 * lr_default, 2.0 * lr_default]

    logger = utils.Logger(os.path.join(save_path, 'log.txt'))
    for arg in vars(args):
        logger.write('{:<20}: {}'.format(arg, getattr(args, arg)))
    best_eval_score = 0
    model.train()
    start_time = time.time()
    best_cnt = 0

    print('Training ... ')
    for epoch in range(args.epochs):

        total_loss = 0
        count = 0
        train_score = 0
        t = time.time()
        train_iter = iter(train_loader)

        # TODO: get learning rate
        # lr = adjust_learning_rate(optim, epoch, args.lr)
        if epoch < 4:
            optim.param_groups[0]['lr'] = gradual_warmup_steps[epoch]
            lr = optim.param_groups[0]['lr']
        elif epoch in lr_decay_epochs:
            optim.param_groups[0]['lr'] *= lr_decay_rate
            lr =  optim.param_groups[0]['lr']
        else:
            lr = optim.param_groups[0]['lr']
        iter_step = 0
        for i in range(len(train_loader)):
            average_loss_tmp = 0
            count_tmp = 0
            train_data = next(train_iter)
            image, image_id, history, question, answer, answerT, ans_len, ans_idx, ques_ori, opt, opt_len, opt_idx = train_data
            batch_size = question.size(0)
            image = image.view(image.size(0), -1, args.img_feat_size)
            img_input = Variable(image).cuda()
            for rnd in range(10):
                ques = question[:, rnd, :]
                his = history[:, :rnd + 1, :].clone().view(-1, args.his_length)
                ans = answer[:, rnd, :]
                tans = answerT[:, rnd, :]
                opt_ans = opt[:, rnd, :].clone().view(-1, args.ans_length)

                ques = Variable(ques).cuda().long()
                his = Variable(his).cuda().long()
                ans = Variable(ans).cuda().long()
                tans = Variable(tans).cuda().long()
                opt_ans = Variable(opt_ans).cuda().long()

                pred = model(img_input, ques, his, ans, tans, rnd + 1)
                loss = LangMCriterion(pred.view(-1, args.vocab_size), tans)
                loss = loss / torch.sum(tans.data.gt(0))

                loss.backward()
                nn.utils.clip_grad_norm(model.parameters(), 0.25)
                optim.step()
                model.zero_grad()

                average_loss_tmp += loss.data[0]
                total_loss += loss.data[0]
                count += 1
                count_tmp += 1
            sys.stdout.write('Training: Epoch {:d} Step {:d}/{:d}  \r'.format(epoch + 1, i + 1, len(train_loader)))
            if (i+1) % 50 == 0:
                average_loss_tmp /= count_tmp
                print("step {} / {} (epoch {}), g_loss {:.3f}, lr = {:.6f}".format(i + 1, len(train_loader), epoch + 1, average_loss_tmp, lr))
        iter_step += 1
        total_loss /= count
        logger.write('Epoch %d : learningRate %4f train loss %4f Time: %3f' % (epoch + 1, lr, total_loss, time.time() - start_time))
        model.eval()
        print('Evaluating ... ')
        start_time = time.time()
        rank_all = evaluate(model, eval_loader, args)
        R1 = np.sum(np.array(rank_all) == 1) / float(len(rank_all))
        R5 = np.sum(np.array(rank_all) <= 5) / float(len(rank_all))
        R10 = np.sum(np.array(rank_all) <= 10) / float(len(rank_all))
        ave = np.sum(np.array(rank_all)) / float(len(rank_all))
        mrr = np.sum(1 / (np.array(rank_all, dtype='float'))) / float(len(rank_all))
        logger.write('Epoch %d: mrr: %f R1: %f R5 %f R10 %f Mean %f time: %.2f' % (epoch + 1, mrr, R1, R5, R10, ave, time.time()-start_time))

        eval_score = mrr

        model_path = os.path.join(save_path, 'model_epoch_%d.pth' % (epoch + 1))
        torch.save({'epoch': epoch,
                    'args': args,
                    'model': model.state_dict()}, model_path)

        if eval_score > best_eval_score:
            model_path = os.path.join(save_path, 'best_model.pth')
            torch.save({'epoch': epoch,
                        'args': args,
                        'model': model.state_dict()}, model_path)
            best_eval_score = eval_score
            best_cnt = 0
        else:
            best_cnt = best_cnt + 1
            if best_cnt > 10:
                break
    return model

def evaluate(model, eval_loader, args, Eval=False):
    rank_all_tmp = []
    eval_iter = iter(eval_loader)
    step = 0
    for i in range(len(eval_loader)):
        eval_data = next(eval_iter)
        image, image_id, history, question, answer, answerT, questionL, opt_answer, \
        opt_answerT, answer_ids, answerLen, opt_answerLen = eval_data

        image = image.view(image.size(0), -1, args.img_feat_size)
        img_input = Variable(image).cuda()
        batch_size = question.size(0)
        for rnd in range(10):
            ques, tans = question[:, rnd, :], opt_answerT[:, rnd, :].clone().view(-1, args.ans_length)
            his = history[:, :rnd + 1, :].clone().view(-1, args.his_length)
            ans = opt_answer[:, rnd, :, :].clone().view(-1, args.ans_length)
            gt_id = answer_ids[:, rnd]

            ques = Variable(ques).cuda().long()
            tans = Variable(tans).cuda().long()
            his = Variable(his).cuda().long()
            ans = Variable(ans).cuda().long()
            gt_index = Variable(gt_id).cuda().long()

            pred = model(img_input, ques, his, ans, tans, rnd + 1, Training=False)
            logprob = - pred.permute(1, 0, 2).contiguous().view(-1, args.vocab_size)

            logprob_select = torch.gather(logprob, 1, tans.t().contiguous().view(-1, 1))

            mask = tans.t().data.eq(0)  # generate the mask
            if isinstance(logprob, Variable):
                mask = Variable(mask, volatile=logprob.volatile)
            logprob_select.masked_fill_(mask.view_as(logprob_select), 0)

            prob = logprob_select.view(args.ans_length, -1, 100).sum(0).view(-1, 100)

            for b in range(batch_size):
                gt_index.data[b] = gt_index.data[b] + b * 100

            gt_score = prob.view(-1).index_select(0, gt_index)
            sort_score, sort_idx = torch.sort(prob, 1)

            count = sort_score.lt(gt_score.view(-1, 1).expand_as(sort_score))
            rank = count.sum(1) + 1

            rank_all_tmp += list(rank.view(-1).data.cpu().numpy())
        step += 1
        if Eval:
            sys.stdout.write('Evaluating: {:d}/{:d}  \r'.format(i, len(eval_loader)))
            if (i+1) % 50 == 0:
                R1 = np.sum(np.array(rank_all_tmp) == 1) / float(len(rank_all_tmp))
                R5 = np.sum(np.array(rank_all_tmp) <= 5) / float(len(rank_all_tmp))
                R10 = np.sum(np.array(rank_all_tmp) <= 10) / float(len(rank_all_tmp))
                ave = np.sum(np.array(rank_all_tmp)) / float(len(rank_all_tmp))
                mrr = np.sum(1 / (np.array(rank_all_tmp, dtype='float'))) / float(len(rank_all_tmp))
                print('%d/%d: mrr: %f R1: %f R5 %f R10 %f Mean %f' % (i+1, len(eval_loader), mrr, R1, R5, R10, ave))

    return rank_all_tmp

ext: py

max_stars_repo_name: Unmesh-Kumar/DMRM

max_issues_repo_name: Unmesh-Kumar/DMRM

max_forks_repo_name: Unmesh-Kumar/DMRM

reason: import numpy


text:
from n2v import node2vec
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import networkx as nx
import pandas as pd
import numpy as np
from graphwave import graphwave
from graphwave.utils import utils
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans

def _get_components(networkXGraph):
    return pd.DataFrame([{"node": k, "component": v} for k, v in nx.get_node_attributes(
        G=networkXGraph, name='component').items()]).groupby('component')['node'].apply(list).to_dict()

def divide_chunks(l, n):

    # looping till length l
    for i in range(0, len(l), n):
        yield l[i:i + n]

def get_structural_signatures(networkXGraph, vocab_size=1, params={'num_kmeans_clusters': 4, "num_pca_components": 6, "num_batch":500}):
    """
    Get structural embeddings using GraphWave.

    Learns structural signatures of each node using GraphWave, and adds these learned structures to the node attributes on the NetworkX graph.

    Example:

    Attributes:
        networkXGraph (networkx.classes.graph.Graph): A graph containing nodes, edges, and their attributes. Nodes must include an attribute called 'component'.

    Returns: a NetworkX graph where structural embeddings are added as `structure` node attributes.

    """
    nb_clust = params['num_kmeans_clusters']
    n_components = params['num_pca_components']
    batch_size = params['num_batch']

    components = list(divide_chunks(list(_get_components(networkXGraph).values()),batch_size))

    heat_signatures = []
    nodes_list = []

    for n in components:
        #flatten list of lists
        nodes =  [item for sublist in n for item in sublist]
        subgraph = networkXGraph.subgraph(nodes)
        chi, heat_print, taus = graphwave.graphwave_alg(
            subgraph, np.linspace(0, 10, 20), taus='auto', verbose=True)
        if len(subgraph.nodes) < n_components:
            print("Omitting graph " + str(subgraph_id) + " with node count: " +
                str(len(subgraph.nodes)) + " < " + str(n_components))
        else:
            heat_signatures += chi.tolist()
            nodes_list += nodes
    print("finished graphwave_alg batches")
    pca = PCA(n_components = n_components)
    trans_data_all = pca.fit_transform(StandardScaler().fit_transform(np.array(heat_signatures)))
    km = KMeans(n_clusters = nb_clust).fit(trans_data_all)

    labels_pred = km.labels_

    out = pd.DataFrame(labels_pred.astype(int), index=nodes_list)
    out[0] += vocab_size
    structure_labels = out[0].to_dict()
    nx.set_node_attributes(G=networkXGraph, values=structure_labels, name='structure')

    return networkXGraph, pca, km

def walk_as_string(networkXGraph, componentLabels, params={'num_walks': 20, 'walk_length': 30}):
    """
    Generate random walks over a graph.

    Create a language to define nodes/edges and their characteristics, and generate random walks in this language.

    Example:

    Attributes:
        networkXGraph (networkx.classes.graph.Graph): A graph containing nodes, edges, and their attributes. Nodes must include an attribute called 'component'.
        componentLabels (dict): A dictionary mapping graph components to their y-values.

    Returns: a DataFrame containing each walk, and the associated graph label.
    """

    # TODO: Allow custom mapping for features e.g. pass a dict for node labels to convert them to chemical names

    # graphComponentLabels = nx.get_node_attributes(G=networkXGraph, name='component')

    num_walks = params['num_walks']
    walk_length = params['walk_length']

    nodeFeatures = list(
        set([z for x in list(networkXGraph.nodes(data=True)) for z in x[1].keys()]))

    nodeFeatures.remove('component')

    edgeFeatures = list(
        set([z for x in list(networkXGraph.edges(data=True)) for z in x[2].keys()]))

    # Remove 'cheating' features (e.g. component)
    if "component" in nodeFeatures:
        nodeFeatures.remove("component")

    # Make sure edges have weights
    if "weight" in edgeFeatures:
        edgeFeatures.remove("weight")
    else:
        nx.set_edge_attributes(
            G=networkXGraph, values=1, name='weight')

    n2vG = node2vec.Graph(nx_G=networkXGraph, is_directed=False, p=1, q=.7)
    n2vG.preprocess_transition_probs()
    walks = n2vG.simulate_walks(num_walks, walk_length)

    node_words = {}
    edge_words = {}

    def express_node(node_idx):
        if node_idx not in node_words:
            node = networkXGraph.nodes[node_idx]
            # node_words[node_idx] = " ".join([str(attribute)+"_"+str(node[attribute])
            #                 for attribute in nodeFeatures if attribute in node])
            node_words[node_idx] = " ".join([str(node[attribute])
                                             for attribute in nodeFeatures if attribute in node])

        return node_words[node_idx]

    def express_edge(src_node, dst_node):
        if (src_node, dst_node) not in edge_words:
            edge = networkXGraph.edges[src_node, dst_node]
            # edge_words[src_node, dst_node] = " ".join([str(attribute)+"_"+str(edge[attribute])
            #                 for attribute in edgeFeatures if attribute in edge])
            edge_words[src_node, dst_node] = " ".join(
                [str(edge[attribute]) for attribute in edgeFeatures if attribute in edge])
        return edge_words[src_node, dst_node]

    sorted_walks = pd.DataFrame(walks).sort_values(0).as_matrix()

    print(sorted_walks[0])
    print(sorted_walks[1])
    print(sorted_walks[2])

    walks = [list(a) for a in sorted_walks]

    walks_as_words = [express_node(walk[0]) + " " + " ".join([express_edge(walk[step], walk[step+1]) + " " +
                       express_node(walk[step+1]) for step in range(len(walk) - 1)]) for walk in walks]

    result = pd.DataFrame({"walk": walks_as_words, "start_node": np.array(walks)[:,0]})

    result['component'] = result['start_node'].map(nx.get_node_attributes(networkXGraph, name='component'))
    result['label'] = result['component'].map(componentLabels)

    result = result[['walk', 'label', 'start_node', 'component']]

    return result

max_stars_repo_name: dtylor/WalkAsString

max_issues_repo_name: dtylor/WalkAsString

max_forks_repo_name: dtylor/WalkAsString

reason: import numpy,import networkx


text:
"""Data utilities."""
import torch
from torch.autograd import Variable
import operator
import json
import pickle
import numpy as np
import datetime

def read_config(file_path):
	"""Read JSON config."""
	json_object = json.load(open(file_path, 'r'))
	return json_object

def read_vocab(file_path):
	vocab = [word for word in pickle.load(open(file_path, 'rb'))]
	word2id = {}
	id2word = {}
	for ind, word in enumerate(vocab):
		word2id[word] = ind
		id2word[ind] = word

	return word2id, id2word

def read_unigram(file_path):
	return pickle.load(open(file_path, 'rb'))

def read_bucket_data(word2id, id2word, src, config):
	"""Read data from files."""
	src_lines = []
	with open(src, 'r') as f:
		for ind, line in enumerate(f):
			src_lines.append(line.strip())

	src = {'data': src_lines, 'word2id': word2id, 'id2word': id2word}
	del src_lines

	return src

#from https://www.geeksforgeeks.org/python-find-maximum-length-sub-list-in-a-nested-list/
def findMaxLength(lst):
    #maxList = max(lst, key = len)
    #print("MAX LENGTH",lst)
    maxLength = max(map(len, lst))
    return maxLength

def get_story_batch_helper(lines, line_index, t_index, end_token):

	batch = []
	for i in range(line_index,len(lines)+1):
		if end_token in lines[i]: break
		batch.append(lines[i].strip().split("|||"))

	max_story_length = findMaxLength(batch)
	if t_index == max_story_length-1:
		return None, None, i+1, 0
		#if t_index is the last one, line_index = i+1

	src_lines = []
	trg_lines = []
	for story in batch:
		if len(story)>t_index:
			src_lines.append(story[t_index].split())
		else:
			src_lines.append(["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"])
		if len(story)> t_index+1:
			trg_lines.append(story[t_index+1].split())
		else:
			trg_lines.append(["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"])

	return src_lines, trg_lines, line_index, t_index+1

def get_story_batch_fit_helper(lines, line_index, t_index, end_token, max_size, min_size):

	batch = []
	curr_size = 0
	for i in range(line_index,len(lines)):
		if end_token in lines[i]:
			if curr_size <= max_size and curr_size >= min_size:
				break
			continue
		if curr_size >= max_size:
			i-=1
			break
		batch.append(lines[i].strip().split("|||"))
		curr_size+=1

	max_story_length = findMaxLength(batch)
	if t_index == max_story_length-1:
		return None, None, i+1, 0
		#if t_index is the last one, line_index = i+1

	src_lines = []
	trg_lines = []
	for story in batch:
		if len(story)>t_index:
			src_lines.append(story[t_index].split())
		else:
			src_lines.append(["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"])
		if len(story)> t_index+1:
			trg_lines.append(story[t_index+1].split())
		else:
			trg_lines.append(["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"])

	#print(i, t_index)
	#print("SRC",src_lines)
	#print("TRG", trg_lines)
	return src_lines, trg_lines, line_index, t_index+1

def get_story_batch_fit_helper_bigram(lines, line_index, t_index, end_token, max_size, min_size):

	batch = []
	curr_size = 0
	for i in range(line_index,len(lines)):
		if end_token in lines[i]:
			if curr_size <= max_size and curr_size >= min_size:
				break
			continue
		if curr_size > max_size:
			i-=1
			break
		batch.append(lines[i].strip().split("|||"))
		curr_size+=1

	max_story_length = findMaxLength(batch)
	if t_index == max_story_length-2:
		return None, None, i+1, 0
		#if t_index is the last one, line_index = i+1

	src_lines = []
	trg_lines = []
	for story in batch:
		src = []
		if len(story)>t_index:
			src = story[t_index].split()
		else:
			src = ["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"]
		if len(story)> t_index+1:
			src += story[t_index+1].split()
		else:
			src+= ["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"]
		src_lines.append(src)
		if len(story)> t_index+2:
			trg_lines.append(story[t_index+2].split())
		else:
			trg_lines.append(["<pad>", "<pad>", "<pad>", "<pad>", "<pad>"])

	#print(i, t_index)
	#print("SRC",src_lines)
	#print("TRG", trg_lines)
	return src_lines, trg_lines, line_index, t_index+1

def getAllBatches(lines, end_token, max_size, min_size, hidden = True):
	j = 0
	t = 0
	count = 0
	batches = {}
	restart = True
	while j < len(lines):
		if t == 0 or hidden == False: restart = True
		else: restart = False
		src_lines, trg_lines, j, t = get_story_batch_fit_helper(lines, j, t, end_token, max_size, min_size)
		#print(src_lines, trg_lines)
		if src_lines !=None:
			batches[count] = {'restart': restart, 'src':src_lines, 'trg':trg_lines}
			count+=1
	return batches

def getAllBigramBatches(lines, end_token, max_size, min_size, hidden = True):
	j = 0
	t = 0
	count = 0
	batches = {}
	restart = True
	while j < len(lines):
		if t == 0 or hidden == False: restart = True
		else: restart = False
		src_lines, trg_lines, j, t = get_story_batch_fit_helper_bigram(lines, j, t, end_token, max_size, min_size)
		#print(src_lines, trg_lines)
		if src_lines !=None:
			batches[count] = {'restart': restart, 'src':src_lines, 'trg':trg_lines}
			count+=1
	return batches

def fillStartEnd(lines, add_start, add_end):
	if add_start and add_end:
		lines = [
			['<s>'] + line + ['</s>']
			for line in lines
		]
	elif add_start and not add_end:
		lines = [
			['<s>'] + line
			for line in lines
		]
	elif not add_start and add_end:
		lines = [
			line + ['</s>']
			for line in lines
		]
	elif not add_start and not add_end:
		lines = [
			line
			for line in lines
		]
	#print("PADDED",lines)
	return lines

def process_batch(
	src_lines, trg_lines, word2ind, add_start=True, add_end=True
):
	"""Find story at start_index until end_token"""
	#print("FIRST STORY",lines[line_index])

	#src_lines, trg_lines, line_index, t_index = get_story_batch_helper(lines, line_index, t_index, end_token)

	#print("SOURCE",src_lines)
	#print("TARGET",trg_lines)

	#if src_lines == None: return None, None, None, None, line_index, t_index

	src_lines = fillStartEnd(src_lines, add_start, add_end)
	trg_lines = fillStartEnd(trg_lines, add_start, add_end)

	#pad events, convert to indices, and create output lines for trg

	max_src_event_length = findMaxLength(src_lines)
	max_trg_event_length = findMaxLength(trg_lines)
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	trg_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_trg_event_length - len(line))
		for line in trg_lines
	]

	trg_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_trg_event_length - len(line))
		for line in trg_lines
	]

	src_input_lines = Variable(torch.LongTensor(src_input_lines)).cuda()
	trg_input_lines = Variable(torch.LongTensor(trg_input_lines)).cuda()
	src_output_lines = Variable(torch.LongTensor(src_output_lines)).cuda()
	trg_output_lines = Variable(torch.LongTensor(trg_output_lines)).cuda()

	return src_input_lines, src_output_lines, trg_input_lines, trg_output_lines

def get_story_batch(
	lines, word2ind, line_index, t_index, end_token, add_start=True, add_end=True
):
	"""Find story at start_index until end_token"""
	#print("FIRST STORY",lines[line_index])

	src_lines, trg_lines, line_index, t_index = get_story_batch_helper(lines, line_index, t_index, end_token)

	#print("SOURCE",src_lines)
	#print("TARGET",trg_lines)

	if src_lines == None: return None, None, None, None, line_index, t_index

	src_lines = fillStartEnd(src_lines, add_start, add_end)
	trg_lines = fillStartEnd(trg_lines, add_start, add_end)

	#pad events, convert to indices, and create output lines for trg
	max_src_event_length = findMaxLength(src_lines)
	max_trg_event_length = findMaxLength(trg_lines)
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	trg_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_trg_event_length - len(line))
		for line in trg_lines
	]

	trg_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_trg_event_length - len(line))
		for line in trg_lines
	]

	src_input_lines = Variable(torch.LongTensor(src_input_lines)).cuda()
	trg_input_lines = Variable(torch.LongTensor(trg_input_lines)).cuda()
	src_output_lines = Variable(torch.LongTensor(src_output_lines)).cuda()
	trg_output_lines = Variable(torch.LongTensor(trg_output_lines)).cuda()

	return src_input_lines, src_output_lines, trg_input_lines, trg_output_lines, line_index, t_index

def read_data_pipeline(data, t, word2id, add_start=True, add_end=True):
	#data is a single story
	split = data.split("|||")
	if len(split) > t+1:
		#print(split[t], "|||", split[t+1])
		src_line = fillStartEnd([split[t].split()], add_start, add_end)
		trg_line = fillStartEnd([split[t+1].split()], add_start, add_end)

		src_line = [[word2id[w] if w in word2id else word2id['<unk>'] for w in src_line[0]]]
		trg_line = [[word2id[w] if w in word2id else word2id['<unk>'] for w in trg_line[0]]]
		src_line = torch.LongTensor(src_line).cuda()
		t = t+1
		"""
		elif len(data) == t+1:
			split = data.split("|||")
			src_line = [word2id[w] if w in word2id else word2id['<unk>'] for w in split[t].split()]
			trg_line = ["<pad>" for w in split[t].split()]
		"""
	elif len(split) == t+1:
		#print(split[t])
		src_line = fillStartEnd([split[t].split()], add_start, add_end)
		src_line = [[word2id[w] if w in word2id else word2id['<unk>'] for w in src_line[0]]]
		src_line = torch.LongTensor(src_line).cuda()
		trg_line = None
		t = 0
	else:
		#print("NONE")
		src_line = None
		trg_line = None
		t = 0
	return src_line, trg_line, t

def pipeline_srcOnly(src_lines, word2ind, add_start=True, add_end=True):
	#data is a single story
	src_lines = fillStartEnd(src_lines, add_start, add_end)
	max_src_event_length = findMaxLength(src_lines)
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_input_lines = Variable(torch.LongTensor(src_input_lines)).cuda()
	src_output_lines = Variable(torch.LongTensor(src_output_lines)).cuda()

	return src_input_lines, src_output_lines

def pipeline_srcOnly_duplicate(src_lines, word2ind, num_dups ,add_start=True, add_end=True):
	#data is a single EVENT
	if len(src_lines) > 1: raise ValueError("Data should only be a single input event")
	src_lines = fillStartEnd(src_lines, add_start, add_end)
	max_src_event_length = findMaxLength(src_lines)
	src_lines = src_lines[0]
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in src_lines[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(src_lines))
		for n in range(0, num_dups)
	]

	src_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in src_lines[1:]] +
		[word2ind['<pad>']] * (max_src_event_length - len(src_lines))
		for n in range(0, num_dups)
	]

	src_input_lines = Variable(torch.LongTensor(src_input_lines)).cuda()
	src_output_lines = Variable(torch.LongTensor(src_output_lines)).cuda()

	return src_input_lines, src_output_lines

def process_batch_pipeline(
	src_lines, trg_lines, word2ind, add_start=True, add_end=True
):

	src_lines = fillStartEnd(src_lines, add_start, add_end)
	trg_lines = fillStartEnd(trg_lines, add_start, add_end)

	#pad events, convert to indices, and create output lines for trg
	max_src_event_length = findMaxLength(src_lines)
	max_trg_event_length = findMaxLength(trg_lines)
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	trg_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_trg_event_length - len(line))
		for line in trg_lines
	]

	src_input_lines = torch.LongTensor(src_input_lines).cuda()
	trg_input_lines = torch.LongTensor(trg_input_lines).cuda()

	return src_input_lines, trg_input_lines

def init_state(input, num_layers, bidirectional, src_hidden_dim):
    """Get cell states and hidden states."""
    #print(input.size(0))
    #print(input.size(1))
    batch_size = input.size(0) #\
        #if self.encoder.batch_first else input.size(1)

    num_directions = 2 if bidirectional else 1
    src_hidden_dim = src_hidden_dim // 2 if bidirectional else src_hidden_dim

    h0_encoder = Variable(torch.zeros(
        num_layers * num_directions,
        batch_size,
        src_hidden_dim
    ), requires_grad=False)

    c0_encoder = Variable(torch.zeros(
        num_layers * num_directions,
        batch_size,
        src_hidden_dim
    ), requires_grad=False)

    return h0_encoder.cuda(), c0_encoder.cuda()

def read_policy_data_set(encoder_inputs, outputs, buckets, batch_size):
    data_set = [[] for _ in buckets]
    for i in range(0,batch_size):
        cur_vector = [encoder_inputs[i][3-j] for j in xrange(5)]
        target_vector = [outputs[j][i] for j in xrange(5)] #this was already 5, should it be 6?
        cur_vector = np.array(cur_vector)
        target_vector = np.array(target_vector)
        data_set[0].append([cur_vector,target_vector])
    return data_set

def process_batch_srcOnly(
	src_lines, word2ind, add_start=True, add_end=True
):
	print("SRC_LINES",src_lines)
	src_lines = fillStartEnd(src_lines, add_start, add_end)

	#pad events, convert to indices, and create output lines for trg

	max_src_event_length = findMaxLength(src_lines)
	print("SRC_LINES, filled",src_lines)
	src_input_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[:-1]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_output_lines = [
		[word2ind[w] if w in word2ind else word2ind['<unk>'] for w in line[1:]] +
		[word2ind['<pad>']] * (max_src_event_length - len(line))
		for line in src_lines
	]

	src_input_lines = Variable(torch.LongTensor(src_input_lines)).cuda()
	src_output_lines = Variable(torch.LongTensor(src_output_lines)).cuda()

	return src_input_lines, src_output_lines

max_stars_repo_name: lara-martin/ASTER-X

max_issues_repo_name: lara-martin/ASTER-X

max_forks_repo_name: lara-martin/ASTER-X


text:
# Copyright (c) Meta Platforms, Inc. and affiliates.
#
# This source code is licensed under the MIT license found in the
# LICENSE file in the root directory of this source tree.

import warnings

import beanmachine.ppl as bm
import numpy as np
import pytest
import torch
import torch.distributions as dist
from beanmachine.ppl.inference.proposer.hmc_utils import (
    DualAverageAdapter,
    MassMatrixAdapter,
    WelfordCovariance,
    WindowScheme,
    RealSpaceTransform,
)
from beanmachine.ppl.world import World

class SampleModel:
    @bm.random_variable
    def foo(self):
        return dist.Uniform(0.0, 1.0)

    @bm.random_variable
    def bar(self):
        return dist.Normal(self.foo(), 1.0)

def test_dual_average_adapter():
    adapter = DualAverageAdapter(0.1)
    epsilon1 = adapter.step(torch.tensor(1.0))
    epsilon2 = adapter.step(torch.tensor(0.0))
    assert epsilon2 < adapter.finalize() < epsilon1

def test_dual_average_with_different_delta():
    adapter1 = DualAverageAdapter(1.0, delta=0.8)
    adapter2 = DualAverageAdapter(1.0, delta=0.2)
    prob = torch.tensor(0.5)
    # prob > delta means we can increase the step size, wherease prob < delta means
    # we need to decrease the step size
    epsilon1 = adapter1.step(prob)
    epsilon2 = adapter2.step(prob)
    assert epsilon1 < epsilon2

def test_small_window_scheme():
    num_adaptive_samples = 10
    scheme = WindowScheme(num_adaptive_samples)
    for _ in range(num_adaptive_samples):
        # no window should be created if num_adaptive_samples is too small
        assert not scheme.is_in_window
        scheme.step()

def test_middle_window_scheme():
    num_adaptive_samples = 125
    scheme = WindowScheme(num_adaptive_samples)
    num_windows = 0
    for i in range(num_adaptive_samples):
        if scheme.is_in_window:
            # there should be a margin at the beginning and the end of a window
            assert i > 0
            if scheme.is_end_window:
                num_windows += 1
                assert i < num_adaptive_samples
        scheme.step()
    # there should only be a single window
    assert num_windows == 1

@pytest.mark.parametrize("num_adaptive_samples", [175, 300, 399, 543])
def test_large_window_scheme(num_adaptive_samples):
    scheme = WindowScheme(num_adaptive_samples)
    window_sizes = []
    for _ in range(num_adaptive_samples):
        if scheme.is_end_window:
            window_sizes.append(scheme._window_size)
        scheme.step()
    # size of windows should be monotonically increasing
    sorted_window_sizes = sorted(window_sizes)
    assert window_sizes == sorted_window_sizes
    for win1, win2 in zip(window_sizes[:-1], window_sizes[1:-1]):
        # except for last window, window size should keep doubling
        assert win2 == win1 * 2

def test_mass_matrix_adapter():
    model = SampleModel()
    world = World()
    world.call(model.bar())
    positions = RealSpaceTransform(world, world.latent_nodes)(dict(world))
    mass_matrix_adapter = MassMatrixAdapter()
    momentums = mass_matrix_adapter.initialize_momentums(positions)
    for node, z in positions.items():
        assert node in momentums
        assert isinstance(momentums[node], torch.Tensor)
        assert len(momentums[node]) == z.numel()
        # after seeing the nodes for the first time, the adapter should've initialized
        # the mass matrix and the distribution to generate momentum
        assert node in mass_matrix_adapter.mass_inv
        assert node in mass_matrix_adapter.momentum_dists
    mass_inv_old = mass_matrix_adapter.mass_inv.copy()
    mass_matrix_adapter.step(positions)

    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        mass_matrix_adapter.finalize()

    # mass matrix adapter has seen less than 2 samples, so mass_inv is not updated
    for node in mass_inv_old:
        assert torch.allclose(mass_inv_old[node], mass_matrix_adapter.mass_inv[node])

def test_diagonal_welford_covariance():
    samples = dist.MultivariateNormal(
        loc=torch.rand(5), scale_tril=torch.randn(5, 5).tril()
    ).sample((1000,))
    welford = WelfordCovariance(diagonal=True)
    for sample in samples:
        welford.step(sample)
    sample_var = torch.var(samples, dim=0)
    estimated_var = welford.finalize(regularize=False)
    assert torch.allclose(estimated_var, sample_var)
    regularized_var = welford.finalize(regularize=True)
    assert (torch.argsort(regularized_var) == torch.argsort(estimated_var)).all()

def test_dense_welford_covariance():
    samples = dist.MultivariateNormal(
        loc=torch.rand(5), scale_tril=torch.randn(5, 5).tril()
    ).sample((1000,))
    welford = WelfordCovariance(diagonal=False)
    for sample in samples:
        welford.step(sample)
    sample_cov = torch.from_numpy(np.cov(samples.T.numpy())).to(samples.dtype)
    estimated_cov = welford.finalize(regularize=False)
    assert torch.allclose(estimated_cov, sample_cov)
    regularized_cov = welford.finalize(regularize=True)
    assert (torch.argsort(regularized_cov) == torch.argsort(estimated_cov)).all()

def test_welford_exception():
    welford = WelfordCovariance()
    welford.step(torch.rand(5))
    with pytest.raises(RuntimeError):  # number of samples is too small
        welford.finalize()


text:
#Problem Number 119
#Solution Number 1
import numpy
class Solution:
    def getRow(self, rowIndex: int) -> List[int]:
        rows, cols=(rowIndex+1,rowIndex+1)
        arr= [[0 for i in range(cols)] for j in range(rows)]

        for i in range(rowIndex+1):
            for j in range(rowIndex+1):
                if j==0 or j==rowIndex or i==0:
                    if j==0 or j==rowIndex:
                        arr[i][j]=1
                    else:
                        arr[i][j]=0
                else:
                    arr[i][j]=arr[i-1][j] + arr[i-1][j-1]
        return arr[rowIndex]

#############################################################
############################################# Solution Number 2 ######################################################

class Solution:
    def getRow(self, rowIndex: int) -> List[int]:
        if rowIndex==0:
            return [1]
        prev_row= self.getRow(rowIndex-1)
        middle=[]

        for i in range(len(prev_row)-1):
            middle.append(prev_row[i] + prev_row[i+1])

        return [1] + middle + [1]


text:
import argparse

import keras.backend as K

from csbdeep.models import BaseConfig
from csbdeep.utils import _raise, axes_check_and_normalize, axes_dict, backend_channels_last

from six import string_types

import numpy as np
from logging.config import BaseConfigurator

# This class is a adapted version of csbdeep.models.config.py.
class N2VConfig(argparse.Namespace):
    """Default configuration for a N2V trainable CARE model.

    This class is meant to be used with :class:`N2V`.

    Parameters
    ----------
    X      : array(float)
             The training data 'X', with dimensions 'SZYXC' or 'SYXC'
    kwargs : dict
             Overwrite (or add) configuration attributes (see below).

    Example
    -------
    >>> n2v_config = N2VConfig(X, unet_n_depth=3)

    Attributes
    ----------
    unet_residual : bool
        Parameter `residual` of :func:`n2v_old.nets.common_unet`. Default: ``n_channel_in == n_channel_out``
    unet_n_depth : int
        Parameter `n_depth` of :func:`n2v_old.nets.common_unet`. Default: ``2``
    unet_kern_size : int
        Parameter `kern_size` of :func:`n2v_old.nets.common_unet`. Default: ``5 if n_dim==2 else 3``
    unet_n_first : int
        Parameter `n_first` of :func:`n2v_old.nets.common_unet`. Default: ``32``
    batch_norm : bool
        Activate batch norm
    unet_last_activation : str
        Parameter `last_activation` of :func:`n2v_old.nets.common_unet`. Default: ``linear``
    train_loss : str
        Name of training loss. Default: ``'laplace' if probabilistic else 'mae'``
    train_epochs : int
        Number of training epochs. Default: ``100``
    train_steps_per_epoch : int
        Number of parameter update steps per epoch. Default: ``400``
    train_learning_rate : float
        Learning rate for training. Default: ``0.0004``
    train_batch_size : int
        Batch size for training. Default: ``16``
    train_tensorboard : bool
        Enable TensorBoard for monitoring training progress. Default: ``True``
    train_checkpoint : str
        Name of checkpoint file for model weights (only best are saved); set to ``None`` to disable. Default: ``weights_best.h5``
    train_reduce_lr : dict
        Parameter :class:`dict` of ReduceLROnPlateau_ callback; set to ``None`` to disable. Default: ``{'factor': 0.5, 'patience': 10}``
    n2v_perc_pix : float
        Percentage of pixel to manipulate per patch. Default: ``1.5``
    n2v_patch_shape : tuple
        Random patches of this shape are extracted from the given training data. Default: ``(64, 64) if n_dim==2 else (64, 64, 64)``
    n2v_manipulator : str
        Noise2Void pixel value manipulator. Default: ``uniform_withCP``
    n2v_neighborhood_radius : int
        Neighborhood radius for n2v_old manipulator. Default: ``5``

        .. _ReduceLROnPlateau: https://keras.io/callbacks/#reducelronplateau
    """

    def __init__(self, X, **kwargs):

        # X is empty if config is None
        if (X.size != 0):

            assert len(X.shape) == 4 or len(X.shape) == 5, "Only 'SZYXC' or 'SYXC' as dimensions is supported."

            n_dim = len(X.shape) - 2
            n_channel_in = X.shape[-1]
            n_channel_out = n_channel_in

            means, stds = [], []
            for i in range(n_channel_in):
                means.append(np.mean(X[...,i]))
                stds.append(np.std(X[...,i]))

            if n_dim == 2:
                axes = 'SYXC'
            elif n_dim == 3:
                axes = 'SZYXC'

            # parse and check axes
            axes = axes_check_and_normalize(axes)
            ax = axes_dict(axes)
            ax = {a: (ax[a] is not None) for a in ax}

            (ax['X'] and ax['Y']) or _raise(ValueError('lateral axes X and Y must be present.'))
            not (ax['Z'] and ax['T']) or _raise(ValueError('using Z and T axes together not supported.'))

            axes.startswith('S') or (not ax['S']) or _raise(ValueError('sample axis S must be first.'))
            axes = axes.replace('S','') # remove sample axis if it exists

            if backend_channels_last():
                if ax['C']:
                    axes[-1] == 'C' or _raise(ValueError('channel axis must be last for backend (%s).' % K.backend()))
                else:
                    axes += 'C'
            else:
                if ax['C']:
                    axes[0] == 'C' or _raise(ValueError('channel axis must be first for backend (%s).' % K.backend()))
                else:
                    axes = 'C'+axes

            # normalization parameters
            self.means                 = [str(el) for el in means]
            self.stds                  = [str(el) for el in stds]
            # directly set by parameters
            self.n_dim                 = n_dim
            self.axes                  = axes
            self.n_channel_in          = int(n_channel_in)
            self.n_channel_out         = int(n_channel_out)

            # default config (can be overwritten by kwargs below)
            self.unet_residual         = False
            self.unet_n_depth          = 2
            self.unet_kern_size        = 5 if self.n_dim==2 else 3
            self.unet_n_first          = 32
            self.unet_last_activation  = 'linear'
            if backend_channels_last():
                self.unet_input_shape  = self.n_dim*(None,) + (self.n_channel_in,)
            else:
                self.unet_input_shape  = (self.n_channel_in,) + self.n_dim*(None,)

            self.train_loss            = 'mae'
            self.train_epochs          = 100
            self.train_steps_per_epoch = 400
            self.train_learning_rate   = 0.0004
            self.train_batch_size      = 16
            self.train_tensorboard     = True
            self.train_checkpoint      = 'weights_best.h5'
            self.train_reduce_lr       = {'factor': 0.5, 'patience': 10}
            self.batch_norm            = True
            self.n2v_perc_pix           = 1.5
            self.n2v_patch_shape       = (64, 64) if self.n_dim==2 else (64, 64, 64)
            self.n2v_manipulator       = 'uniform_withCP'
            self.n2v_neighborhood_radius = 5

            # disallow setting 'n_dim' manually
            try:
                del kwargs['n_dim']
                # warnings.warn("ignoring parameter 'n_dim'")
            except:
                pass

        self.probabilistic         = False

        for k in kwargs:
            setattr(self, k, kwargs[k])

    def is_valid(self, return_invalid=False):
        """Check if configuration is valid.

        Returns
        -------
        bool
            Flag that indicates whether the current configuration values are valid.
        """
        def _is_int(v,low=None,high=None):
            return (
                isinstance(v,int) and
                (True if low is None else low <= v) and
                (True if high is None else v <= high)
            )

        ok = {}
        ok['means'] = True
        for mean in self.means:
            ok['means'] &= np.isscalar(float(mean))
        ok['stds'] = True
        for std in self.stds:
            ok['stds'] &= np.isscalar(float(std))  and float(std) > 0.0
        ok['n_dim'] = self.n_dim in (2,3)
        try:
            axes_check_and_normalize(self.axes,self.n_dim+1,disallowed='S')
            ok['axes'] = True
        except:
            ok['axes'] = False
        ok['n_channel_in']  = _is_int(self.n_channel_in,1)
        ok['n_channel_out'] = _is_int(self.n_channel_out,1)

        ok['unet_residual'] = (
            isinstance(self.unet_residual,bool) and
            (not self.unet_residual or (self.n_channel_in==self.n_channel_out))
        )
        ok['unet_n_depth']         = _is_int(self.unet_n_depth,1)
        ok['unet_kern_size']       = _is_int(self.unet_kern_size,1)
        ok['unet_n_first']         = _is_int(self.unet_n_first,1)
        ok['unet_last_activation'] = self.unet_last_activation in ('linear','relu')
        ok['unet_input_shape'] = (
            isinstance(self.unet_input_shape,(list,tuple)) and
            len(self.unet_input_shape) == self.n_dim+1 and
            self.unet_input_shape[-1] == self.n_channel_in and
            all((d is None or (_is_int(d) and d%(2**self.unet_n_depth)==0) for d in self.unet_input_shape[:-1]))
        )
        ok['train_loss'] = (
            (self.train_loss in ('mse','mae'))
        )
        ok['train_epochs']          = _is_int(self.train_epochs,1)
        ok['train_steps_per_epoch'] = _is_int(self.train_steps_per_epoch,1)
        ok['train_learning_rate']   = np.isscalar(self.train_learning_rate) and self.train_learning_rate > 0
        ok['train_batch_size']      = _is_int(self.train_batch_size,1)
        ok['train_tensorboard']     = isinstance(self.train_tensorboard,bool)
        ok['train_checkpoint']      = self.train_checkpoint is None or isinstance(self.train_checkpoint,string_types)
        ok['train_reduce_lr']       = self.train_reduce_lr  is None or isinstance(self.train_reduce_lr,dict)
        ok['batch_norm']            = isinstance(self.batch_norm, bool)
        ok['n2v_perc_pix']           = self.n2v_perc_pix > 0 and self.n2v_perc_pix <= 100
        ok['n2v_patch_shape']       = (
            isinstance(self.n2v_patch_shape, (list,tuple)) and
            len(self.n2v_patch_shape) == self.n_dim and
            all(d > 0 for d in self.n2v_patch_shape)
        )
        ok['n2v_manipulator']       = self.n2v_manipulator in ['normal_withoutCP', 'uniform_withCP', 'normal_additive',
                                                               'normal_fitted', 'identity']

        if return_invalid:
            return all(ok.values()), tuple(k for (k,v) in ok.items() if not v)
        else:
            return all(ok.values())

    def update_parameters(self, allow_new=True, **kwargs):
        if not allow_new:
            attr_new = []
            for k in kwargs:
                try:
                    getattr(self, k)
                except AttributeError:
                    attr_new.append(k)
            if len(attr_new) > 0:
                raise AttributeError("Not allowed to add new parameters (%s)" % ', '.join(attr_new))
        for k in kwargs:
            setattr(self, k, kwargs[k])


text:
#!/usr/bin/env python

import sys
import rospy
import copy
import moveit_commander
import moveit_msgs.msg
import geometry_msgs.msg

from math import pi
from std_msgs.msg import String
from moveit_commander.conversions import pose_to_list

#initialise moveit commander and rospy nodes
moveit_commander.roscpp_initialize(sys.argv)
rospy.init_node('move_group_python_interface_tutorial',anonymous=True)

# instantiate Robotcommander object. This object is the outer inerface to the robot
robot = moveit_commander.RobotCommander()

# instantiate planningsceneinterfacae object. This object is an interface to the world surrounding the robot
scene=moveit_commander.PlanningSceneInterface()

#instantiate MoveGroupCommander for one of the planning groups of robot
group_name = "manipulator" #arm
group = moveit_commander.MoveGroupCommander(group_name)

#create a display trajectory publisher which is used later to publish trajectories for rvis to visulalize
display_trajectory_publisher=rospy.Publisher('/move_group/display_planned_path', moveit_msgs.msg.DisplayTrajectory, queue_size=20)

# # We can get the name of the reference frame for this robot:
# planning_frame = group.get_planning_frame()
# print "============ Reference frame: %s" % planning_frame

# # We can also print the name of the end-effector link for this group:
# eef_link = group.get_end_effector_link()
# print "============ End effector: %s" % eef_link

# # We can get a list of all the groups in the robot:
# group_names = robot.get_group_names()
# print "============ Robot Groups:", robot.get_group_names()

# # Sometimes for debugging it is useful to print the entire state of the
# # robot:
# print "============ Printing robot state"
print(robot.get_current_state())
# print ""

# move robot to non singular position by adjusting joint vlaues in group
joint_goal = group.get_current_joint_values()
print(joint_goal)
joint_goal[0] = 0
joint_goal[1] = 0
joint_goal[2] =  -1
joint_goal[3] =  1

# #use go command to execute the joint angles
group.go(joint_goal,wait=True)

# # call stop to avoid residual movements
group.stop()

import numpy as np
jac = group.get_jacobian_matrix(joint_goal)
value = np.linalg.det(jac)
print(value)


text:
#!/usr/bin/env python3
import pandas as pd
import numpy as np

# retrieved from https://github.com/BuysDB/IWSS
def sparseDistance( X, minPresence=1, minMeasurementsPerCell=1, weight=True ):
    """Calculate a distance matrix based on a boolean sparse cells/observations matrix
    Parameters
    ----------
    X : pandas dataframe
        sparse boolean cell by observation matrix, contains 1s 0s and NaNs for missing data
    minPresence : int
        mimumum amount of 1s or zeros present per feature, lower sums are pruned
    minMeasurementsPerCell : int
        minimum of total observations per cell
    Returns
    -------
    keptX : pandas dataframe
        Cell by observation matrix with only the cells and features which have been kept
    jointMatrix : numpy matrix
        Cell by cell distance matrix where the Similarity and Difference have been combined
    simMatrix : numpy matrix
        Similarity matrix
    distanceMatrix: numpy matrix
        Difference matrix
    normalisationFactor : numpy vector
        Weight per feature
    """

    prev= None
    selectedCells=None
    keptX = X.copy()
    while prev is None or prev[0]!=sum(selectedCells) or prev[1]!=sum(selectedColumns):
        if selectedCells is not None:
            prev = (sum(selectedCells), sum(selectedColumns))
        selectedCells = (keptX>=0).sum(axis=1) >= minMeasurementsPerCell # Select only cells/rows with at least one measurement
        # selectedColumns = ((keptX[selectedCells]>=0).sum(axis=0)>=minPresence)
        selectedColumns = ((keptX[selectedCells]==1).sum(axis=0)>=minPresence) & ((keptX[selectedCells]==0).sum(axis=0)>0)

        #print( f'We have {sum(selectedCells)} rows and {sum(selectedColumns)} X left' )
        keptX = keptX[selectedCells].loc[:,selectedColumns]

    if keptX.shape[0]<2:
        raise ValueError('Not enough data')

    pOnes = []
    pZeros = []
    for feature in keptX.columns:
        # Weights:
        column = keptX[feature]

        pOnes.append( -np.log2( ( np.sum(column==1)/len(column) )**2 ) ) #probability of two cells both having feature
        pZeros.append( -np.log2( ( np.sum(column==0)/len(column) )**2 ) )#probability of two cells not having feature,  (and we know it)

    pOnes = np.array(pOnes)
    pZeros = np.array(pZeros)
    #print( np.sum(np.isnan(pZeros)), np.sum(np.isnan(pOnes)),  np.min(pZeros), np.min(pOnes))
    #print(keptX.shape)
    rawMatrix = keptX.values

    #return rawMatrix
    iteration = 0

    # Similarity: how much do cells look alike?
    simMatrix =  np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )
    # What is the difference between the cells?
    distanceMatrix =  np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )

    #jointMatrix = np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )

    mv = int((len(simMatrix)*(len(simMatrix)-1))/2)

    for cai in range(rawMatrix.shape[0]):
        a = rawMatrix[cai,:]

        for cbi in range(rawMatrix.shape[0]):
            b = rawMatrix[cbi,:]

            # Un normalized distance
            pairwiseUnnormalizedDistance = np.logical_and( a==1, b==0 ) * \
            (pOnes + pZeros) + \
             np.logical_and( a==0, b==1 ) * (pZeros + pOnes) # For different batches the pOnes/pZeros is batch depended

            # Normalize the distance:
            normalisationFactor = np.sum( pOnes* (a==1)) + np.sum( pZeros*(a==0)) + \
                    np.sum( pOnes* (b==1)) + np.sum( pZeros* (b==0))

            pairwiseNormalizedDistance = np.sum(pairwiseUnnormalizedDistance) / (
                    normalisationFactor )
            distanceMatrix[cai, cbi] = pairwiseNormalizedDistance

            # Similarity calculation:
            sim = np.sum( (pOnes+pOnes) * np.logical_and( a==1, b==1 )) + \
                  np.sum( (pZeros+pZeros)*np.logical_and( a==0, b==0 ))

            normalisedSim = sim/normalisationFactor
            simMatrix[cai, cbi] = normalisedSim

            #joinedDistance =  (pairwiseNormalizedDistance*0.5+ (1.0-normalisedSim )*0.5)
            #jointMatrix[cai,cbi]= joinedDistance

            if cai==cbi:
                break

    for i in range(simMatrix.shape[0]):
        for j in range(simMatrix.shape[0]):
            simMatrix[j,i]=simMatrix[i,j]
            distanceMatrix[j,i]=distanceMatrix[i,j]
            if i==j:
                break

    simMatrix = pd.DataFrame(np.clip(simMatrix,0,1))
    simMatrix.index = keptX.index
    simMatrix.columns = keptX.index

    distanceMatrix = pd.DataFrame(np.clip(distanceMatrix,0,1))
    distanceMatrix.index = keptX.index
    distanceMatrix.columns = keptX.index

    jointMatrix = pd.DataFrame( (1-simMatrix)*0.5 + distanceMatrix*0.5 )
#    jointMatrix.index = keptX.index
    #jointMatrix.columns = keptX.index

    return keptX, jointMatrix, simMatrix, distanceMatrix, normalisationFactor

def sparseDistance_no_square( X, minPresence=1, minMeasurementsPerCell=1, weight=True ):
    """Calculate a distance matrix based on a boolean sparse cells/observations matrix
    Parameters
    ----------
    X : pandas dataframe
        sparse boolean cell by observation matrix, contains 1s 0s and NaNs for missing data
    minPresence : int
        mimumum amount of 1s or zeros present per feature, lower sums are pruned
    minMeasurementsPerCell : int
        minimum of total observations per cell
    Returns
    -------
    keptX : pandas dataframe
        Cell by observation matrix with only the cells and features which have been kept
    jointMatrix : numpy matrix
        Cell by cell distance matrix where the Similarity and Difference have been combined
    simMatrix : numpy matrix
        Similarity matrix
    distanceMatrix: numpy matrix
        Difference matrix
    normalisationFactor : numpy vector
        Weight per feature
    """

    prev= None
    selectedCells=None
    keptX = X.copy()
    while prev is None or prev[0]!=sum(selectedCells) or prev[1]!=sum(selectedColumns):
        if selectedCells is not None:
            prev = (sum(selectedCells), sum(selectedColumns))
        selectedCells = (keptX>=0).sum(axis=1) >= minMeasurementsPerCell # Select only cells/rows with at least one measurement
        selectedColumns = ((keptX[selectedCells]==1).sum(axis=0)>=minPresence) & ((keptX[selectedCells]==0).sum(axis=0)>0)
        print( f'We have {sum(selectedCells)} rows and {sum(selectedColumns)} X left' )
        keptX = keptX[selectedCells].loc[:,selectedColumns]

    if keptX.shape[0]<2:
        raise ValueError('Not enough data')

    pOnes = []
    pZeros = []
    for feature in keptX.columns:
        # Weights:
        column = keptX[feature]
        if weight:
            pOnes.append( -np.log2( ( np.sum(column==1)/len(column) ) ) ) #probability of two cells both having feature
            pZeros.append( -np.log2( ( np.sum(column==0)/len(column) ) ) )#probability of two cells not having feature,  (and we know it)
        else:
            pOnes.append( -np.log2( ( 0.5 ) ) ) #probability of two cells both having feature
            pZeros.append( -np.log2( (0.5) ))

    pOnes = np.array(pOnes)
    pZeros = np.array(pZeros)
    #print( np.sum(np.isnan(pZeros)), np.sum(np.isnan(pOnes)),  np.min(pZeros), np.min(pOnes))
    #print(keptX.shape)
    rawMatrix = keptX.values

    #return rawMatrix
    iteration = 0

    # Similarity: how much do cells look alike?
    simMatrix =  np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )
    # What is the difference between the cells?
    distanceMatrix =  np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )

    #jointMatrix = np.zeros( (rawMatrix.shape[0], rawMatrix.shape[0]) )

    mv = int((len(simMatrix)*(len(simMatrix)-1))/2)

    for cai in range(rawMatrix.shape[0]):
        a = rawMatrix[cai,:]

        for cbi in range(rawMatrix.shape[0]):
            b = rawMatrix[cbi,:]

            # Un normalized distance
            pairwiseUnnormalizedDistance = np.logical_and( a==1, b==0 ) * \
            (pOnes + pZeros) + \
             np.logical_and( a==0, b==1 ) * (pZeros + pOnes) # For different batches the pOnes/pZeros is batch depended

            # Normalize the distance:
            normalisationFactor = np.sum( pOnes* (a==1)) + np.sum( pZeros*(a==0)) + \
                    np.sum( pOnes* (b==1)) + np.sum( pZeros* (b==0))

            pairwiseNormalizedDistance = np.sum(pairwiseUnnormalizedDistance) / (
                    normalisationFactor )
            distanceMatrix[cai, cbi] = pairwiseNormalizedDistance

            # Similarity calculation:
            sim = np.sum( (pOnes+pOnes) * np.logical_and( a==1, b==1 )) + \
                  np.sum( (pZeros+pZeros)*np.logical_and( a==0, b==0 ))

            normalisedSim = sim/normalisationFactor
            simMatrix[cai, cbi] = normalisedSim

            #joinedDistance =  (pairwiseNormalizedDistance*0.5+ (1.0-normalisedSim )*0.5)
            #jointMatrix[cai,cbi]= joinedDistance

            if cai==cbi:
                break

    # Copy half of the matrix to the other side
    for i in range(simMatrix.shape[0]):
        for j in range(simMatrix.shape[0]):
            simMatrix[j,i]=simMatrix[i,j]
            distanceMatrix[j,i]=distanceMatrix[i,j]
            if i==j:
                break

    simMatrix = pd.DataFrame(np.clip(simMatrix,0,1))
    simMatrix.index = keptX.index
    simMatrix.columns = keptX.index

    distanceMatrix = pd.DataFrame(np.clip(distanceMatrix,0,1))
    distanceMatrix.index = keptX.index
    distanceMatrix.columns = keptX.index

    jointMatrix = pd.DataFrame( (1-simMatrix)*0.5 + distanceMatrix*0.5 )
#    jointMatrix.index = keptX.index
    #jointMatrix.columns = keptX.index

    return keptX, jointMatrix, simMatrix, distanceMatrix, normalisationFactor

max_stars_repo_name: BuysDB/siCloneFitIO

max_issues_repo_name: BuysDB/siCloneFitIO

max_forks_repo_name: BuysDB/siCloneFitIO