text:
from tkinter import W
import numpy as np
import pandas as pd
from TACT.computation.adjustments import Adjustments

def perform_match_input(inputdata):
    """
    correct the TI inputs separately
    """
    results = pd.DataFrame(
        columns=[
            "sensor",
            "height",
            "adjustment",
            "m",
            "c",
            "rsquared",
            "difference",
            "mse",
            "rmse",
        ]
    )
    inputdata_train = inputdata[inputdata["split"] == True].copy()
    inputdata_test = inputdata[inputdata["split"] == False].copy()
    adj = Adjustments()

    if inputdata.empty or len(inputdata) < 2:
        results = adj.post_adjustment_stats([None], results, "Ref_TI", "adjTI_RSD_TI")
        if "Ane_TI_Ht1" in inputdata.columns and "RSD_TI_Ht1" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
            )
        if "Ane_TI_Ht2" in inputdata.columns and "RSD_TI_Ht2" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
            )
        if "Ane_TI_Ht3" in inputdata.columns and "RSD_TI_Ht3" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
            )
        if "Ane_TI_Ht4" in inputdata.columns and "RSD_TI_Ht4" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
            )
        m = np.NaN
        c = np.NaN
    else:
        full = pd.DataFrame()
        full["Ref_WS"] = inputdata_test["Ref_WS"]
        full["RSD_WS"] = inputdata_test["RSD_WS"]
        full = full.dropna()
        if len(full) < 10:
            results = adj.post_adjustment_stats(
                [None], results, "Ref_TI", "adjTI_RSD_TI"
            )
            m = np.NaN
            c = np.NaN
        else:
            WS_output = hist_match(inputdata_train, inputdata_test, "Ref_WS", "RSD_WS")
            SD_output = hist_match(inputdata_train, inputdata_test, "Ref_SD", "RSD_SD")
            inputdata_test["adjTI_RSD_TI"] = SD_output / WS_output
            results = adj.post_adjustment_stats(
                inputdata_test, results, "Ref_TI", "adjTI_RSD_TI"
            )
        if (
            "Ane_WS_Ht1" in inputdata.columns
            and "Ane_SD_Ht1" in inputdata.columns
            and "RSD_WS_Ht1" in inputdata.columns
            and "RSD_SD_Ht1" in inputdata.columns
        ):
            full = pd.DataFrame()
            full["Ref_WS"] = inputdata_test["Ane_WS_Ht1"]
            full["RSD_WS"] = inputdata_test["RSD_WS_Ht1"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
                )
                m = np.NaN
                c = np.NaN
            else:
                WS_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_WS_Ht1", "RSD_WS_Ht1"
                )
                SD_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_SD_Ht1", "RSD_SD_Ht1"
                )
                inputdata_test["adjTI_RSD_TI_Ht1"] = SD_output / WS_output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
                )
        if (
            "Ane_WS_Ht2" in inputdata.columns
            and "Ane_SD_Ht2" in inputdata.columns
            and "RSD_WS_Ht2" in inputdata.columns
            and "RSD_SD_Ht2" in inputdata.columns
        ):
            full = pd.DataFrame()
            full["Ref_WS"] = inputdata_test["Ane_WS_Ht2"]
            full["RSD_WS"] = inputdata_test["RSD_WS_Ht2"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
                )
                m = np.NaN
                c = np.NaN
            else:
                WS_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_WS_Ht2", "RSD_WS_Ht2"
                )
                SD_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_SD_Ht2", "RSD_SD_Ht2"
                )
                inputdata_test["adjTI_RSD_TI_Ht2"] = SD_output / WS_output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
                )
        if (
            "Ane_WS_Ht3" in inputdata.columns
            and "Ane_SD_Ht3" in inputdata.columns
            and "RSD_WS_Ht3" in inputdata.columns
            and "RSD_SD_Ht3" in inputdata.columns
        ):
            full = pd.DataFrame()
            full["Ref_WS"] = inputdata_test["Ane_WS_Ht3"]
            full["RSD_WS"] = inputdata_test["RSD_WS_Ht3"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
                )
                m = np.NaN
                c = np.NaN
            else:
                WS_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_WS_Ht3", "RSD_WS_Ht3"
                )
                SD_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_SD_Ht3", "RSD_SD_Ht3"
                )
                inputdata_test["adjTI_RSD_TI_Ht3"] = SD_output / WS_output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
                )
        if (
            "Ane_WS_Ht4" in inputdata.columns
            and "Ane_SD_Ht4" in inputdata.columns
            and "RSD_WS_Ht4" in inputdata.columns
            and "RSD_SD_Ht4" in inputdata.columns
        ):
            full = pd.DataFrame()
            full["Ref_WS"] = inputdata_test["Ane_WS_Ht4"]
            full["RSD_WS"] = inputdata_test["RSD_WS_Ht4"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
                )
                m = np.NaN
                c = np.NaN
            else:
                WS_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_WS_Ht4", "RSD_WS_Ht4"
                )
                SD_output = hist_match(
                    inputdata_train, inputdata_test, "Ane_SD_Ht4", "RSD_SD_Ht4"
                )
                inputdata_test["adjTI_RSD_TI_Ht4"] = SD_output / WS_output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
                )

    results["adjustment"] = ["SS-Match_input"] * len(results)
    results = results.drop(columns=["sensor", "height"])
    return inputdata_test, results

def perform_match(inputdata):
    # manipulate histogram to match template histogram (ref) - virtually a look-up table.
    results = pd.DataFrame(
        columns=[
            "sensor",
            "height",
            "adjustment",
            "m",
            "c",
            "rsquared",
            "difference",
            "mse",
            "rmse",
        ]
    )
    inputdata_train = inputdata[inputdata["split"] == True].copy()
    inputdata_test = inputdata[inputdata["split"] == False].copy()
    adj = Adjustments()

    if inputdata.empty or len(inputdata) < 2:
        results = adj.post_adjustment_stats([None], results, "Ref_TI", "adjTI_RSD_TI")
        if "Ane_TI_Ht1" in inputdata.columns and "RSD_TI_Ht1" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
            )
        if "Ane_TI_Ht2" in inputdata.columns and "RSD_TI_Ht2" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
            )
        if "Ane_TI_Ht3" in inputdata.columns and "RSD_TI_Ht3" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
            )
        if "Ane_TI_Ht4" in inputdata.columns and "RSD_TI_Ht4" in inputdata.columns:
            results = adj.post_adjustment_stats(
                [None], results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
            )
        m = np.NaN
        c = np.NaN
    else:
        full = pd.DataFrame()
        full["Ref_TI"] = inputdata_test["Ref_TI"]
        full["RSD_TI"] = inputdata_test["RSD_TI"]
        full = full.dropna()
        if len(full) < 10:
            results = adj.post_adjustment_stats(
                [None], results, "Ref_TI", "adjTI_RSD_TI"
            )
            m = np.NaN
            c = np.NaN
        else:
            output = hist_match(inputdata_train, inputdata_test, "Ref_TI", "RSD_TI")
            inputdata_test["adjTI_RSD_TI"] = output
            results = adj.post_adjustment_stats(
                inputdata_test, results, "Ref_TI", "adjTI_RSD_TI"
            )
        if "Ane_TI_Ht1" in inputdata.columns and "RSD_TI_Ht1" in inputdata.columns:
            full = pd.DataFrame()
            full["Ref_TI"] = inputdata_test["Ane_TI_Ht1"]
            full["RSD_TI"] = inputdata_test["RSD_TI_Ht1"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
                )
                m = np.NaN
                c = np.NaN
            else:
                output = hist_match(
                    inputdata_train, inputdata_test, "Ane_TI_Ht1", "RSD_TI_Ht1"
                )
                inputdata_test["adjTI_RSD_TI_Ht1"] = output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht1", "adjTI_RSD_TI_Ht1"
                )
        if "Ane_TI_Ht2" in inputdata.columns and "RSD_TI_Ht2" in inputdata.columns:
            full = pd.DataFrame()
            full["Ref_TI"] = inputdata_test["Ane_TI_Ht2"]
            full["RSD_TI"] = inputdata_test["RSD_TI_Ht2"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
                )
                m = np.NaN
                c = np.NaN
            else:
                output = hist_match(
                    inputdata_train, inputdata_test, "Ane_TI_Ht2", "RSD_TI_Ht2"
                )
                inputdata_test["adjTI_RSD_TI_Ht2"] = output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht2", "adjTI_RSD_TI_Ht2"
                )
        if "Ane_TI_Ht3" in inputdata.columns and "RSD_TI_Ht3" in inputdata.columns:
            full = pd.DataFrame()
            full["Ref_TI"] = inputdata_test["Ane_TI_Ht3"]
            full["RSD_TI"] = inputdata_test["RSD_TI_Ht3"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
                )
                m = np.NaN
                c = np.NaN
            else:
                output = hist_match(
                    inputdata_train, inputdata_test, "Ane_TI_Ht3", "RSD_TI_Ht3"
                )
                inputdata_test["adjTI_RSD_TI_Ht3"] = output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht3", "adjTI_RSD_TI_Ht3"
                )
        if "Ane_TI_Ht4" in inputdata.columns and "RSD_TI_Ht4" in inputdata.columns:
            full = pd.DataFrame()
            full["Ref_TI"] = inputdata_test["Ane_TI_Ht4"]
            full["RSD_TI"] = inputdata_test["RSD_TI_Ht4"]
            full = full.dropna()
            if len(full) < 10:
                results = adj.post_adjustment_stats(
                    [None], results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
                )
                m = np.NaN
                c = np.NaN
            else:
                output = hist_match(
                    inputdata_train, inputdata_test, "Ane_TI_Ht4", "RSD_TI_Ht4"
                )
                inputdata_test["adjTI_RSD_TI_Ht4"] = output
                results = adj.post_adjustment_stats(
                    inputdata_test, results, "Ane_TI_Ht4", "adjTI_RSD_TI_Ht4"
                )

    import matplotlib.pyplot as plt

    #    plt.plot(inputdata_test['Ref_TI'])
    #    plt.plot(inputdata_test['RSD_TI'])
    #    plt.plot(inputdata_test['adjTI_RSD_TI'])

    #    plt.scatter(inputdata_test['Ref_TI'], inputdata_test['RSD_TI'], label='RefvsRSD')
    #    plt.scatter(inputdata_test['Ref_TI'], inputdata_test['adjTI_RSD_TI'], label='RefvsCorrectedRSD')
    #    plt.scatter(inputdata_test['Ref_TI'], inputdata_test['Ane2_TI'], label='RefvsRedundant')
    #    plt.legend()
    #    plt.show()
    results["adjustment"] = ["SS-Match"] * len(results)
    results = results.drop(columns=["sensor", "height"])
    return inputdata_test, results

def hist_match(inputdata_train, inputdata_test, refCol, testCol):

    test1 = inputdata_test[refCol].copy
    source = inputdata_test[refCol].copy().dropna()
    template = inputdata_train[testCol].copy().dropna()

    # get the set of unique pixel values and their corresponding indices and
    # counts
    s_values, bin_idx, s_counts = np.unique(
        source, return_inverse=True, return_counts=True
    )
    t_values, t_counts = np.unique(template, return_counts=True)
    n_bins = 200
    # take the cumsum of the counts and normalize by the number of pixels to
    # get the empirical cumulative distribution functions for the source and
    # template images (maps pixel value --> quantile)
    s_quantiles = np.cumsum(s_counts).astype(np.float64)
    s_quantiles /= s_quantiles[-1]
    t_quantiles = np.cumsum(t_counts).astype(np.float64)
    t_quantiles /= t_quantiles[-1]

    #    import matplotlib.pyplot as plt

    #    plt.plot(s_quantiles, label='source')
    #    plt.plot(t_quantiles, label='template')
    #    plt.legend()
    #    plt.show()
    # interpolate linearly to find the pixel values in the template image
    # that correspond most closely to the quantiles in the source image
    interp_t_values = np.interp(s_quantiles, t_quantiles, t_values)
    output = interp_t_values[bin_idx]

    # test number 2
    imhist_source, bins_source = np.histogram(source, n_bins, density=True)
    cdf_source = imhist_source.cumsum()  # cumulative distribution function
    cdf_source = n_bins * cdf_source / cdf_source[-1]  # normalize

    imhist_template, bins_template = np.histogram(template, n_bins, density=True)
    cdf_template = imhist_template.cumsum()  # cumulative distribution function
    cdf_template = n_bins * cdf_template / cdf_template[-1]  # normalize

    im2 = np.interp(source, bins_template[:-1], cdf_source)
    output = np.interp(im2, cdf_template, bins_template[:-1])

    #    plt.plot(cdf_source,label='source')
    #    plt.plot(cdf_template,label='template')

    #    plt.legend()
    #    plt.show()

    output_df = source
    output_df = output_df.to_frame()
    output_df["output"] = output
    not_outs = output_df.columns.to_list()
    not_outs = [i for i in not_outs if i != "output"]
    output_df = output_df.drop(columns=not_outs)

    res = inputdata_test.join(output_df, how="left")
    output_res = res["output"].values

    return output_res

max_stars_repo_name: CFARS/TACT

max_issues_repo_name: CFARS/TACT

max_forks_repo_name: CFARS/TACT


text:
# PyVot Python Variational Optimal Transportation
# Author: Liang Mi <icemiliang@gmail.com>
# Date: April 28th 2020
# Licence: MIT

import os
import sys
import torch
import numpy as np
import matplotlib
from mpl_toolkits import mplot3d
matplotlib.use('Agg')
import matplotlib.pyplot as plt
sys.path.append(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
from vot_torch import VWB
import imageio

np.random.seed(19)

# K = [16, 64, 256]
K = [16, 80, 336]
dot_size = 0.5
dot_size_scale = 50
alpha = 0.3

# -------------- load data ------------------ #

x1 = np.loadtxt("data/mountain1_array.csv", delimiter=',')
x2 = np.loadtxt("data/mountain2_array.csv", delimiter=',')
x3 = np.loadtxt("data/canyon_array.csv", delimiter=',')
x_all = np.loadtxt("data/all_array.csv", delimiter=',')

C1 = np.loadtxt("data/mountain1_C.csv", delimiter=',')
C2 = np.loadtxt("data/mountain2_C.csv", delimiter=',')
C3 = np.loadtxt("data/canyon_C.csv", delimiter=',')
C_all = np.loadtxt("data/all_C.csv", delimiter=',')

idx1 = np.loadtxt("data/mountain1_idx.csv", delimiter=',', dtype=int)
idx2 = np.loadtxt("data/mountain2_idx.csv", delimiter=',', dtype=int)
idx3 = np.loadtxt("data/canyon_idx.csv", delimiter=',', dtype=int)
idx_all = np.loadtxt("data/all_idx.csv", delimiter=',', dtype=int)

C1_16 = C1[:K[0], :]
C2_16 = C2[:K[0], :]
C3_16 = C3[:K[0], :]

idx1_16 = idx1[:, 0] - 1
idx2_16 = idx2[:, 0] - 1
idx3_16 = idx3[:, 0] - 1

C_all_16 = C_all[:K[0], :]
idx_all_16 = idx_all[:, 0] - 1

fig1 = plt.figure(figsize=(8, 8))
ax1 = fig1.add_subplot(111, projection='3d')
ax1.scatter(x1[:, 0], x1[:, 1], x1[:, 2], s=dot_size, color=x1, alpha=alpha)
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.set_xlabel('R')
ax1.set_ylabel('G')
ax1.set_zlabel('B')
# fig1.savefig("x1_histogram.svg", bbox_inches='tight')
fig1.savefig("x1_histogram.png", dpi=300, bbox_inches='tight')

fig2 = plt.figure(figsize=(8, 8))
ax2 = fig2.add_subplot(111, projection='3d')
ax2.scatter(x2[:, 0], x2[:, 1], x2[:, 2], s=dot_size, color=x2, alpha=alpha)
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.set_xlabel('R')
ax2.set_ylabel('G')
ax2.set_zlabel('B')
# fig2.savefig("x2_histogram.svg", bbox_inches='tight')
fig2.savefig("x2_histogram.png", dpi=300, bbox_inches='tight')

fig3 = plt.figure(figsize=(8, 8))
ax3 = fig3.add_subplot(111, projection='3d')
ax3.scatter(x3[:, 0], x3[:, 1], x3[:, 2], s=dot_size, color=x3, alpha=alpha)
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.set_xlabel('R')
ax3.set_ylabel('G')
ax3.set_zlabel('B')
# fig3.savefig("x3_histogram.svg", bbox_inches='tight')
fig3.savefig("x3_histogram.png", dpi=300, bbox_inches='tight')

# -------------- kmeans ------------------ #

x1_kmeans_16 = C1_16[idx1_16, :]
x2_kmeans_16 = C2_16[idx2_16, :]
x3_kmeans_16 = C3_16[idx3_16, :]

x1_kmeans_16 = np.transpose(np.reshape(x1_kmeans_16 * 255, (128, 128, 3)), (1, 0, 2))
x2_kmeans_16 = np.transpose(np.reshape(x2_kmeans_16 * 255, (128, 128, 3)), (1, 0, 2))
x3_kmeans_16 = np.transpose(np.reshape(x3_kmeans_16 * 255, (128, 128, 3)), (1, 0, 2))

imageio.imwrite("x1_kmeans_16.png", x1_kmeans_16)
imageio.imwrite("x2_kmeans_16.png", x2_kmeans_16)
imageio.imwrite("x3_kmeans_16.png", x3_kmeans_16)

fig4 = plt.figure(figsize=(8, 8))
ax4 = fig4.add_subplot(111, projection='3d')
ce1 = C1[idx1_16]
ax4.scatter(x1[:, 0], x1[:, 1], x1[:, 2], s=dot_size, color=ce1, alpha=alpha)
ax4.xaxis.pane.fill = False
ax4.yaxis.pane.fill = False
ax4.zaxis.pane.fill = False
ax4.scatter(C1_16[:, 0], C1_16[:, 1], C1_16[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax4.set_xlabel('R')
ax4.set_ylabel('G')
ax4.set_zlabel('B')
# fig4.savefig("x1_histogram_kmeans1.svg", bbox_inches='tight')
fig4.savefig("x1_histogram_kmeans1.png", dpi=300, bbox_inches='tight')

fig5 = plt.figure(figsize=(8, 8))
ax5 = fig5.add_subplot(111, projection='3d')
ce2 = C2[idx2_16]
ax5.scatter(x2[:, 0], x2[:, 1], x2[:, 2], s=dot_size, color=ce2, alpha=alpha)
ax5.xaxis.pane.fill = False
ax5.yaxis.pane.fill = False
ax5.zaxis.pane.fill = False
ax5.scatter(C2_16[:, 0], C2_16[:, 1], C2_16[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax5.set_xlabel('R')
ax5.set_ylabel('G')
ax5.set_zlabel('B')
# fig5.savefig("x2_histogram_kmeans1.svg", bbox_inches='tight')
fig5.savefig("x2_histogram_kmeans1.png", dpi=300, bbox_inches='tight')

fig6 = plt.figure(figsize=(8, 8))
ax6 = fig6.add_subplot(111, projection='3d')
ce3 = C3[idx3_16]
ax6.scatter(x3[:, 0], x3[:, 1], x3[:, 2], s=dot_size, color=ce3, alpha=alpha)
ax6.xaxis.pane.fill = False
ax6.yaxis.pane.fill = False
ax6.zaxis.pane.fill = False
ax6.scatter(C3_16[:, 0], C3_16[:, 1], C3_16[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax6.set_xlabel('R')
ax6.set_ylabel('G')
ax6.set_zlabel('B')
# fig6.savefig("x3_histogram_kmeans1.svg", bbox_inches='tight')
fig6.savefig("x3_histogram_kmeans1.png", dpi=300, bbox_inches='tight')

# -------------- VWB ------------------ #

use_gpu = False
if use_gpu and torch.cuda.is_available():
    device = 'cuda:0'
else:
    device = 'cpu'

x1 = torch.from_numpy(x1.clip(-0.99, 0.99))
x2 = torch.from_numpy(x2.clip(-0.99, 0.99))
x3 = torch.from_numpy(x3.clip(-0.99, 0.99))

x = torch.from_numpy(C1_16.clip(-0.99, 0.99))
vwb = VWB(x, [x1], device=device, verbose=False)
output = vwb.cluster(lr=1, max_iter_h=5000, max_iter_p=1)
e_idx = output['idx']
x1_vwb_16 = vwb.data_p[e_idx[0], :].cpu().numpy()
x1_vwb_16 = np.transpose(np.reshape(x1_vwb_16 * 255, (128, 128, 3)), (1, 0, 2))
imageio.imwrite("x1_vwb_16.png", x1_vwb_16)
fig1 = plt.figure(figsize=(8, 8))
ax1 = fig1.add_subplot(111, projection='3d')
ce2 = vwb.data_p[e_idx].numpy()
p11 = vwb.data_p.numpy()
ax1.scatter(x1[:, 0], x1[:, 1], x1[:, 2], s=dot_size, color=ce2, alpha=alpha)
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.scatter(p11[:, 0], p11[:, 1], p11[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax1.set_xlabel('R')
ax1.set_ylabel('G')
ax1.set_zlabel('B')
# fig1.savefig("x1_histogram_vwb_16.svg", bbox_inches='tight')
fig1.savefig("x1_histogram_vwb_16.png", dpi=300, bbox_inches='tight')

x = torch.from_numpy(C2_16.clip(-0.99, 0.99))
vwb = VWB(x, [x2], device=device, verbose=False)
output = vwb.cluster(lr=1, max_iter_h=5000, max_iter_p=1)
e_idx = output['idx']
x2_vwb_16 = vwb.data_p[e_idx[0], :].cpu().numpy()
x2_vwb_16 = np.transpose(np.reshape(x2_vwb_16 * 255, (128, 128, 3)), (1, 0, 2))
imageio.imwrite("x2_vwb_16.png", x2_vwb_16)
fig2 = plt.figure(figsize=(8, 8))
ax2 = fig2.add_subplot(111, projection='3d')
ce2 = vwb.data_p[e_idx].numpy()
p21 = vwb.data_p.numpy()
ax2.scatter(x2[:, 0], x2[:, 1], x2[:, 2], s=dot_size, color=ce2, alpha=alpha)
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.scatter(p21[:, 0], p21[:, 1], p21[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax2.set_xlabel('R')
ax2.set_ylabel('G')
ax2.set_zlabel('B')
# fig2.savefig("x2_histogram_vwb_16.svg", bbox_inches='tight')
fig2.savefig("x2_histogram_vwb_16.png", dpi=300, bbox_inches='tight')

x = torch.from_numpy(C3_16.clip(-0.99, 0.99))
vwb = VWB(x, [x3], device=device, verbose=False)
output = vwb.cluster(lr=1, max_iter_h=5000, max_iter_p=1)
e_idx = output['idx']
x3_vwb_16 = vwb.data_p[e_idx[0], :].cpu().numpy()
x3_vwb_16 = np.transpose(np.reshape(x3_vwb_16 * 255, (128, 128, 3)), (1, 0, 2))
imageio.imwrite("x3_vwb_16.png", x3_vwb_16)
fig3 = plt.figure(figsize=(8, 8))
ax3 = fig3.add_subplot(111, projection='3d')
ce3 = vwb.data_p[e_idx].numpy()
p31 = vwb.data_p.numpy()
ax3.scatter(x3[:, 0], x3[:, 1], x3[:, 2], s=dot_size, color=ce3, alpha=alpha)
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.scatter(p31[:, 0], p31[:, 1], p31[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax3.set_xlabel('R')
ax3.set_ylabel('G')
ax3.set_zlabel('B')
# fig3.savefig("x3_histogram_vwb_16.svg", bbox_inches='tight')
fig3.savefig("x3_histogram_vwb_16.png", dpi=300, bbox_inches='tight')

# -------------- VWB ALL ------------------ #

x = torch.from_numpy(C_all_16.clip(-0.99, 0.99))
vwb = VWB(x, [x1, x2, x3], device=device, verbose=False)
output = vwb.cluster(lr=1, max_iter_h=5000, max_iter_p=1)
e_idx = output['idx']

x1_vwb_all = vwb.data_p[e_idx[0], :].cpu().numpy()
x2_vwb_all = vwb.data_p[e_idx[1], :].cpu().numpy()
x3_vwb_all = vwb.data_p[e_idx[2], :].cpu().numpy()
x1_vwb_all = np.transpose(np.reshape(x1_vwb_all*255, (128, 128, 3)), (1, 0, 2))
x2_vwb_all = np.transpose(np.reshape(x2_vwb_all*255, (128, 128, 3)), (1, 0, 2))
x3_vwb_all = np.transpose(np.reshape(x3_vwb_all*255, (128, 128, 3)), (1, 0, 2))
imageio.imwrite("x1_vwb_all.png", x1_vwb_all)
imageio.imwrite("x2_vwb_all.png", x2_vwb_all)
imageio.imwrite("x3_vwb_all.png", x3_vwb_all)

fig1 = plt.figure(figsize=(8, 8))
ax1 = fig1.add_subplot(111, projection='3d')
ce1 = vwb.data_p[e_idx[0]].numpy()
p1 = vwb.data_p.numpy()
ax1.scatter(x1[:, 0], x1[:, 1], x1[:, 2], s=dot_size, color=ce1, alpha=alpha)
ax1.xaxis.pane.fill = False
ax1.yaxis.pane.fill = False
ax1.zaxis.pane.fill = False
ax1.scatter(p1[:, 0], p1[:, 1], p1[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax1.set_xlabel('R')
ax1.set_ylabel('G')
ax1.set_zlabel('B')
# fig1.savefig("x1_histogram_vwb1all.svg", bbox_inches='tight')
fig1.savefig("x1_histogram_vwb1all.png", dpi=300, bbox_inches='tight')

fig2 = plt.figure(figsize=(8, 8))
ax2 = fig2.add_subplot(111, projection='3d')
ce2 = vwb.data_p[e_idx[1]].numpy()
p2 = vwb.data_p.numpy()
ax2.scatter(x2[:, 0], x2[:, 1], x2[:, 2], s=dot_size, color=ce2, alpha=alpha)
ax2.xaxis.pane.fill = False
ax2.yaxis.pane.fill = False
ax2.zaxis.pane.fill = False
ax2.scatter(p2[:, 0], p2[:, 1], p2[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax2.set_xlabel('R')
ax2.set_ylabel('G')
ax2.set_zlabel('B')
# fig2.savefig("x2_histogram_vwb1all.svg", bbox_inches='tight')
fig2.savefig("x2_histogram_vwb1all.png", dpi=300, bbox_inches='tight')

fig3 = plt.figure(figsize=(8, 8))
ax3 = fig3.add_subplot(111, projection='3d')
ce3 = vwb.data_p[e_idx[2]].numpy()
p3 = vwb.data_p.numpy()
ax3.scatter(x3[:, 0], x3[:, 1], x3[:, 2], s=dot_size, color=ce3, alpha=alpha)
ax3.xaxis.pane.fill = False
ax3.yaxis.pane.fill = False
ax3.zaxis.pane.fill = False
ax3.scatter(p3[:, 0], p3[:, 1], p3[:, 2], s=dot_size*dot_size_scale, color='k', zorder=1)
ax3.set_xlabel('R')
ax3.set_ylabel('G')
ax3.set_zlabel('B')
# fig3.savefig("x3_histogram_vwb1all.svg", bbox_inches='tight')
fig3.savefig("x3_histogram_vwb1all.png", dpi=300, bbox_inches='tight')

ext: py

max_stars_repo_name: icemiliang/pyvot

max_issues_repo_name: icemiliang/pyvot

max_forks_repo_name: icemiliang/pyvot

reason: import numpy


text:
import datetime

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

def clean_pair(pair):
    if len(pair) > 2:
        pair = [pair[0], "=".join(pair[1:])]
        return pair
    elif len(pair) <= 1:
        raise TypeError(
            "Pair must be an iterable of length 2. Is there a blank row in your .env file?")
    return pair

def load_env(path):
    """Load environment variables from .env."""
    with open(path) as f:
        lines = f.readlines()
    pairs = [x.split("=") for x in lines]
    pairs = [clean_pair(p) for p in pairs]
    env = {k: v.replace("\n", "") for k, v in pairs}
    try:
        return env['DATA_DIR'], env['IMAGE_DIR'], env['SAVE_DATA_DIR']
    except KeyError as e:
        print("Environment variable missing:")
        print(e)

_DATA_DIR, _IMAGE_DIR, _SAVE_DATA_DIR = load_env("../.env")

def load_population():
    """Load and clean population data."""
    population = pd.read_csv(
        _DATA_DIR + "population_2018.csv", encoding="ISO-8859-1")
    # fix for Namibia (country code NA)
    population["country_code"].fillna("NA", inplace=True)
    return population

def load_country_codes():
    """Load and clean UN Country Codes."""
    un_codes = pd.read_csv(
        _DATA_DIR + "un_countrycode_new.csv", encoding="ISO-8859-1")
    un_codes = un_codes.fillna("NA")
    return un_codes

def load_members():
    """Load and clean members approval data."""
    members = pd.read_csv(_DATA_DIR + "approval.csv")
    # clean up data column names
    members.columns = (
        members.columns.to_series()
        .str.strip()
        .str.lower()
        .str.replace(" ", "_", regex=False)
        .str.replace("(", "", regex=False)
        .str.replace(")", "", regex=False)
    )
    # format dates for sorting
    members["application_date"] = pd.to_datetime(members["application_date"])
    members["approval_date"] = pd.to_datetime(members["approval_date"])
    members["acceptance_date"] = pd.to_datetime(members["acceptance_date"])
    # convert ids to int, set id NA values to 0
    members["member_id"].fillna("0", inplace=True)
    members["member_id"] = members.member_id.astype(int)
    # fix for Namibia (country code NA)
    members["country_code"].fillna("NA", inplace=True)
    return members

def process_apps_by_date(members):
    """Build applicants by date aggregates."""
    # group and sort applicants and approved members by date
    apps_by_date = (
        members.groupby(["application_date"])
        .application_date.count()
        .reset_index(name="applicants")
        .sort_values(["application_date"])
    )
    return apps_by_date

def process_exp(members):
    """Build member counts with health experience."""
    approved = members[members.approval_date.notnull()]
    exp = (
        approved.groupby(["health_experience"])
        .health_experience.count()
        .reset_index(name="count")
    )

    # total experience
    total_exp = exp["count"].sum()
    # add column for percent values
    exp["percent"] = exp["count"] * 100 / total_exp
    exp.percent = exp.percent.round().astype(int)
    exp.columns = ["Experience", "Members", "%"]
    return exp

def process_members(members, population, un_codes):
    """Build merged dataframes split by members/no members."""
    approved = members.loc[members["user_status"] == "Approved"].copy()
    app_country = (
        approved.groupby(["country", "country_code", "who_region"])
        .country.count()
        .reset_index(name="applicants")
        .sort_values(["country"])
    )
    # merge country and population tables
    app_country_pop = app_country.merge(
        population.drop("country", axis=1), how="left", on="country_code"
    )

    # merge with 3 letter UN country codes
    mem_countries = app_country_pop.merge(
        un_codes.drop("country", axis=1), how="left", on="country_code"
    )

    # include members density (responders per 1 million population) and sort by density
    mem_countries["member_density"] = (
        1000 * mem_countries["applicants"] / mem_countries["population"]
    ).round(2)
    mem_countries = mem_countries.sort_values(["member_density"])

    # get countries with no members
    all_country = un_codes.merge(
        app_country.drop("country", axis=1), how="left", on="country_code"
    )
    all_country.sort_values(["applicants"], inplace=True)
    mask = all_country["applicants"].isnull()
    no_mem_countries = all_country.loc[mask].copy()
    no_mem_countries = no_mem_countries[["country", "un_country_code"]]
    no_mem_countries = no_mem_countries[no_mem_countries.un_country_code.notnull(
    )]
    return mem_countries, no_mem_countries

def process_new_region(members):
    """Build members by WHO region dataframe."""
    approved = members.loc[members["user_status"] == "Approved"].copy()
    # group by region, count members in regions, and sort by region
    app_region = (
        approved.groupby(["who_region"])
        .who_region.count()
        .reset_index(name="applicants")
        .sort_values(["who_region"])
    )

    # re-arrange regions
    central_america = app_region[
        app_region["who_region"].str.lower().str.contains("central america")
    ]
    central_america_total = central_america["applicants"].sum()
    app_region.at[5, "who_region"] = "Central-South America"
    app_region.at[5, "applicants"] = int(central_america_total)
    new_region = app_region.drop([3, 7])
    new_region.rename(columns={"who_region": "Region"}, inplace=True)
    return new_region

def get_masks(df, dt_col):
    """Get boolean masks for filtering by timeframe."""
    last_month, year = get_last_month_year()

    last_month = (df[dt_col].dt.month == last_month) & (
        df[dt_col].dt.year == year)
    ytd = df[dt_col].dt.year == year
    last_year = df[dt_col].dt.year == (year - 1)

    return last_month, ytd, last_year

def get_last_month_year():
    """Get last full month and associated year."""
    today = datetime.date.today()
    first_day_of_month = today.replace(day=1)
    last_day_last_month = first_day_of_month - datetime.timedelta(days=1)
    last_month = last_day_last_month.month
    year = last_day_last_month.year
    return last_month, year

def get_total_applicants_month(apps_by_date):
    """Get total count of applicants for last full month."""
    last_month, _, _ = get_masks(apps_by_date, "application_date")
    applicants_month = apps_by_date.loc[last_month].copy()
    total_apps = applicants_month["applicants"].sum()
    return total_apps

def get_approved_totals(members):
    """Get total count of approvals for all time and last month."""
    approved_by_date = (
        members.groupby(["approval_date"])
        .approval_date.count()
        .reset_index(name="approved")
        .sort_values(["approval_date"])
    )
    last_month, _, _ = get_masks(approved_by_date, "approval_date")
    approved_month = approved_by_date.loc[last_month].copy()

    tot_appr = approved_by_date["approved"].sum()
    tot_appr_month = approved_month["approved"].sum()
    return tot_appr, tot_appr_month

def get_comma_separated(codes):
    """Takes list of codes, returns comma separated string."""
    return ", ".join(codes)

def plot_experience(exp):
    """Save chart of members by health experience."""
    # get total human experience, filtered by "human"
    human_exp = exp[exp["Experience"].str.lower().str.contains("human")]
    total_human_exp = human_exp["%"].sum()
    # get total animal experience, filtered by "animal:
    animal_exp = exp[exp["Experience"].str.lower().str.contains("animal")]
    total_animal_exp = animal_exp["%"].sum()
    # get total environmental experience, filtered by "environmental"
    environmental_exp = exp[exp["Experience"].str.lower(
    ).str.contains("environmental")]
    total_environmental_exp = environmental_exp["%"].sum()
    # plot expertise bar chart
    fig2 = plt.figure()
    x = np.arange(3)
    y = [total_human_exp, total_animal_exp, total_environmental_exp]
    hm, an, en = plt.bar(x, y, align="center", edgecolor="none")
    hm.set_facecolor("#015163")
    an.set_facecolor("#3fc8c8")
    en.set_facecolor("#2e3335")
    plt.xticks(x, ("Human", "Animal", "Environmental"))
    plt.tick_params(
        top=False,
        bottom=False,
        left=False,
        right=False,
        labelleft=False,
        labelbottom=True,
        labelsize=14,
    )
    plt.box(False)
    # add value labels to bars
    for a, b in zip(x, y):
        plt.text(a, b + 2, str(b) + "%", ha="center",
                 color="#2e3335", fontsize=18)
    fig2.savefig(_IMAGE_DIR + "health_expertise.svg",
                 bbox_inches="tight", format="svg")

def plot_applicants(apps_by_date, year):
    """Save plot of applicants by date."""
    fig = plt.figure()
    plt.plot(
        apps_by_date.application_date,
        apps_by_date.applicants,
        color="#3fc8c8",
    )
    plt.tick_params(
        top=False,
        bottom=False,
        left=False,
        right=False,
        labelleft=True,
        labelbottom=True,
        labelsize=12,
        labelcolor="#878c8d",
    )
    plt.xticks(rotation=20)
    plt.gca().spines["right"].set_color("none")
    plt.gca().spines["top"].set_color("none")
    plt.gca().spines["left"].set_color("#d5d7d8")
    plt.gca().spines["bottom"].set_color("#d5d7d8")
    plt.xlim([datetime.date((year - 1), 1, 1), datetime.datetime.now()])
    plt.ylim([0, 10])
    fig.savefig(_IMAGE_DIR + "applicants_year.svg",
                format="svg", bbox_inches="tight")

def plot_applicants_last_month(apps_by_date):
    """Save plot of applicants by date for last month."""
    last_month, year = get_last_month_year()
    last_month_mask, _, _ = get_masks(apps_by_date, "application_date")
    applicants_plot = apps_by_date.loc[last_month_mask].copy()
    # generate plot
    fig = plt.figure()
    plt.plot(
        applicants_plot.application_date,
        applicants_plot.applicants,
        color="#3fc8c8",
    )
    plt.tick_params(
        top=False,
        bottom=False,
        left=False,
        right=False,
        labelleft=True,
        labelbottom=True,
        labelsize=12,
        labelcolor="#878c8d",
    )
    plt.xticks(rotation=20)
    plt.gca().spines["right"].set_color("none")
    plt.gca().spines["top"].set_color("none")
    plt.gca().spines["left"].set_color("#d5d7d8")
    plt.gca().spines["bottom"].set_color("#d5d7d8")
    plt.xlim([datetime.date(year, last_month, 1), datetime.datetime.now()])
    plt.ylim([0, 10])
    fig.savefig(_IMAGE_DIR + "applicants_month.svg",
                format="svg", bbox_inches="tight")

def create_experience_table(exp):
    """Export experience report."""
    exp.to_html(_SAVE_DATA_DIR + "experience_table.html", index=False)

def create_app_no_training_table(members):
    """Export accepted applicants w/ no training report."""
    today = datetime.datetime.now()
    two_months_ago = today - datetime.timedelta(days=60)
    mask = (
        (members["user_status"] == "Accepted")
        & (members["course_type"].isnull())
        & (members["acceptance_date"] > two_months_ago)
    )
    app_no_training = members.loc[mask].copy()
    app_no_training.rename(
        columns={"acceptance_date": "Acceptance Date",
                 "member_id": "members ID"},
        inplace=True,
    )
    app_no_training = app_no_training[["Acceptance Date", "members ID"]]
    app_no_training.to_html(
        _SAVE_DATA_DIR + "app_no_training_table.html", index=False)

def create_heard_about_table(members):
    """Export heard about Epicore report."""
    last_month_mask, _, _ = get_masks(members, "application_date")
    members_month = members.loc[
        last_month_mask,
    ].copy()
    members_month.rename(
        columns={"heard_about_epicore_by": "Source"}, inplace=True)
    heard_about = (
        members_month.groupby(["Source"])
        .Source.count()
        .reset_index(name="Applicants")
        .sort_values(["Source"])
    )
    heard_about.to_html(_SAVE_DATA_DIR + "heard_about_table.html", index=False)

def create_new_applicants_table(members):
    """Export new applicants report."""
    members = members[
        ["application_date", "approval_date", "country_code", "country"]
    ].copy()

    new_by_country = (
        members.groupby(["country", "approval_date"])
        .country.count()
        .reset_index(name="Members")
        .sort_values(["country"])
    )
    last_month_mask, _, _ = get_masks(new_by_country, "approval_date")
    new_last_month = new_by_country.loc[last_month_mask].copy()

    report = new_last_month.groupby(["country"]).sum().reset_index()
    report.to_html(_SAVE_DATA_DIR + "new_applicants_table.html", index=False)

def create_no_member_countries_table(no_mem_countries):
    """Export countries with no members report."""
    no_mem_countries.to_html(
        _SAVE_DATA_DIR + "no_applicants_countries_table.html", index=False
    )

def create_no_members_region_table(members, no_mem_countries, un_codes):
    """Export regions with no members report."""
    approved = members.loc[members["user_status"] == "Approved"].copy()
    # merge no members countries with regions
    no_member_regions = no_mem_countries.merge(
        un_codes, how="left", on="un_country_code"
    )
    no_member_regions = no_member_regions[
        ["country_x", "country_code", "un_country_code", "who_region"]
    ]
    no_member_regions.rename(
        columns={"country_x": "country", "who_region": "Region"}, inplace=True
    )

    # group by region
    no_members_reg = no_member_regions.groupby("Region").agg(
        {"un_country_code": get_comma_separated}
    )

    # Merge current approval data to the new regions list
    no_member_merged = approved.merge(un_codes, how="left", on="country_code")

    # Pick few columns for ease
    no_member_merged = no_member_merged[
        [
            "approval_date",
            "country_code",
            "un_country_code",
            "who_region_x",
            "who_region_y",
        ]
    ]

    # Rename Columns
    no_member_merged.rename(
        columns={
            "approval_date": "approval_date",
            "country_code": "country_code",
            "un_country_code": "un_country_code",
            "who_region_x": "who_region_old",
            "who_region_y": "Region",
        },
        inplace=True,
    )

    # Group by Region and sum up the total applicants from that region
    po_grouped = (
        no_member_merged.groupby(["Region"])
        .Region.count()
        .reset_index(name="applicants")
        .sort_values(["Region"])
    )

    # Finally Merge the above output to the list of NO members countries in the region
    report = po_grouped.merge(no_members_reg, how="left", on="Region")

    total_region_applicants = report.applicants.sum()

    report["Members(%)"] = (
        (100 * report["applicants"] /
         total_region_applicants).round().astype(int)
    )
    report.rename(
        columns={
            "applicants": "Members",
            "un_country_code": "Countries with No Members",
        },
        inplace=True,
    )
    report.fillna("", inplace=True)

    report.to_html(_SAVE_DATA_DIR +
                   "no_members_region_table.html", index=False)
    return report

def create_country_table(mem_countries):
    """Export applicants by country report."""
    app_country_density = mem_countries[
        ["country", "applicants", "member_density"]
    ].copy()
    app_country_density.rename(
        columns={"country": "Country", "applicants": "n",
                 "member_density": "Density"},
        inplace=True,
    )
    app_country_density.to_html(
        _SAVE_DATA_DIR + "country_table.html", index=False)

def create_region_table(new_region):
    """Export region table report."""
    # data frame for members by region merged with counties with no members
    # calculate %
    total_region_applicants = new_region.applicants.sum()
    new_region["Members(%)"] = (
        (100 * new_region["applicants"].astype(float) /
         float(total_region_applicants)).round(1)
    )
    new_region.rename(columns={"applicants": "Members"}, inplace=True)
    new_region.to_html(_SAVE_DATA_DIR + "region_table.html", index=False)

def create_membership_table(
    tot_appr,
    total_apps,
    tot_appr_month,
    total_member_countries,
):
    """Export membership summary report."""
    # creat data frame for memebership summary
    data = [
        ["Responders", tot_appr],
        ["New Applicants", str(total_apps)],
        ["New Responders", str(tot_appr_month)],
        ["Countries Represented", str(total_member_countries)],
    ]

    membership = pd.DataFrame(data, columns=["Metric", "Value"])
    membership.to_csv(_SAVE_DATA_DIR + "membership.csv", sep="|", index=False)
    membership.to_html(_SAVE_DATA_DIR + "membership.html", index=False)

def create_members_regions_table(total_member_countries, total_no_member_countries):
    """Export regions with/without members report."""
    data_region = [
        [
            "Countries and territories included: " +
            str(total_member_countries),
            "Countries and territories missing: " +
            str(total_no_member_countries),
        ]
    ]

    data_region = pd.DataFrame(data_region, columns=["", ""])
    data_region.to_html(_SAVE_DATA_DIR + "members_regions.html", index=False)

def main():
    """Main function for generating plots/reports."""
    population = load_population()
    un_codes = load_country_codes()
    members = load_members()

    exp = process_exp(members)
    apps_by_date = process_apps_by_date(members)
    mem_countries, no_mem_countries = process_members(
        members, population, un_codes)
    new_region = process_new_region(members)

    _, year = get_last_month_year()
    total_apps = get_total_applicants_month(apps_by_date)
    tot_appr, tot_appr_month = get_approved_totals(members)

    plot_experience(exp)
    plot_applicants(apps_by_date, year)
    plot_applicants_last_month(apps_by_date)

    create_experience_table(exp)
    create_app_no_training_table(members)
    create_heard_about_table(members)
    create_new_applicants_table(members)
    create_no_members_region_table(members, no_mem_countries, un_codes)
    create_country_table(mem_countries)
    create_region_table(new_region)
    create_membership_table(
        tot_appr,
        total_apps,
        tot_appr_month,
        mem_countries.shape[0],
    )
    create_members_regions_table(
        mem_countries.shape[0], no_mem_countries.shape[0])

if __name__ == "__main__":
    main()

max_stars_repo_name: healthmap/epicore-open

max_issues_repo_name: healthmap/epicore-open

max_forks_repo_name: healthmap/epicore-open


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 28 10:10:25 2020

@author: bressler
"""

import SBCcode as sbc
import numpy as np
import scipy
import runlistscatalogue as rlc
import os
import pulse_integrator as pi
import matplotlib.pyplot as plt
plt.style.use('seaborn-darkgrid')
import gc
from gaincalc import get_gain

CONVERSION_TO_CHARGE = (125.0/128)*(1/50.0)*(1/1000.0)*(1/(1.602e-19))
datapath = '/bluearc/storage/SBC-17-data/'
"""
run='20170628_9'
runreconpath = "/pnfs/coupp/persistent/grid_output/SBC-17/output/%s/"%run
timingname = runreconpath+"TimingAnalysis_%s.bin"%run
timing = sbc.DataHandling.ReadBinary.ReadBlock(timingname)
print(timing.keys())

#runs = rlc.CfCombinedMultiTemp

with open("/nashome/b/bressler/sbcoutput/20170628_9_merged.txt","r") as fin:
    data = fin.readlines()

bubdata = {}

headers = data[0].split()

runind = headers.index("run")
evind = headers.index("event")
xind = headers.index("x")
yind = headers.index("y")
zind = headers.index("z")
pmtt0ind = headers.index("PMTt0")
pmttraceind = headers.index("PMTt0_ind")
lagind = headers.index("lag")
spectind = headers.index("PMTphe")
blockedind = headers.index('isBlocked')
print(headers)
x=[]
y=[]
z=[]
lag = []
spect = []
theirspect = []
goodsrcbubbles = 0
i=0
for line in data:
    if i>0:
        split_line = line.split()
        if (not np.isnan(float(split_line[xind]))) and float(split_line[blockedind])==0:
            goodsrcbubbles += 1
            evid = split_line[runind]+"-"+split_line[evind]
            print(evid)
            print(i)
            x.append(float(split_line[xind]))
            y.append(float(split_line[yind]))
            z.append(float(split_line[zind]))
            lag.append(float(split_line[lagind]))
            spect.append(float(split_line[spectind]))
            theirspect.append(timing["PMTmatch_nphe"][i-1][0])
            bubdata[evid]=[float(split_line[spectind]),float(split_line[pmttraceind])]

    i+=1
plt.figure()
plt.scatter(spect,theirspect,18)
plt.xlabel('My pulse integration',fontsize=18)
plt.ylabel('PMTmatch_nphe',fontsize=18)
plt.show

#print(bubdata)
"""

run='20170804_2'
runreconpath = "/pnfs/coupp/persistent/grid_output/SBC-17/output/%s/"%run
runpath = datapath+run
g = get_gain(datapath, run)
print("gain: %e"%g)
events = [evnt for evnt in os.listdir(runpath) if not os.path.isfile(os.path.join(runpath,evnt))]
Nevents = len(events)

pmtfilename = runreconpath+"PMTpulseAnalysis_%s.bin"%run
pmt = sbc.DataHandling.ReadBinary.ReadBlock(pmtfilename)

evdone = []
myphe = [0 for i in range(len(pmt["PMT_pulse_area"]))]
my_Vmax = [0 for i in range(len(pmt["PMT_pulse_area"]))]

wNIM = [0 for i in range(len(pmt["PMT_pulse_area"]))]

theirphe = [0 for i in range(len(pmt["PMT_pulse_area"]))]
ratio = [0 for i in range(len(pmt["PMT_pulse_area"]))]
j=0
for n in range(1):
    print("event %d"%n)
    e = sbc.DataHandling.GetSBCEvent.GetEvent(runpath, n)
    print(len(e["PMTtraces"]["traces"]))
    for k in range(len(e["PMTtraces"]["traces"])):
        trace = np.fabs(e["PMTtraces"]["traces"][k][0])
        trace = trace[:-50]

        #if ch0 saturated, stitch in low res channel:
        if max(trace) == 128:
            trace = pi.stitchTraces(trace, np.fabs(e["PMTtraces"]["traces"][k][1]))
        dt = e["PMTtraces"]["dt"][k][0]
        baseline = np.mean(trace[:50])
        trace -= baseline
        my_Vmax[j] = max(trace)
        #integrate and convert to phe:
        [a,npeaks,totInt,pktimes] = pi.SBC_pulse_integrator_bressler(trace,dt)

        caught = False
        if min(e["PMTtraces"]["traces"][k][1]) < -0.5:
            caught = True

        myphe[j]=a/g

        if npeaks >= 1 and a/6 < 2 and not caught:
            print("event %d trace %d (had %d peaks)"%(n, k, npeaks))
            baseline_std = np.std(trace[0:50])
            tPMT = np.arange(len(trace))*dt
            integral_t0_index = np.argmax(np.diff(trace)>4)
            integral_t0 = tPMT[integral_t0_index]
            p,t = pi.get_pulse(trace,tPMT,dt, 500e-9, integral_t0,baseline_std)
            ret = scipy.integrate.trapz(p,t)*CONVERSION_TO_CHARGE
            print("ret: %f"%ret)
            print("a: %f"%a)

            plt.figure()
            plt.scatter(t,p, 10, c='r')
            plt.plot(tPMT,trace, c='b', lw=1)
            plt.plot([integral_t0, integral_t0], [0, max(trace)], 'g')
            plt.xlim([0, 7.5e-7])
            plt.title('%e phe'%(a/g))
            plt.plot()

        j+= 1
        gc.collect()

#ft = np.polyfit(myphe,ratio,0)
#print(ft)
plt.figure()
plt.scatter(myphe, my_Vmax)
plt.title('my phe vs vmax')
plt.xlabel('photoelectrons')
plt.ylabel('trace maximum (ADC)')
plt.show()

plt.figure()
plt.hist2d(myphe, my_Vmax,150)
plt.title('my phe vs vmax')
plt.xlabel('photoelectrons')
plt.ylabel('trace maximum (ADC)')
plt.show()

"""
run='20170628_9'
runreconpath = "/pnfs/coupp/persistent/grid_output/SBC-17/output/%s/"%run
runpath = datapath+run
events = [evnt for evnt in os.listdir(runpath) if not os.path.isfile(os.path.join(runpath,evnt))]
Nevents = len(events)

pmtfilename = runreconpath+"PMTpulseAnalysis_%s.bin"%run
pmt = sbc.DataHandling.ReadBinary.ReadBlock(pmtfilename)
#print(len(pmt["PMT_pulse_area"]))
print(pmt.keys())
#evid = pmt["runid"]+"-"+pmt["ev"]
#area = pmt["PMT_pulse_area"]
#bubdata[]
evdone = []
myphe = [0 for i in range(len(pmt["PMT_pulse_area"]))]
my_Vmax = [0 for i in range(len(pmt["PMT_pulse_area"]))]
theirphe = [0 for i in range(len(pmt["PMT_pulse_area"]))]
ratio = [0 for i in range(len(pmt["PMT_pulse_area"]))]
j=0
for n in range(1):
    print("event %d"%n)
    indices = [i for i,x in enumerate(pmt["ev"]) if int(x)==n]
    #print(indices)
    areas=pmt["PMT_pulse_area"][indices]
    #print(len(areas))
    e = sbc.DataHandling.GetSBCEvent.GetEvent(runpath,n)
    print(len(e["PMTtraces"]["traces"]))
    for k in range(len(e["PMTtraces"]["traces"])):
        trace = np.fabs(e["PMTtraces"]["traces"][k][0])
        #if ch0 saturated, stitch in low res channel:
        if max(trace) == 128:
            trace = pi.stitchTraces(trace,np.fabs(e["PMTtraces"]["traces"][k][1]))
        dt = e["PMTtraces"]["dt"][k][0]
        my_Vmax[k] = max(trace)/e["PMTtraces"]["v_scale"][k,0,None]
        #integrate and convert to phe:
        [a,n,totInt,pktimes] = pi.SBC_pulse_integrator_bressler(trace,dt)
        myphe[j]=a
        theirphe[j]=-areas[k]*CONVERSION_TO_CHARGE/e["PMTtraces"]["v_scale"][k,0,None]
        ratio[j] = myphe[j]/theirphe[j]

        j+= 1
        gc.collect()

print(len(myphe))
print(j)
#ft = np.polyfit(myphe,ratio,0)
#print(ft)

plt.figure()
plt.scatter(myphe, my_Vmax)
plt.title('my phe vs vmax')
plt.xlabel('photoelectrons')
plt.ylabel('trace maximum (V)')
plt.show()

plt.figure()
plt.scatter(myphe,theirphe)
plt.xlabel("My Area [elecrons]")
plt.ylabel("Their area [electrons]")
plt.xlim([min(myphe),max(myphe)])
plt.ylim([min(theirphe),max(theirphe)])
plt.show()

plt.figure()
plt.scatter(myphe,ratio)
#plt.plot([min(myphe),max(myphe)],[ft[0],ft[0]])
plt.show()

plt.figure()
plt.hist(theirphe)
plt.show()

"""

max_stars_repo_name: cericdahl/SBCcode

max_issues_repo_name: SBC-Collaboration/SBC-Analysis

max_forks_repo_name: SBC-Collaboration/SBC-Analysis

reason: import numpy,import scipy


text:
import numpy as np
import pytest
import xarray as xr
from xarray.testing import assert_equal

from climpred.comparisons import (
    __ALL_COMPARISONS__ as all_comparisons,
    PM_COMPARISONS,
    PROBABILISTIC_PM_COMPARISONS,
    Comparison,
    __e2c,
    __m2c,
    __m2e,
    __m2m,
)
from climpred.metrics import PM_METRICS, __mse as metric
from climpred.prediction import compute_perfect_model
from climpred.utils import get_comparison_class, get_metric_class

def test_e2c(PM_da_initialized_1d):
    """Test ensemble_mean-to-control (which can be any other one member) (e2c)
    comparison basic functionality.

    Clean comparison: Remove one control member from ensemble to use as reference.
    Take the remaining member mean as forecasts."""
    ds = PM_da_initialized_1d
    aforecast, areference = __e2c.function(ds, metric=metric)

    control_member = ds.member.values[0]
    reference = ds.sel(member=control_member, drop=True)
    # drop the member being reference
    ds = ds.drop_sel(member=control_member)
    forecast = ds.mean("member")

    eforecast, ereference = forecast, reference
    # very weak testing on shape
    assert eforecast.size == aforecast.size
    assert ereference.size == areference.size

    assert_equal(eforecast, aforecast)
    assert_equal(ereference, areference)

def test_m2c(PM_da_initialized_1d):
    """Test many-to-control (which can be any other one member) (m2c) comparison basic
    functionality.

    Clean comparison: Remove one control member from ensemble to use as reference.
    Take the remaining members as forecasts."""
    ds = PM_da_initialized_1d
    aforecast, areference = __m2c.function(ds, metric=metric)

    control_member = ds.member.values[0]
    reference = ds.sel(member=control_member, drop=True)
    # drop the member being reference
    ds_dropped = ds.drop_sel(member=control_member)
    forecast, reference = xr.broadcast(ds_dropped, reference)

    eforecast, ereference = forecast, reference
    # very weak testing on shape
    assert eforecast.size == aforecast.size
    assert ereference.size == areference.size

    assert_equal(eforecast, aforecast)
    assert_equal(ereference, areference)

def test_m2e(PM_da_initialized_1d):
    """Test many-to-ensemble-mean (m2e) comparison basic functionality.

    Clean comparison: Remove one member from ensemble to use as reference.
    Take the remaining members as forecasts."""
    ds = PM_da_initialized_1d
    aforecast, areference = __m2e.function(ds, metric=metric)

    reference_list = []
    forecast_list = []
    for m in ds.member.values:
        forecast = ds.drop_sel(member=m).mean("member")
        reference = ds.sel(member=m, drop=True)
        forecast, reference = xr.broadcast(forecast, reference)
        forecast_list.append(forecast)
        reference_list.append(reference)
    reference = xr.concat(reference_list, "member")
    forecast = xr.concat(forecast_list, "member")
    forecast["member"] = np.arange(forecast.member.size)
    reference["member"] = np.arange(reference.member.size)

    eforecast, ereference = forecast, reference
    # very weak testing on shape
    assert eforecast.size == aforecast.size
    assert ereference.size == areference.size

    assert_equal(eforecast, aforecast)
    assert_equal(ereference, areference)

def test_m2m(PM_da_initialized_1d):
    """Test many-to-many (m2m) comparison basic functionality.

    Clean comparison: Remove one member from ensemble to use as reference. Take the
    remaining members as forecasts."""
    ds = PM_da_initialized_1d
    aforecast, areference = __m2m.function(ds, metric=metric)

    reference_list = []
    forecast_list = []
    for m in ds.member.values:
        forecast = ds.drop_sel(member=m)
        forecast["member"] = np.arange(1, 1 + forecast.member.size)
        reference = ds.sel(member=m, drop=True)
        forecast, reference = xr.broadcast(forecast, reference)
        reference_list.append(reference)
        forecast_list.append(forecast)
    supervector_dim = "forecast_member"
    reference = xr.concat(reference_list, supervector_dim)
    forecast = xr.concat(forecast_list, supervector_dim)
    reference[supervector_dim] = np.arange(reference[supervector_dim].size)
    forecast[supervector_dim] = np.arange(forecast[supervector_dim].size)
    eforecast, ereference = forecast, reference
    # very weak testing here
    assert eforecast.size == aforecast.size
    assert ereference.size == areference.size

@pytest.mark.parametrize("metric", ["crps", "mse"])
@pytest.mark.parametrize("comparison", PM_COMPARISONS)
def test_all(PM_da_initialized_1d, comparison, metric):
    metric = get_metric_class(metric, PM_METRICS)
    ds = PM_da_initialized_1d
    comparison = get_comparison_class(comparison, PM_COMPARISONS)
    forecast, obs = comparison.function(ds, metric=metric)
    assert not forecast.isnull().any()
    assert not obs.isnull().any()
    if not metric.probabilistic:
        # same dimensions for deterministic metrics
        assert forecast.dims == obs.dims
    else:
        if comparison.name in PROBABILISTIC_PM_COMPARISONS:
            # same but member dim for probabilistic
            assert set(forecast.dims) - set(["member"]) == set(obs.dims)

def my_m2me_comparison(ds, metric=None):
    """Identical to m2e but median."""
    reference_list = []
    forecast_list = []
    supervector_dim = "member"
    for m in ds.member.values:
        forecast = ds.drop_sel(member=m).median("member")
        reference = ds.sel(member=m).squeeze()
        forecast_list.append(forecast)
        reference_list.append(reference)
    reference = xr.concat(reference_list, supervector_dim)
    forecast = xr.concat(forecast_list, supervector_dim)
    forecast[supervector_dim] = np.arange(forecast[supervector_dim].size)
    reference[supervector_dim] = np.arange(reference[supervector_dim].size)
    return forecast, reference

my_m2me_comparison = Comparison(
    name="m2me", function=my_m2me_comparison, probabilistic=False, hindcast=False,
)

@pytest.mark.parametrize("metric", ("rmse", "pearson_r"))
def test_new_comparison_passed_to_compute(
    PM_da_initialized_1d, PM_da_control_1d, metric
):
    actual = compute_perfect_model(
        PM_da_initialized_1d,
        PM_da_control_1d,
        comparison=my_m2me_comparison,
        metric=metric,
    )

    expected = compute_perfect_model(
        PM_da_initialized_1d, PM_da_control_1d, comparison="m2e", metric="mse"
    )

    assert (actual - expected).mean() != 0

def test_Comparison_display():
    summary = __m2m.__repr__()
    assert "Kind: deterministic and probabilistic" in summary.split("\n")[2]

def test_no_repeating_comparison_aliases():
    """Tests that there are no repeating aliases for comparison, which would overwrite
    the earlier defined comparison."""
    COMPARISONS = []
    for c in all_comparisons:
        if c.aliases is not None:
            for a in c.aliases:
                COMPARISONS.append(a)
    duplicates = set([x for x in COMPARISONS if COMPARISONS.count(x) > 1])
    print(f"Duplicate comparisons: {duplicates}")
    assert len(duplicates) == 0

max_stars_repo_name: bradyrx/climpred

max_issues_repo_name: bradyrx/climpred

max_forks_repo_name: bradyrx/climpred


text:
import numpy as np
import mpnum as mp
from itertools import repeat

def get_number_ground_state(dims):
    """
        Embeds ground states: (1 0 ... 0).T (of possibly different dimensions) onto a chain as mparrays
    :param dims: List of Hilbert space dimensions of all ground states
    :return: MPArray which contains psi_0[0], psi_0[1], ... psi_0[L] embedded in a chain of length L
             where psi_0[i] is a ground state of Hilbert space dimension dims[i]. Returns None if dims is empty
    """
    psi = []
    if len(dims) == 0:
        return None
    for dim in dims:
        next_psi = np.zeros(dim)
        next_psi[0] = 1
        psi.append(next_psi)
    return mp.MPArray.from_kron(psi)

def broadcast_number_ground_state(dim, L):
    """
        Broadcasts a ground state: (1 0 ... 0).T (dimension dim) onto a chain of length L as mparrays
    :param dim: Hilbert space dimension of all ground states
    :param L: Size of the chain to broadcast state onto
    :return: MPArray with L copies of the ground state specified above, ranks 1 between them. Returns None if L is <= 0
    """
    if L <= 0:
        return None
    psi = np.zeros(dim)
    psi[0] = 1
    return mp.MPArray.from_kron(repeat(psi, L))

max_stars_repo_name: fhoeb/py-tmps

max_issues_repo_name: fhoeb/py-tmps

max_forks_repo_name: fhoeb/py-tmps


text:
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
import io
import base64
import sys
import argparse

# See https://matplotlib.org/3.1.0/users/dflt_style_changes.html
plt.style.use('seaborn-ticks')
mpl.rcParams['grid.color'] = 'grey'
mpl.rcParams['grid.linestyle'] = ':'
mpl.rcParams['grid.linewidth'] = 0.5
mpl.rcParams['axes.spines.right'] = False
mpl.rcParams['axes.spines.top'] = False
mpl.rcParams['font.size'] = 12
mpl.rcParams['legend.fontsize'] = 'medium'
mpl.rcParams['figure.titlesize'] = 'large'

def build_graph(f, export):
    if export:
        f.savefig(export, format='png')
        return
    img = io.BytesIO()
    f.set_size_inches(11.7, 8.27)
    f.savefig(img, format='png')
    img.seek(0)
    graph_url = base64.b64encode(img.getvalue()).decode()
    return graph_url
    # return 'data:image/png;base64,{}'.format(graph_url)

def compute_graph(df):
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(16, 12), sharex=False)

    ax1.fill_between(x=df['time'], y1=df['control_msgs'], y2=0, alpha=0.75, label="Control Plane")
    ax1.plot(df['time'], df['control_msgs'], marker='D', markersize=2, alpha=1, linewidth=1)

    ax1.fill_between(x=df['time'], y1=df['sr_read'], y2=0, alpha=0.75, label="SharedR READ")
    ax1.plot(df['time'], df['sr_read'], marker='^', markersize=2, alpha=1, linewidth=1)

    ax1.fill_between(x=df['time'], y1=df['sr_write'], y2=0, alpha=0.75, label="SharedR WRITE")
    ax1.plot(df['time'], df['sr_write'], marker='v', markersize=2, alpha=1, linewidth=1)

    ax1.fill_between(x=df['time'], y1=df['health_checks'], y2=0, alpha=0.75, label="Health Checks")
    ax1.plot(df['time'], df['health_checks'], marker='s', markersize=2, alpha=1, linewidth=1)

    ax1.fill_between(x=df['time'], y1=df['data_msgs'], y2=0, alpha=0.65, label="Data Plane")
    ax1.plot(df['time'], df['data_msgs'], marker='o', markersize=2, alpha=1, linewidth=1)

    ax1.legend()
    ax1.set_title('Messages Transmitted')
    ax1.set_ylabel('# Msgs')
    ax1.grid()

    ax2.plot(df['time'], df['control_msgs'].cumsum(), alpha=1, label="Control Plane")
    ax2.plot(df['time'], df['sr_read'].cumsum(), alpha=1, label="SharedR Read")
    ax2.plot(df['time'], df['sr_write'].cumsum(), alpha=1, label="SharedR Write")
    ax2.plot(df['time'], df['health_checks'].cumsum(), alpha=1, label="Health Checks")
    ax2.plot(df['time'], df['data_msgs'].cumsum(), alpha=1, label="Data Plane")
    ax2.legend()
    ax2.grid()
    ax2.set_title('CDF of messages transmitted')
    ax2.set_xlabel('Simulation Step')
    return fig

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument("--export", help="Export the plot", default=None)
    args = parser.parse_args()

    df = pd.read_csv(sys.stdin)

    fig = compute_graph(df)

    build_graph(fig, export=args.export)


text:
import sympy, types, math, inspect
import tensorflow as tf
import numpy as np
from collections import defaultdict

def split(desc):
    elements = [""]
    nested_parentheses = 0
    for c in desc:
        if nested_parentheses == 0 and c == "(":
            nested_parentheses += 1
            elements[-1] += c
        elif nested_parentheses > 0 and c == ")":
            elements[-1] += c
            nested_parentheses -= 1
        elif nested_parentheses == 0 and c == " ":
            elements = elements + [""]
        else:
            elements[-1] += c
        if nested_parentheses < 0:
            raise ValueError("Invalid parentheses")
    if nested_parentheses != 0:
        raise ValueError("Invalid parentheses")

    return [e.strip() for e in elements if len(e) > 0]

def inferred_value(n):
    if isinstance(n, int):
        return n
    elif "__dict__" in dir(n) and "_inferred_value" in vars(n) and not n._inferred_value is None and not n._inferred_value[0] is None:
        return int(n._inferred_value[0])
    else:
        return None

class Value:
    def __init__(self, name, value):
        self.name = name
        self.value = value

    def __str__(self):
        value = inferred_value(self.value)
        if not value is None:
            return str(value)
        else:
            return self.name

    def get_ordered_values(self):
        return [self]

    def get_value(self):
        value = inferred_value(self.value)
        return self.value if value is None else self.value

class Variable:
    def __init__(self, name, context, ellipses=None):
        if name.isnumeric():
            self.value = int(name)
            self.name = f"__constantdim{context.constant_name_counter}({self.value})"
            context.constant_name_counter += 1
        else:
            self.value = None
            self.name = name
        self.basename = self.name
        self.ellipses = ellipses

    def __str__(self):
        return self.name

    def get_all(self):
        return [self]

    def sympy_variable_num_occurrence(self, name):
        return 1 if name == self.name else 0

    def sympy_total_num_occurrence(self):
        return 1

    def sympy_element_num(self):
        return 1

    def expand(self, context, values, ellipses=[]):
        return [EllipsisVariable(self, ellipses, context) if len(ellipses) > 0 else self]

    def fill(self, values):
        return Value(self.name, values[self.name])

class EllipsisVariable:
    def __init__(self, variable, ellipses, context):
        self.basename = variable.name
        self.value = variable.value
        self.name = self.basename + "-" + "-".join([str(i) for _, _, i in ellipses])
        self.ellipses = ellipses

    def __str__(self):
        return self.name

    def get_all(self):
        return [self]

    def constant_shape(self):
        return np.asarray([d for _, d, _ in self.ellipses])

    def fill(self, values):
        return Value(self.name, values[self.name])

class Ellipsis:
    def __init__(self, sub, context):
        self.sub = sub
        self.sympy_expansion_count = sympy.Symbol(f"__ellipsis{context.ellipsis_name_counter}", integer=True)
        context.ellipsis_name_counter += 1

    def __str__(self):
        return str(self.sub) + "..."

    def get_all(self):
        return [self] + self.sub.get_all()

    def sympy_variable_num_occurrence(self, name):
        return self.sympy_expansion_count * self.sub.sympy_variable_num_occurrence(name)

    def sympy_total_num_occurrence(self):
        return self.sympy_expansion_count * self.sub.sympy_total_num_occurrence()

    def sympy_element_num(self):
        return self.sympy_expansion_count

    def expand(self, context, values, ellipses=[]):
        expansion_count = values[self.sympy_expansion_count.name]
        return [x for i in range(expansion_count) for x in self.sub.expand(context, values, ellipses + [(self, expansion_count, i)])]

class Group:
    def __init__(self, children):
        self.children = children

    def __str__(self):
        return "(" + " ".join([str(c) for c in self.children]) + ")"

    def get_all(self):
        return [self] + [x for c in self.children for x in c.get_all()]

    def sympy_variable_num_occurrence(self, name):
        return sum(c.sympy_variable_num_occurrence(name) for c in self.children)

    def sympy_total_num_occurrence(self):
        return sum(c.sympy_total_num_occurrence() for c in self.children)

    def sympy_element_num(self):
        return 1

    def expand(self, context, values, ellipses=[]):
        children = []
        for c in self.children:
            children.extend(c.expand(context, values, ellipses))
        return [Group(children)]

    def fill(self, values):
        return Group([c.fill(values) for c in self.children])

    def get_ordered_values(self):
        return [x for c in self.children for x in c.get_ordered_values()]

    def get_value(self):
        return math.prod([c.value for c in self.children])

class Root(Group):
    def __init__(self, children):
        Group.__init__(self, children)

    def __str__(self):
        return " ".join([str(c) for c in self.children])

    def sympy_root_children_num(self):
        return sum(c.sympy_element_num() for c in self.children)

    def expand(self, context, values):
        children = []
        for c in self.children:
            children.extend(c.expand(context, values))
        return Root(children)

    def fill(self, values):
        return Root([c.fill(values) for c in self.children])

def parse_desc(desc, context, ellipses=[], root=True):
    desc = split(desc)
    if len(desc) == 1:
        desc = desc[0]
        if desc.endswith("..."):
            ellipsis = Ellipsis.__new__(Ellipsis)
            sub = parse_desc(desc[:-3], context=context, ellipses=ellipses + [ellipsis], root=False)
            assert len(sub) == 1
            ellipsis.__init__(sub[0], context=context)
            result = [ellipsis]
        elif desc[0] == "(":
            assert desc[-1] == ")"
            result = [Group(parse_desc(desc[1:-1], context=context, ellipses=ellipses, root=False))]
        else:
            assert not "(" in desc and not ")" in desc
            result = [Variable(name=desc, context=context, ellipses=ellipses)]
    else:
        result = [x for d in desc for x in parse_desc(d, context=context, ellipses=ellipses, root=False)]

    if root:
        result = Root(result)
    return result

cache = {}
def apply(desc, *tensors, reduction=None, output_shape=None, output_ndims=None, keepdims=False, **constants):
    tensors = list(tensors)
    input_shapes = []
    for tensor in tensors:
        tf_input_shape = None
        input_shape = []
        for j in range(len(tensor.shape)):
            static_dim = tensor.shape[j]
            if not static_dim is None:
                input_shape.append(static_dim)
            else:
                if tf_input_shape is None:
                    tf_input_shape = tf.shape(tensor)
                input_shape.append(tf_input_shape[j])
        input_shapes.append(input_shape)

    output_shape_len = None if output_shape is None else (output_shape.shape[0] if "shape" in dir(output_shape) else len(output_shape))
    def get_constants(descs_in, desc_out):
        result = []

        name_to_variables = defaultdict(list)
        all_variables = []
        for v in [v for l in ([desc_in.get_all() for desc_in in descs_in] + [desc_out.get_all()]) for v in l]:
            if isinstance(v, EllipsisVariable) or isinstance(v, Variable):
                name_to_variables[v.basename].append(v)
                all_variables.append(v)

        # Input shape
        for i, desc_in in enumerate(descs_in):
            for j, c in enumerate(desc_in.children):
                result.append((f"__input{i}shape{j}", c, input_shapes[i][j]))

        # Output shape if given
        if not output_shape is None:
            for i, c in enumerate(desc_out.children):
                result.append((f"__outputshape{i}", c, output_shape[i]))

        # Constant arguments if given
        for name, value in constants.items():
            value_broadcasted = None
            for variable in name_to_variables[name]:
                if isinstance(variable, EllipsisVariable):
                    dest_shape = variable.constant_shape()
                    if value_broadcasted is None:
                        value_broadcasted = tf.broadcast_to(value, dest_shape)
                    index = [int(i) for i in variable.name.split("-")[-dest_shape.shape[0]:]]

                    result.append((f"__constantarg-{variable.name}", variable, value_broadcasted[index]))
                else:
                    result.append((f"__constantarg-{variable.name}", variable, value))

        # Constant dimensions if given
        for variable in all_variables:
            if not variable.value is None:
                result.append((f"__constantdim-{variable.name}", variable, variable.value))

        return result

    cache_key = (
        tuple(len(input_shape) for input_shape in input_shapes),
        desc,
        reduction,
        output_shape_len if not output_shape is None else None,
        output_ndims,
        keepdims,
        tuple((k, tuple(tf.convert_to_tensor(constants[k]).shape.as_list())) for k in sorted(list(constants.keys()))),
    )
    if not cache_key in cache:
        # ########################### Parse description ###########################
        desc = desc.split("->")
        if len(desc) != 2:
            raise ValueError("Description must contain exactly one '->'")
        descs_in, desc_out = desc
        descs_in = descs_in.split(",")
        if len(descs_in) != len(tensors):
            raise ValueError(f"Expected {len(descs_in)} input tensors, got {len(tensors)}")
        context = types.SimpleNamespace(
            ellipsis_name_counter=0,
            constant_name_counter=0,
        )
        descs_in = [parse_desc(desc_in, context=context) for desc_in in descs_in]
        desc_out = parse_desc(desc_out, context=context)

        # ########################### Expand ellipses ###########################
        descs = descs_in + [desc_out]
        variable_names = [set(v.name for v in desc.get_all() if isinstance(v, Variable)) for desc in descs]
        # variable_names_in = set(v.name for v in desc_in.get_all() if isinstance(v, Variable))
        # variable_names_out = set(v.name for v in desc_out.get_all() if isinstance(v, Variable))
        # sympy_variable_in_num_occurrences = {name: sympy.Symbol(f"__in-{name}", integer=True) for name in variable_names_in}
        # sympy_variable_out_num_occurrences = {name: sympy.Symbol(f"__out-{name}", integer=True) for name in variable_names_out}
        sympy_variable_num_occurrences = [{name: sympy.Symbol(f"__in-{name}", integer=True) for name in vn} for vn in variable_names]
        all_nodes = [x for desc in descs for x in desc.get_all()]

        equations = []

        # Inputs/ outputs from desc
        for desc, vn, svno in zip(descs, variable_names, sympy_variable_num_occurrences):
            for name in vn:
                equations.append(sympy.Eq(svno[name], desc.sympy_variable_num_occurrence(name)))
        for desc1, vn1, svno1 in zip(descs, variable_names, sympy_variable_num_occurrences):
            for desc2, vn2, svno2 in zip(descs, variable_names, sympy_variable_num_occurrences):
                if id(desc1) != id(desc2):
                    for name in vn1.intersection(vn2):
                        equations.append(sympy.Eq(svno1[name], svno2[name]))

        # Inputs/ outputs from tensors and parameters
        for i, desc_in in enumerate(descs_in):
            equations.append(sympy.Eq(desc_in.sympy_root_children_num(), len(input_shapes[i])))
        if not output_shape is None:
            equations.append(sympy.Eq(desc_out.sympy_root_children_num(), output_shape_len))
        if not output_ndims is None:
            equations.append(sympy.Eq(desc_out.sympy_root_children_num(), output_ndims))

        # Constants
        name_to_any_variable = {v.name: v for v in all_nodes if isinstance(v, Variable)}
        for name, value in constants.items():
            ellipses = name_to_any_variable[name].ellipses
            expected_rank = len(name_to_any_variable[name].ellipses)
            def get_shape(value):
                if tf.is_tensor(value) or isinstance(value, np.ndarray):
                    return [int(i) for i in value.shape]
                elif isinstance(value, list) or isinstance(value, tuple):
                    if len(set(tuple(get_shape(v)) for v in value)) != 1:
                        raise ValueError(f"Invalid value for parameter {name}")
                    return [len(value)] + get_shape(value[0])
                else:
                    return []
            got_shape = get_shape(value)
            if expected_rank == len(got_shape):
                for i in range(len(got_shape)):
                    equations.append(sympy.Eq(ellipses[i].sympy_expansion_count, got_shape[i]))

        ellipses = [e for e in all_nodes if isinstance(e, Ellipsis)]
        variables = [e.sympy_expansion_count for e in ellipses] + [v for svno in sympy_variable_num_occurrences for v in svno.values()]
        equations2 = []
        for e in equations:
            try_equations = equations2 + [e]
            result = sympy.solve(try_equations, variables)
            if len(result) > 0:
                equations2 = try_equations
        equations = equations2

        expansion_values = sympy.solve(equations, variables)
        for k, v in expansion_values.items():
            if not v.is_number:
                k = str(k)
                if k.startswith("__ellipsis"):
                    raise ValueError(f"Failed to deduce expansion count of ellipsis {k[10:]}")
                else:
                    raise ValueError(f"Failed to deduce variable occurrence count of {k}")
        expansion_values = {str(k): int(v) for k, v in expansion_values.items()}

        descs = [desc.expand(context, expansion_values) for desc in descs]
        descs_in = descs[:-1]
        desc_out = descs[-1]

        # ########################### Find dimension values ###########################
        variable_names = [set(v.name for v in desc.get_all() if isinstance(v, Variable) or isinstance(v, EllipsisVariable)) for desc in descs]
        variable_names_in = variable_names[:-1]
        variable_names_out = variable_names[-1]
        sympy_variable_values = {name: sympy.Symbol(name, integer=True) for name in set(v for vn in variable_names for v in vn)}
        def sympy_value(n):
            if isinstance(n, Root) or isinstance(n, Group):
                return math.prod([sympy_value(c) for c in n.children])
            elif isinstance(n, EllipsisVariable) or isinstance(n, Variable):
                return sympy_variable_values[n.name]
            else:
                assert False

        equations = []
        sympy_constants = {}
        constants_values = {}
        # remaining_assertions = []
        for sympy_constant_name, node, value in get_constants(descs_in, desc_out):
            sympy_computed_value = sympy_value(node)
            sympy_constant = sympy.Symbol(sympy_constant_name, integer=True)
            equations_try = equations + [sympy.Eq(sympy_constant, sympy_computed_value)]
            if len(sympy.solve(equations_try, list(sympy_variable_values.values()))) > 0:
                equations = equations_try
                sympy_constants[sympy_constant_name] = sympy_constant
                constants_values[sympy_constant_name] = value
            # else:
            #     remaining_assertions.append((sympy_constant, sympy_computed_value, value))

        # Solve symbolic
        sympy_solved_variable_values = sympy.solve(equations, list(sympy_variable_values.values()))

        sympy_constants_names = list(sympy_constants.keys())
        sympy_constants = [sympy_constants[k] for k in sympy_constants_names]
        def to_func(expr, sympy_values):
            func = sympy.lambdify(sympy_values, expr)
            source = inspect.getsource(func)
            source = source.replace("/", "//")
            d = {}
            exec(source, d, d)
            func = d["_lambdifygenerated"]
            return func
        variable_funcs = {str(k): to_func(v, sympy_constants) for k, v in sympy_solved_variable_values.items()}

        cache[cache_key] = (tuple(descs), sympy_constants_names, variable_funcs, variable_names)
    else:
        descs, sympy_constants_names, variable_funcs, variable_names = cache[cache_key]
        descs_in = descs[:-1]
        desc_out = descs[-1]
        variable_names_in = variable_names[:-1]
        variable_names_out = variable_names[-1]

        constants_values = {}
        for sympy_constant_name, node, value in get_constants(descs_in, desc_out):
            constants_values[sympy_constant_name] = value

    # Compute shapes in tensorflow from solved equations
    constants_values = [constants_values[k] for k in sympy_constants_names]
    variable_values = {k: func(*constants_values) for k, func in variable_funcs.items()}

    descs = [desc.fill(variable_values) for desc in descs]
    descs_in = descs[:-1]
    desc_out = descs[-1]

    # TODO: make assertions work in keras and tf.function
    # sympy_variable_names = list(sympy_variable_values.keys())
    # sympy_variables = [sympy_variable_values[k] for k in sympy_variable_names]
    # variables = [variable_values[k] for k in sympy_variable_names]
    # for sympy_constant, sympy_computed_value, value in remaining_assertions:
    #     got = evaluate(sympy_computed_value, sympy_variables, variables)
    #     tf.debugging.assert_equal(value, got, f"Assertion failed for {sympy_constant.name}. Expected {value}, got {got}")

    # ########################### Transform tensor ###########################
    ordered_names = [[v.name for v in desc.get_ordered_values()] for desc in descs]
    ordered_values = [[v.value for v in desc.get_ordered_values()] for desc in descs]
    in_ordered_names = ordered_names[:-1]
    in_ordered_values = ordered_values[:-1]
    out_ordered_names = ordered_names[-1]
    out_ordered_values = ordered_values[-1]

    # Reshape nested input to flat input
    for i in range(len(tensors)):
        shape = in_ordered_values[i]
        if len(shape) != len(tensors[i].shape):
            # TODO: possibly two consecutive reshapes here and in reduction
            tensors[i] = tf.reshape(tensors[i], shape)

    # Reduce and squeeze input dimensions
    for i in range(len(tensors)):
        reduced_names = [name for name in variable_names[i] if not any(name in variable_names[j] for j in range(len(descs)) if i != j)]
        reduced_len1_names = set(name for name in reduced_names if inferred_value(variable_values[name]) == 1)
        if len(reduced_names) - len(reduced_len1_names) > 0 and reduction is None:
            raise ValueError("Expected reduction argument")
        if len(reduced_names) == 0 and not reduction is None:
            raise ValueError("No dimensions are reduced, but got reduction argument")
        if len(reduced_names) > 0:
            shape = [v for v, name in zip(in_ordered_values[i], in_ordered_names[i]) if not name in reduced_len1_names]
            if len(shape) != len(tensors[i].shape):
                tensors[i] = tf.reshape(tensors[i], shape)

            axes = sorted([in_ordered_names[i].index(name) for name in reduced_names if not name in reduced_len1_names])
            if len(axes) > 0:
                tf_reduce_name = f"reduce_{reduction}"
                if tf_reduce_name in vars(tf):
                    tensors[i] = vars(tf)[tf_reduce_name](tensors[i], axis=axes)
                else:
                    raise ValueError(f"Invalid reduction argument {reduction}")
    # Transform to flat output
    if len(descs_in) == 1:
        # Use tf.transpose for single inputs
        in_ordered_names_union = [n for n in in_ordered_names[0] if n in out_ordered_names]
        out_ordered_names_union = [n for n in out_ordered_names if n in in_ordered_names[0]]
        assert len(in_ordered_names_union) == len(tensors[0].shape)
        assert len(out_ordered_names_union) == len(tensors[0].shape)

        perm = []
        for out_dim in range(len(out_ordered_names_union)):
            out_name = out_ordered_names_union[out_dim]
            in_index = in_ordered_names_union.index(out_name)
            perm.append(in_index)
        if perm != list(range(len(perm))):
            tensors[0] = tf.transpose(tensors[0], perm)

        output_tensor = tensors[0]
    else:
        # Use tf.einsum for multiple inputs
        einsum_str = ""
        einsum_variables = {}
        def get_einsum_variable(name):
            if name in einsum_variables:
                return einsum_variables[name]
            else:
                v = chr(ord("a") + len(einsum_variables))
                if ord(v) > ord("z"):
                    raise ValueError(f"Only supports up to {ord('z') - ord('a') + 1} input tensors")
                einsum_variables[name] = v
                return v
        for i in range(len(ordered_names)):
            ordered_names_union = [name for name in ordered_names[i] if any(name in ordered_names[j] for j in range(len(ordered_names)) if i != j)]
            for name in ordered_names_union:
                einsum_str += get_einsum_variable(name) + " "
            if i < len(ordered_names) - 2:
                einsum_str += ", "
            elif i == len(ordered_names) - 2:
                einsum_str += " -> "

        output_tensor = tf.einsum(einsum_str, *tensors, optimize="optimal")

    # Expand and broadcast missing output dimensions
    broadcasted_names = [n for n in variable_names_out if not any(n in vn for vn in variable_names_in)]
    if len(broadcasted_names) > 0:
        slices = [(tf.newaxis if name in broadcasted_names else slice(None)) for name in out_ordered_names]
        output_tensor = output_tensor[slices]
        if not all(inferred_value(value) == 1 for name, value in zip(out_ordered_names, out_ordered_values) if name in broadcasted_names):
            output_tensor = tf.broadcast_to(output_tensor, out_ordered_values)

    # Reshape flat output to nested output
    shape = [v.get_value() for v in desc_out.children]
    if len(shape) != len(out_ordered_values):
        output_tensor = tf.reshape(output_tensor, shape)

    return output_tensor

max_stars_repo_name: fferflo/tfcv

max_issues_repo_name: fferflo/tfcv

max_forks_repo_name: fferflo/tfcv

reason: import numpy,import sympy


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from seaborn.palettes import color_palette

'''
def plotArchCapability():
    plt.figure()
    plt.style.use('classic')
    df = pd.DataFrame({'DGEMM':cats, 'GFLOPs/s':data})
    ax = sns.barplot(x="DGEMM", y="GFLOPs/s", data=df)
    fig = ax.get_figure()
    pltName = "dgemm_performance_variations.png"
    path = results + pltName

    #attainableGFs = 4800 #P100=4800, V100=7065
    #ax.axhline(25,linewidth=3, color='red',label="Exascale Power (MW) Target",linestyle='--')
    yX = 8000 #P100=5500
    plt.ylim(0, yX)
    plt.xlabel('DGEMM', weight='bold')
    plt.ylabel('GFLOPs/s', weight='bold')
    plt.legend(loc='upper center',prop={'size': 9})
    plt.grid(True)
    fig.savefig(path)

def grouped_DEVbarplot(df, cat, subcat, val, path ):

    u = df[cat].unique()
    x = np.arange(len(u))
    subx = df[subcat].unique()
    offsets = (np.arange(len(subx))-np.arange(len(subx)).mean())/(len(subx)+1.)
    width= np.diff(offsets).mean()

    plt.figure()
    plt.style.use('classic')

    for i,gr in enumerate(subx):
        dfg = df[df[subcat] == gr]
        plt.bar(x+offsets[i], dfg[val].values, width=width,label="{}".format(gr))
    plt.xlabel(cat, weight='bold')
    plt.ylabel(val, weight='bold')
    #plt.ylim(0,  300)
    plt.xticks(x, u)

    plt.axhline(20,linewidth=3, color='red',label="Exascale Power (MW) Target",linestyle='--')

    plt.legend(loc='center left',prop={'size': 9})
    plt.grid(True)
        #plt.show()
    plt.savefig(path)

def grouped_MWbarplot(df, cat, subcat, val, path ):

    u = df[cat].unique()
    x = np.arange(len(u))
    subx = df[subcat].unique()
    offsets = (np.arange(len(subx))-np.arange(len(subx)).mean())/(len(subx)+1.)
    width= np.diff(offsets).mean()

    plt.figure()
    plt.style.use('classic')

    for i,gr in enumerate(subx):
        dfg = df[df[subcat] == gr]
        plt.bar(x+offsets[i], dfg[val].values, width=width,label="{}".format(gr))
    plt.xlabel(cat, weight='bold')
    plt.ylabel(val, weight='bold')
    #plt.ylim(0,  300)
    plt.xticks(x, u)

    plt.axhline(20,linewidth=3, color='red',label="Exascale Power (MW) Target",linestyle='--')

    plt.legend(loc='center left',prop={'size': 9})
    plt.grid(True)
        #plt.show()
    plt.savefig(path)

def plotPwrConsumption():
    plt.figure()
    plt.style.use('classic')
    df = pd.DataFrame({'DGEMM':cats, 'GFLOPs/s':data})
    ax = sns.barplot(x="DGEMM", y="GFLOPs/s", data=df)
    fig = ax.get_figure()
    pltName = "dgemm_performance_variations.png"
    path = results + pltName

    #attainableGFs = 4800 #P100=4800, V100=7065
    ax.axhline(attainableGFs,linewidth=3, color='black',label="Attainable GFLOPs/s",linestyle='--')
    yX = 8000 #P100=5500
    plt.ylim(0, yX)
    plt.xlabel('DGEMM', weight='bold')
    plt.ylabel('GFLOPs/s', weight='bold')
    plt.legend(loc='upper center',prop={'size': 9})
    plt.grid(True)
    fig.savefig(path)

pltPath = "C:/rf/lbnl/plots/examples/"
df = pd.DataFrame(columns=['Total Power Per Architecture (MW)', 'cat', 'Power (MW)'])
apps = ['GP100','GV100','GA100']
cats = ['Power (MW)','Power (MW)','Power (MW)']
data = [56.6,38.46,41.24]

for n in range(len(data)):
    df.loc[n] = [apps[n]] + [cats[n]] + [data[n]]

cat = "Total Power Per Architecture (MW)"
subcat = "cat"
val = "Power (MW)"
pltName = "exa-pwr.png"
results = pltPath + pltName
grouped_MWbarplot(df, cat, subcat, val, results)

# DEVICES
df1 = pd.DataFrame(columns=['No. of GPU Devices Required Per Architecture', 'cat', 'Devices'])
apps = ['GP100','GV100','GA100']
cats = ['GPU Devices','GPU Devices','GPU Devices']
data = [188680,128205,103093]

for n in range(len(data)):
    df1.loc[n] = [apps[n]] + [cats[n]] + [data[n]]

cat = "No. of GPU Devices Required Per Architecture"
subcat = "cat"
val = "GPU Devices"
pltName = "exa-dev.png"
results = pltPath + pltName
grouped_DEVbarplot(df1, cat, subcat, val, results)

def GPUCharacteristic(pltPath):
    fig, axs = plt.subplots(2, 3,sharex=True)

    #DVFS for STREAM on P100
    x = ['GP100','GV100','GA100']
    y0 = [56.6,38.46,41.24]
    axs[0, 0].plot(x, y0,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[0,0].set_ylabel("Exaflop System Power Consumption (MW)",weight='bold',fontsize=10)

    # Power capping for STREAM on P100
    y1 = [5.3, 7.8, 9.7]
    axs[0, 1].plot(x, y1,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[0,1].set_ylabel("Peak FP64 TFLOPS",weight='bold',fontsize=10)

    y2 = [188680,128205,103093]
    axs[0, 2].plot(x, y2,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[1,0].set_ylabel("Total GPU Devices",weight='bold',fontsize=10)

    #DVFS for STREAM on V100
    y0 = [300,300,400]
    axs[1, 0].plot(x, y0,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[1,0].set_ylabel("TDP (W)",weight='bold',fontsize=10)

    # Power capping for STREAM on V100
    y1 = [56, 80, 108]
    axs[1, 1].plot(x, y1,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[1,1].set_ylabel("Total SM Count",weight='bold',fontsize=10)

    y2 = [15.3, 21.1, 54.2]
    axs[1, 2].plot(x, y2,c='red',label="X",linestyle='dashed', marker='^', markersize=2, kind='bar', grid=True)
    axs[1,1].set_ylabel("Total Transistors Count (Billions)",weight='bold',fontsize=10)
    # axs[0,0].set_ylim([0, 65])
    # axs[0,1].set_ylim([0, 65])
    # axs[1,0].set_ylim([0, 65])
    # axs[1,1].set_ylim([0, 65])

    plt.suptitle("Power Consumption for Exaflop System Using GPUs",weight='bold',fontsize=12)
    #,bbox_inches='tight'
    fig.savefig(pltPath+'exaflop-data.png',bbox_inches='tight')

pltPath = "C:/rf/lbnl/plots/examples/"
GPUCharacteristic(pltPath)
'''
'''
def GPUCharacteristic(pltPath):

    plt.figure(figsize=(13,8))
    plt.subplot(2, 3,1)

    #DVFS for STREAM on P100
    x = ['GP100','GV100','GA100']
    y0 = [56.6,38.46,41.24]
    p1 =plt.bar(x, y0,color='red')
    plt.ylabel("Exaflop Power (MW)",weight='bold',fontsize=8)
    plt.axhline(20,linewidth=3, color='green',label="Exascale Power (MW)",linestyle='--')
    plt.legend(loc=0)
    plt.xticks([],[])

    # Power capping for STREAM on P100
    y1 = [5.3, 7.8, 9.7]
    plt.subplot(2, 3,2)
    p2=plt.bar(x, y1,label="X",color='orange')
    plt.ylabel("Peak FP64 TFLOPS",weight='bold',fontsize=8)
    plt.xticks([],[])

    y2 = [189,128,103]
    plt.subplot(2, 3,3)
    p3=plt.bar(x, y2,label="X",color='brown')
    plt.ylabel("GPU Devices (xK)",weight='bold',fontsize=8)
    plt.xticks([],[])

    #DVFS for STREAM on V100
    y0 = [300,300,400]
    plt.subplot(2, 3,4)
    p4=plt.bar(x, y0,label="X",color='royalblue')
    plt.ylabel("TDP (W)",weight='bold',fontsize=8)
    plt.xticks(x)

    # Power capping for STREAM on V100
    y1 = [56, 80, 108]
    plt.subplot(2, 3,5)
    p5=plt.bar(x, y1,color='olive')
    plt.ylabel("Total SM Count",weight='bold',fontsize=8)
    plt.xticks(x)

    y2 = [15.3, 21.1, 54.2]
    plt.subplot(2, 3,6)
    p6=plt.bar(x, y2,color='limegreen')
    plt.ylabel("Total Transistors Count (Billions)",weight='bold',fontsize=8)
    plt.xticks(x)
    # axs[0,0].set_ylim([0, 65])
    # axs[0,1].set_ylim([0, 65])
    # axs[1,0].set_ylim([0, 65])
    # axs[1,1].set_ylim([0, 65])
    #plt.grid(True)
    plt.suptitle("Power Consumption for Exaflop System Using GPUs",weight='bold',fontsize=12)
    #,bbox_inches='tight'
    plt.savefig(pltPath+'exaflop-data.png',bbox_inches='tight')
'''
'''
def GPUCharacteristic1(pltPath):
    plt.figure()
    plt.style.use('classic')
    #DVFS for STREAM on P100
    x = ['GP100','GV100','GA100']
    # y0 = [56.6,38.46,41.24]
    # p1 =plt.bar(x, y0,color='red')
    # plt.ylabel("Exaflop Power (MW)",weight='bold',fontsize=8)
    # plt.axhline(20,linewidth=3, color='green',label="Exascale Power (MW)",linestyle='--')
    # plt.legend(loc=0)
    # plt.xticks([],[])

    # Power capping for STREAM on P100
    #y1 = [5.3, 7.8, 9.7]
    y1 = [5300,7800,9700]
    #plt.subplot(2, 3,2)
    plt.plot(x, y1,label="Peak FP64 (GFLOPS)",color='orange')
    #plt.ylabel("Peak FP64 TFLOPS",weight='bold',fontsize=8)
    #plt.xticks([],[])

    #y2 = [189,128,103]
    #plt.subplot(2, 3,3)
    #plt.plot(x, y2,label="GPU Devices (xK)",color='brown')
    # plt.ylabel("GPU Devices (xK)",weight='bold',fontsize=8)
    # plt.xticks([],[])

    #DVFS for STREAM on V100
    #y0 = [300,300,400]
    # plt.subplot(2, 3,4)
    #p4=plt.plot(x, y0,label="TDP (W)",color='royalblue')
    # plt.ylabel("TDP (W)",weight='bold',fontsize=8)
    # plt.xticks(x)

    # # Power capping for STREAM on V100
    # y1 = [56, 80, 108]
    # plt.subplot(2, 3,5)
    # p5=plt.bar(x, y1,color='olive')
    # plt.ylabel("Total SM Count",weight='bold',fontsize=8)
    # plt.xticks(x)

    y2 = [15.3, 21.1, 54.2]
    #y2 = [15300000000000,21100000000000,54200000000000]
    #plt.subplot(2, 3,6)
    plt.plot(x, y2,label = 'Transistors Count (Billions)',color='limegreen')
    #plt.ylabel("Transistors Count (Billions)",weight='bold',fontsize=8)
    plt.xticks(x)
    plt.yscale("log")

    # axs[0,0].set_ylim([0, 65])
    # axs[0,1].set_ylim([0, 65])
    # axs[1,0].set_ylim([0, 65])
    # axs[1,1].set_ylim([0, 65])
    #plt.grid(True)
    #plt.suptitle("Power Consumption for Exaflop System Using GPUs",weight='bold',fontsize=12)
    #,bbox_inches='tight'
    plt.legend(loc=0)
    plt.grid(True)
    plt.savefig(pltPath+'exaflop-data.png')
'''

def GPUCharacteristic2(pltPath):
    plt.figure()
    plt.style.use('classic')
    #DVFS for STREAM on P100
    x = ['GP100','GV100','GA100']
    y0 = [56.6,38.46,41.24]
    # p1 =plt.bar(x, y0,color='red')
    plt.plot(x, y0,color='red',label='Exaflop Power (MW)')
    # plt.ylabel("Exaflop Power (MW)",weight='bold',fontsize=8)
    plt.axhline(20,linewidth=3, color='green',label="Exascale Power Limit (MW)",linestyle='--')
    # plt.legend(loc=0)
    # plt.xticks([],[])

    # Power capping for STREAM on P100
    #y1 = [5.3, 7.8, 9.7]
    #y1 = [5300,7800,9700]
    #plt.subplot(2, 3,2)
    #plt.plot(x, y1,label="Peak FP64 (GFLOPS)",color='orange')
    #plt.ylabel("Peak FP64 TFLOPS",weight='bold',fontsize=8)
    #plt.xticks([],[])

    #y2 = [189,128,103]
    #plt.subplot(2, 3,3)
    #plt.plot(x, y2,label="GPU Devices (xK)",color='brown')
    # plt.ylabel("GPU Devices (xK)",weight='bold',fontsize=8)
    # plt.xticks([],[])

    #DVFS for STREAM on V100
    y0 = [300,300,400]
    # plt.subplot(2, 3,4)
    plt.plot(x, y0,label="TDP (W)",color='royalblue')
    # plt.ylabel("TDP (W)",weight='bold',fontsize=8)
    plt.xticks(x)

    # # Power capping for STREAM on V100
    # y1 = [56, 80, 108]
    # plt.subplot(2, 3,5)
    # p5=plt.bar(x, y1,color='olive')
    # plt.ylabel("Total SM Count",weight='bold',fontsize=8)
    # plt.xticks(x)

    #y2 = [15.3, 21.1, 54.2]
    #y2 = [15300000000000,21100000000000,54200000000000]
    #plt.subplot(2, 3,6)
    #plt.plot(x, y2,label = 'Transistors Count (Billions)',color='limegreen')
    #plt.ylabel("Transistors Count (Billions)",weight='bold',fontsize=8)
    #plt.xticks(x)
    plt.yscale("log")

    # axs[0,0].set_ylim([0, 65])
    # axs[0,1].set_ylim([0, 65])
    # axs[1,0].set_ylim([0, 65])
    # axs[1,1].set_ylim([0, 65])
    #plt.grid(True)
    #plt.suptitle("Power Consumption for Exaflop System Using GPUs",weight='bold',fontsize=12)
    #,bbox_inches='tight'
    plt.legend(loc=0)
    plt.grid(True)
    plt.savefig(pltPath+'exaflop-pwr-usage.png')

pltPath = "C:/rf/lbnl/plots/examples/"
GPUCharacteristic2(pltPath)