text:
# import important libraries
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import seaborn as sns
import os
import argparse

# import machine learning libraries
from sklearn.linear_model import Ridge, LogisticRegression
from sklearn.svm import SVR, SVC
from sklearn.metrics import f1_score
from sklearn.model_selection import cross_val_predict, GridSearchCV, KFold
from sklearn.metrics import confusion_matrix, classification_report,\
                                              precision_recall_fscore_support

# import custom functions for vectorizing & visualizing data
import utils
from process_data import get_split
plt.rcParams.update({'font.size': 12})

all_props = ['bulk_modulus',
             'thermal_conductivity',
             'shear_modulus',
             'band_gap',
             'debye_temperature',
             'thermal_expansion']

symbols = ['B', '$\\kappa$', 'G', 'E$_{\\mathrm{g}}$', 'D', '$\\alpha$']

prop_labels = ['Bulk Modulus (GPa)',
               'Log$_{10}$ Thermal Conductivity $\\left(\\dfrac{\\mathrm{W}}' +
               '{\\mathrm{m}\\cdot \\mathrm{K}}\\right)$',
               'Shear Modulus (GPa)',
               'Band Gap (eV)',
               'Log$_{10}$ Debye Temperature (K)',
               'Log$_{10}$ Thermal Expansion $(\\mathrm{K}^{-1})$']

arg2prop = {'bulk_modulus': 'ael_bulk_modulus_vrh',
            'thermal_conductivity': 'agl_log10_thermal_conductivity_300K',
            'shear_modulus': 'ael_shear_modulus_vrh',
            'band_gap': 'Egap',
            'debye_temperature': 'ael_log10_debye_temperature',
            'thermal_expansion': 'agl_log10_thermal_expansion_300K'}

prop2label = dict([[v, k] for k, v
                  in zip(prop_labels, arg2prop.values())])

parser_desc = 'Reproduce the results of this work'
parser = argparse.ArgumentParser(description=parser_desc)
group = parser.add_mutually_exclusive_group(required=True)

group.add_argument('--properties',
                   type=str,
                   nargs='+',
                   metavar='Property to reproduce',
                   choices=all_props,
                   help=('example:\n\t' +
                         'python reproduce.py --properties bulk_modulus\n\t'))

group.add_argument('--all',
                   action='store_true',
                   help='Run through each property one at a time '
                   'and generate results and figures.')

args = parser.parse_args()

if not args.all:
    mat_props = []
    for j in args.properties:
        mat_props.append(arg2prop[j])
else:
    mat_props = list(map(lambda p: arg2prop[p], all_props))
print('Reproducing results for the following data:', mat_props)

def optimize_threshold(y_train_labeled, y_train_pred):
    """Given a DataFrame of labels and predictions, return the
     optimal threshold for a high F1 score"""
    y_train_ = y_train_labeled.copy()
    y_train_pred_ = pd.Series(y_train_pred).copy()
    f1score_max = 0
    for threshold in np.arange(0.1, 1, 0.1):
        diff = (max(y_train_pred) - min(y_train_pred))
        threshold = min(y_train_pred) + threshold * diff
        y_train_pred_[y_train_pred < threshold] = 0
        y_train_pred_[y_train_pred >= threshold] = 1
        f1score = f1_score(y_train_, y_train_pred_)
        if f1score > f1score_max:
            f1score_max = f1score
            opt_thresh = threshold
    return opt_thresh

def get_performance(mat_props, seed):
    metrics_dict = {}
    for mat_prop in mat_props:
        os.makedirs('figures/'+mat_prop, exist_ok=True)
        data = get_split(mat_prop, elem_prop='oliynyk', seed_num=seed)
        X_train_scaled, X_test_scaled = data[0:2]
        y_train, y_test = data[2:4]
        y_train_labeled, y_test_labeled = data[4:6]
        formula_train, formula_test = data[6:8]

        test_threshold = y_test.iloc[-y_test_labeled.sum().astype(int)]
        train_threshold = y_train.iloc[-y_train_labeled.sum().astype(int)]

        y = pd.concat([y_train, y_test])

        plt.figure(1, figsize=(7, 7))
        ax = sns.distplot(y, bins=50, kde=False)

        rect1 = patches.Rectangle((test_threshold, 0),
                                  ax.get_xlim()[1]-test_threshold,
                                  ax.get_ylim()[1], linewidth=1,
                                  edgecolor='k',
                                  facecolor='g',
                                  alpha=0.2)
        ax.add_patch(rect1)

        text_size = 18

        ax.text(.1,
                .5,
                'Ordinary\nCompounds',
                size=text_size,
                horizontalalignment='left',
                verticalalignment='center',
                transform=ax.transAxes)
        ax.text(.98,
                .15,
                'Extraordinary\nCompounds',
                size=text_size,
                horizontalalignment='right',
                verticalalignment='center',
                transform=ax.transAxes)
        ax.tick_params(direction='in',
                       length=5,
                       bottom=True,
                       top=True,
                       left=True,
                       right=True,
                       labelsize=text_size)
        ax.set_xlabel(prop2label[mat_prop], size=text_size)
        ax.set_ylabel('number of occurrences'.title(), size=text_size)
        plt.savefig('figures/' + mat_prop + '/distplot',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()
        # ## Learn with a Ridge Regression (linear model)
        # define ridge regression object
        rr = Ridge()
        # define k-folds
        cv = KFold(n_splits=5, shuffle=True, random_state=1)
        # choose search space
        parameter_candidates = {'alpha': np.logspace(-5, 2, 10)}

        # define the grid search
        grid = GridSearchCV(estimator=rr,
                            param_grid=parameter_candidates,
                            cv=cv)
        # run grid search
        grid.fit(X_train_scaled, y_train)

        # plot grid search to ensure good values)
        utils.plot_1d_grid_search(grid, midpoint=0.75)
        print('best parameters:', grid.best_params_)
        plt.savefig('figures/' + mat_prop + '/rr_1d_search',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()
        best_params_rr = grid.best_params_

    #    best_params_rr = {'alpha': 0.0021544346900318843}
        rr = Ridge(**best_params_rr)
        rr.fit(X_train_scaled, y_train)
        y_test_predicted_rr = rr.predict(X_test_scaled)
        y_train_predicted_rr = rr.predict(X_train_scaled)
        # plot the data
        plt.figure(figsize=(6, 6))
        plt.plot(y_test,
                 y_test_predicted_rr,
                 marker='o',
                 mfc='none',
                 color='#0077be',
                 linestyle='none',
                 label='test')
        plt.plot(y_train,
                 y_train_predicted_rr,
                 marker='o',
                 mfc='none',
                 color='#e34234',
                 linestyle='none',
                 label='train')
        max_val = max(y_test.max(), y_test_predicted_rr.max())
        min_val = min(y_test.min(), y_test_predicted_rr.min())
        plt.plot([min_val, max_val], [min_val, max_val], 'k--')
        limits = [min_val, max_val]
        plt.xlim(limits)
        plt.ylim(limits)
        plt.xlabel('actual')
        plt.ylabel('predicted')
        plt.legend(loc=4)
        plt.tick_params(direction='in',
                        length=5,
                        bottom=True,
                        top=True,
                        left=True,
                        right=True)
        plt.savefig('figures/' + mat_prop + '/rr_act_vs_pred',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        # ## Learn with a support vector regression (non-linear model)

        # to speed up the grid search, optimize on a subsample of data
        X_train_scaled_sampled = X_train_scaled.sample(500, random_state=1)
        y_train_sampled = y_train.loc[X_train_scaled_sampled.index.values]

        # define support vector regression object (default to rbf kernel)
        svr = SVR()
        # define k-folds
        cv = KFold(n_splits=5, shuffle=True, random_state=1)
        # choose search space
        parameter_candidates = {'C': np.logspace(2, 4, 8),
                                'gamma': np.logspace(-3, 1, 8)}

        # define the grid search
        grid = GridSearchCV(estimator=svr,
                            param_grid=parameter_candidates,
                            cv=cv)
        # run grid search
        grid.fit(X_train_scaled_sampled, y_train_sampled)

        # plot grid search to ensure good values
        utils.plot_2d_grid_search(grid, midpoint=0.7)
        plt.savefig('figures/' + mat_prop + '/svr_2d_search',
                    dpi=300, bbox_inches='tight')
        plt.clf()
        print('best parameters:', grid.best_params_)
        best_params_svr = grid.best_params_

        svr = SVR(**best_params_svr)
        svr.fit(X_train_scaled, y_train)

        y_test_predicted_svr = svr.predict(X_test_scaled)
        y_train_predicted_svr = svr.predict(X_train_scaled)

        # plot the data
        plt.figure(figsize=(6, 6))
        plt.plot(y_test,
                 y_test_predicted_svr,
                 marker='o',
                 mfc='none',
                 color='#0077be',
                 linestyle='none',
                 label='test')
        plt.plot(y_train,
                 y_train_predicted_svr,
                 marker='o',
                 mfc='none',
                 color='#e34234',
                 linestyle='none',
                 label='train')

        max_val = max(y_test.max(), y_test_predicted_svr.max())
        min_val = min(y_test.min(), y_test_predicted_svr.min())
        plt.plot([min_val, max_val], [min_val, max_val], 'k--')
        limits = [min_val, max_val]
        plt.xlim(limits)
        plt.ylim(limits)

        plt.xlabel('actual')
        plt.ylabel('predicted')
        plt.legend(loc=4)
        plt.tick_params(direction='in',
                        length=5,
                        bottom=True,
                        top=True,
                        left=True,
                        right=True)
        plt.savefig('figures/' + mat_prop + '/svr_act_vs_pred',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        # # Approach the problem as a classification task
        # ## Learn with a logistic regression (linear classification)
        # define logistic regression object
        lr = LogisticRegression(solver='lbfgs')
        # define k-folds
        cv = KFold(n_splits=5, shuffle=True, random_state=1)

        # choose search space
        class_1_weight = [{0: 1, 1: weight} for weight in
                          np.linspace(1, 50, 5)]
        parameter_candidates = {'C': np.logspace(-1, 4, 5),
                                'class_weight': class_1_weight}

        # define the grid search. We use log-loss to decide which
        # parameters to use.
        grid = GridSearchCV(estimator=lr,
                            param_grid=parameter_candidates,
                            scoring='neg_log_loss',
                            cv=cv)

        # run grid search
        grid.fit(X_train_scaled, y_train_labeled)

        # plot grid search to ensure good values
        utils.plot_2d_grid_search(grid, midpoint=-0.05, vmin=-0.13, vmax=0)
        plt.savefig('figures/' + mat_prop + '/lr_2d_search',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()
        print('best parameters:', grid.best_params_)
        best_params_lr = grid.best_params_

        lr = LogisticRegression(solver='lbfgs', penalty='l2', **best_params_lr)
        lr.fit(X_train_scaled, y_train_labeled)

        # define k-folds
        cv = KFold(n_splits=5, shuffle=True, random_state=1)

        y_pred_train_lr = cross_val_predict(lr,
                                            X_train_scaled,
                                            y_train_labeled,
                                            cv=cv)
        y_prob_train_lr = cross_val_predict(lr,
                                            X_train_scaled,
                                            y_train_labeled,
                                            cv=cv,
                                            method='predict_proba')
        y_probability_train_lr = [probability[1] for probability in
                                  y_prob_train_lr]

        y_prob_test_lr = lr.predict_proba(X_test_scaled)
        y_probability_test_lr = [probability[1] for probability in
                                 y_prob_test_lr]

        df_cm = pd.DataFrame(confusion_matrix(y_train_labeled,
                                              y_pred_train_lr))

        ax = sns.heatmap(df_cm,
                         square=True,
                         annot=True,
                         annot_kws={"size": 18},
                         cbar=False,
                         linewidths=.5,
                         cmap="YlGnBu",
                         center=-10000000)

        ax.set_ylabel('actual')
        ax.set_xlabel('predicted')
        ax.xaxis.tick_top()
        plt.savefig('figures/' + mat_prop + '/lr_cm',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        threshold = 0.5
        utils.plot_prob(threshold,
                        y_train,
                        y_probability_train_lr,
                        threshold_x=train_threshold,
                        mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/lr_train_prob_thresh={:0.2f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        # ### Check our perfromance on the test set!

        utils.plot_prob(threshold,
                        y_test,
                        y_probability_test_lr,
                        threshold_x=test_threshold,
                        mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/lr_test_prob_thresh={:0.2f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        # ### Compare this performance to regression models
        #
        # **For the same recall, we are three times more likely that predicted
        #   compound is not actually extraordinary.**

        threshold = optimize_threshold(y_train_labeled, y_train_predicted_rr)
        utils.plot_regression(threshold,
                              y_train,
                              y_train_predicted_rr,
                              threshold_x=train_threshold,
                              mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/rr_train_reg_thresh={:0.2f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        utils.plot_regression(threshold,
                              y_test,
                              y_test_predicted_rr,
                              threshold_x=test_threshold,
                              mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/rr_test_reg_thresh={:0.2f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        threshold = optimize_threshold(y_train_labeled, y_train_predicted_svr)
        utils.plot_regression(threshold,
                              y_train,
                              y_train_predicted_svr,
                              threshold_x=train_threshold,
                              mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/svr_train_reg_thresh={:0.02f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        utils.plot_regression(threshold,
                              y_test,
                              y_test_predicted_svr,
                              threshold_x=test_threshold,
                              mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/svr_test_reg_thresh={:0.02f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        # ## Learn with a support vector classification (non-linear)
        # to speed up the grid search, optimize on a subsample of data
        index_location = X_train_scaled_sampled.index.values
        y_train_labeled_sampled = y_train_labeled.loc[index_location]

        # define suppor vector classification object
        # (need to set probability to True)
        svc = SVC(probability=True)
        # define k-folds
        cv = KFold(n_splits=5, shuffle=True, random_state=1)

        # choose search space (we will start with class_weight=1
        # as that was optimal for svc)
        parameter_candidates = {'C': np.logspace(-1, 4, 5),
                                'gamma': np.logspace(-2, 2, 5)}

        # define the grid search. We use log-loss to decide
        # which parameters to use.
        grid = GridSearchCV(estimator=svc,
                            param_grid=parameter_candidates,
                            scoring='neg_log_loss',
                            cv=cv)

        # run grid search
        grid.fit(X_train_scaled_sampled, y_train_labeled_sampled)

        # plot grid search to ensure good values
        utils.plot_2d_grid_search(grid, midpoint=-0.04, vmin=-0.13, vmax=0)
        plt.savefig('figures/' + mat_prop +
                    '/svc_2d_search.png',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()
        print('best parameters:', grid.best_params_)
        best_params_svc = grid.best_params_

        svc = SVC(probability=True, **best_params_svc)
        svc.fit(X_train_scaled, y_train_labeled)

        cv = KFold(n_splits=5, shuffle=True, random_state=1)

        y_pred_train_svc = cross_val_predict(svc,
                                             X_train_scaled,
                                             y_train_labeled,
                                             cv=cv)
        y_prob_train_svc = cross_val_predict(svc,
                                             X_train_scaled,
                                             y_train_labeled,
                                             cv=cv,
                                             method='predict_proba')
        y_probability_train_svc = [probability[1] for probability in
                                   y_prob_train_svc]

        y_prob_test_svc = svc.predict_proba(X_test_scaled)
        y_probability_test_svc = [probability[1] for probability in
                                  y_prob_test_svc]

        metrics = precision_recall_fscore_support(y_train_labeled,
                                                  y_pred_train_svc)
        precision, recall, fscore, support = metrics
        print('precision: {:0.2f}\nrecall: {:0.2f}'.format(precision[1],
              recall[1]))
        df_cm = pd.DataFrame(confusion_matrix(y_train_labeled,
                                              y_pred_train_svc))

        ax = sns.heatmap(df_cm,
                         square=True,
                         annot=True,
                         annot_kws={"size": 18},
                         cbar=False,
                         linewidths=0.5,
                         cmap="YlGnBu",
                         center=-10000000)
        ax.set_ylabel('actual')
        ax.set_xlabel('predicted')
        ax.xaxis.tick_top()
        plt.savefig('figures/' + mat_prop +
                    '/svc_cm',
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        threshold = 0.5
        utils.plot_prob(threshold,
                        y_train,
                        y_probability_train_svc,
                        threshold_x=train_threshold,
                        mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/svc_train_prob_thresh={:0.02f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        utils.plot_prob(threshold,
                        y_test,
                        y_probability_test_svc,
                        threshold_x=test_threshold,
                        mat_prop=prop2label[mat_prop])
        plt.savefig('figures/' + mat_prop +
                    '/svc_test_prob_thresh={:0.2f}.png'.format(threshold),
                    dpi=300,
                    bbox_inches='tight')
        plt.clf()

        metrics_dict[mat_prop] = {'precision': [], 'recall': [], 'f1': []}
        threshold = optimize_threshold(y_train_labeled, y_train_predicted_rr)
        y_pred_rr = [1 if x >= threshold else 0 for x in y_test_predicted_rr]
        threshold = optimize_threshold(y_train_labeled, y_train_predicted_svr)
        y_pred_svr = [1 if x >= threshold else 0 for x in y_test_predicted_svr]
        threshold = 0.5
        y_pred_lr = [1 if x >= threshold else 0 for x in y_probability_test_lr]
        threshold = 0.5
        y_pred_svc = [1 if x >= threshold else 0 for x in
                      y_probability_test_svc]

        predictions = [y_pred_rr,
                       y_pred_svr,
                       y_pred_lr,
                       y_pred_svc]

        for prediction in predictions:
            print(classification_report(y_test_labeled, prediction))
            metrics = precision_recall_fscore_support(y_test_labeled,
                                                      prediction)
            precision, recall, f1, support = metrics
            if precision[1] == 0:
                if precision == 10:
                    pass
            metrics_dict[mat_prop]['precision'].append(precision[1])
            metrics_dict[mat_prop]['recall'].append(recall[1])
            metrics_dict[mat_prop]['f1'].append(f1[1])
    return metrics_dict

def build_metrics():
    for seed in [1]:
        metrics = get_performance(mat_props, seed)
        for prop in metrics:
            metric_csv = prop+'_metrics_seed_{:0.0f}.csv'.format(seed)
            computed_metrics = os.listdir('data/metrics/')
            if metric_csv in computed_metrics:
                continue
            else:
                df_prop_metric = pd.DataFrame(metrics[prop],
                                              index=['rr', 'svr', 'lr', 'svc'])
                df_prop_metric.to_csv('data/metrics/'+metric_csv)

def plot_metrics():
    metric_mean = {}
    metric_std = {}
    metric_mean[0] = {}
    metric_mean[1] = {}
    metric_mean[2] = {}
    metric_std[0] = {}
    metric_std[1] = {}
    metric_std[2] = {}
    for prop in mat_props:
        rr = []
        svr = []
        lr = []
        svc = []
        for seed in [1, 2, 3, 4, 5]:
            metric_csv = prop+'_metrics_seed_{:0.0f}.csv'.format(seed)
            df_prop_metric = pd.read_csv('data/metrics/'+metric_csv)
            rr.append(df_prop_metric.iloc[0, 1:].tolist())
            svr.append(df_prop_metric.iloc[1, 1:].tolist())
            lr.append(df_prop_metric.iloc[2, 1:].tolist())
            svc.append(df_prop_metric.iloc[3, 1:].tolist())
        for i in [0, 1, 2]:
            metric_mean[i][prop] = [pd.DataFrame(rr).mean()[i],
                                    pd.DataFrame(svr).mean()[i],
                                    pd.DataFrame(lr).mean()[i],
                                    pd.DataFrame(svc).mean()[i]]
            metric_std[i][prop] = [pd.DataFrame(rr).std()[i],
                                   pd.DataFrame(svr).std()[i],
                                   pd.DataFrame(lr).std()[i],
                                   pd.DataFrame(svc).std()[i]]

    df_p_mean = pd.DataFrame(metric_mean[0], index=['rr', 'svr', 'lr', 'svc'])
    df_p_std = pd.DataFrame(metric_std[0], index=['rr', 'svr', 'lr', 'svc'])
    df_r_mean = pd.DataFrame(metric_mean[1], index=['rr', 'svr', 'lr', 'svc'])
    df_r_std = pd.DataFrame(metric_std[1], index=['rr', 'svr', 'lr', 'svc'])
    df_f_mean = pd.DataFrame(metric_mean[2], index=['rr', 'svr', 'lr', 'svc'])
    df_f_std = pd.DataFrame(metric_std[2], index=['rr', 'svr', 'lr', 'svc'])

    plt.rcParams.update({'font.size': 12})
    means = [df_p_mean, df_r_mean, df_f_mean]
    stds = [df_p_std, df_r_std, df_f_std]
    metric_type = ['Precision', 'Recall', 'F1-Score']
    i = 0
    f, (ax1, ax2, ax3) = plt.subplots(3, 1, sharey=False, figsize=(7, 9))
    f.subplots_adjust(hspace=.05, wspace=1)
    axes = [ax1, ax2, ax3]
    colors = ['#d7191c', '#fdae61', '#abdda4', '#2b83ba']

    prop_loc = [1, 2, 3, 4, 5, 6]
    for df_mean, df_std, ax in zip(means, stds, axes):
        alpha = 0.15

        ax.fill_between(prop_loc,
                        df_mean.loc['rr']+df_std.loc['rr'],
                        df_mean.loc['rr']-df_std.loc['rr'],
                        color=colors[0],
                        alpha=alpha)
        ax.fill_between(prop_loc,
                        df_mean.loc['svr']+df_std.loc['svr'],
                        df_mean.loc['svr']-df_std.loc['svr'],
                        color=colors[3],
                        alpha=alpha)
        ax.fill_between(prop_loc,
                        df_mean.loc['lr']+df_std.loc['lr'],
                        df_mean.loc['lr']-df_std.loc['lr'],
                        color=colors[2],
                        alpha=alpha)
        ax.fill_between(prop_loc,
                        df_mean.loc['svc']+df_std.loc['svc'],
                        df_mean.loc['svc']-df_std.loc['svc'],
                        color=colors[1],
                        alpha=alpha)

        ax.plot(prop_loc, df_mean.loc['rr'], '-x', color=colors[0],
                linewidth=2, label='Ridge')
        ax.plot(prop_loc, df_mean.loc['svr'], '-s', color=colors[3],
                linewidth=2, label='SVR')
        ax.plot(prop_loc, df_mean.loc['lr'], '--d', color=colors[2],
                linewidth=3, label='Logistic')
        ax.plot(prop_loc, df_mean.loc['svc'], '--*', color=colors[1],
                linewidth=3, label='SVC')
        ax.set_ylabel(metric_type[i])
        ax.tick_params(top=True, right=True, direction='in', length=6)

        if i == 0:
            ax.set_ylim(0.25, 1)
            ax.xaxis.set_ticklabels([])
        if i == 1:
            ax.set_ylim(0.25, 1)
            ax.xaxis.set_ticklabels([])
            ax.legend(loc=4, fancybox=True, framealpha=0.3)
        if i == 2:
            ax.set_ylim(0.25, 1)
            ax.xaxis.set_ticklabels([])
        i += 1

    plt.xticks(np.array(prop_loc), labels=[prop for prop in symbols],
               rotation=0, ha='center', rotation_mode='anchor')

    plt.savefig('figures/metrics.png', dpi=300, bbox_inches='tight')

if __name__ == '__main__':
    build_metrics()


text:
# -*- coding: utf-8 -*-
# Copyright 2020 Minh Nguyen (@dathudeptrai)
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""Extract durations based-on tacotron-2 alignments for FastSpeech."""

import argparse
import logging
import os
from numba import jit
import sys

sys.path.append(".")

import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import yaml
from tqdm import tqdm

from examples.tacotron2.tacotron_dataset import CharactorMelDataset
from tensorflow_tts.configs import Tacotron2Config
from tensorflow_tts.models import TFTacotron2

@jit(nopython=True)
def get_duration_from_alignment(alignment):
    D = np.array([0 for _ in range(np.shape(alignment)[0])])

    for i in range(np.shape(alignment)[1]):
        max_index = list(alignment[:, i]).index(alignment[:, i].max())
        D[max_index] = D[max_index] + 1

    return D

def main():
    """Running extract tacotron-2 durations."""
    parser = argparse.ArgumentParser(
        description="Extract durations from charactor with trained Tacotron-2 "
        "(See detail in tensorflow_tts/example/tacotron-2/extract_duration.py)."
    )
    parser.add_argument(
        "--rootdir",
        default=None,
        type=str,
        required=True,
        help="directory including ids/durations files.",
    )
    parser.add_argument(
        "--outdir", type=str, required=True, help="directory to save generated speech."
    )
    parser.add_argument(
        "--checkpoint", type=str, required=True, help="checkpoint file to be loaded."
    )
    parser.add_argument(
        "--use-norm", default=1, type=int, help="usr norm-mels for train or raw."
    )
    parser.add_argument("--batch-size", default=8, type=int, help="batch size.")
    parser.add_argument("--win-front", default=2, type=int, help="win-front.")
    parser.add_argument("--win-back", default=2, type=int, help="win-front.")
    parser.add_argument(
        "--use-window-mask", default=1, type=int, help="toggle window masking."
    )
    parser.add_argument("--save-alignment", default=0, type=int, help="save-alignment.")
    parser.add_argument(
        "--config",
        default=None,
        type=str,
        required=True,
        help="yaml format configuration file. if not explicitly provided, "
        "it will be searched in the checkpoint directory. (default=None)",
    )
    parser.add_argument(
        "--verbose",
        type=int,
        default=1,
        help="logging level. higher is more logging. (default=1)",
    )
    args = parser.parse_args()

    # set logger
    if args.verbose > 1:
        logging.basicConfig(
            level=logging.DEBUG,
            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
        )
    elif args.verbose > 0:
        logging.basicConfig(
            level=logging.INFO,
            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
        )
    else:
        logging.basicConfig(
            level=logging.WARN,
            format="%(asctime)s (%(module)s:%(lineno)d) %(levelname)s: %(message)s",
        )
        logging.warning("Skip DEBUG/INFO messages")

    # check directory existence
    if not os.path.exists(args.outdir):
        os.makedirs(args.outdir)

    # load config
    with open(args.config) as f:
        config = yaml.load(f, Loader=yaml.Loader)
    config.update(vars(args))

    if config["format"] == "npy":
        char_query = "*-ids.npy"
        mel_query = "*-raw-feats.npy" if args.use_norm is False else "*-norm-feats.npy"
        char_load_fn = np.load
        mel_load_fn = np.load
    else:
        raise ValueError("Only npy is supported.")

    # define data-loader
    dataset = CharactorMelDataset(
        dataset=config["tacotron2_params"]["dataset"],
        root_dir=args.rootdir,
        charactor_query=char_query,
        mel_query=mel_query,
        charactor_load_fn=char_load_fn,
        mel_load_fn=mel_load_fn,
        reduction_factor=config["tacotron2_params"]["reduction_factor"],
        use_fixed_shapes=True,
    )
    dataset = dataset.create(allow_cache=True, batch_size=args.batch_size, drop_remainder=False)

    # define model and load checkpoint
    tacotron2 = TFTacotron2(
        config=Tacotron2Config(**config["tacotron2_params"]),
        name="tacotron2",
    )
    tacotron2._build()  # build model to be able load_weights.
    tacotron2.load_weights(args.checkpoint)

    # apply tf.function for tacotron2.
    tacotron2 = tf.function(tacotron2, experimental_relax_shapes=True)

    for data in tqdm(dataset, desc="[Extract Duration]"):
        utt_ids = data["utt_ids"]
        input_lengths = data["input_lengths"]
        mel_lengths = data["mel_lengths"]
        utt_ids = utt_ids.numpy()
        real_mel_lengths = data["real_mel_lengths"]
        del data["real_mel_lengths"]

        # tacotron2 inference.
        mel_outputs, post_mel_outputs, stop_outputs, alignment_historys = tacotron2(
            **data,
            use_window_mask=args.use_window_mask,
            win_front=args.win_front,
            win_back=args.win_back,
            training=True,
        )

        # convert to numpy
        alignment_historys = alignment_historys.numpy()

        for i, alignment in enumerate(alignment_historys):
            real_char_length = input_lengths[i].numpy()
            real_mel_length = real_mel_lengths[i].numpy()
            alignment_mel_length = int(
                np.ceil(
                    real_mel_length / config["tacotron2_params"]["reduction_factor"]
                )
            )
            alignment = alignment[:real_char_length, :alignment_mel_length]
            d = get_duration_from_alignment(alignment)  # [max_char_len]

            d = d * config["tacotron2_params"]["reduction_factor"]
            assert (
                np.sum(d) >= real_mel_length
            ), f"{d}, {np.sum(d)}, {alignment_mel_length}, {real_mel_length}"
            if np.sum(d) > real_mel_length:
                rest = np.sum(d) - real_mel_length
                # print(d, np.sum(d), real_mel_length)
                if d[-1] > rest:
                    d[-1] -= rest
                elif d[0] > rest:
                    d[0] -= rest
                else:
                    d[-1] -= rest // 2
                    d[0] -= rest - rest // 2

                assert d[-1] >= 0 and d[0] >= 0, f"{d}, {np.sum(d)}, {real_mel_length}"

            saved_name = utt_ids[i].decode("utf-8")

            # check a length compatible
            assert (
                len(d) == real_char_length
            ), f"different between len_char and len_durations, {len(d)} and {real_char_length}"

            assert (
                np.sum(d) == real_mel_length
            ), f"different between sum_durations and len_mel, {np.sum(d)} and {real_mel_length}"

            # save D to folder.
            np.save(
                os.path.join(args.outdir, f"{saved_name}-durations.npy"),
                d.astype(np.int32),
                allow_pickle=False,
            )

            # save alignment to debug.
            if args.save_alignment == 1:
                figname = os.path.join(args.outdir, f"{saved_name}_alignment.png")
                fig = plt.figure(figsize=(8, 6))
                ax = fig.add_subplot(111)
                ax.set_title(f"Alignment of {saved_name}")
                im = ax.imshow(
                    alignment, aspect="auto", origin="lower", interpolation="none"
                )
                fig.colorbar(im, ax=ax)
                xlabel = "Decoder timestep"
                plt.xlabel(xlabel)
                plt.ylabel("Encoder timestep")
                plt.tight_layout()
                plt.savefig(figname)
                plt.close()

if __name__ == "__main__":
    main()

ext: py

max_stars_repo_name: kewlbear/TensorFlowTTS

max_issues_repo_name: kewlbear/TensorFlowTTS

max_forks_repo_name: kewlbear/TensorFlowTTS

reason: import numpy,from numba


text:
from ddd_subplots import subplots
import numpy as np
import os
from sklearn import datasets
from sklearn.decomposition import PCA

def test_subplots():
    X, y = datasets.load_iris(return_X_y=True)
    np.random.shuffle(X)
    X_reduced = PCA(n_components=3).fit_transform(X)
    colors = np.array(["red", "green", "blue"])[y]
    fig, axes = subplots(1, 3, figsize=(15, 5))
    for axis in axes.flatten():
        axis.scatter(*X_reduced.T, depthshade=False, c=colors, marker='o', s=20)
    fig.tight_layout()
    fig.savefig("test.png")

reason: import numpy


text:
import numpy as np
from psecas import Solver, FourierGrid, IO
from psecas.systems.kh_uniform import KelvinHelmholtzUniform
import time
from mpi4py import MPI

comm = MPI.COMM_WORLD

"""
    This example illustrates how to use mpi4py to run a parameter study.

    The physics explored is that of the KH instability with anisotropic
    viscosity and a constant magnetic field in the x-direction.
    The equilibrium is assumed to have constant density, temperature and
    pressure. The velocity profile varies smoothly and the setup is periodic.

    More details about this specific setup can be found in

    Berlok, T. & Pfrommer, C. (2019). *On the Kelvin-Helmholtz instability
    with smooth initial conditions – Linear theory and simulations*, MNRAS,
    485, 908

    Another reference for the KHI with anisotropic viscosity is

    Suzuki, K., Ogawa, T., Matsumoto, Y., & Matsumoto, R. (2013).
    Magnetohydrodynamic simulations of the formation of cold fronts in
    clusters of galaxies: Effects of anisotropic viscosity. Astrophysical
    Journal, 768(2). https://doi.org/10.1088/0004-637X/768/2/175
"""

directory = './data/'
kx_global = np.linspace(3, 4, 5)
kx_local = kx_global[comm.rank :: comm.size]

grid = FourierGrid(N=64, zmin=0, zmax=2)
system = KelvinHelmholtzUniform(grid, beta=1e4, nu=1e-2, kx=0)
io = IO(system, directory, __file__, len(kx_global))

solver = Solver(grid, system)

for i in range(len(kx_local)):
    t1 = time.time()
    system.kx = kx_local[i]
    omega, v = solver.solve()
    io.save_system(i)
    io.log(i, time.time() - t1, 'kx = {:1.4e}'.format(system.kx))

io.finished()

max_stars_repo_name: tberlok/psecas

max_issues_repo_name: tberlok/psecas

max_forks_repo_name: tberlok/psecas


text:
# coding=utf-8
# Copyright (c) DIRECT Contributors
import abc
import ast
import functools
import importlib
import logging
import os
import pathlib
import random
import subprocess
import sys
from collections import OrderedDict
from typing import Any, Callable, Dict, KeysView, List, Optional, Tuple, Union

import numpy as np
import torch

logger = logging.getLogger(__name__)

def is_complex_data(data: torch.Tensor, complex_last: bool = True) -> bool:
    """Returns True if data is a complex tensor, i.e. has a complex axis of dimension 2, and False otherwise.

    Parameters
    ----------
    data: torch.Tensor
        For 2D data the shape is assumed ([batch], [coil], height, width, [complex])
        or ([batch], [coil], [complex], height, width).
        For 3D data the shape is assumed ([batch], [coil], slice, height, width, [complex])
        or ([batch], [coil], [complex], slice, height, width).
    complex_last: bool
        If true, will require complex axis to be at the last axis. Default: True.

    Returns
    -------
    bool
    """
    if 2 not in data.shape:
        return False
    if complex_last:
        if data.size(-1) != 2:
            return False
    else:
        if data.ndim == 6:
            if data.size(2) != 2 and data.size(-1) != 2:  # (B, C, 2, S, H, 2) or (B, C, S, H, W, 2)
                return False

        elif data.ndim == 5:
            # (B, 2, S, H, W) or (B, C, 2, H, W) or (B, S, H, W, 2) or (B, C, H, W, 2)
            if data.size(1) != 2 and data.size(2) != 2 and data.size(-1) != 2:
                return False

        elif data.ndim == 4:
            if data.size(1) != 2 and data.size(-1) != 2:  # (B, 2, H, W) or (B, H, W, 2) or (S, H, W, 2)
                return False

        elif data.ndim == 3:
            if data.size(-1) != 2:  # (H, W, 2)
                return False

        else:
            raise ValueError(f"Not compatible number of dimensions for complex data. Got {data.ndim}.")

    return True

def is_power_of_two(number: int) -> bool:
    """Check if input is a power of 2.

    Parameters
    ----------
    number: int

    Returns
    -------
    bool
    """
    return number != 0 and ((number & (number - 1)) == 0)

def ensure_list(data: Any) -> List:
    """Ensure input is a list.

    Parameters
    ----------
    data: object

    Returns
    -------
    list
    """
    if data is None:
        return []

    if not isinstance(data, (list, tuple)):
        return [data]

    return list(data)

def cast_as_path(data: Optional[Union[pathlib.Path, str]]) -> Optional[pathlib.Path]:
    """Ensure the the input is a path.

    Parameters
    ----------
    data: str or pathlib.Path

    Returns
    -------
    pathlib.Path
    """
    if data is None:
        return None

    return pathlib.Path(data)

def str_to_class(module_name: str, function_name: str) -> Callable:
    """Convert a string to a class Base on: https://stackoverflow.com/a/1176180/576363.

    Also support function arguments, e.g. ifft(dim=2) will be parsed as a partial and return ifft where dim has been
    set to 2.

    Examples
    --------
    >>> def mult(f, mul=2):
    >>>    return f*mul

    >>> str_to_class(".", "mult(mul=4)")
    >>> str_to_class(".", "mult(mul=4)")
    will return a function which multiplies the input times 4, while

    >>> str_to_class(".", "mult")
    just returns the function itself.

    Parameters
    ----------
    module_name: str
        e.g. direct.data.transforms
    function_name: str
        e.g. Identity
    Returns
    -------
    object
    """
    tree = ast.parse(function_name)
    func_call = tree.body[0].value  # type: ignore
    args = [ast.literal_eval(arg) for arg in func_call.args] if hasattr(func_call, "args") else []
    kwargs = (
        {arg.arg: ast.literal_eval(arg.value) for arg in func_call.keywords} if hasattr(func_call, "keywords") else {}
    )

    # Load the module, will raise ModuleNotFoundError if module cannot be loaded.
    module = importlib.import_module(module_name)

    if not args and not kwargs:
        return getattr(module, function_name)
    return functools.partial(getattr(module, func_call.func.id), *args, **kwargs)

def dict_to_device(
    data: Dict[str, torch.Tensor],
    device: Union[torch.device, str, None],
    keys: Union[List, Tuple, KeysView, None] = None,
) -> Dict:
    """Copy tensor-valued dictionary to device. Only torch.Tensor is copied.

    Parameters
    ----------
    data: Dict[str, torch.Tensor]
    device: torch.device, str
    keys: List, Tuple
        Subselection of keys to copy.

    Returns
    -------
    Dictionary at the new device.
    """
    if keys is None:
        keys = data.keys()
    return {k: v.to(device) if isinstance(v, torch.Tensor) else v for k, v in data.items() if k in keys}

def detach_dict(data: Dict[str, torch.Tensor], keys: Optional[Union[List, Tuple, KeysView]] = None) -> Dict:
    """Return a detached copy of a dictionary. Only torch.Tensor's are detached.

    Parameters
    ----------
    data: Dict[str, torch.Tensor]
    keys: List, Tuple
        Subselection of keys to detach

    Returns
    -------
    Dictionary at the new device.
    """
    if keys is None:
        keys = data.keys()
    return {k: v.detach() for k, v in data.items() if k in keys if isinstance(v, torch.Tensor)}

def reduce_list_of_dicts(data: List[Dict[str, torch.Tensor]], mode="average", divisor=None) -> Dict[str, torch.Tensor]:
    """Average a list of dictionary mapping keys to Tensors.

    Parameters
    ----------
    data: List[Dict[str, torch.Tensor]])
    mode: str
        Which reduction mode, average reduces the dictionary, sum just adds while average computes the average.
    divisor: None or int
        If given values are divided by this factor.

    Returns
    -------
    Dict[str, torch.Tensor]: Reduced dictionary.
    """
    if not data:
        return {}

    if mode not in ["average", "sum"]:
        raise ValueError("Reduction can only be `sum` or `average`.")

    if not divisor:
        divisor = 1.0

    result_dict = {k: torch.zeros_like(v) for k, v in data[0].items()}

    for elem in data:
        result_dict = {k: result_dict[k] + v for k, v in elem.items()}

    if mode == "average":
        divisor *= len(data)

    return {k: v / divisor for k, v in result_dict.items()}

def merge_list_of_dicts(list_of_dicts: List[Dict]) -> Dict:
    """A list of dictionaries is merged into one dictionary.

    Parameters
    ----------
    list_of_dicts: List[Dict]

    Returns
    -------
    Dict
    """
    if not list_of_dicts:
        return {}

    return functools.reduce(lambda a, b: {**dict(a), **dict(b)}, list_of_dicts)

def evaluate_dict(
    fns_dict: Dict[str, Callable], source: torch.Tensor, target: torch.Tensor, reduction: str = "mean"
) -> Dict:
    """Evaluate a dictionary of functions.

    Examples
    --------
    > evaluate_dict({'l1_loss: F.l1_loss, 'l2_loss': F.l2_loss}, a, b)

    Will return
    > {'l1_loss', F.l1_loss(a, b, reduction='mean'), 'l2_loss': F.l2_loss(a, b, reduction='mean')

    Parameters
    ----------
    fns_dict: Dict[str, Callable]
    source: torch.Tensor
    target: torch.Tensor
    reduction: str

    Returns
    -------
    Dict[str, torch.Tensor]
        Evaluated dictionary.
    """
    return {k: fns_dict[k](source, target, reduction=reduction) for k, v in fns_dict.items()}

def prefix_dict_keys(data: Dict[str, Any], prefix: str) -> Dict[str, Any]:
    """Append a prefix to a dictionary keys.

    Parameters
    ----------
    data: Dict[str, Any]
    prefix: str

    Returns
    -------
    Dict[str, Any]
    """
    return {prefix + k: v for k, v in data.items()}

def git_hash() -> str:
    """Returns the current git hash.

    Returns
    -------
    _git_hash: str
        The current git hash.
    """
    try:
        _git_hash = subprocess.check_output(["git", "rev-parse", "HEAD"], stderr=subprocess.PIPE).decode().strip()
    except FileNotFoundError:
        _git_hash = "git not installed."
    except subprocess.CalledProcessError as e:
        exit_code = e.returncode
        stdout = e.output.decode(sys.getfilesystemencoding())
        stderr = e.stderr.decode(sys.getfilesystemencoding())
        _git_hash = f"cannot get git hash: git returned {exit_code}\n" f"stdout: {stdout}.\n" f"stderr: {stderr}."

    return _git_hash

def normalize_image(image: torch.Tensor, eps: float = 0.00001) -> torch.Tensor:
    r"""Normalize image to range [0,1] for visualization.

    Given image :math:`x` and :math:`\epsilon`, it returns:

    .. math::
        \frac{x - \min{x}}{\max{x} + \epsilon}.

    Parameters
    ----------
    image: torch.Tensor
        Image to scale.
    eps: float

    Returns
    -------
    image: torch.Tensor
        Scaled image.
    """

    image = image - image.min()
    image = image / (image.max() + eps)
    return image

def multiply_function(multiplier: float, func: Callable) -> Callable:
    """Create a function which multiplier another one with a multiplier.

    Parameters
    ----------
    multiplier: float
        Number to multiply with.
    func: callable
        Function to multiply.

    Returns
    -------
    return_func: Callable
    """

    def return_func(*args, **kwargs):
        return multiplier * func(*args, **kwargs)

    return return_func

class DirectTransform:
    """Direct transform class.

    Defines :meth:`__repr__` method for Direct transforms.
    """

    def __init__(self):
        """Inits DirectTransform."""
        super().__init__()

    def __repr__(self):
        """Representation of DirectTransform."""
        repr_string = self.__class__.__name__ + "("
        for k, v in self.__dict__.items():
            if k == "logger":
                continue
            repr_string += f"{k}="
            if callable(v):
                if hasattr(v, "__class__"):
                    repr_string += type(v).__name__ + ", "
                else:
                    # TODO(jt): better way to log functions
                    repr_string += str(v) + ", "
            elif isinstance(v, (dict, OrderedDict)):
                repr_string += f"{k}=dict(len={len(v)}), "
            elif isinstance(v, list):
                repr_string = f"{k}=list(len={len(v)}), "
            elif isinstance(v, tuple):
                repr_string = f"{k}=tuple(len={len(v)}), "
            else:
                repr_string += str(v) + ", "

        if repr_string[-2:] == ", ":
            repr_string = repr_string[:-2]
        return repr_string + ")"

class DirectModule(torch.nn.Module, DirectTransform, abc.ABC):
    @abc.abstractmethod
    def __init__(self):
        super().__init__()

    def forward(self, sample: Dict):
        pass  # This comment passes "Function/method with an empty body PTC-W0049" error.

def count_parameters(models: Dict) -> None:
    """Count the number of parameters of a dictionary of models.

    Parameters
    ----------
    models: Dict
        Dictionary mapping model name to model.
    """
    total_number_of_parameters = 0
    for model_name in models:
        n_params = sum(p.numel() for p in models[model_name].parameters())
        logger.info(f"Number of parameters model {model_name}: {n_params} ({n_params / 10.0 ** 3:.2f}k).")
        logger.debug(models[model_name])
        total_number_of_parameters += n_params
    logger.info(
        f"Total number of parameters model: {total_number_of_parameters} "
        f"({total_number_of_parameters / 10.0 ** 3:.2f}k)."
    )

def _select_random_seed(min_seed_value: int = 1, max_seed_value: int = 2**32) -> int:
    """Selects random seed.

    Parameters
    ----------
    min_seed_value: int
        Minimum seed value. Default: 1.
    max_seed_value: int
        Maximum seed value. Default: 2**32.

    Returns
    -------
    seed: int
        Random integer in range(min_seed_value, max_seed_value).
    """
    return random.randint(min_seed_value, max_seed_value)  # nosec

def set_all_seeds(seed: int) -> None:
    """Sets seed for deterministic runs.

    Parameters
    ----------
    seed:  int
        Seed for random module.

    Returns
    -------
    """
    # Global seed.
    random.seed(seed)

    # Set individual seeds
    torch.manual_seed(_select_random_seed())
    torch.cuda.manual_seed(_select_random_seed())
    np.random.seed(_select_random_seed())
    os.environ["PYTHONHASHSEED"] = str(_select_random_seed())
    os.environ["PL_GLOBAL_SEED"] = str(_select_random_seed())

def chunks(list_to_chunk: List, number_of_chunks: int):
    """Yield `number_of_chunks number` of sequential chunks from `list_to_chunk`. Adapted from [1]_.

    Parameters
    ----------
    list_to_chunk: List
    number_of_chunks: int

    References
    ----------

    .. [1] https://stackoverflow.com/a/54802737
    """
    d, r = divmod(len(list_to_chunk), number_of_chunks)
    for idx in range(number_of_chunks):
        si = (d + 1) * (idx if idx < r else r) + d * (0 if idx < r else idx - r)
        yield list_to_chunk[si : si + (d + 1 if idx < r else d)]

def remove_keys(input_dict: Dict, keys: Union[str, List[str], Tuple[str]]) -> Dict:
    """Removes `keys` from `input_dict`.

    Parameters
    ----------
    input_dict: Dict
    keys: Union[str, List[str], Tuple[str]]

    Returns
    -------
    Dict
    """
    input_dict = dict(input_dict).copy()
    if not isinstance(keys, (list, tuple)):
        keys = [keys]
    for key in keys:
        if key not in input_dict:
            continue
        del input_dict[key]
    return input_dict

max_stars_repo_name: cnmy-ro/direct-custom

max_issues_repo_name: cnmy-ro/direct-custom

max_forks_repo_name: cnmy-ro/direct-custom


text:
from pathlib import Path
from io import BytesIO
import numpy as np
from PIL import Image
from cairosvg import svg2png
import cv2
import argparse, os
from transforms import transforms

parser = argparse.ArgumentParser()
parser.add_argument('--path', help= 'paste path to biog.txt file')
args = parser.parse_args()

svgdir = Path(args.path) # to change directory to argument passed for '--path'
outputdir = Path(os.getcwd() + '/output')
# make sure output dir exists
if outputdir.exists() == False:
     outputdir.mkdir()

for svg in svgdir.iterdir():
    if svg.is_dir():
        pass
    else:
        png = svg2png(bytestring=svg.read_bytes())
        pil_img = Image.open(BytesIO(png)).convert('RGBA')
        #pil_img.save(outputdir.stem + '/' + svg.name.split('.')[0] +'-pil.png')
        cv_img = cv2.cvtColor(np.array(pil_img), cv2.COLOR_RGBA2BGRA)
        transformer = transforms()
        manipulated_img = transformer.filterHologram(cv_img)
        cv2.imwrite(outputdir.stem + '/' + svg.name.split('.')[0] + '.png', manipulated_img)

max_stars_repo_name: knkp/svg-transform-tool

max_issues_repo_name: knkp/svg-transform-tool

max_forks_repo_name: knkp/svg-transform-tool


text:
"""
kkpy.util
========================

Utility functions for my research

.. currentmodule:: util

Winds
-------
.. autosummary::
    kkpy.util.wind2uv
    kkpy.util.uv2wind
    kkpy.util.ms2knot
    kkpy.util.knot2ms

Radars
--------
.. autosummary::
    kkpy.util.dbzmean

Microphysics
---------------
.. autosummary::
    kkpy.util.vel_atlas
    kkpy.util.calc_dsdmoments

Maps
-------
.. autosummary::
    kkpy.util.proj_dfs
    kkpy.util.proj_icepop
    kkpy.util.dist_bearing

Spatial calculations
----------------------
.. autosummary::
    kkpy.util.nn_idx_2d
    kkpy.util.cross_section_2d
    kkpy.util.to_lower_resolution

ICE-POP 2018
--------------
.. autosummary::
    kkpy.util.icepop_events
    kkpy.util.icepop_sites
    kkpy.util.icepop_extent

Miscellaneous
---------------
.. autosummary::
    kkpy.util.std2d
    kkpy.util.nanstd2d
    kkpy.util.nanconvolve2d
    kkpy.util.derivative
    kkpy.util.stats
    kkpy.util.summary
"""
import numpy as np

def wind2uv(wd=None, ws=None, radians=False, knots=False):
    """
    Convert wind direction and speed to u and v wind components.

    Examples
    ---------
    >>> u, v = kkpy.util.wind2uv(wd, ws)

    Parameters
    ----------
    wd : array_like
        Array containing wind direction in **degree**. It should be **meteorological** direction, not mathematical.
    ws : array_like
        Array containing wind speed in **m/s**.
    radians : boolean, optional
        If this is set to True, the unit of *wd* is **radian**. The default is False (i.e. **degree**).
    knots : boolean, optional
        If this is set to True, the unit of *ws* is **knots**. The default is False (i.e. **m/s**).

    Returns
    ---------
    u : ndarray
        Return u component of wind in **m/s**.
    v : ndarray
        Return v component of wind in **m/s**.
    """
    if not radians:
        wd = np.radians(wd)

    if knots:
        ws = knot2ms(ws)

    u = -ws * np.sin(wd)
    v = -ws * np.cos(wd)

    return u, v

def uv2wind(u, v):
    """
    Convert u and v wind components to wind direction and speed.

    Examples
    ---------
    >>> wd, ws = kkpy.util.uv2wind(u, v)

    Parameters
    ----------
    u : array_like
        Array containing u component of wind in **m/s**.
    v : array_like
        Array containing v component of wind in **m/s**.

    Returns
    ---------
    wd : ndarray
        Return wind direction in **degree**.
    ws : ndarray
        Return wind speed in **m/s**.
    """

    ws = np.hypot(u, v)
    wd = 270 - np.rad2deg(np.arctan2(v, u))
    wd = wd % 360

    return wd, ws

def ms2knot(ws_ms):
    """
    Convert unit of wind speed from meter per seconds to knots.

    Examples
    ---------
    >>> ws_knot = kkpy.util.ms2knot(ws_ms)

    Parameters
    ----------
    ws_ms : array_like
        Array containing wind speed in **m/s**.

    Returns
    ---------
    ws_knot : ndarray
        Return wind speed in **knots**.
    """

    ws_knot = ws_ms * 1. / 0.5144444

    return ws_knot

def knot2ms(ws_knot):
    """
    Convert unit of wind speed from knots to meter per seconds.

    Examples
    ---------
    >>> ws_ms = kkpy.util.knot2ms(ws_knot)

    Parameters
    ----------
    ws_knot : array_like
        Array containing wind speed in **knots**.

    Returns
    ---------
    ws_ms : ndarray
        Return wind speed in **m/s**.
    """

    ws_ms = ws_knot * 0.5144444

    return ws_ms

def nn_idx_2d(xtarget, ytarget, x2d, y2d):
    """
    Get nearest indices of the x and y values of target if 2D array is given.

    Examples
    ---------
    >>> x2d, y2d = np.meshgrid(np.arange(128,130,0.1), np.arange(36,38,0.1))
    >>> idx2d = kkpy.util.nn_idx_2d(128.35, 37.65, x2d, y2d)
    >>> print(x2d[idx2d], y2d[idx2d])

    Parameters
    ----------
    xtarget : float
        Target value of x.
    ytarget : float
        Target value of y.
    x2d : 2D array
        Numpy 2D array containing x value of each grid point. The shape of **x2d** and **y2d** should be same.
    y2d : 2D array
        Numpy 2D array containing y value of each grid point. The shape of **x2d** and **y2d** should be same.

    Returns
    ---------
    idx2d : 1D array
        Return an array of nearest x and y indicies.
    """
    import bottleneck as bn
    diff = (x2d - xtarget)**2 + (y2d - ytarget)**2
    idx2d = np.float_(np.unravel_index(bn.nanargmin(diff), x2d.shape))
    return np.int_(idx2d)

def cross_section_2d(xy0, xy1, x2d, y2d, value2d, linewidth=1, along='x', reduce_func=np.mean):
    """
    Get values along the transect of two points (lon/lat) for 2D array.

    Examples
    ---------
    >>> dem, lon, lat, proj = kkpy.io.read_dem(area='pyeongchang')
    >>> start = [127.52, 37.3]
    >>> end = [129, 36.52]
    >>> xaxis, vcross = kkpy.util.cross_section_2d(start, end, lon, lat, dem)
    >>> plt.plot(xaxis, vcross)
    >>> plt.show()

    Parameters
    ----------
    xy0 : array_like
        Array of x and y values of starting point.
    xy1 : array_like
        Array of x and y values of ending point.
    x2d : 2D array
        Numpy 2D array containing x value of each grid point. The shape of **x2d**, **y2d**, and **value2d** should be same.
    y2d : 2D array
        Numpy 2D array containing y value of each grid point. The shape of **x2d**, **y2d**, and **value2d** should be same.
    value2d : 2D array
        Numpy 2D array containing data value of each grid point. The shape of **x2d**, **y2d**, and **value2d** should be same.
    linewidth : int, optional
        The linewidth used in average over perpendicular direction along the transect. Default is 1 (i.e. no average).
    along : str, optional
        Set 'x' to return **xaxis** in x axis, otherwise return it in y axis. Default is 'x'.
    reduce_func : callable, optional
        Function used to calculate the aggregation of pixel values perpendicular to the profile_line direction when linewidth > 1. See skimage.measure.profile_line for detail.

    Returns
    ---------
    xaxis : 1D array
        Return xaxis of cross-section in longitude or latitude unit. The unit is determined by along keyword.
    vcross : 1D array
        Return averaged value along the cross-section.
    """
    import skimage.measure
    minx, maxx = np.nanmin(x2d), np.nanmax(x2d)
    miny, maxy = np.nanmin(y2d), np.nanmax(y2d)

    if not minx < xy0[0] < maxx:
        raise ValueError(f'xy0[0]({xy0[0]}) is out of range ({minx} ~ {maxx})')
    if not miny < xy0[1] < maxy:
        raise ValueError(f'xy0[1]({xy0[1]}) is out of range ({miny} ~ {maxy})')
    if not minx < xy1[0] < maxx:
        raise ValueError(f'xy1[0]({xy1[0]}) is out of range ({minx} ~ {maxx})')
    if not miny < xy1[1] < maxy:
        raise ValueError(f'xy1[1]({xy1[1]}) is out of range ({miny} ~ {maxy})')
    assert along in ['x', 'y'], "Check along keyword (should be 'x' or 'y')"

    # find indices
    istart = nn_idx_2d(xy0[0], xy0[1], x2d, y2d)
    iend = nn_idx_2d(xy1[0], xy1[1], x2d, y2d)

    # get xaxis
    ixlen, iylen = iend - istart
    ixsize = int(np.ceil(np.hypot(ixlen, iylen) + 1))
    if along in 'x':
        xaxis = np.linspace(xy0[0], xy1[0], ixsize)
    else:
        xaxis = np.linspace(xy0[1], xy1[1], ixsize)

    # get cross-section
    vcross = skimage.measure.profile_line(
        value2d,
        istart, iend,
        linewidth=linewidth,
        reduce_func=reduce_func,
        order=0,
        mode='constant',
        cval=np.nan
    )

    return xaxis, vcross

def proj_dfs():
    """
    Return a lambert conformal conic projection of DFS (Digital Forecasting System) used in KMA.

    Examples
    ---------
    >>> import cartopy.crs as ccrs
    >>> ax = plt.subplot(proj=kkpy.util.proj_dfs())
    >>> ax.scatter([126], [38], transform=ccrs.PlateCarree())
    >>> plt.show()

    Returns
    ---------
    proj : cartopy projection
        Return a map projection of DFS.
    """
    import cartopy.crs as ccrs

    globe = ccrs.Globe(ellipse=None,
                       semimajor_axis=6371008.77,
                       semiminor_axis=6371008.77)

    proj = ccrs.LambertConformal(central_longitude=126,
                                 central_latitude=38,
                                 standard_parallels=(30,60),
                                 false_easting=400000,
                                 false_northing=789000,
                                 globe=globe)
    return proj

def proj_icepop():
    """
    Return a lambert conformal conic projection for ICE-POP 2018 domain.

    Examples
    ---------
    >>> import cartopy.crs as ccrs
    >>> ax = plt.subplot(proj=kkpy.util.proj_icepop())
    >>> ax.scatter([129], [37.7], transform=ccrs.PlateCarree())
    >>> ax.set_extent(kkpy.util.icepop_extent(), crs=ccrs.PlateCarree())
    >>> plt.show()

    Returns
    ---------
    proj : cartopy projection
        Return a map projection of DFS.
    """
    import cartopy.crs as ccrs

    globe = ccrs.Globe(ellipse=None,
                       semimajor_axis=6371008.77,
                       semiminor_axis=6371008.77)

    proj = ccrs.LambertConformal(central_longitude=128.75,
                                 central_latitude=37.65,
                                 standard_parallels=(30,60),
                                 false_easting=75000,
                                 false_northing=75000,
                                 globe=globe)
    return proj

def dist_bearing(lonlat0, lonlat1, radians=False):
    """
    Get distance [km] and bearing [deg] between two lon/lat points.

    Examples
    ---------
    >>> dist_km, bearing_deg = kkpy.util.dist_bearing([127.5,36], [130,37])

    Parameters
    ----------
    lonlat0 : 1D Array
        Array containing longitude and latitude of the first point. Longitude (latitude) should be at the first (second) element.
    lonlat1 : 1D Array
        Array containing longitude and latitude of the second point. Longitude (latitude) should be at the first (second) element.
    radians : boolean, optional
        If this is set to True, the unit of *bearing* is **radian**. The default is False (i.e. **degree**).

    Returns
    ---------
    distance : float
        Return distance between two points in **km**.
    bearing : ndarray
        Return bearing of two points in **degree**. If radians is True, the unit is **radians**.
    """
    from haversine import haversine
    from math import sin, cos, atan2
    lon0 = lonlat0[0]
    lat0 = lonlat0[1]
    lon1 = lonlat1[0]
    lat1 = lonlat1[1]

    dist = haversine((lat0,lon0),(lat1,lon1)) # km

    rlat0, rlon0, rlat1, rlon1 = np.radians((lat0, lon0, lat1, lon1))
    coslt0,  sinlt0  = cos(rlat0), sin(rlat0)
    coslt1,  sinlt1  = cos(rlat1), sin(rlat1)
    cosl0l1, sinl0l1 = cos(rlon1-rlon0), sin(rlon1-rlon0)
    cosc = sinlt0*sinlt1 + coslt0*coslt1*cosl0l1
    sinc = np.sqrt(1.0 - cosc**2)
    cosaz = (coslt0*sinlt1 - sinlt0*coslt1*cosl0l1) / sinc
    sinaz = sinl0l1*coslt1/sinc

    bearing = np.arctan(sinaz/cosaz)
    if not radians:
        bearing = np.degrees(bearing)

    return dist, bearing

def nanconvolve2d(slab, kernel, max_missing=0.99):
    """
    Get 2D convolution with missings ignored.

    Examples
    ---------
    Moving 2D standard deviation

    >>> import astropy.convolution
    >>> kernel = np.array(astropy.convolution.Box2DKernel(5))
    >>> c1 = kkpy.util.nanconvolve2d(arr2d, kernel)
    >>> c2 = kkpy.util.nanconvolve2d(arr2d*arr2d, kernel)
    >>> stddev2d = np.sqrt(c2 - c1*c1)

    Moving 2D average

    >>> import astropy.convolution
    >>> kernel = np.array(astropy.convolution.Box2DKernel(5))
    >>> avg2d = kkpy.util.nanconvolve2d(arr2d, kernel)

    Parameters
    ----------
    slab : 2D Array
        Input array to convolve. Can have numpy.nan or masked values.
    kernel : 1D Array
        Convolution kernel, must have sizes as odd numbers.
    max_missing : float, optional
        Float in (0,1), max percentage of missing in each convolution window is tolerated before a missing is placed in the result.

    Returns
    ---------
    result : 2D Array
        Return convolution result. Missings are represented as numpy.nans if they are in slab, or masked if they are masked in slab.

    Notes
    ---------
    This code is from Stack Overflow answer (https://stackoverflow.com/a/40416633/12272819), written by Jason (https://stackoverflow.com/users/2005415/jason).
    This is licensed under the Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA 3.0).
    Modified by Kwonil Kim in November 2020: modify docstring format, remove verbose argument, modify default value of max_missing, change numpy to np
    """
    from scipy.ndimage import convolve as sciconvolve

    assert np.ndim(slab)==2, "<slab> needs to be 2D."
    assert np.ndim(kernel)==2, "<kernel> needs to be 2D."
    assert kernel.shape[0]%2==1 and kernel.shape[1]%2==1, "<kernel> shape needs to be an odd number."
    assert max_missing > 0 and max_missing < 1, "<max_missing> needs to be a float in (0,1)."

    #--------------Get mask for missings--------------
    if not hasattr(slab,'mask') and np.any(np.isnan(slab))==False:
        has_missing=False
        slab2=slab.copy()

    elif not hasattr(slab,'mask') and np.any(np.isnan(slab)):
        has_missing=True
        slabmask=np.where(np.isnan(slab),1,0)
        slab2=slab.copy()
        missing_as='nan'

    elif (slab.mask.size==1 and slab.mask==False) or np.any(slab.mask)==False:
        has_missing=False
        slab2=slab.copy()

    elif not (slab.mask.size==1 and slab.mask==False) and np.any(slab.mask):
        has_missing=True
        slabmask=np.where(slab.mask,1,0)
        slab2=np.where(slabmask==1,np.nan,slab)
        missing_as='mask'

    else:
        has_missing=False
        slab2=slab.copy()

    #--------------------No missing--------------------
    if not has_missing:
        result=sciconvolve(slab2,kernel,mode='constant',cval=0.)
    else:
        H,W=slab.shape
        hh=int((kernel.shape[0]-1)/2)  # half height
        hw=int((kernel.shape[1]-1)/2)  # half width
        min_valid=(1-max_missing)*kernel.shape[0]*kernel.shape[1]

        # dont forget to flip the kernel
        kernel_flip=kernel[::-1,::-1]

        result=sciconvolve(slab2,kernel,mode='constant',cval=0.)
        slab2=np.where(slabmask==1,0,slab2)

        #------------------Get nan holes------------------
        miss_idx=zip(*np.where(slabmask==1))

        if missing_as=='mask':
            mask=np.zeros([H,W])

        for yii,xii in miss_idx:

            #-------Recompute at each new nan in result-------
            hole_ys=range(max(0,yii-hh),min(H,yii+hh+1))
            hole_xs=range(max(0,xii-hw),min(W,xii+hw+1))

            for hi in hole_ys:
                for hj in hole_xs:
                    hi1=max(0,hi-hh)
                    hi2=min(H,hi+hh+1)
                    hj1=max(0,hj-hw)
                    hj2=min(W,hj+hw+1)

                    slab_window=slab2[hi1:hi2,hj1:hj2]
                    mask_window=slabmask[hi1:hi2,hj1:hj2]
                    kernel_ij=kernel_flip[max(0,hh-hi):min(hh*2+1,hh+H-hi),
                                     max(0,hw-hj):min(hw*2+1,hw+W-hj)]
                    kernel_ij=np.where(mask_window==1,0,kernel_ij)

                    #----Fill with missing if not enough valid data----
                    ksum=np.sum(kernel_ij)
                    if ksum<min_valid:
                        if missing_as=='nan':
                            result[hi,hj]=np.nan
                        elif missing_as=='mask':
                            result[hi,hj]=0.
                            mask[hi,hj]=True
                    else:
                        result[hi,hj]=np.sum(slab_window*kernel_ij)

        if missing_as=='mask':
            result=np.ma.array(result)
            result.mask=mask

    return result

def nanstd2d(X, window_size):
    """
    Get 2D standard deviation of 2D array efficiently with missings ignored.

    Examples
    ---------
    >>> std2d = kkpy.util.nanstd2d(arr2d, 3)

    Parameters
    ----------
    X : 2D Array
        Array containing the data.
    window_size : float
        Window size of x and y. Window sizes of x and y should be same.

    Returns
    ---------
    std2d : 2D Array
        Return 2D standard deviation.
    """
    import astropy.convolution
    kernel = np.array(astropy.convolution.Box2DKernel(window_size))
    c1 = nanconvolve2d(X, kernel)
    c2 = nanconvolve2d(X*X, kernel)
    return np.sqrt(c2 - c1*c1)

def std2d(X, window_size):
    """
    Get 2D standard deviation of 2D array efficiently.

    Examples
    ---------
    >>> std2d = kkpy.util.std2d(arr2d, 3)

    Parameters
    ----------
    X : 2D Array
        Array containing the data.
    window_size : float or 1D array
        Window size. If array of two elements, window sizes of x and y will be window_size[0] and window_size[1], respectively.

    Returns
    ---------
    std2d : 2D Array
        Return 2D standard deviation.

    Notes
    ---------
    This code is from https://nickc1.github.io/python,/matlab/2016/05/17/Standard-Deviation-(Filters)-in-Matlab-and-Python.html, written by Nick Cortale.
    Modified by Kwonil Kim in November 2020: add docstring, modify function name
    """
    from scipy.ndimage.filters import uniform_filter

    r,c = X.shape
    X+=np.random.rand(r,c)*1e-6
    c1 = uniform_filter(X, window_size, mode='reflect')
    c2 = uniform_filter(X*X, window_size, mode='reflect')
    return np.sqrt(c2 - c1*c1)

def dbzmean(dbz_arr, outside_radar=-9999., noprecip=-9998., qced=-9997., axis=None):
    """
    Get linear-scale average of reflectivity (dBZ).

    Examples
    ---------
    >>> dbz_avg = kkpy.util.dbzmean(dbz_arr)

    Parameters
    ----------
    X : array_like
        Array containing the reflectivity in dBZ.
    outside_radar : float, optional
        Value indicating the grid point is located outside radar observable range. This will not be included in the count during the average.
    noprecip : float, optional
        Value indicating the grid point is clear echo. This will be considered as 0.0 mm6 m-3.
    qced : float, optional
        Value indicating the grid point is masked by QC process. This will not be included in the count during the average.
    axis : None or int or tuple of ints, optional
        Axis or axes along which the means are computed.

        .. versionadded:: 0.2.0

    Returns
    ---------
    dbz_avg : array_like
        Return averaged reflectivity in dBZ.

    """

    dbz_arr[dbz_arr == outside_radar] = np.nan
    dbz_arr[dbz_arr == qced] = np.nan

    dbz_lin = 10.**(dbz_arr/10.)
    dbz_lin[dbz_arr == noprecip] = 0.

    return 10*np.log10(np.nanmean(dbz_lin, axis=axis))

def icepop_events(KST=False):
    """
    Get datetimes for nine major evnets in ICE-POP 2018.

    Examples
    ---------
    >>> dts = kkpy.util.icepop_events()
    >>> for eventno in np.arange(9)+1:
    >>>     dt = dts[eventno]
    >>>     print(f'eventno = {eventno}, start = {dt[0]}, end = {dt[1]}')

    Parameters
    ----------
    KST : boolean, optional
        True if return the datetime in KST, rather than UTC.

    Returns
    ---------
    dict_datetimes : dictionary
        Return a dictionary of array containing start and finish datetime for each event.

    """
    import datetime

    if not KST:
        events = {
            1: [datetime.datetime(2017,11,24,20), datetime.datetime(2017,11,25,19)],
            2: [datetime.datetime(2017,12,23,20), datetime.datetime(2017,12,25, 3)],
            3: [datetime.datetime(2018, 1,22, 3), datetime.datetime(2018, 1,23, 0)],
            4: [datetime.datetime(2018, 2,27,23), datetime.datetime(2018, 3, 1, 3)],
            5: [datetime.datetime(2018, 3, 4, 8), datetime.datetime(2018, 3, 5,10)],
            6: [datetime.datetime(2018, 3, 7, 5), datetime.datetime(2018, 3, 9, 3)],
            7: [datetime.datetime(2018, 3,14,19), datetime.datetime(2018, 3,16,12)],
            8: [datetime.datetime(2018, 3,18,10), datetime.datetime(2018, 3,19, 9)],
            9: [datetime.datetime(2018, 3,19,21), datetime.datetime(2018, 3,21,14)],
        }
    else:
        events = {
            1: [datetime.datetime(2017,11,25, 5), datetime.datetime(2017,11,26, 4)],
            2: [datetime.datetime(2017,12,24, 5), datetime.datetime(2017,12,25,12)],
            3: [datetime.datetime(2018, 1,22,12), datetime.datetime(2018, 1,23, 9)],
            4: [datetime.datetime(2018, 2,28, 8), datetime.datetime(2018, 3, 1,12)],
            5: [datetime.datetime(2018, 3, 4,17), datetime.datetime(2018, 3, 5,19)],
            6: [datetime.datetime(2018, 3, 7,14), datetime.datetime(2018, 3, 9,12)],
            7: [datetime.datetime(2018, 3,15, 4), datetime.datetime(2018, 3,16,21)],
            8: [datetime.datetime(2018, 3,18,19), datetime.datetime(2018, 3,19,18)],
            9: [datetime.datetime(2018, 3,20, 6), datetime.datetime(2018, 3,21,23)],
        }

    return events

def icepop_sites():
    """
    Get metadata of ICE-POP 2018 sites.

    Examples
    ---------
    >>> # Plotting in Lon/Lat coordinates
    >>> for lon, lat, hgt, site in kkpy.util.icepop_sites():
    >>>     plt.plot(lon, lat, marker='o', color='red', markersize=3, alpha=0.5, transform=ccrs.PlateCarree())
    >>>     plt.text(lon+0.01, lat, site, verticalalignment='center', fontsize=10, transform=ccrs.PlateCarree())

    >>> # Plotting in Lon/Hgt coordinates
    >>> for lon, lat, hgt, site in kkpy.util.icepop_sites():
    >>>     plt.plot(lon, hgt/1e3, marker='o', color='red', markersize=3, alpha=0.5)
    >>>     plt.text(lon+0.01, hgt/1e3, site, verticalalignment='center', fontsize=10)

    Returns
    ---------
    tuple_sites : tuple
        Return a tuple containing longitude, latitude, height, and sitename for each ICE-POP 2018 supersite.

    """

    sites = (
        (128.866858, 37.770897,  36, 'GWU'),
        (128.805847, 37.738157, 175, 'BKC'),
        (128.758636, 37.686953, 855, 'CPO'),
        (128.718825, 37.677331, 773, 'DGW'),
        (128.699611, 37.665208, 789, 'MHS'),
        (128.670494, 37.643342, 772, 'YPO'),
        (128.564700, 38.250900,  18, 'SCW'),
        (128.629700, 38.087200,   4, 'YYO'),
        (128.540700, 38.007400, 146, 'YDO'),
        (128.821100, 37.898300,  10, 'JMO'),
        (129.028900, 37.613500,  58, 'OGO'),
        (129.124300, 37.507100,  40, 'DHW'),
        (128.377600, 37.562000, 532, 'MOO'),
        (128.394600, 37.377900, 303, 'PCO')
    )

    return sites

def icepop_extent():
    """
    Get extent of ICE-POP 2018 domain.

    Examples
    ---------
    >>> ax.set_extent(kkpy.util.icepop_extent(), crs=ccrs.PlateCarree())

    Returns
    ---------
    list_extent : list
        Return a list containing lon0, lon1, lat0, and lat1 of default ICE-POP 2018 domain.

    """
    return [128, 129.5, 37, 38.3]

def vel_atlas(D):
    """
    Get Atlas et al. (1973) fall velocity.

    Examples
    ---------
    >>> V_arr = kkpy.util.vel_atlas(D_arr)

    Parameters
    ----------
    D : array_like
        Array containing the diameter in mm.

    Returns
    ---------
    V : array_like
        Return a fall velocity of raindrop for a given diameter.
    """
    return 9.65 - 10.3 * np.exp(-0.6 * D)

def calc_dsdmoments(ND, D, dD, VD=None, Vatlas=False):
    """
    Get DSD moments from drop size distribution.

    Examples
    ---------
    >>> tdvd = kkpy.util.calc_dsdmoments(ND_arr, D_arr, dD_arr, VD=VD_arr)

    >>> poss = kkpy.util.calc_dsdmoments(ND_arr, D_arr, dD_arr, Vatlas=True)

    Parameters
    ----------
    ND : array_like
        Array containing the number concentration in mm-1 m-3. If 2D array, the diameter should be at the first axis (i.e. axis = 0).
    D : array_like
        Array containing the median of each diameter channel in mm.
    dD : array_like
        Array containing the spread of each diameter channel in mm.
    VD : array_like, optional
        Array containing the mean velocity of each diameter channel in m s-1. Vatlas should be False if VD contains data.
    Vatlas : boolean, optional
        True if VD is replaced by the velocity-diameter relationship of Atlas et al. (1973).

    Returns
    ---------
    dict_moments : dictionary
        Return a dictionary containing the DSD moments: z[mm6 m-3], Z[dBZ], and R[mm hr-1].
    """

    if Vatlas:
        VD = vel_atlas(D)

    dict_ = {}
    if ND.ndim == 2:
        dict_['z'] = np.sum(ND*(D**6)*dD, axis=1) # mm6 m-3
        dict_['R'] = np.pi*6*1e-4*np.sum(ND*VD*(D**3)*dD, axis=1) # mm hr-1
    else:
        dict_['z'] = np.sum(ND*(D**6)*dD) # mm6 m-3
        dict_['R'] = np.pi*6*1e-4*np.sum(ND*VD*(D**3)*dD) # mm hr-1

    dict_['Z'] = 10*np.log10(dict_['z']) # dBZ

    return dict_

def stats(x, y,
          fmtbias='.3f', fmtrmse='.3f',
          fmtstd='.3f', fmtcorr='.3f'):
    """
    Get performance matrix BIAS, RMSE, STD, and CORR values and its f-string.

    Examples
    ---------
    >>> stats, str_stat = kkpy.util.stats(xarr, yarr)
    >>> print(stats)
    >>> plt.text(3, 5, str_stat)

    >>> stats, str_stat = kkpy.util.stats(xarr, yarr, fmtbias='.1f', fmtcorr='.2f')

    Parameters
    ----------
    x : array_like
        Array containing multiple variables and observations.
    y : array_like
        Array containing multiple variables and observations. The shape should be same as **x**.
    fmtbias : str, optional
        String format for BIAS. Default is '.3f'.
    fmtrmse : str, optional
        String format for RMSE. Default is '.3f'.
    fmtstd : str, optional
        String format for STD. Default is '.3f'.
    fmtcorr : str, optional
        String format for CORR. Default is '.3f'.

    Returns
    ---------
    dict_moments : dictionary
        Return a dictionary containing the DSD moments: z[mm6 m-3], Z[dBZ], and R[mm hr-1].
    """
    import pandas as pd

    bias = np.nanmean(y - x)
    rmse = np.nanmean((y - x) ** 2)**0.5
    std = np.nanstd(y - x)
    corr = pd.Series(x).corr(pd.Series(y))

    stats = [bias, rmse, std, corr]
    str_stat = f'BIAS={bias:.3f}\nRMSE={rmse:.3f}\nSTD={std:.3f}\nCORR={corr:.3f}'

    return stats, str_stat

def summary(arr, cnt_print=10):
    """
    Get summary of the array.

    Parameters
    ----------
    arr : array_like
        Array containing multiple variables and observations.
    cnt_print : int, optional
        The number of elements to print.
    """
    import pandas as pd
    import matplotlib.pyplot as plt

    print(pd.Series(arr).describe())
    cnt = arr.size
    count_nan = np.count_nonzero(np.isnan(arr))
    print(f'\nNans   {count_nan} ({count_nan/arr.size:.1f} %)\n')

    if cnt > cnt_print:
        print(f'First {cnt_print} values:\n {arr[:cnt_print]}\n')
        print(f'Last {cnt_print} values:\n {arr[-cnt_print:]}\n')

        arr = arr[~np.isnan(arr)]
        arr_sorted = np.sort(arr)
        print(f'Lowest {cnt_print} values:\n {arr_sorted[:cnt_print]}\n')
        print(f'Highest {cnt_print} values:\n {arr_sorted[-cnt_print:]}\n')
        arr_unique = np.unique(arr)
        print(f'Lowest {cnt_print} unique values:\n {arr_unique[:cnt_print]}\n')
        print(f'Highest {cnt_print} unique values:\n {arr_unique[-cnt_print:]}\n')

        hist = plt.hist(arr)
        print('Histograms')
        for y, x in zip(hist[0], hist[1]):
            print(f'{x:.3f} -- {y:.0f}')

    return

def derivative(var, cnt_x, pixelsize=1, axis=0):
    """
    Get derivative of 1D or 2D array.

    Examples
    ---------
    >>> # Derivative of 1D array with a window size of 5.
    >>> # Note that the horizontal resolution of xarr is 0.1.
    >>> xarr = np.arange(100)*0.1
    >>> yarr = np.exp(xarr)
    >>> dydx = kkpy.util.derivative(yarr, 5, pixelsize=0.1)

    >>> # This example gives the derivative of arr2d along x axis with a window size of 20.
    >>> # Assume that you have an arr2d[iy,ix] with horizontal resolution of 3 km.
    >>> dzdx = kkpy.util.derivative(arr2d, 20, pixelsize=3000, axis=1)

    Parameters
    ----------
    var : array_like
        Array containing multiple variables and observations.
    cnt_x : int
        The number of the pixel for derivative. This should be odd number. The derivative will be (Y[i+N/2]-Y[i-N/2]) / (N-1).
    pixelsize : float, optional
        The size of one pixel of x axis. If default (1), the derivative means the rate of change on an interval of **one pixel** with resolution of 1.
    axis : int, optional
        Axis along which the derivatives are computed. Default is 0.

    Returns
    ---------
    derivative : array_like
        Return a derivative for a given array.
    """
    import astropy.convolution

    if cnt_x % 2 != 1:
        raise ValueError('cnt_x must be odd.')

    if var.ndim > 2:
        raise ValueError('dimension size of var should be <= 2')

    kernel = np.zeros(cnt_x)
    kernel[0] = 1
    kernel[-1] = -1

    if var.ndim == 2:
        if axis == 0:
            kernel = kernel[:,np.newaxis]
        if axis == 1:
            kernel = kernel[np.newaxis,:]

    result = astropy.convolution.convolve(
            var,
            kernel,
            boundary='extend',
            normalize_kernel=False,
            nan_treatment='fill',
            fill_value=np.nan,
            preserve_nan=True
        ) / (cnt_x-1)

    if var.ndim == 1:
        for j in np.arange(int(cnt_x/2)+1):
            result[j] = (var[j+int(cnt_x/2)]-var[0]) / (j+int(cnt_x/2))
        for j in np.arange(var.size-1-int(cnt_x/2), var.size):
            result[j] = (var[-1]-var[j-int(cnt_x/2)]) / (var.size-j+int(cnt_x/2)-1)

    if var.ndim == 2:
        if axis == 0:
            for j in np.arange(int(cnt_x/2)+1):
                result[j,:] = (var[j+int(cnt_x/2),:]-var[0,:]) / (j+int(cnt_x/2))
            for j in np.arange(var.shape[0]-1-int(cnt_x/2), var.shape[0]):
                result[j,:] = (var[-1,:]-var[j-int(cnt_x/2),:]) / (var.shape[0]-j+int(cnt_x/2)-1)
        if axis == 1:
            for j in np.arange(int(cnt_x/2)+1):
                result[:,j] = (var[:,j+int(cnt_x/2)]-var[:,0]) / (j+int(cnt_x/2))
            for j in np.arange(var.shape[1]-1-int(cnt_x/2), var.shape[1]):
                result[:,j] = (var[:,-1]-var[:,j-int(cnt_x/2)]) / (var.shape[1]-j+int(cnt_x/2)-1)

    derivative = result / pixelsize

    return derivative

def to_lower_resolution(arr_in, ratio_x, ratio_y, dBZ=False):
    """
    Get 2D array with lower resolution.

    Examples
    ---------
    >>> arr1 = np.random.rand(250,150) # shape: (250,150)
    >>> arr2 = to_lower_resolution(arr1,10,10) # shape: (25,15)
    >>> arr3 = to_lower_resolution(arr1,5,10) # shape: (50,15)

    Parameters
    ----------
    arr_in : 2D array
        Array you want to lower the resolution.
    ratio_x : int
        The reduction ratio over x axis. The size of **arr_in** xaxis should be able to be divided by **ratio_x**.
    ratio_y : int
        The reduction ratio over y axis. The size of **arr_in** yaxis should be able to be divided by **ratio_y**.
    dBZ : boolean, optional
        True if return a linear-scale average of reflectivity (dBZ).

    Returns
    ---------
    arr_out : 2D array
        Return an array with lower resolution.

    Notes
    ---------
    This code is from Stack Overflow answer (https://stackoverflow.com/a/14916963/12272819), written by Jaime (https://stackoverflow.com/users/110026/jaime).
    This is licensed under the Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA 3.0).
    """

    if arr_in.ndim != 2:
        raise ValueError("arr_in must be a 2D array")
    if ratio_x <= 1 or ratio_y <= 1:
        raise ValueError("ratio_x and ratio_y must be greater than 1")
    if not isinstance(ratio_x, int) or not isinstance(ratio_y, int):
        raise ValueError("ratio_x and ratio_y must be an integer")

    nx_in, ny_in = arr_in.shape
    nx_out, ny_out = nx_in/ratio_x, ny_in/ratio_y

    if not nx_out.is_integer():
        raise ValueError(f"check ratio_x: nx_out={nx_out}")
    if not ny_out.is_integer():
        raise ValueError(f"check ratio_y: ny_out={ny_out}")

    arr_out = arr_in.reshape(
        [int(nx_out), ratio_x, int(ny_out), ratio_y])
    if not dBZ:
        arr_out = arr_out.mean(axis=3).mean(axis=1)
    else:
        arr_out = dbzmean(arr_out, axis=3)
        arr_out = dbzmean(arr_out, axis=1)

    return arr_out

max_stars_repo_name: Kwonil-Kim/kkpy

max_issues_repo_name: Kwonil-Kim/kkpy

max_forks_repo_name: Kwonil-Kim/kkpy

reason: import numpy,from scipy,import astropy


text:
"""
Prepare vocabulary and initial word vectors.
"""
import json
import msgpack
import pickle
import argparse
import numpy as np
from collections import Counter

from utils import vocab, constant, helper

def parse_args():
    parser = argparse.ArgumentParser(description='Prepare vocab for relation extraction.')
    parser.add_argument('data_dir', help='TACRED directory.')
    parser.add_argument('squad_dir', help='SQuAD directory.')
    parser.add_argument('vocab_dir', help='Output vocab directory.')
    parser.add_argument('--glove_dir', default='dataset/glove', help='GloVe directory.')
    parser.add_argument('--wv_file', default='glove.840B.300d.txt', help='GloVe vector file.')
    parser.add_argument('--wv_dim', type=int, default=300, help='GloVe vector dimension.')
    parser.add_argument('--min_freq', type=int, default=0, help='If > 0, use min_freq as the cutoff.')
    parser.add_argument('--lower', action='store_true', help='If specified, lowercase all words.')

    args = parser.parse_args()
    return args

def process_squad(squad_msgpack):
    train, dev = squad_msgpack
    train_tokens = []
    dev_tokens = []
    for row in train:
        train_tokens += row[1]  # context
        train_tokens += row[5]  # question
    for row in dev:
        dev_tokens += row[1]
        dev_tokens += row[5]
    return train_tokens, dev_tokens

def main():
    args = parse_args()

    # input files
    train_file = args.data_dir + '/train.json'
    dev_file = args.data_dir + '/dev.json'
    test_file = args.data_dir + '/test.json'
    wv_file = args.glove_dir + '/' + args.wv_file
    wv_dim = args.wv_dim

    # output files
    helper.ensure_dir(args.vocab_dir)
    vocab_file = args.vocab_dir + '/vocab.pkl'
    emb_file = args.vocab_dir + '/embedding.npy'

    # load files
    print("loading files...")
    train_tokens = load_tokens(train_file)
    dev_tokens = load_tokens(dev_file)
    test_tokens = load_tokens(test_file)
    # processing squad intermediate files
    with open(args.squad_dir + '/intermediate.msgpack', 'rb') as squad_file:
        squad_msgpack = msgpack.load(squad_file, encoding='utf-8')
    squad_train, squad_dev = squad_msgpack
    squad_train_tokens, squad_dev_tokens = process_squad(squad_msgpack)
    if args.lower:
        train_tokens, dev_tokens, test_tokens = [[t.lower() for t in tokens] for tokens in \
                                                 (train_tokens, dev_tokens, test_tokens)]

        squad_train_tokens, squad_dev_tokens = [[t.lower() for t in tokens] for tokens in \
                                                (squad_train_tokens, squad_dev_tokens)]
    # load glove
    print("loading glove...")
    glove_vocab = vocab.load_glove_vocab(wv_file, wv_dim)
    print("{} words loaded from glove.".format(len(glove_vocab)))

    print("building vocab...")
    v = build_vocab(train_tokens + squad_train_tokens, glove_vocab, args.min_freq)

    print("calculating oov...")
    datasets = {'train': train_tokens, 'dev': dev_tokens, 'test': test_tokens}
    for dname, d in datasets.items():
        total, oov = count_oov(d, v)
        print("{} oov: {}/{} ({:.2f}%)".format(dname, oov, total, oov * 100.0 / total))

    print("building embeddings...")
    embedding = vocab.build_embedding(wv_file, v, wv_dim)
    print("embedding size: {} x {}".format(*embedding.shape))

    print("dumping to files...")
    with open(vocab_file, 'wb') as outfile:
        pickle.dump(v, outfile)
    np.save(emb_file, embedding)
    print("all done.")

    print('converting SQuAD dataset to ids')

    id2word = v
    word2id = dict([(id2word[idx], idx) for idx in range(len(id2word))])

    def to_id(row, unk_id=constant.UNK_ID):
        context_tokens = row[1]
        context_features = row[2]
        context_tags = row[3]
        context_ents = row[4]
        question_tokens = row[5]
        question_ids = [word2id[w] if w in word2id else unk_id for w in question_tokens]
        context_ids = [word2id[w] if w in word2id else unk_id for w in context_tokens]
        tag_ids = [constant.POS_TO_ID[w] if w in constant.POS_TO_ID else unk_id for w in context_tags]
        ent_ids = [constant.NER_TO_ID[w] if w in constant.NER_TO_ID else unk_id for w in context_ents]
        return [row[0], context_ids, context_features, tag_ids, ent_ids, question_ids] + row[6:]

    squad_train = list(map(to_id, squad_train))
    squad_dev = list(map(to_id, squad_dev))
    result = {
        'train': squad_train,
        'dev': squad_dev
    }
    # train: id, context_id, context_features, tag_id, ent_id,
    #        question_id, context, context_token_span, answer_start, answer_end
    # dev:   id, context_id, context_features, tag_id, ent_id,
    #        question_id, context, context_token_span, answer
    with open('dataset/SQuAD/data.msgpack', 'wb') as f:
        msgpack.dump(result, f)

def load_tokens(filename):
    with open(filename) as infile:
        data = json.load(infile)
        tokens = []
        for d in data:
            tokens += d['tokens']
    print("{} tokens from {} examples loaded from {}.".format(len(tokens), len(data), filename))
    return tokens

def build_vocab(tokens, glove_vocab, min_freq):
    """ build vocab from tokens and glove words. """
    counter = Counter(t for t in tokens)
    # if min_freq > 0, use min_freq, otherwise keep all glove words
    if min_freq > 0:
        v = sorted([t for t in counter if counter.get(t) >= min_freq], key=counter.get, reverse=True)
    else:
        v = sorted([t for t in counter if t in glove_vocab], key=counter.get, reverse=True)
    # add special tokens and entity mask tokens
    v = constant.VOCAB_PREFIX + entity_masks() + v
    print("vocab built with {}/{} words.".format(len(v), len(counter)))
    return v

def count_oov(tokens, vocab):
    c = Counter(t for t in tokens)
    total = sum(c.values())
    matched = sum(c[t] for t in vocab)
    return total, total - matched

def entity_masks():
    """ Get all entity mask tokens as a list. """
    masks = []
    subj_entities = list(constant.SUBJ_NER_TO_ID.keys())[2:]
    obj_entities = list(constant.OBJ_NER_TO_ID.keys())[2:]
    masks += ["SUBJ-" + e for e in subj_entities]
    masks += ["OBJ-" + e for e in obj_entities]
    return masks

if __name__ == '__main__':
    main()


text:
""" ---------------------------------------------------------------"""
"""          Neural networks 1 - Multilayer perceptron             """
""" ---------------------------------------------------------------"""
"""https://machinelearningmastery.com/handwritten-digit-recognition-using-convolutional-neural-networks-python-keras/"""

import numpy
from keras.datasets import mnist
from keras.models import Sequential
from keras.layers import Dense
from keras.layers import Dropout
from keras.utils import np_utils
from sklearn.model_selection import train_test_split

# fix random seed for reproducibility
seed = 7
numpy.random.seed(seed)

# load data
(X_train, y_train), (X_test, y_test) = mnist.load_data()

# flatten 28*28 images to a 784 vector for each image
num_pixels = X_train.shape[1] * X_train.shape[2]
X_train = X_train.reshape(X_train.shape[0], num_pixels).astype('float32')
X_test = X_test.reshape(X_test.shape[0], num_pixels).astype('float32')

# train-validation split
X_train, X_validation, y_train, y_validation = train_test_split(X_train, y_train, test_size=0.1, random_state=42)

# normalize inputs from 0-255 to 0-1
X_train = X_train / 255
X_validation = X_validation / 255
X_test = X_test / 255

# one hot encode outputs
y_train = np_utils.to_categorical(y_train)
y_validation = np_utils.to_categorical(y_validation)
y_test = np_utils.to_categorical(y_test)
num_classes = y_test.shape[1]

# define baseline model
def baseline_model():
    # create model
    model = Sequential()
    model.add(Dense(num_pixels, input_dim=num_pixels, kernel_initializer='normal', activation='relu'))
    model.add(Dense(num_classes, kernel_initializer='normal', activation='softmax'))
    # Compile model
    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])
    return model

# build the model
model = baseline_model()
# Fit the model
model.fit(X_train, y_train, validation_data=(X_validation, y_validation), epochs=10, batch_size=200, verbose=2)

# Final evaluation of the model
print("Final evaluation of the model")
scores = model.evaluate(X_test, y_test, verbose=1)
print("Baseline Error: %.2f%%" % (100 - scores[1] * 100))