text:
import os
import json
import sys

from matplotlib import pyplot as plt
import numpy as np
import math
import argparse

dir_path = os.path.dirname(os.path.realpath(__file__))

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('-p', '--plots', nargs='+', help='''
                        One or more plots to be shown.
                        Available: latency, throughput, overhead, messages, wants, tcp.
                        ''')
    parser.add_argument('-o', '--outputs', nargs='+', help='''
                        One or more outputs to be shown.
                        Available: data
                        ''')
    parser.add_argument('-dir', '--dir', type=str, help='''
                        Result directory to process
                        ''')

    return parser.parse_args()

# def groupby(arr, key):
#     res = {}
#     for l in arr:
#         if !l[]
#         if len(res[""])
def process_result_line(l):
    l = json.loads(l)
    name = l["name"].split('/')
    value = (l["measures"])["value"]
    item = {}
    for attr in name:
        attr = attr.split(":")
        item[attr[0]] = attr[1]
    item["value"] = value
    return item

def aggregate_results(results_dir):
    res = []
    for subdir, _, files in os.walk(results_dir):
        for filename in files:
            filepath = subdir + os.sep + filename
            if filepath.split("/")[-1] == "results.out":
                # print (filepath)
                resultFile = open(filepath, 'r')
                for l in resultFile.readlines():
                    res.append(process_result_line(l))
    return res, len(os.listdir(results_dir))

def groupBy(agg, metric):
    res = {}
    for item in agg:
        if not item[metric] in res:
            res[item[metric]] = []
        res[item[metric]].append(item)
    return res

def plot_latency(byLatency, byBandwidth, byFileSize):

    plt.figure()

    p1, p2 = len(byLatency), len(byBandwidth)
    pindex = 1
    x = []
    y = {}
    tc = {}
    for l in byLatency:

        for b in byBandwidth:
            ax =plt.subplot(p1, p2, pindex)
            ax.set_title("latency: "+l + " bandwidth: " + b)
            ax.set_xlabel('File Size (MB)')
            ax.set_ylabel('time_to_fetch (ms)')

            for f in byFileSize:

                x.append(int(f)/1e6)

                y[f] = []
                tc[f] = []
                for i in byFileSize[f]:
                    if i["latencyMS"] == l and i["bandwidthMB"] == b and\
                            i["nodeType"]=="Leech":
                            if i["name"] == "time_to_fetch":
                                y[f].append(i["value"])
                            if i["name"] == "tcp_fetch":
                                tc[f].append(i["value"])

                avg = []
                for i in y:
                    scaled_y = [i/1e6 for i in y[i]]
                    ax.scatter([int(i)/1e6]*len(y[i]), scaled_y, marker="+")
                    avg.append(sum(scaled_y)/len(scaled_y))
                avg_tc = []
                for i in tc:
                    scaled_tc = [i/1e6 for i in tc[i]]
                    ax.scatter([int(i)/1e6]*len(tc[i]), scaled_tc, marker="*")
                    avg_tc.append(sum(scaled_tc)/len(scaled_tc))

            #print(y)
            ax.plot(x, avg, label="Bitswap fetch")
            ax.plot(x, avg_tc, label="TCP fetch")

            ax.legend()

            pindex+=1
            x = []
            y = {}
            tc = {}

def plot_tcp_latency(byLatency, byBandwidth, byFileSize):

    plt.figure()

    p1, p2 = len(byLatency), len(byBandwidth)
    pindex = 1
    x = []
    tc = {}
    for l in byLatency:

        for b in byBandwidth:
            ax =plt.subplot(p1, p2, pindex)
            ax.set_title("latency: "+l + " bandwidth: " + b)
            ax.set_xlabel('File Size (MB)')
            ax.set_ylabel('time_to_fetch (ms)')

            for f in byFileSize:

                x.append(int(f)/1e6)

                tc[f] = []
                for i in byFileSize[f]:
                    if i["latencyMS"] == l and i["bandwidthMB"] == b and\
                            i["nodeType"]=="Leech":
                            if i["name"] == "tcp_fetch":
                                tc[f].append(i["value"])

                avg_tc = []
                for i in tc:
                    scaled_tc = [i/1e6 for i in tc[i]]
                    ax.scatter([int(i)/1e6]*len(tc[i]), scaled_tc, marker="*")
                    avg_tc.append(sum(scaled_tc)/len(scaled_tc))

            # print(x, tc)
            ax.plot(x, avg_tc, label="TCP fetch")
            ax.legend()

            pindex+=1
            x = []
            tc = {}

def autolabel(ax, rects):
    """Attach a text label above each bar in *rects*, displaying its height."""
    for rect in rects:
        height = rect.get_height()
        ax.annotate('{}'.format(height),
                    xy=(rect.get_x() + rect.get_width() / 2, height),
                    xytext=(0, 3),  # 3 points vertical offset
                    textcoords="offset points",
                    ha='center', va='bottom')

def plot_messages(byFileSize, byTopology):

    # plt.figure()

    for t in byTopology:
        labels = []
        arr_blks_sent = []
        arr_blks_rcvd = []
        arr_dup_blks_rcvd = []
        arr_msgs_rcvd = []
        for f in byFileSize:
            labels.append(int(f)/1e6)
            x = np.arange(len(labels))  # the label locations
            blks_sent = blks_rcvd = dup_blks_rcvd = msgs_rcvd = 0
            blks_sent_n = blks_rcvd_n = dup_blks_rcvd_n = msgs_rcvd_n = 0
            width = 1/4

            for i in byFileSize[f]:
                if i["topology"]==t:
                    if i["name"] == "blks_sent":
                        blks_sent += i["value"]
                        blks_sent_n += 1
                    if i["name"] == "blks_rcvd":
                        blks_rcvd += i["value"]
                        blks_rcvd_n += 1
                    if i["name"] == "dup_blks_rcvd":
                        dup_blks_rcvd += i["value"]
                        dup_blks_rcvd_n += 1
                    if i["name"] == "msgs_rcvd":
                        msgs_rcvd += i["value"]
                        msgs_rcvd_n += 1

            # Computing averages
            # Remove the division if you want to see total values
            arr_blks_rcvd.append(round(blks_rcvd/blks_rcvd_n,1))
            arr_blks_sent.append(round(blks_sent/blks_sent_n,1))
            arr_dup_blks_rcvd.append(round(dup_blks_rcvd/dup_blks_rcvd_n,1))
            arr_msgs_rcvd.append(round(msgs_rcvd/msgs_rcvd_n,1))
            blks_sent = blks_rcvd = dup_blks_rcvd = msgs_rcvd = 0
            blks_sent_n = blks_rcvd_n = dup_blks_rcvd_n = msgs_rcvd_n = 0

        fig, ax = plt.subplots()
        bar1 = ax.bar(x-(3/2)*width, arr_msgs_rcvd, width, label="Msgs Received")
        bar2 = ax.bar(x-width/2, arr_blks_rcvd, width, label="Blocks Rcv")
        bar3 = ax.bar(x+width/2, arr_blks_sent, width, label="Blocks Sent")
        bar4 = ax.bar(x+(3/2)*width, arr_dup_blks_rcvd, width, label="Duplicates blocks")

        autolabel(ax, bar1)
        autolabel(ax, bar2)
        autolabel(ax, bar3)
        autolabel(ax, bar4)

        ax.set_ylabel('Number of messages')
        ax.set_ylabel('File Size (MB)')
        ax.set_title('Average number of messages exchanged ' + t)
        ax.set_xticks(x)
        ax.set_xticklabels(labels)
        ax.legend()
        fig.tight_layout()

def plot_want_messages(byFileSize, byTopology):

    # plt.figure()

    for t in byTopology:
        # t_aux = t.replace("(","").replace(")","").split(",")
        # instances = int(t_aux[0]) + int(t_aux[1])

        labels = []
        # arr_wants_avg = []
        arr_wants_max = []
        arr_wants_total = []
        arr_wants_avg_single = []
        arr_want_haves = []
        arr_want_blocks = []

        for f in byFileSize:
            labels.append(int(f)/1e6)
            x = np.arange(len(labels))  # the label locations
            wants = want_haves = want_blocks = 0
            wants_n = want_max = 0
            width = 1/4

            for i in byFileSize[f]:
                if i["topology"]==t:
                    if i["name"] == "wants_rcvd":
                        wants += i["value"]
                        wants_n += 1
                        if want_max < i["value"]:
                            want_max = i["value"]
                    if i["name"] == "want_blocks_rcvd":
                        want_blocks += i["value"]
                    if i["name"] == "want_haves_rcvd":
                        want_haves += i["value"]

            # Computing averages
            # Remove the division if you want to see total values
            # arr_wants_avg.append(round(wants/instances/1000,1))
            arr_want_haves.append(round(want_haves/wants_n,1))
            arr_want_blocks.append(round(want_blocks/wants_n,1))
            arr_wants_avg_single.append(round(wants/wants_n,1))
            arr_wants_max.append(want_max)
            arr_wants_total.append(wants/1000)

            wants  = 0
            wants_n = want_max = 0

        fig, ax = plt.subplots()
        bar1a = ax.bar(x-(3/2)*width, arr_want_haves, width, label="Average want-haves")
        bar1b = ax.bar(x-(3/2)*width, arr_want_blocks, width, label="Average want-blocks")
        bar2 = ax.bar(x-width/2, arr_wants_avg_single, width, label="Average wants per node in single file")
        bar3 = ax.bar(x+width/2, arr_wants_max, width, label="Max wants received by node in single file")
        bar4 = ax.bar(x+(3/2)*width, arr_wants_total, width, label="Total want messages exchanged in test (KMessages)")

        autolabel(ax, bar1a)
        autolabel(ax, bar1b)
        autolabel(ax, bar2)
        autolabel(ax, bar3)
        autolabel(ax, bar4)

        ax.set_ylabel('Number of messages')
        ax.set_title('Average number of WANTs exchanged ' + t)
        ax.set_xticks(x)
        ax.set_xticklabels(labels)
        ax.legend()
        fig.tight_layout()

def plot_bw_overhead(byFileSize, byTopology):

    # plt.figure()

    for t in byTopology:
        labels = []
        arr_control_rcvd = []
        arr_block_data_rcvd = []
        arr_dup_data_rcvd = []
        arr_overhead = []
        for f in byFileSize:
            #TODO: Considering a 5.5% overhead of TPC
            leechCount = t.replace("(", "").replace(")", "").split("-")[1]
            TCP_BASELINE = int(leechCount)*1.055*int(f)
            labels.append(int(f)/1e6)
            x = np.arange(len(labels))  # the label locations
            data_rcvd = block_data_rcvd = dup_data_rcvd = overhead = 0
            data_rcvd_n = block_data_rcvd_n = dup_data_rcvd_n = overhead_n = 0
            width = 1/4

            for i in byFileSize[f]:
                # We are only interested in leechers so we don't duplicate measurements.
                if i["nodeType"] == "Leech" and i["topology"]==t:
                    if i["name"] == "data_rcvd":
                        data_rcvd += i["value"]
                        data_rcvd_n += 1
                        overhead = (data_rcvd-TCP_BASELINE)*100/TCP_BASELINE
                        overhead_n += 1
                    if i["name"] == "block_data_rcvd":
                        block_data_rcvd += i["value"]
                        block_data_rcvd_n += 1
                    if i["name"] == "dup_data_rcvd":
                        dup_data_rcvd += i["value"]
                        dup_data_rcvd_n += 1

            control_rcvd = data_rcvd - block_data_rcvd
            # Computing averages
            # Remove the division if you want to see total values
            arr_control_rcvd.append(round(control_rcvd/data_rcvd_n/1e6,3))
            arr_block_data_rcvd.append(round(block_data_rcvd/block_data_rcvd_n/1e6,3))
            arr_dup_data_rcvd.append(round(dup_data_rcvd/dup_data_rcvd_n/1e6,3))
            arr_overhead.append(round(overhead/overhead_n,3))
            control_rcvd = data_rcvd = block_data_rcvd = dup_data_rcvd = overhead = 0
            data_rcvd_n = block_data_rcvd_n = dup_data_rcvd_n = overhead_n = 0

        fig, ax = plt.subplots()
        bar1 = ax.bar(x-(3/2)*width, arr_control_rcvd, width, label="Control data received (MB)")
        bar2 = ax.bar(x-width/2, arr_block_data_rcvd, width, label="Total data received from blocks (MB)")
        bar3 = ax.bar(x+width/2, arr_dup_data_rcvd, width, label="Total data received from duplicates (MB)")
        bar4 = ax.bar(x+(3/2)*width, arr_overhead, width, label="Bandwidth overhead (%)")

        autolabel(ax, bar1)
        autolabel(ax, bar2)
        autolabel(ax, bar3)
        autolabel(ax, bar4)

        ax.set_ylabel('Number of messages')
        # ax.set_ylabel('File Size (MB)')
        ax.set_title('Data received '+ t)
        ax.set_xticks(x)
        ax.set_xticklabels(labels)
        ax.legend()
        fig.tight_layout()

def output_avg_stream_data_bitswap(byFileSize, byTopology):

    for t in byTopology:
        labels = []
        arr_data_stream_sent = []

        for f in byFileSize:
            labels.append(int(f)/1e6)

            stream_data_sent = 0
            stream_data_sent_n = 0

            for i in byFileSize[f]:
                # We are only interested in leechers so we don't duplicate measurements.
                if i["nodeType"] == "Seed" and i["topology"]==t:
                    if i["name"] == "stream_data_sent":
                        stream_data_sent += i["value"]
                        stream_data_sent_n += 1

            # Computing averages
            # Remove the division if you want to see total values
            arr_data_stream_sent.append(round(stream_data_sent/stream_data_sent_n/1e6,3))
            stream_data_sent = 0
            stream_data_sent_n = 0

        print("=== Stream Data Sent ===")
        print("[*] Topology: ", t)
        i = 0
        for x in labels:
            print("Filesize: %s MB -- Avg. Stream Data Sent Seeders: %s MB" % (x, arr_data_stream_sent[i]) )
            i+=1

def output_latency(byFileSize, byTopology):

    for t in byTopology:
        labels = []
        arr_time_to_fetch = []

        for f in byFileSize:
            labels.append(int(f)/1e6)

            time_to_fetch = 0
            time_to_fetch_n = 0

            for i in byFileSize[f]:
                # We are only interested in leechers so we don't duplicate measurements.
                if i["nodeType"] == "Leech" and i["topology"]==t:
                    if i["name"] == "time_to_fetch":
                        time_to_fetch += i["value"]
                        time_to_fetch_n += 1

            # Computing averages
            # Remove the division if you want to see total values
            arr_time_to_fetch.append(round(time_to_fetch/1e6,3))

            time_to_fetch = 0
            time_to_fetch_n = 0

        print("=== Time to fetch ====")
        print("[*] Topology: ", t)
        i = 0
        for x in labels:
            print("[*]Filesize: %s MB" % x)
            print("Avg. Time to Fetch: %s ms" % (arr_time_to_fetch[i]) )
            i+=1

def output_avg_data(byFileSize, byTopology, nodeType):

    for t in byTopology:
        labels = []
        arr_total_data_in = []
        arr_total_data_out = []
        arr_total_rate_in = []
        arr_total_rate_out = []

        for f in byFileSize:
            labels.append(int(f)/1e6)

            total_data_in = total_data_out = total_rate_in = total_rate_out = 0
            total_data_in_n = total_data_out_n = total_rate_in_n = total_rate_out_n = 0

            for i in byFileSize[f]:
                # We are only interested in leechers so we don't duplicate measurements.
                if i["nodeType"] == nodeType and i["topology"]==t:
                    if i["name"] == "total_in":
                        total_data_in += i["value"]
                        total_data_in_n += 1
                    elif i["name"] == "total_out":
                        total_data_out += i["value"]
                        total_data_out_n += 1
                    elif i["name"] == "rate_in":
                        total_rate_in += i["value"]
                        total_rate_in_n += 1
                    elif i["name"] == "rate_out":
                        total_rate_out += i["value"]
                        total_rate_out_n += 1

            # Computing averages
            # Remove the division if you want to see total values
            arr_total_data_in.append(round(total_data_in/total_data_in_n/1e6,3))
            arr_total_data_out.append(round(total_data_out/total_data_out_n/1e6,3))
            arr_total_rate_in.append(round(total_rate_in/total_rate_in_n/1e6,3))
            arr_total_rate_out.append(round(total_rate_out/total_rate_out_n/1e6,3))

            total_data_in = total_data_out = total_rate_in = total_rate_out = 0
            total_data_in_n = total_data_out_n = total_rate_in_n = total_rate_out_n = 0

        print("=== Data Exchanges for %s ===" % nodeType)
        print("[*] Topology: ", t)
        i = 0
        for x in labels:
            print("[*]Filesize: %s MB" % x)
            print("Avg. Data In: %s MB" % (arr_total_data_in[i]) )
            print("Avg. Data Out: %s MB" % (arr_total_data_out[i]) )
            print("Avg. Rate In: %s MBps" % (arr_total_rate_in[i]) )
            print("Avg. Rate Out: %s MBps" % (arr_total_rate_out[i]) )
            i+=1

def plot_througput(byLatency, byBandwidth, byFileSize, byTopology, testcases):

    plt.figure()

    p1, p2 = 2, math.ceil(testcases/2)
    pindex = 1
    x = []
    y = {}
    toPlot = False

    for t in byTopology:
        for l in byLatency:

            for b in byBandwidth:
                ax =plt.subplot(p1, p2, pindex)
                ax.set_title("Average Throughput - latency: "+l + "ms bandwidth: " + b + "MB topology: " + t)
                ax.set_xlabel('File Size (MB)')
                ax.set_ylabel('throughput (Mbps)')

                for f in byFileSize:
                    time_to_fetch = block_data_rcvd = dup_data_rcvd = 0
                    time_to_fetch_n = block_data_rcvd_n = dup_data_rcvd_n = 0
                    x.append(int(f)/1e6)

                    y[f] = []
                    for i in byFileSize[f]:
                        if i["latencyMS"] == l and i["bandwidthMB"] == b and\
                            i["topology"] == t and i["nodeType"]=="Leech":
                                if i["name"] == "time_to_fetch":
                                    time_to_fetch += i["value"]/1e6  #Get in ms
                                    time_to_fetch_n +=1
                                if i["name"] == "block_data_rcvd":
                                    block_data_rcvd += i["value"]
                                    block_data_rcvd_n +=1
                                if i["name"] == "dup_data_rcvd":
                                    dup_data_rcvd += i["value"]
                                    dup_data_rcvd_n +=1

                    if time_to_fetch_n != 0:
                        avg_time_to_fetch = time_to_fetch / time_to_fetch_n / 1e3   # Use it in s
                        avg_data = ((block_data_rcvd/block_data_rcvd_n)-(dup_data_rcvd/dup_data_rcvd_n)) / 1e6   # IN MB
                        y[f].append(avg_data / avg_time_to_fetch)

                        time_to_fetch = block_data_rcvd = dup_data_rcvd = 0
                        time_to_fetch_n = block_data_rcvd_n = dup_data_rcvd_n = 0
                        toPlot = True

                avg = []
                if toPlot:
                    for i in y:
                        ax.scatter([int(i)/1e6]*len(y[i]), y[i], marker="+")
                        avg.append(sum(y[i])/len(y[i]))

                    ax.plot(x, avg)
                    pindex+=1
                    toPlot = False

                x = []
                y = {}

if __name__ == "__main__":
    args = parse_args()

    results_dir = dir_path + '/results'
    if args.dir:
        results_dir = args.dir

    agg, testcases = aggregate_results(results_dir)
    byLatency = groupBy(agg, "latencyMS")
    byNodeType = groupBy(agg, "nodeType")
    byFileSize = groupBy(agg, "fileSize")
    byBandwidth = groupBy(agg, "bandwidthMB")
    byTopology = groupBy(agg, "topology")
    byConnectionRate = groupBy(agg, "maxConnectionRate")

    if args.plots is None and args.outputs is None:
        print("[!!] No plots or outputs provided...")
        sys.exit()

    if args.plots is not None:
        if "latency" in args.plots:
            plot_latency(byLatency, byBandwidth, byFileSize)
        if "messages" in args.plots:
            plot_messages(byFileSize, byTopology)
        if "overhead" in args.plots:
            plot_bw_overhead(byFileSize, byTopology)
        if "throughput" in args.plots:
            plot_througput(byLatency, byBandwidth, byFileSize, byTopology, testcases)
        if "tcp" in args.plots:
            plot_tcp_latency(byLatency, byBandwidth, byFileSize)
        if "wants" in args.plots:
            plot_want_messages(byFileSize, byTopology)

    if args.outputs is not None:
            if "bitswap-data" in args.outputs:
                output_avg_stream_data_bitswap(byFileSize, byTopology)
            if "latency" in args.outputs:
                output_latency(byFileSize, byTopology)
            if "data" in args.outputs:
                output_avg_data(byFileSize, byTopology, "Seed")
                output_avg_data(byFileSize, byTopology, "Leech")

    plt.show()

max_stars_repo_name: obsidiansystems/beyond-bitswap

max_issues_repo_name: obsidiansystems/beyond-bitswap

max_forks_repo_name: obsidiansystems/beyond-bitswap


text:
import sympy
import torch
from torch.distributions import kl_divergence

from ..utils import get_dict_values
from .losses import Divergence

def KullbackLeibler(p, q, dim=None, analytical=True, sample_shape=torch.Size([1])):
    r"""
    Kullback-Leibler divergence (analytical or Monte Carlo Apploximation).

    .. math::

        D_{KL}[p||q] &= \mathbb{E}_{p(x)}\left[\log \frac{p(x)}{q(x)}\right] \qquad \text{(analytical)}\\
        &\approx \frac{1}{L}\sum_{l=1}^L \log\frac{p(x_l)}{q(x_l)},
         \quad \text{where} \quad  x_l \sim p(x) \quad \text{(MC approximation)}.

    Examples
    --------
    >>> import torch
    >>> from pixyz.distributions import Normal, Beta
    >>> p = Normal(loc=torch.tensor(0.), scale=torch.tensor(1.), var=["z"], features_shape=[64], name="p")
    >>> q = Normal(loc=torch.tensor(1.), scale=torch.tensor(1.), var=["z"], features_shape=[64], name="q")
    >>> loss_cls = KullbackLeibler(p,q,analytical=True)
    >>> print(loss_cls)
    D_{KL} \left[p(z)||q(z) \right]
    >>> loss_cls.eval()
    tensor([32.])
    >>> loss_cls = KullbackLeibler(p,q,analytical=False,sample_shape=[64])
    >>> print(loss_cls)
    \mathbb{E}_{p(z)} \left[\log p(z) - \log q(z) \right]
    >>> loss_cls.eval() # doctest: +SKIP
    tensor([31.4713])
    """
    if analytical:
        loss = AnalyticalKullbackLeibler(p, q, dim)
    else:
        loss = (p.log_prob() - q.log_prob()).expectation(p, sample_shape=sample_shape)
    return loss

class AnalyticalKullbackLeibler(Divergence):
    def __init__(self, p, q, dim=None):
        self.dim = dim
        super().__init__(p, q)

    @property
    def _symbol(self):
        return sympy.Symbol("D_{{KL}} \\left[{}||{} \\right]".format(self.p.prob_text, self.q.prob_text))

    def forward(self, x_dict, **kwargs):
        if (not hasattr(self.p, 'distribution_torch_class')) or (not hasattr(self.q, 'distribution_torch_class')):
            raise ValueError("Divergence between these two distributions cannot be evaluated, "
                             "got %s and %s." % (self.p.distribution_name, self.q.distribution_name))

        input_dict = get_dict_values(x_dict, self.p.input_var, True)
        self.p.set_dist(input_dict)

        input_dict = get_dict_values(x_dict, self.q.input_var, True)
        self.q.set_dist(input_dict)

        divergence = kl_divergence(self.p.dist, self.q.dist)

        if self.dim:
            divergence = torch.sum(divergence, dim=self.dim)
            return divergence, {}

        dim_list = list(torch.arange(divergence.dim()))
        divergence = torch.sum(divergence, dim=dim_list[1:])
        return divergence, {}

        """
        if (self._p1.distribution_name == "vonMisesFisher" and \
            self._p2.distribution_name == "HypersphericalUniform"):
            inputs = get_dict_values(x, self._p1.input_var, True)
            params1 = self._p1.get_params(inputs, **kwargs)

            hyu_dim = self._p2.dim
            return vmf_hyu_kl(params1["loc"], params1["scale"],
                              hyu_dim, self.device), x

        raise Exception("You cannot use these distributions, "
                        "got %s and %s." % (self._p1.distribution_name,
                                            self._p2.distribution_name))

        #inputs = get_dict_values(x, self._p2.input_var, True)
        #self._p2.set_dist(inputs)

        #divergence = kl_divergence(self._p1.dist, self._p2.dist)

        if self.dim:
            _kl = torch.sum(divergence, dim=self.dim)
            return divergence, x
        """

"""
def vmf_hyu_kl(vmf_loc, vmf_scale, hyu_dim, device):
    __m = vmf_loc.shape[-1]
    vmf_entropy = vmf_scale * ive(__m/2, vmf_scale) / ive((__m/2)-1, vmf_scale)
    vmf_log_norm = ((__m / 2 - 1) * torch.log(vmf_scale) - (__m / 2) * math.log(2 * math.pi) - (
        vmf_scale + torch.log(ive(__m / 2 - 1, vmf_scale))))
    vmf_log_norm = vmf_log_norm.view(*(vmf_log_norm.shape[:-1]))
    vmf_entropy = vmf_entropy.view(*(vmf_entropy.shape[:-1])) + vmf_log_norm

    hyu_entropy = math.log(2) + ((hyu_dim + 1) / 2) * math.log(math.pi) - torch.lgamma(
        torch.Tensor([(hyu_dim + 1) / 2])).to(device)
    return - vmf_entropy + hyu_entropy
"""

ext: py

max_stars_repo_name: masa-su/pixyz

max_issues_repo_name: masa-su/pixyz

max_forks_repo_name: masa-su/pixyz

reason: import sympy


text:
import cv2
import argparse
import numpy as np
import zipfile
import random
from tensorpack import RNGDataFlow, MapDataComponent, dftools

"""
DBDIR=/graphics/projects/scratch/datasets/mscoco/
python data_sampler.py --lmdb ${DBDIR}val2017.lmdb --input ${DBDIR}val2017.zip --create
python data_sampler.py --lmdb ${DBDIR}train2017.lmdb --input ${DBDIR}train2017.zip --create
"""

class ImageDataFromZIPFile(RNGDataFlow):
    """ Produce images read from a list of zip files. """
    def __init__(self, zipFile, shuffle=False, maxFiles=None):
        """
        Args:
            zipFile (list): list of zip file paths.
            channel (int): 1 or 3. Will convert grayscale to RGB images if channel==3.
            resize (tuple): int or (h, w) tuple. If given, resize the image.
        """
        assert len(zipFile), "No zip zipFile given to ImageFromFile!"
        self.shuffle = shuffle
        self.max = maxFiles
        self.archivefiles = []
        # print zipFile
        try:
            archive = zipfile.ZipFile(zipFile)
            imagesInArchive = archive.namelist()
            for img_name in imagesInArchive:
                if img_name.endswith('.jpg'):
                    self.archivefiles.append((archive, img_name))
        except zipfile.BadZipfile:
            print "Couldn't read " + zipFile
        if self.max is None:
            self.max = self.size()

    def size(self):
        return len(self.archivefiles)

    def get_data(self):
        if self.shuffle:
            self.rng.shuffle(self.archivefiles)
        # return <self.max> shuffled images  from zip
        self.archivefiles = random.sample(self.archivefiles, self.max)
        for archive in self.archivefiles:
            # print archive
            im_data = archive[0].read(archive[1])
            yield [im_data]

class ImageEncode(MapDataComponent):
    def __init__(self, ds, mode='.jpg', dtype=np.uint8, index=0):
        def func(img):
            img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
            return np.asarray(bytearray(cv2.imencode(mode, img)[1].tostring()), dtype=dtype)
        super(ImageEncode, self).__init__(ds, func, index=index)

class ImageDecode(MapDataComponent):
    def __init__(self, ds, mode='.jpg', dtype=np.uint8, index=0):
        def func(im_data):
            img = cv2.imdecode(np.asarray(bytearray(im_data), dtype=dtype), cv2.IMREAD_COLOR)
            return img
        super(ImageDecode, self).__init__(ds, func, index=index)

class RejectTooSmallImages(MapDataComponent):
    def __init__(self, ds, thresh=384, index=0):
        def func(img):
            h, w, _ = img.shape
            if (h < thresh) or (w < thresh):
                return None
            else:
                return img
        super(RejectTooSmallImages, self).__init__(ds, func, index=index)

class CenterSquareResize(MapDataComponent):
    def __init__(self, ds, index=0):
        """See section 5.3
        """
        def func(img):
            try:
                h, w, _ = img.shape
                if h > w:
                    off = (h - w) // 2
                    if off > 0:
                        img = img[off:-off, :, :]
                if w > h:
                    off = (w - h) // 2
                    if off > 0:
                        img = img[:, off:-off, :]

                img = cv2.resize(img, (256, 256))
                return img
            except:  # noqa
                return None
        super(CenterSquareResize, self).__init__(ds, func, index=index)

# Testcode for encode/decode.
if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--create', action='store_true', help='create lmdb')
    parser.add_argument('--debug', action='store_true', help='debug images')
    parser.add_argument('--input', type=str, help='path to lmdb', required=True)
    parser.add_argument('--lmdb', type=str, help='path to lmdb', required=True)
    parser.add_argument('--size', type=int, help='number of images', default=1000)
    parser.add_argument('--proc', type=int, help='number of processes', default=1)
    args = parser.parse_args()

    ds = ImageDataFromZIPFile(args.input)
    ds = ImageDecode(ds, index=0)
    ds = RejectTooSmallImages(ds, index=0)
    ds = CenterSquareResize(ds, index=0)
    if args.create:
        ds = ImageEncode(ds, index=0)
        dftools.dump_dataflow_to_lmdb(ds, args.lmdb)
    if args.debug:
        ds.reset_state()
        for i in ds.get_data():
            cv2.imshow('example', i[0])
            cv2.waitKey(0)

max_stars_repo_name: PatWie/tensorpack-recipes

max_issues_repo_name: AmeerAnsari/tensorflow-recipes

max_forks_repo_name: AmeerAnsari/tensorflow-recipes

reason: import numpy


text:
# Copyright 2017 The TensorFlow Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
# Copyright 2021 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from npu_bridge.npu_init import *
import time
import numpy as np

G_TIMERS = {}
G_TM_BEGINS = {}

def TM_CLEAR():
    G_TIMERS = {}

def TM_DECLEAR(name):
    if name not in G_TIMERS:
        G_TIMERS[name] = []

def TM_BEGIN(name):
    G_TM_BEGINS[name] = time.time()

def TM_END(name):
    t = time.time()
    TM_DECLEAR(name)
    G_TIMERS[name].append(t - G_TM_BEGINS.pop(name, t))

def TM_DISPLAY(remove_first_one=True):
    for key, values in G_TIMERS.items():
        startPos = 0
        if remove_first_one:
            startPos = 1
            if len(values) <= 1:
                continue
        print("{}:{} times, cost time avg:{}".format(key, len(values), np.average(values[startPos:])))

def TM_PICK(name):
    if name not in G_TIMERS or len(G_TIMERS[name])==0:
        return None
    else:
        return G_TIMERS[name][-1]

max_stars_repo_name: Ascend/modelzoo

max_issues_repo_name: Ascend/modelzoo

max_forks_repo_name: Ascend/modelzoo


text:
import numpy
import numpy.random

import csb.test as test

from csb.statistics import Cumulative
from csb.statistics import kurtosis, skewness, autocorrelation

@test.functional
class TestStatFunction(test.Case):

    def testCumulative(self):
        from scipy.stats import norm

        x = numpy.linspace(-5., 5., 200)
        samples = numpy.random.normal(size=100000)
        cumula = Cumulative(samples)
        c = cumula(x)

        cx = norm.cdf(x)
        for i in range(199):
            self.assertAlmostEqual(cx[i], c[i], delta=1e-2)

    def testKurtosis(self):
        samples = numpy.random.normal(size=100000)
        self.assertAlmostEqual(kurtosis(samples), 0., delta=1e-1)

        samples = numpy.random.uniform(-2., 2., size=100000)
        self.assertAlmostEqual(kurtosis(samples), -1.2, delta=1e-1)

    def testSkewness(self):
        samples = numpy.random.gamma(2., 0.5, size=100000)
        self.assertAlmostEqual(skewness(samples), 2. / numpy.sqrt(2.), delta=1e-1)

    def testAutorcorrelation(self):
        x = numpy.random.normal(size=1000) + numpy.sin(numpy.linspace(0., 2 * numpy.pi, 1000))
        n = 10
        ac = autocorrelation(x, n)
        self.assertAlmostEqual(ac[0], 1., delta=1e-1)

    def testEntropy(self):
        pass

    def testCircvar(self):
        pass

    def testCircmean(self):
        pass

max_stars_repo_name: ujjwalsh/CSB

max_issues_repo_name: ujjwalsh/CSB

max_forks_repo_name: ujjwalsh/CSB

reason: import numpy,from scipy


text:
# -*- coding: utf-8 -*-
import sys,os,json
import sys,os

father_path = os.path.join(os.getcwd())
print(father_path, "==father path==")

def find_bert(father_path):
	if father_path.split("/")[-1] == "BERT":
		return father_path

	output_path = ""
	for fi in os.listdir(father_path):
		if fi == "BERT":
			output_path = os.path.join(father_path, fi)
			break
		else:
			if os.path.isdir(os.path.join(father_path, fi)):
				find_bert(os.path.join(father_path, fi))
			else:
				continue
	return output_path

bert_path = find_bert(father_path)
t2t_bert_path = os.path.join(bert_path, "t2t_bert")
sys.path.extend([bert_path, t2t_bert_path])

import numpy as np
import tensorflow as tf
from bunch import Bunch
from data_generator import tokenization
import json

flags = tf.flags

FLAGS = flags.FLAGS

tf.logging.set_verbosity(tf.logging.INFO)

flags.DEFINE_string("buckets", "", "oss buckets")

## Required parameters
flags.DEFINE_string(
	"train_file", None,
	"Input TF example files (can be a glob or comma separated).")

flags.DEFINE_string(
	"output_file", None,
	"Input TF example files (can be a glob or comma separated).")

flags.DEFINE_string(
	"word_piece_model", None,
	"Input TF example files (can be a glob or comma separated).")

flags.DEFINE_string(
	"tokenizer_type", None,
	"Input TF example files (can be a glob or comma separated).")

def main(_):

	input_config = {
		"word_piece_model":FLAGS.word_piece_model
	}
	if FLAGS.tokenizer_type == 'spm':
		tokenizer = tokenization.SPM(input_config)
		tokenizer.load_model()
	elif FLAGS.tokenizer_type == 'jieba':
		tokenizer = jieba
	fwobj = open(FLAGS.output_file, "w")
	with open(FLAGS.train_file, "r") as frobj:
		for line in frobj:
			content = line.strip()
			if len(content) >= 1:
				token_lst = my_spm.tokenize(content)
				fwobj.write(" ".join(token_lst)+"\n")
			else:
				fwobj.write("\n")
	fwobj.close()

if __name__ == "__main__":
	tf.app.run()

max_stars_repo_name: yyht/bert

max_issues_repo_name: yyht/bert

max_forks_repo_name: yyht/bert


text:
from __future__ import print_function, division
import sys,os
# line 4 and line 5 below are for development purposes and can be removed
qspin_path = os.path.join(os.getcwd(),"../../")
sys.path.insert(0,qspin_path)
#
from quspin.operators import hamiltonian # Hamiltonians and operators
from quspin.basis import spinless_fermion_basis_1d, tensor_basis # Hilbert space fermion and tensor bases
import numpy as np # generic math functions
##### define model parameters #####
L=4 # system size
J=1.0 # hopping
U=np.sqrt(2.0) # interaction
mu=0.0 # chemical potential
##### construct Fermi-Hubbard Hamiltonian #####
# define boson basis with 3 states per site L bosons in the lattice
N_up = L//2 + L % 2 # number of fermions with spin up
N_down = L//2 # number of fermions with spin down
basis_up=spinless_fermion_basis_1d(L,Nf=N_up)
basis_down=spinless_fermion_basis_1d(L,Nf=N_down)
basis = tensor_basis(basis_up,basis_down) # spinful fermions
print(basis)
# define site-coupling lists
hop_right=[[-J,i,(i+1)%L] for i in range(L)] #PBC
hop_left= [[+J,i,(i+1)%L] for i in range(L)] #PBC
pot=[[-mu,i] for i in range(L)] # -\mu \sum_j n_{j \sigma}
interact=[[U,i,i] for i in range(L)] # U/2 \sum_j n_{j,up} n_{j,down}
# define static and dynamic lists
static=[
		['+-|',hop_left],  # up hops left
		['-+|',hop_right], # up hops right
		['|+-',hop_left],  # down hops left
		['|-+',hop_right], # down hops right
		['n|',pot],		   # up on-site potention
		['|n',pot],		   # down on-site potention
		['n|n',interact]   # up-down interaction
							]
dynamic=[]
# build Hamiltonian
no_checks = dict(check_pcon=False,check_symm=False,check_herm=False)
H=hamiltonian(static,dynamic,basis=basis,dtype=np.float64,**no_checks)

max_stars_repo_name: anton-buyskikh/QuSpin

max_issues_repo_name: cileeky/QuSpin

max_forks_repo_name: cileeky/QuSpin


text:
import cv2
import numpy as np
from PIL import ImageGrab
import time
from pyautogui import press
import keyboard

ver = 2

# Ekran:
# https://stackoverflow.com/questions/49479552/imagegrab-grabbbox-and-image-getpixel-used-together
# https://www.pyimagesearch.com/2021/01/20/opencv-getting-and-setting-pixels/

# Klawiatura:
# https://stackoverflow.com/posts/43567521/revisions

# Tylko Red dla sprawdzenia czy jest galaz
galaz = 56
galaz = 161 # bez filtra (Blue)
darkGalaz = 136 # bez filtra (Blue)

keyboardSleep = 0.0406 # Ver 3
keyboardSleep = 0.06
stepSleep = 0.107 # Ver 3
stepSleep = 0.12

nextStep = ['l','r']

# Poprawne koordynaty dla ekranu 1680x1050 na windowsie
# bbox=((1680/2)-100,(1050/2)-25,(1680/2)+100,(1050/2)-20)

# Zaczekaj na focus na okno
if ver>0:
    img = ImageGrab.grab(bbox=((1680/2)-100, 50,(1680/2)+100,(1050/2)+100))
    img_np = np.array(img)
    cv2.imshow("test", img_np)
    cv2.waitKey(1)
    print("Czekam na focus...")
    time.sleep(7)
    print("3...")
    time.sleep(1)
    print("2...")
    time.sleep(1)
    print("1...")
    time.sleep(1)

# Poczatkowo dwa razy gdziekolwiek
#press("left")
#press("left")
while ver == 1:
    # Pobierz obraz
    time.sleep(0.1)
    img = ImageGrab.grab(bbox=((1680/2)-100,(1050/2)-25,(1680/2)+100,(1050/2)-20))
    img_np = np.array(img)
    img_new = cv2.cvtColor(img_np, cv2.COLOR_BGR2RGB)
    #cv2.imshow("test", img_new)
    #cv2.waitKey(1)
    #print("Left: " + str(img_new[0,0+20]) + " Right: " + str(img_new[0,200-1-20]))
    # Sprawdz strone
    if(img_new[0,0+20][0] == galaz):
        print("Galaz po lewej stronie")
        nextStep[0] = nextStep[1]
        nextStep[1] = 'r'
    else:
        print("Galaz po prawej stronie")
        nextStep[0] = nextStep[1]
        nextStep[1] = 'l'
    # Wydaj polecenie
    if nextStep[0] == 'l':
        #press("left")
        #press("left")
        keyboard.press_and_release('left')
        time.sleep(0.05)
        keyboard.press_and_release('left')
        time.sleep(0.05)
    else:
        #press("right")
        #press("right")
        keyboard.press_and_release('right')
        time.sleep(0.05)
        keyboard.press_and_release('right')
        time.sleep(0.05)

while ver == 2:
    # Pobierz obraz
    time.sleep(stepSleep)
    img = ImageGrab.grab(bbox=((1680/2)-100, 50,(1680/2)+100,(1050/2)+100))
    img_new = np.array(img)
    #img_new = cv2.cvtColor(img_np, cv2.COLOR_BGR2RGB)
    #cv2.imshow("test", img_new)
    #cv2.waitKey(1)
    # Skanowanie
    galezie = [img_new[50-1, 155-1][0],
               img_new[150-1, 155-1][0],
               img_new[250-1, 155-1][0],
               img_new[350-1, 155-1][0],
               img_new[450-1, 155-1][0],
               img_new[550-1, 155-1][0]]
    for x in galezie[::-1]:
        if x == galaz:
            #print("Galaz po prawej")
            #press("left")
            #press("left")
            keyboard.press_and_release('left')
            time.sleep(keyboardSleep)
            keyboard.press_and_release('left')
            time.sleep(keyboardSleep)
        else:
            #print("Galaz po lewej")
            #press("right")
            #press("right")
            keyboard.press_and_release('right')
            time.sleep(keyboardSleep)
            keyboard.press_and_release('right')
            time.sleep(keyboardSleep)

while ver == 3:
    # ~~Srednia wartosc z boxow na drzewach~~
    # Znajdz galaz jeszcze podczas animacji opadania
    # Pobierz obraz
    time.sleep(stepSleep)
    img = ImageGrab.grab(bbox=((1680/2)-100, 50,(1680/2)+100,(1050/2)+100))
    img_new = np.array(img)
    #img_new = cv2.cvtColor(img_np, cv2.COLOR_BGR2RGB)

    cv2.imshow("test", img_new)
    #cv2.waitKey(10000)
    cv2.waitKey(1)
    # Skanowanie
    base = 47
    offset = 0
    for x in range(0,80):
        if img_new[base-x-1, 155-1][0] == galaz or img_new[base-x-1, 155-1][0] == darkGalaz or img_new[base-x-1, 40-1][0] == galaz or img_new[base-x-1, 40-1][0] == darkGalaz:
            offset = x
            break

    galezie = [img_new[base-offset-1, 155-1][0],
               img_new[100+base-offset-1, 155-1][0],
               img_new[200+base-offset-1, 155-1][0],
               img_new[300+base-offset-1, 155-1][0],
               img_new[400+base-offset-1, 155-1][0],
               img_new[500+base-offset-1, 155-1][0]]

    # DEBUG - Pokaz wartosc podczas wywalenia sie na przelocie
    #print("DEBUG: " + str(img_new[47-1, 40-1][0]) + " " + str(img_new[47-1, 155-1][0]))
    print("DEBUG: " + str(offset))

    for x in galezie[::-1]:
        if x == galaz or x == darkGalaz:
            #print("Galaz po prawej")
            #press("left")
            #press("left")
            keyboard.press_and_release('left')
            time.sleep(keyboardSleep)
            keyboard.press_and_release('left')
            time.sleep(keyboardSleep)
        else:
            #print("Galaz po lewej")
            #press("right")
            #press("right")
            keyboard.press_and_release('right')
            time.sleep(keyboardSleep)
            keyboard.press_and_release('right')
            time.sleep(keyboardSleep)

while ver == -1:
    # Centrowanie
    time.sleep(0.06)
    img = ImageGrab.grab(bbox=((1680/2)-100, 50,(1680/2)+100,(1050/2)+100))
    img_new = np.array(img)
    #img_new = cv2.cvtColor(img_np, cv2.COLOR_BGR2RGB)
    cv2.imshow("test", img_new)
    cv2.waitKey(1)
    # Skanowanie
    galezie = [img_new[50-1, 155-1][0],
               img_new[150-1, 155-1][0],
               img_new[250-1, 155-1][0],
               img_new[350-1, 155-1][0],
               img_new[450-1, 155-1][0],
               img_new[550-1, 155-1][0]]
    for x in galezie:
        print(x)
        if x == galaz:
            print("Galaz po prawej")
        else:
            print("Galaz po lewej")
    time.sleep(60)


text:
import numpy as np

#flatten an 3D np array
def flatten(X, pos = -1):
    flattened_X = np.empty((X.shape[0], X.shape[2]))  # sample x features array.
    if pos == -1:
        pos = X.shape[1]-1
    for i in range(X.shape[0]):
        flattened_X[i] = X[i, pos, :]
    return(flattened_X)

# Scale 3D array. X = 3D array, scalar = scale object from sklearn. Output = scaled 3D array.
def scale(X, scaler):
    for i in range(X.shape[0]):
        X[i, :, :] = scaler.transform(X[i, :, :])
    return X

# Symbolize a 3D array. X = 3D array, scalar = SAX symbolizer object. Output = symbolic 3D string array.
def symbolize(X, scaler):
    X_s = scaler.transform(X)
    #X_s = X.astype('U13')
    #for i in range(X.shape[0]):
        #X_s[i, :, :][:,0] = scaler.transform(np.array([X[i, :, :][:,0]]))
        #X_t.append(' '.join(X_s[i, :, :][:,0]))
    return X_s

# maps a string from vocab to [-1,1] in uniform steps
def trans(val, vocab) -> float:
    for i in range(len(vocab)):
        if val == vocab[i]:
            halfSize = (len(vocab)-1)/2
            return (i - halfSize) / halfSize

# Transform the data based on the given already fit SAX scaler and maps the symbols to [-1,1] via the trans method
def symbolizeTrans(X, scaler, bins = 5):
    vocab = scaler._check_params(bins)
    X_s = scaler.transform(X)
    #X_s = X.astype(str)
    for i in range(X.shape[0]):
        X = X.astype(float)

        #z1 = scaler.transform(np.array([X[i, :, :][:,0]]))

        for j in range(X.shape[1]):
            X[i][j] = trans(X_s[i][j], vocab)
    return X

#saves a object
def save_obj(obj, name ):
    with open(name + '.pkl', 'wb') as f:
        np.save(f, obj)

#loades a saved object
def load_obj(name ):
    with open(name + '.pkl', 'rb') as f:
        return np.load(f, allow_pickle=True)

#truncate a value to 3 zeros
def truncate(n):
    return int(n * 1000) / 1000

max_stars_repo_name: cslab-hub/LocalTSMHAInterpretability

max_issues_repo_name: cslab-hub/LocalTSMHAInterpretability

max_forks_repo_name: cslab-hub/LocalTSMHAInterpretability


text:
"""Code for FineNet in paper "Robust Minutiae Extractor: Integrating Deep Networks and Fingerprint Domain Knowledge" at ICB 2018
  https://arxiv.org/pdf/1712.09401.pdf

  If you use whole or partial function in this code, please cite paper:

  @inproceedings{Nguyen_MinutiaeNet,
    author    = {Dinh-Luan Nguyen and Kai Cao and Anil K. Jain},
    title     = {Robust Minutiae Extractor: Integrating Deep Networks and Fingerprint Domain Knowledge},
    booktitle = {The 11th International Conference on Biometrics, 2018},
    year      = {2018},
    }
"""

from __future__ import absolute_import
from __future__ import division

from keras.models import Model
from keras.layers import Activation, AveragePooling2D, BatchNormalization, Concatenate, Conv2D, Dense, \
    GlobalAveragePooling2D
from keras.layers import Input, Lambda, MaxPooling2D
from keras_applications.imagenet_utils import _obtain_input_shape
from tensorflow.keras import backend as K
import numpy as np
from CoarseNet.CoarseNet_utils import *

def orientation_loss(y_true, y_pred, lamb=1.):
    # clip
    y_pred = K.tf.clip_by_value(y_pred, K.epsilon(), 1 - K.epsilon())
    # get ROI
    label_seg = K.sum(y_true, axis=-1, keepdims=True)
    label_seg = K.tf.cast(K.tf.greater(label_seg, 0), K.tf.float32)
    # weighted cross entropy loss
    lamb_pos, lamb_neg = 1., 1.
    logloss = lamb_pos * y_true * K.log(y_pred) + lamb_neg * (1 - y_true) * K.log(1 - y_pred)
    logloss = logloss * label_seg  # apply ROI
    logloss = -K.sum(logloss) / (K.sum(label_seg) + K.epsilon())

    # coherence loss, nearby ori should be as near as possible
    # Oritentation coherence loss

    # 3x3 ones kernel
    mean_kernal = np.reshape(np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]], dtype=np.float32) / 8, [3, 3, 1, 1])

    sin2angle_ori, cos2angle_ori, modulus_ori = ori2angle(y_pred)
    sin2angle = K.conv2d(sin2angle_ori, mean_kernal, padding='same')
    cos2angle = K.conv2d(cos2angle_ori, mean_kernal, padding='same')
    modulus = K.conv2d(modulus_ori, mean_kernal, padding='same')

    coherence = K.sqrt(K.square(sin2angle) + K.square(cos2angle)) / (modulus + K.epsilon())
    coherenceloss = K.sum(label_seg) / (K.sum(coherence * label_seg) + K.epsilon()) - 1
    loss = logloss + lamb * coherenceloss
    return loss

def orientation_output_loss(y_true, y_pred):
    # clip
    y_pred = tf.clip_by_value(y_pred, K.epsilon(), 1 - K.epsilon())
    # get ROI
    label_seg = K.sum(y_true, axis=-1, keepdims=True)
    label_seg = tf.cast(tf.math.greater(label_seg, 0), tf.float32)
    # weighted cross entropy loss
    lamb_pos, lamb_neg = 1., 1.
    logloss = lamb_pos * y_true * K.log(y_pred) + lamb_neg * (1 - y_true) * K.log(1 - y_pred)
    logloss = logloss * label_seg  # apply ROI
    logloss = -K.sum(logloss) / (K.sum(label_seg) + K.epsilon())
    return logloss

def segmentation_loss(y_true, y_pred, lamb=1.):
    # clip
    y_pred = tf.clip_by_value(y_pred, K.epsilon(), 1 - K.epsilon())
    # weighted cross entropy loss
    total_elements = K.sum(tf.ones_like(y_true))
    label_pos = tf.cast(tf.math.greater(y_true, 0.0), tf.float32)
    lamb_pos = 0.5 * total_elements / K.sum(label_pos)
    lamb_neg = 1 / (2 - 1 / lamb_pos)
    logloss = lamb_pos * y_true * K.log(y_pred) + lamb_neg * (1 - y_true) * K.log(1 - y_pred)
    logloss = -K.mean(K.sum(logloss, axis=-1))
    # smooth loss
    smooth_kernal = np.reshape(np.array([[-1, -1, -1], [-1, 8, -1], [-1, -1, -1]], dtype=np.float32) / 8, [3, 3, 1, 1])
    smoothloss = K.mean(K.abs(K.conv2d(y_pred, smooth_kernal)))
    loss = logloss + lamb * smoothloss
    return loss

def minutiae_score_loss(y_true, y_pred):
    # clip
    y_pred = tf.clip_by_value(y_pred, K.epsilon(), 1 - K.epsilon())
    # get ROI
    label_seg = tf.cast(tf.math.not_equal(y_true, 0.0), tf.float32)
    y_true = tf.where(tf.math.less(y_true, 0.0), tf.zeros_like(y_true), y_true)  # set -1 -> 0
    # weighted cross entropy loss
    total_elements = K.sum(label_seg) + K.epsilon()
    lamb_pos, lamb_neg = 10., .5
    logloss = lamb_pos * y_true * K.log(y_pred) + lamb_neg * (1 - y_true) * K.log(1 - y_pred)
    # apply ROI
    logloss = logloss * label_seg
    logloss = -K.sum(logloss) / total_elements
    return logloss


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# Load the training results
train_accs = np.load('train_accs.npy')
test_accs = np.load('test_accs.npy')
train_losses = np.load('train_losses.npy')

plot_data = []
for acc_type, acc_data in {'Training Dataset': train_accs, 'Test Dataset': test_accs}.items():
    print(acc_data.shape)
    for epoch in range(acc_data.shape[0]):
        for subj in range(acc_data.shape[1]):
            for task in range(acc_data.shape[2]):
                plot_data.append({
                    'epoch': epoch,
                    'subj': subj,
                    'task': task,
                    'type': acc_type,
                    'value': acc_data[epoch, subj, task]
                })

plot_df = pd.DataFrame(plot_data)

# Plot
sns.set(style='whitegrid', palette='colorblind')

plt.subplot(211)
sns.lineplot(np.arange(len(train_losses)), train_losses)
plt.ylabel('Mean Squared Error')
plt.ylim(bottom=0)

plt.subplot(212)
sns.lineplot(x='epoch', y='value', hue='type', data=plot_df)

handles, labels = plt.gca().get_legend_handles_labels()
plt.gca().legend(
    handles=handles[1:], labels=labels[1:], bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
    ncol=2, mode="expand", borderaxespad=0., frameon=False)

plt.xlabel('Epochs')
plt.ylabel('Predictive Accuracy')

plt.suptitle('Autoencoder')
plt.tight_layout()
plt.subplots_adjust(top=0.88)
plt.savefig('train_eval_plot.pdf')
plt.show()

max_stars_repo_name: nriesterer/iccm-neural-bound

max_issues_repo_name: nriesterer/iccm-neural-bound

max_forks_repo_name: nriesterer/iccm-neural-bound


text:
import numpy as np
import pytest

from ...exceptions import NeuroLangNotImplementedError
from ..distributions import Distribution, TableDistribution
from ..exceptions import UncomparableDistributionsError

dog_cat_table = {"cat": 0.2, "dog": 0.8}

beach_jellyfish_table = {
    "chironex_fleckeri": 0.01,
    "portuguese_man_o_war": 0.05,
    "pelagia_noctiluca": 0.1,
    "lion_mane": 0.14,
    "chrysaora_hysoscella": 0.2,
    "moon_jelly": 0.2,
    "aequorea_victoria": 0.3,
}

jellyfish_dangerosity = {
    "chironex_fleckeri": 1.0,
    "portuguese_man_o_war": 0.8,
    "pelagia_noctiluca": 0.4,
    "lion_mane": 0.4,
    "chrysaora_hysoscella": 0.1,
    "moon_jelly": 0.05,
    "aequorea_victoria": 0.0,
}

def test_abstract_distribution_class():
    with pytest.raises(NeuroLangNotImplementedError):
        _ = Distribution().probability(42)
    with pytest.raises(NeuroLangNotImplementedError):
        _ = Distribution().support
    with pytest.raises(NeuroLangNotImplementedError):
        _ = Distribution().expectation(lambda x: x)

def test_table_distribution():
    distrib = TableDistribution(dog_cat_table)
    assert distrib.probability("cat") == 0.2
    assert distrib.probability("dog") == 0.8
    assert distrib.support == frozenset({"cat", "dog"})

def test_expectation_table_distribution():
    beach_jellyfish_dist = TableDistribution(beach_jellyfish_table)
    beach_dangerosity = beach_jellyfish_dist.expectation(
        lambda v: jellyfish_dangerosity[v]
    )
    expected_value = sum(
        jellyfish_dangerosity[k] * beach_jellyfish_table[k]
        for k in jellyfish_dangerosity.keys()
    )
    assert beach_dangerosity == expected_value

def test_conditional_table_distribution():
    observed_jellyfishes = {"chironex_fleckeri", "moon_jelly"}
    condition = lambda value: value in observed_jellyfishes
    beach_jellyfish_dist = TableDistribution(beach_jellyfish_table)
    beach_jellyfish_dist = beach_jellyfish_dist.conditioned_on(condition)
    beach_dangerosity = beach_jellyfish_dist.expectation(
        lambda v: jellyfish_dangerosity[v]
    )
    assert np.isclose(beach_dangerosity, 0.09523809523809525)

def test_compare_table_distribs():
    td1 = TableDistribution(beach_jellyfish_table)
    assert td1 == td1
    td2_table = beach_jellyfish_table.copy()
    td2_table["pelagia_noctiluca"] = 0.14
    td2_table["lion_mane"] = 0.1
    td2 = TableDistribution(td2_table)
    assert td1 != td2

def test_compare_table_distrib_with_non_table_distrib():
    td = TableDistribution(beach_jellyfish_table)
    with pytest.raises(UncomparableDistributionsError):
        if td == Distribution():
            pass

max_stars_repo_name: hndgzkn/NeuroLang

max_issues_repo_name: hndgzkn/NeuroLang

max_forks_repo_name: hndgzkn/NeuroLang


text:
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 21 15:42:18 2010

Author: josef-pktd
"""

import numpy as np
from scipy import linalg

def tiny2zero(x, eps = 1e-15):
    '''replace abs values smaller than eps by zero, makes copy
    '''
    mask = np.abs(x.copy()) <  eps
    x[mask] = 0
    return x

nobs = 5
autocov = 0.8**np.arange(nobs)
#from statsmodels.tsa import arima_process as ap
#autocov = ap.arma_acf([1, -0.8, 0.2], [1])[:10]
autocov = np.array([ 3.,  2.,  1.,  0.4,  0.12,  0.016, -0.0112,
        0.016    , -0.0112   , -0.01216  , -0.007488 , -0.0035584])/3.
autocov = autocov[:nobs]
sigma = linalg.toeplitz(autocov)
sigmainv = linalg.inv(sigma)

c = linalg.cholesky(sigma, lower=True)
ci = linalg.cholesky(sigmainv, lower=True)

print sigma
print tiny2zero(ci/ci.max())

"this is the text book transformation"
print 'coefficient for first observation', np.sqrt(1-autocov[1]**2)
ci2 = ci[::-1,::-1].T
print tiny2zero(ci2/ci2.max())

print np.dot(ci/ci.max(), np.ones(nobs))

print np.dot(ci2/ci2.max(), np.ones(nobs))

max_stars_repo_name: toobaz/statsmodels

max_issues_repo_name: AnaMP/statsmodels

max_forks_repo_name: AnaMP/statsmodels

reason: import numpy,from scipy,from statsmodels


text:
import numpy as np
import benchmark.museum_of_forks.p0.errors as er
import benchmark.museum_of_forks.p0.core.evolution as ev
from benchmark.museum_of_forks.p0.helpers.matrix import is_unitary

class FixedSystem(object):
    """
    Class that defines one specific instance of a Floquet problem, i.e.
    some system for which the dynamics have to be calculated. Provides the
    methods to compute U and dU.

    Provides the following methods:
    - u: computes u / returns already computed u
    - du: computes du / returns already computed du

    Has the following attributes:
    - hf: the Fourier transformed Hamiltonian (ndarray, square)
    - dhf: its derivative with respect to the controls (ndarray of square ndarrays)
    - params: an instance of FixedSystemParameters

    - omega: The frequency associated with the period T of the control pulse
    - t: Control duration

    The following parameters are inferred from hf and dhf during initialisation:
    - dim: the size of the Hilbert space
    - nc: the number of components in Hf
    - np: the number of control parameters
    """

    def __init__(self, hf, dhf, nz, omega, t, decimals=10, sparse=True):
        self.hf = hf
        self.dhf = dhf

        # Inferred parameters
        dim = hf.shape[1]
        nc = hf.shape[0]
        np = dhf.shape[0]

        self.params = FixedSystemParameters(dim, nz, nc, np, omega, t, decimals)

        self._u = None
        self._du = None
        self._vals, self._vecs, self._phi, self._psi = None, None, None, None

    def __eq__(self, other):
        assert isinstance(other, FixedSystem)
        hf_same = np.array_equal(self.hf, other.hf)
        dhf_same = np.array_equal(self.dhf, other.dhf)
        params_same = (self.params.nz, self.params.omega, self.params.t) \
            == (other.params.nz, other.params.omega, other.params.t)

        return hf_same and dhf_same and params_same

    @property
    def u(self):
        if self._u is not None:
            return self._u
        else:
            self._compute_u()
            return self._u

    @property
    def du(self):
        if self._du is not None:
            return self._du
        else:
            self._compute_du()
            return self._du

    def _compute_u(self):
        # Increase nz until U can be computed,
        # then set U and the intermediary results
        [nz_okay, results] = self._test_nz()
        while nz_okay is False:
            self.params.nz += 2
            [nz_okay, results] = self._test_nz()

        self._u, self._vals, self._vecs, self._phi, self._psi = results

    def _test_nz(self):
        # Try to compute U with the current nz,
        # if an error occurs or U is not unitary
        # return [False, []], else return [u, vecs, vals, phi, psi]
        try:
            results = ev.get_u_and_eigensystem(self.hf, self.params)
            if is_unitary(results[0]):
                return [True, results]
            else:
                return [False, []]
        except er.EigenvalueNumberError:
            return [False, []]

    def _compute_du(self):
        if self._u is None:
            self._compute_u()
        self._du = ev.get_du_from_eigensystem(self.dhf, self._psi,
                                              self._vals, self._vecs, self.params)

class DummyFixedSystem(FixedSystem):
    """
    A dummy FixedSystem that can be initialised with arbitrary dimensions
    without specifying hf and dhf (mainly for testing)
    """

    def __init__(self, **kwargs):
        self.params = FixedSystemParameters.optional(**kwargs)

class FixedSystemParameters(object):
    """
    Hold parameters for a FixedSystem

    - dim: the size of the Hilbert space
    - nz: Number of Fourier components taken into account for K, and of Brillouin
          zones in the resulting set of eigenvalues

    Derived from nz:
        - nz_min/max, the cutoff for the integers labelling the Fourier components
        - k_kim = dim*nz: the size of K

    - nc: number of components in Hf
    - np: number of control parameters
    - omega: The frequency associated with the period T of the control pulse
    - t: Control duration
    - decimals: The number of decimals used for rounding when finding unique eigenvalues
    - sparse: If True, a sparse eigensolver will be used
              -- unless working with very small systems with < 15 zones, this should be True
    """

    def __init__(self, dim, nz, nc, np, omega, t, decimals, sparse=True):
        self._dim = 0
        self._nz = 0
        self._nc = 0

        self.dim = dim
        self.nz = nz

        self.nc = nc
        self.np = np

        self.omega = omega
        self.t = t

        self.decimals = decimals
        self.sparse = sparse

    @classmethod
    def optional(self, dim=0, nz=1, nc=1, np=0, omega=1, t=1, decimals=10, sparse=True):
        """
        Class method to instantiate FixedSystemParameters without specifying
        the full set of parameters -- only needed for testing!
        """
        return FixedSystemParameters(dim, nz, nc, np, omega, t, decimals, sparse)

    @property
    def nz(self):
        return self._nz

    @nz.setter
    def nz(self, value):
        if value % 2 == 0:
            raise er.UsageError("Number of Fourier components in the \
                  extended space (nz) cannot be even.")

        self._nz = value
        self.k_dim = self.dim*value
        self.nz_max = (value-1)/2
        self.nz_min = -(value-1)/2

    @property
    def nc(self):
        return self._nc

    @nc.setter
    def nc(self, value):
        if value % 2 == 0:
            raise er.UsageError("Number of Fourier components of H \
                  cannot be even.")

        self._nc = value
        self.nc_max = (value-1)/2
        self.nc_min = -(value-1)/2

    @property
    def dim(self):
        return self._dim

    @dim.setter
    def dim(self, value):
        self._dim = value
        self.k_dim = self.nz*value

max_stars_repo_name: sirmarcel/floq

max_issues_repo_name: sirmarcel/floq

max_forks_repo_name: sirmarcel/floq


text:
import numpy as np
from open3d import geometry as o3dgm

from urdfpy import (rpy_to_matrix, matrix_to_rpy, matrix_to_xyz_rpy,
                    xyz_rpy_to_matrix)

def test_rpy_to_matrix():
    xr45 = o3dgm.get_rotation_matrix_from_axis_angle([np.pi / 2, 0, 0])
    yr45 = o3dgm.get_rotation_matrix_from_axis_angle([0, np.pi / 2, 0])
    zr45 = o3dgm.get_rotation_matrix_from_axis_angle([0, 0, np.pi / 2])
    c = zr45.dot(yr45.dot(xr45))

    assert np.allclose(rpy_to_matrix([np.pi / 2, 0, 0]), xr45)
    assert np.allclose(rpy_to_matrix([0, np.pi / 2, 0]), yr45)
    assert np.allclose(rpy_to_matrix([0, 0, np.pi / 2]), zr45)
    assert np.allclose(rpy_to_matrix(np.pi / 2 * np.ones(3)), c)

def test_matrix_to_rpy():
    xr45 = o3dgm.get_rotation_matrix_from_axis_angle([np.pi / 2, 0, 0])
    yr45 = o3dgm.get_rotation_matrix_from_axis_angle([0, np.pi / 2, 0])
    zr45 = o3dgm.get_rotation_matrix_from_axis_angle([0, 0, np.pi / 2])

    assert np.allclose([np.pi / 2, 0, 0], matrix_to_rpy(xr45))
    assert np.allclose([0, np.pi / 2, 0], matrix_to_rpy(yr45))
    assert np.allclose([0, 0, np.pi / 2], matrix_to_rpy(zr45))

def test_matrix_to_xyz_rpy():
    xr45 = np.eye(4)
    yr45 = np.eye(4)
    zr45 = np.eye(4)
    xr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([np.pi / 2, 0, 0])
    yr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([0, np.pi / 2, 0])
    zr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([0, 0, np.pi / 2])

    xr45[:3,3] = np.array([1,2,3])
    yr45[:3,3] = np.array([2,3,1])
    zr45[:3,3] = np.array([3,1,2])

    assert np.allclose([1, 2, 3, np.pi / 2, 0, 0], matrix_to_xyz_rpy(xr45))
    assert np.allclose([2, 3, 1, 0, np.pi / 2, 0], matrix_to_xyz_rpy(yr45))
    assert np.allclose([3, 1, 2, 0, 0, np.pi / 2], matrix_to_xyz_rpy(zr45))

def test_xyz_rpy_to_matrix():
    xr45 = np.eye(4)
    yr45 = np.eye(4)
    zr45 = np.eye(4)
    xr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([np.pi / 2, 0, 0])
    yr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([0, np.pi / 2, 0])
    zr45[:3, :3] = o3dgm.get_rotation_matrix_from_axis_angle([0, 0, np.pi / 2])

    xr45[:3,3] = np.array([1,2,3])
    yr45[:3,3] = np.array([2,3,1])
    zr45[:3,3] = np.array([3,1,2])
    assert np.allclose(xyz_rpy_to_matrix([1, 2, 3, np.pi / 2, 0, 0]), xr45)
    assert np.allclose(xyz_rpy_to_matrix([2, 3, 1, 0, np.pi / 2, 0]), yr45)
    assert np.allclose(xyz_rpy_to_matrix([3, 1, 2, 0, 0, np.pi / 2]), zr45)


text:
"""
Density of Sea Water using the Jackett et al. (2006) [1]_ function

Functions:

rho :: computes in-situ density from salinity, potential temperature and
    pressure

rho_s_t :: compute the partial derivatives of in-situ density with
    respect to salinity and potential temperature

rho_p :: compute the partial derivative of in-situ density with
    respect to pressure

Notes:
To make Boussinesq versions of these functions, see
`neutralocean.eos.tools.make_eos_bsq`.

To make vectorized versions of these functions, see
`neutralocean.eos.tools.vectorize_eos`.

.. [1] Jackett, D. R., McDougall, T. J., Feistel, R., Wright, D. G., &
   Griffies, S. M. (2006). Algorithms for Density, Potential Temperature,
   Conservative Temperature, and the Freezing Temperature of Seawater.
   Journal of Atmospheric and Oceanic Technology, 23(12), 1709–1728.
   https://doi.org/10.1175/JTECH1946.1

History
-------
Code adapted from MOM5.1 (Griffies et al) originally written in Fortran
2022 January 14 - David Hutchinson - Translated into python
2022 January 14 - Geoff Stanley    - code optimization, partial derivatives, check vals
"""

import numpy as np
import numba as nb

@nb.njit(nb.f8(nb.f8, nb.f8, nb.f8))
def rho(s, t, p):
    """
    Parameters
    ----------
    s : float
        Practical salinity [PSS-78]
    t : float
        Potential temperature [ITS-90]
    p : float
        Pressure [dbar]

    Returns
    -------
    rho : float
        In-situ density [kg m-3]
    """

    # List of coefficients for the rational function
    # fmt: off
    a0  =  9.9984085444849347e+02
    a1  =  7.3471625860981584e+00
    a2  = -5.3211231792841769e-02
    a3  =  3.6492439109814549e-04
    a4  =  2.5880571023991390e+00
    a5  = -6.7168282786692355e-03
    a6  =  1.9203202055760151e-03
    a7  =  1.1798263740430364e-02
    a8  =  9.8920219266399117e-08
    a9  =  4.6996642771754730e-06
    a10 = -2.5862187075154352e-08
    a11 = -3.2921414007960662e-12

    b0  =  1.0000000000000000e+00
    b1  =  7.2815210113327091e-03
    b2  = -4.4787265461983921e-05
    b3  =  3.3851002965802430e-07
    b4  =  1.3651202389758572e-10
    b5  =  1.7632126669040377e-03
    b6  = -8.8066583251206474e-06
    b7  = -1.8832689434804897e-10
    b8  =  5.7463776745432097e-06
    b9  =  1.4716275472242334e-09
    b10 =  6.7103246285651894e-06
    b11 = -2.4461698007024582e-17
    b12 = -9.1534417604289062e-18
    # fmt: on
    epsln = 1.0e-40

    # Precompute some commonly used terms
    t2 = t * t

    # Rational function for density
    num = (
        a0
        + t * (a1 + t * (a2 + a3 * t))
        + s * (a4 + a5 * t + a6 * s)
        + p * (a7 + a8 * t2 + a9 * s + p * (a10 + a11 * t2))
    )

    inv_den = 1.0 / (
        b0
        + t * (b1 + t * (b2 + t * (b3 + t * b4)))
        + s * (b5 + t * (b6 + b7 * t2) + np.sqrt(s) * (b8 + b9 * t2))
        + p * (b10 + p * t * (b11 * t2 + b12 * p))
        + epsln
    )

    return num * inv_den

@nb.njit
def rho_s_t(s, t, p):
    """
    Parameters
    ----------
    s : float
        Practical salinity [PSS-78]
    t : float
        Potential temperature [ITS-90]
    p : float
        Pressure [dbar]

    Returns
    -------
    rho_s : float
        Partial derivative of in-situ density with respect to salinity [kg m-3 psu-1]

    rho_t : float
        Partial derivative of in-situ density with respect to temperature [kg m-3 degc-1]
    """

    # List of coefficients for the rational function
    # fmt: off
    a0  =  9.9984085444849347e+02
    a1  =  7.3471625860981584e+00
    a2  = -5.3211231792841769e-02
    a3  =  3.6492439109814549e-04
    a4  =  2.5880571023991390e+00
    a5  = -6.7168282786692355e-03
    a6  =  1.9203202055760151e-03
    a7  =  1.1798263740430364e-02
    a8  =  9.8920219266399117e-08
    a9  =  4.6996642771754730e-06
    a10 = -2.5862187075154352e-08
    a11 = -3.2921414007960662e-12

    b0  =  1.0000000000000000e+00
    b1  =  7.2815210113327091e-03
    b2  = -4.4787265461983921e-05
    b3  =  3.3851002965802430e-07
    b4  =  1.3651202389758572e-10
    b5  =  1.7632126669040377e-03
    b6  = -8.8066583251206474e-06
    b7  = -1.8832689434804897e-10
    b8  =  5.7463776745432097e-06
    b9  =  1.4716275472242334e-09
    b10 =  6.7103246285651894e-06
    b11 = -2.4461698007024582e-17
    b12 = -9.1534417604289062e-18
    # fmt: on

    epsln = 1.0e-40

    # Precompute some commonly used terms
    t2 = t * t
    sp5 = np.sqrt(s)
    pt = p * t

    # Rational function for density
    num = (
        a0
        + t * (a1 + t * (a2 + a3 * t))
        + s * (a4 + a5 * t + a6 * s)
        + p * (a7 + a8 * t2 + a9 * s + p * (a10 + a11 * t2))
    )

    inv_den = 1.0 / (
        b0
        + t * (b1 + t * (b2 + t * (b3 + t * b4)))
        + s * (b5 + t * (b6 + b7 * t2) + sp5 * (b8 + b9 * t2))
        + p * (b10 + pt * (b11 * t2 + b12 * p))
        + epsln
    )

    # The density is
    #   rho = num / den
    # Taking the partial derivative w.r.t. s gives
    #   rho_s = (num_s - num * den_s / den ) / den
    # and similarly for rho_t

    num_s = a4 + a5 * t + 2.0 * a6 * s + p * a9
    num_t = (
        a1 + t * (2.0 * a2 + 3.0 * a3 * t) + a5 * s + 2.0 * a8 * pt + 2.0 * a11 * p * pt
    )

    den_s = b5 + t * (b6 + b7 * t2) + sp5 * (1.5 * b8 + 1.5 * b9 * t2)
    den_t = (
        b1
        + t * (2.0 * b2 + t * (3.0 * b3 + 4.0 * b4 * t))
        + s * (b6 + 3.0 * b7 * t2 + 2.0 * b9 * sp5 * t)
        + 3.0 * b11 * pt * pt
        + b12 * p ** 3
    )

    rho_s = (num_s - num * den_s * inv_den) * inv_den
    rho_t = (num_t - num * den_t * inv_den) * inv_den

    return rho_s, rho_t

@nb.njit(nb.f8(nb.f8, nb.f8, nb.f8))
def rho_p(s, t, p):
    """
    Parameters
    ----------
    s : float
        Practical salinity [PSS-78]
    t : float
        Potential temperature [ITS-90]
    p : float
        Pressure [dbar]

    Returns
    -------
    rho_p : float
        Partial derivative of in-situ density with respect to pressure [kg m-3 dbar-1]
    """

    # List of coefficients for the rational function
    # fmt: off
    a0  =  9.9984085444849347e+02
    a1  =  7.3471625860981584e+00
    a2  = -5.3211231792841769e-02
    a3  =  3.6492439109814549e-04
    a4  =  2.5880571023991390e+00
    a5  = -6.7168282786692355e-03
    a6  =  1.9203202055760151e-03
    a7  =  1.1798263740430364e-02
    a8  =  9.8920219266399117e-08
    a9  =  4.6996642771754730e-06
    a10 = -2.5862187075154352e-08
    a11 = -3.2921414007960662e-12

    b0  =  1.0000000000000000e+00
    b1  =  7.2815210113327091e-03
    b2  = -4.4787265461983921e-05
    b3  =  3.3851002965802430e-07
    b4  =  1.3651202389758572e-10
    b5  =  1.7632126669040377e-03
    b6  = -8.8066583251206474e-06
    b7  = -1.8832689434804897e-10
    b8  =  5.7463776745432097e-06
    b9  =  1.4716275472242334e-09
    b10 =  6.7103246285651894e-06
    b11 = -2.4461698007024582e-17
    b12 = -9.1534417604289062e-18
    # fmt: on
    epsln = 1.0e-40

    # Precompute some commonly used terms
    t2 = t * t

    # Rational function for density
    num = (
        a0
        + t * (a1 + t * (a2 + a3 * t))
        + s * (a4 + a5 * t + a6 * s)
        + p * (a7 + a8 * t2 + a9 * s + p * (a10 + a11 * t2))
    )

    inv_den = 1.0 / (
        b0
        + t * (b1 + t * (b2 + t * (b3 + t * b4)))
        + s * (b5 + t * (b6 + b7 * t2) + np.sqrt(s) * (b8 + b9 * t2))
        + p * (b10 + p * t * (b11 * t2 + b12 * p))
        + epsln
    )

    # The density is
    #   rho = num / den
    # Taking the partial derivative w.r.t. p gives
    #   rho_p = (num_p - num * den_p / den ) / den

    num_p = a7 + a8 * t2 + a9 * s + p * (2.0 * a10 + 2.0 * a11 * t2)

    den_p = b10 + p * t * (2.0 * b11 * t2 + 3.0 * b12 * p)

    return (num_p - num * den_p * inv_den) * inv_den

max_stars_repo_name: geoffstanley/neutralocean

max_issues_repo_name: geoffstanley/neutralocean

max_forks_repo_name: geoffstanley/neutralocean

reason: import numpy,import numba


text:
import numpy as np
import tensorflow as tf
from tfoptests.persistor import TensorFlowPersistor
from tfoptests.test_graph import TestGraph

'''
This is meant to be a simple placeholder to implement ONE op at a time in a graph and write out pbs to
test op mapping in nd4j samediff import
Donot check this in to the main branch as this is a one time thing...
Implement your graph with the op in question under def test_simple and run it
Will write pbs and relevant files to $DL4J_TEST_RESOURCES/src/main/resources/tf_graphs/examples/simple_run
Run mvn clean install in dl4j-test-resources and then run TFGraphTestList import test in nd4j with model_name set to simple_run
Make sure to run with both libnd4j executor and samediff executor
eg. below is with one hot op and no placeholders
'''

class SimpleRun(TestGraph):
    def __init__(self, *args, **kwargs):
        super(SimpleRun, self).__init__(*args, **kwargs)
        # self.input_1 = [[0, 2], [1, -1]]

    def list_inputs(self):
        return []

'''
    def get_placeholder_input(self, name):
        if name == "input_1":
            return self.input_1

    def _get_placeholder_shape(self, name):
        if name == "input_1":
            return [2, 2]
'''

def test_simple():
    simple_run = SimpleRun(seed=19)
    in_node_1 = tf.Variable([[0, 2], [1, -1]])
    out_node = tf.one_hot(in_node_1, 3, axis=1, off_value=-2.0, name="output")

    predictions = [out_node]
    # Run and persist
    tfp = TensorFlowPersistor(save_dir="simple_run")
    tfp.set_placeholders([]) \
        .set_output_tensors(predictions) \
        .set_test_data(simple_run.get_test_data()) \
        .build_save_frozen_graph()

if __name__ == '__main__':
    test_simple()

max_stars_repo_name: KonduitAI/ImportTests

max_issues_repo_name: KonduitAI/ImportTests

max_forks_repo_name: KonduitAI/ImportTests


text:
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as colors
import math
from scipy.stats import expon

from matplotlib import cm
from matplotlib import rc

__author__ = 'ernesto'

# if use latex or mathtext
rc('text', usetex=False)
rc('mathtext', fontset='cm')

# auxiliar function for plot ticks of equal length in x and y axis despite its scales.
def convert_display_to_data_coordinates(transData, length=10):
    # create a transform which will take from display to data coordinates
    inv = transData.inverted()
    # transform from display coordinates to data coordinates in x axis
    data_coords = inv.transform([(0, 0), (length, 0)])
    # get the length of the segment in data units
    yticks_len = data_coords[1, 0] - data_coords[0, 0]
    # transform from display coordinates to data coordinates in y axis
    data_coords = inv.transform([(0, 0), (0, length)])
    # get the length of the segment in data units
    xticks_len = data_coords[1, 1] - data_coords[0, 1]
    return xticks_len, yticks_len

#####################################
# PARAMETERS - This can be modified #
#####################################

# exponential PDF
theta = 2
lambd = 0.4
N = 20

#####################
# END OF PARAMETERS #
#####################

# se sortean N muestras de la PDF exponencial con parámetro lambda=theta
np.random.seed(seed=43)
xn = expon.rvs(loc=0, scale=1/theta, size=N)

# abscissa values
tmin = 0
tmax = 7

thetas = np.linspace(tmin, tmax, 300)
priori_pdf = lambd * np.exp(-lambd * thetas)
post_pdf = np.power(thetas, N) * np.exp(-thetas * np.sum(xn))
# post_pdf = np.power(thetas, N) * np.exp(-thetas * 0.37 * N)

# media muestral
sm = np.mean(xn)
post_pdf_sm = (1 / (sm ** N)) * math.exp(-N)
# post_pdf = np.power(1/sm, N) * np.exp(-N)

# axis parameters
dx = 0.4
xmin_ax = tmin - dx
xmax_ax = tmax + dx

dy = 0.1
ymax_ax = np.amax(post_pdf) * (1 + dy)
ymin_ax = -np.amax(post_pdf) * dy

# length of the ticks for all subplot (6 pixels)
display_length = 6  # in pixels
# x ticks labels margin
xtm = -0.15
ytm = -0.2
# font size
fontsize = 14
# colors from coolwarm
cNorm = colors.Normalize(vmin=0, vmax=1)
scalarMap = cm.ScalarMappable(norm=cNorm, cmap=cm.coolwarm)
col10 = scalarMap.to_rgba(0)
col20 = scalarMap.to_rgba(1)

fig = plt.figure(0, figsize=(5, 3), frameon=False)
ax = fig.add_subplot(111)

plt.xlim(xmin_ax, xmax_ax)
plt.ylim(ymin_ax, ymax_ax)

# horizontal and vertical ticks length
xtl, ytl = convert_display_to_data_coordinates(ax.transData, length=display_length)

# axis arrows
plt.annotate("", xytext=(xmin_ax, 0), xycoords='data', xy=(xmax_ax, 0), textcoords='data',
             arrowprops=dict(width=0.1, headwidth=6, headlength=8, facecolor='black', shrink=0.002))
plt.annotate("", xytext=(0, ymin_ax), xycoords='data', xy=(0, ymax_ax), textcoords='data',
             arrowprops=dict(width=0.1, headwidth=6, headlength=8, facecolor='black', shrink=0.002))

plt.plot(thetas, priori_pdf, color=col10, lw=2, label='$p(\\theta)$')
plt.plot(thetas, post_pdf, color=col20, lw=2, label='$p(\mathbf{x|}\\theta)$')

plt.plot([1/sm, 1/sm], [0, post_pdf_sm], 'k--', lw=1)

plt.plot([0, ytl], [lambd, lambd], 'k', lw=1)

# labels
plt.text(xmax_ax, xtm, '$\\theta$', fontsize=fontsize, ha='right', va='baseline')
plt.text(1/sm, xtm, '$1/\\bar{x}$', fontsize=fontsize, ha='center', va='baseline')
plt.text(ytm, lambd, '$\lambda$', fontsize=fontsize, ha='right', va='center')

leg = plt.legend(loc=1, frameon=False, fontsize=fontsize)

plt.axis('off')

# save as pdf image
plt.savefig('example_11_2.pdf', bbox_inches='tight')

plt.show()


text:
import numpy as np
import pandas as pd

from zyplots.data.DataLoader import DataLoader
from zyplots.core.GraphProcessor import GraphProcessor
from zyplots.util.DataUtil import DataUtil

dtype = {"split_id": int, "seq1_time": float, "c0_c1_copy_time": float, "seq2_time": float, "seq_fc_time": float}
usecols = ['split_id', 'seq1_time', 'c0_c1_copy_time', 'seq2_time', 'seq_fc_time']
index_col = False

data_dir_120 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split120_timeline_withheaders.csv"
records_120 = 120
filter_120 = [0, 1, 2, 5, 40, 41, 42, 57, 58, 59, 87, 88, 89, 117, 118, 119]

data_dir_60 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split60_timeline_withheader.csv"
records_60 = 60
filter_60 = [0, 1, 2, 5, 27, 28, 29, 57, 58, 59]

data_dir_40 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split40_timeline_withheaders.csv"
records_40 = 40
filter_40 = [0, 1, 2, 5, 15, 16, 17, 37, 38, 39]

data_dir_30 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split30_timeline_withheaders.csv"
records_30 = 30
filter_30 = [0, 1, 2, 5, 15, 16, 17, 27, 28, 29]

data_dir_24 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split24_timeline_withheaders.csv"
records_24 = 24
filter_24 = [0, 1, 2, 5, 12, 13, 14, 22, 23]

data_dir_20 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split20_timeline_withheaders.csv"
records_20 = 20
filter_20 = [0, 1, 2, 7, 8, 9, 17, 18, 19]

data_dir_15 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split15_timeline_withheaders.csv"
records_15 = 15
filter_15 = [0, 1, 2, 6, 7, 8, 12, 13, 14]

data_dir_12 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split12_timeline_withheaders.csv"
records_12 = 12
filter_12 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]

data_dir_10 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split10_timeline_withheaders.csv"
records_10 = 10
filter_10 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

data_dir_8 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split8_timeline_withheaders.csv"
records_8 = 8
filter_8 = [0, 1, 2, 3, 4, 5, 6, 7]

data_dir_6 = "/home/vibhatha/github/TimelineGraphTool/data/resnet/resnet_split6_timeline_withheaders.csv"
records_6 = 6
filter_6 = [0, 1, 2, 3, 4, 5]

data_dir_list = [data_dir_120, data_dir_60, data_dir_40, data_dir_30, data_dir_24, data_dir_20, data_dir_15,
                 data_dir_12, data_dir_10, data_dir_8, data_dir_6]
record_list = [records_120, records_60, records_40, records_30, records_24, records_20, records_15, records_12,
               records_10, records_8, records_6]
filter_list = [filter_120, filter_60, filter_40, filter_30, filter_24, filter_20, filter_15, filter_12, filter_10,
               filter_8, filter_6]

# id = 0
# data_dir_list = [data_dir_list[id]]
# record_list = [record_list[id]]
# filter_list = [filter_list[id]]

def atomic_analytics():
    for data_dir, record, filter in zip(data_dir_list, record_list, filter_list):
        data_loader = DataLoader(file_path=data_dir)

        pdf = data_loader.read_as_pandas()

        seq1 = pdf[usecols[1]]
        cp = pdf[usecols[2]]
        seq2 = pdf[usecols[3]]
        fc = pdf[usecols[4]]
        new_pdf = pdf.sum()
        print(seq1.sum() + cp.sum() + seq2.sum() + fc.sum())

def summary_analytics():
    dtype = {"split_id": int, "forward_time": float, "copy_time": float, "backward_time": float,
             "optimization_time": float}
    usecols = ['split_id', 'forward_time', 'copy_time', 'backward_time', 'optimization_time']
    data_file = "/home/vibhatha/github/TimelineGraphTool/data/summary/resnet/summary_v1.csv"
    pdf = pd.read_csv(data_file, usecols=usecols, dtype=dtype)
    #print(pdf.columns)
    #print(pdf.head())
    group_pdf = pdf.groupby(['split_id'])
    mean_pdf = group_pdf.mean()

    sum_pdf = mean_pdf.sum(axis=1)
    print(mean_pdf)
    print(sum_pdf)
    mean_pdf['total_time'] = sum_pdf

    print(mean_pdf)

summary_analytics()

max_stars_repo_name: vibhatha/zyplots

max_issues_repo_name: vibhatha/zyplots

max_forks_repo_name: vibhatha/zyplots


text:
"""
RF Prediction
"""

# import
import torch
import pickle
import numpy as np
import pandas as pd
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F

# NN building
class MLP(nn.Module):
    """
    MLP
    """
    def __init__(self, hidden_layers, input_size, output_size, seed=1):
        """
        `hidden_layers`: list, the number of neurons for every layer;
        `input_size`: number of states;
        `output_size`: number of actions;
        `seed`: random seed.
        """
        super().__init__()
        # parameters
        self.seed = torch.manual_seed(seed)
        # NN, adding layers dynamically.
        self.layers = nn.Sequential()
        # ---------------------- input -----------------------
        self.layers.add_module(
            'Linear_inp', nn.Linear(input_size, hidden_layers[0])
        )
        self.layers.add_module('Act_inp', nn.ReLU())
        # ---------------------- hidden ----------------------
        for i in range(1, len(hidden_layers)):
            self.layers.add_module(
                'Linear_{}'.format(i),
                nn.Linear(hidden_layers[i - 1], hidden_layers[i])
            )
            self.layers.add_module('Act_{}'.format(i), nn.ReLU())
        # ----------------------- output ---------------------
        self.layers.add_module(
            'Linear_out', nn.Linear(hidden_layers[-1], output_size)
        )
        self.layers.add_module('Act_out', nn.Softmax(dim=1))

    def forward(self, input_seq):
        """
        `input_seq`: states, torch.FloatTensor.
        """
        # return
        return self.layers(input_seq)

# NN training
def NN_training(hidden_size, X, y, train=True, train_iter=3000):
    """
    train the RF model with the given dataset.
    return a trained model.
    - train_data: with 'label' indicate the classes.
    """
    # model
    NN_model = MLP(
        hidden_layers=hidden_size,
        input_size=X.shape[1],
        output_size=2
    )
    # Train
    if train:
        # optimizer
        NN_optimizer = optim.Adam(
            NN_model.parameters(), lr=0.0001
        )
        loss_hist = {}
        train_iter = 12000
        batch_size = 1000
        print("Training...")
        for iter in range(train_iter):
            if iter % 1000 == 0:
                print("Iteration: {}".format(iter))
            # zero grad
            NN_optimizer.zero_grad()
            # batch training
            sample_ind = np.random.choice(
                range(X.shape[0]), batch_size, False
            )
            # X batch
            X_batch = torch.FloatTensor(X.iloc[sample_ind, :].values)
            # y_batch
            y_true_batch = torch.LongTensor([
                int(y[i]) for i in sample_ind
            ])
            # predict
            y_pred_batch = NN_model(X_batch)
            # loss
            loss = F.cross_entropy(input=y_pred_batch, target=y_true_batch)
            loss_hist[iter] = loss.detach().item()
            # backpropogate
            loss.backward()
            NN_optimizer.step()
        # save model
        torch.save(NN_model.state_dict(), "models/NN_{}.pt".format(
            hidden_size
        ))
    # not train
    else:
        NN_model.load_state_dict(torch.load(
            "models/NN_{}.pt".format(hidden_size)
        ))
        NN_model.eval()
    return NN_model

# temporal test
def temporal_test(horizon, train=True, hidden_size=[100]):
    """
    temporal testing
    """
    # =========================== Train ? =============================
    # load patient ID
    separate_id = pickle.load(open('data/ids/1h_ID.pickle', 'rb'))
    sepsis_id, nonsep_id = separate_id['sepsis'], separate_id['nonsep']
    # load data
    train_data = pd.read_csv(
        "data/feature_data/train_data.csv",
        index_col=False
    )
    # drop lab variables
    train_data = train_data.drop([
        'paO2_FiO2', 'platelets_x_1000', 'total_bilirubin',
        'urinary_creatinine', 'creatinine', 'HCO3', 'pH', 'paCO2',
        'direct_bilirubin', 'excess', 'ast', 'bun', 'calcium',
        'glucose', 'lactate', 'magnesium', 'phosphate', 'potassium',
        'hct', 'hgb', 'ptt', 'wbc', 'fibrinogen', 'troponin',
        'GCS_Score', 'ventilator'
    ], axis=1)
    # split data, ensure balanced train data
    sepsis_data = train_data.loc[
        train_data['patientunitstayid'].isin(sepsis_id)
    ]
    sepsis_data = sepsis_data.drop(['patientunitstayid', 'label'], axis=1)
    nonsep_data = train_data.loc[
        train_data['patientunitstayid'].isin(nonsep_id)
    ]
    nonsep_data = nonsep_data.drop(['patientunitstayid', 'label'], axis=1)
    # combined
    X_tr = sepsis_data.append(nonsep_data)
    y_tr = [1] * sepsis_data.shape[0] + [0] * nonsep_data.shape[0]
    # train model
    model = NN_training(hidden_size, X_tr, y_tr, train)
    # =========================== Test =============================
    print("Tesing...")
    # load patient id
    separate_id = pickle.load(open('data/ids/12h_ID.pickle', 'rb'))
    sepsis_id, nonsep_id = separate_id['sepsis'], separate_id['nonsep']
    all_id = sepsis_id + nonsep_id
    # prediction results
    col_names = ['patientunitstayid', 'label']\
        + list(range(-60 * horizon, -55, 5))
    cohort_pr = pd.DataFrame(columns=col_names)
    # loop
    for p_id in all_id:
        # load data
        patient_data = pd.read_csv(
            "data/feature_data/patient_feature/{}.csv".format(p_id),
            index_col=False
        )
        # 6 h, otherwise delete
        patient_data = patient_data.loc[
            (patient_data['offset'] >= -60 * horizon) &
            (patient_data['offset'] <= -60)
        ]
        # remove p_id, label, offset
        patient_data = patient_data.drop(
            ['patientunitstayid', 'label', 'offset'], axis=1
        )
        # remove lab
        patient_data = patient_data.drop([
            'paO2_FiO2', 'platelets_x_1000', 'total_bilirubin',
            'urinary_creatinine', 'creatinine', 'HCO3', 'pH', 'paCO2',
            'direct_bilirubin', 'excess', 'ast', 'bun', 'calcium',
            'glucose', 'lactate', 'magnesium', 'phosphate', 'potassium',
            'hct', 'hgb', 'ptt', 'wbc', 'fibrinogen', 'troponin',
            'GCS_Score', 'ventilator'
        ], axis=1)
        # make tensor
        patient_data = torch.FloatTensor(patient_data.values)
        # predict
        pred_pr = model(patient_data).detach().numpy()
        pred_pr = [
            pr[1] for pr in pred_pr
        ]
        # true label
        true_label = 1 if p_id in sepsis_id else 0
        # add pr to dataframe
        patient_record = [int(p_id), true_label] + pred_pr
        # dict
        patient_dict = {
            col_names[i]: patient_record[i] for i in range(len(col_names))
        }
        # add pr to dataframe
        cohort_pr = cohort_pr.append(patient_dict, ignore_index=True)
    # write pr table
    cohort_pr.to_csv("data/pr/NN-{}.csv".format(horizon), index=False)
    return

# Main
def main():
    """main"""
    np.random.seed(1)
    horizon = 6
    temporal_test(
        horizon, hidden_size=[100, 500, 100], train=True
    )
    return

if __name__ == "__main__":
    main()