text:
__author__ = 'wagnerca'
import pylab as plt
import numpy as np
from pylab import plot, show, savefig, xlim, figure, \
                hold, ylim, legend, boxplot, setp, axes

import seaborn as sns
import pandas as pd
sns.set_context("poster", font_scale=1.1)

def plot_topk(top_women,top_men,colors, xaxis_values, filename):
    """
    plots the percentage of people who show up in N language editions
    :param men:
    :param women:
    :param colors:
    :return:
    """
    f, axarr  = plt.subplots(5, 2, sharex=True) #sharey='row'
    plt.xlabel("Gender")
    plt.ylabel("Num Languages")

    row = 0
    col = 0
    print xaxis_values
    for i in range(0, len(xaxis_values)):
        alldata = np.concatenate((top_women[i],top_men[i]), axis=0)
        #print alldata

        bp = axarr[row, col].boxplot([top_women[i], top_men[i]], sym="") #notch=0, sym='+', vert=1, whis=1.5
        #sns.boxplot(data=[top_women[i], top_men[i]], sym = "", ax = axarr[row, col])
        #sns.stripplot(data=alldata, size=4, jitter=True, edgecolor="gray", ax = axarr[row, col])

        #plt.setp(bp['boxes'], color='black')
        #plt.setp(bp['whiskers'], color='black')
        #plt.setp(bp['fliers'], color='red', marker='+')
        axarr[row, col].annotate("bottom "+str(xaxis_values[i])+"%", xy=(0.9, 0.7), xycoords='axes fraction', fontsize=10,
                xytext=(-5, 5), textcoords='offset points',
                ha='right', va='bottom')
        #axarr[row, col].set_ylim(np.min(alldata), np.max(alldata))
        axarr[row, col].set_xticks([1, 2], ['women', 'men'])

        print "\n bottom "+str(xaxis_values[i])+"%"
        print (row, col)
        print "min %s, max %s"%(np.min(alldata), np.max(alldata))
        print "mean women %s (median women %s), mean men %s (median men %s)" % (str(np.mean(top_women[i])), str(np.median(top_women[i])), str(np.mean(top_men[i])), str(np.median(top_men[i])))

        row += 1
        if(row >= 5):
            row = 0
            col += 1
    f.savefig(filename, bbox_inches='tight')
    plt.close()

# This function takes an array of numbers and smoothes them out.
# Smoothing is useful for making plots a little easier to read.

def sliding_mean(data_array, window=5):
    data_array = np.array(data_array)
    new_list = []
    for i in range(len(data_array)):
        indices = range(max(i - window + 1, 0),
                        min(i + window + 1, len(data_array)))
        avg = 0
        for j in indices:
            avg += data_array[j]
        avg /= float(len(indices))
        new_list.append(avg)

    return np.array(new_list)

def plot_smooth_lines(years, mean_vals, sem_vals, filename):
    mean_PlyCount = sliding_mean(mean_vals,
                                 window=10)
    sem_PlyCount = sliding_mean(sem_vals.mul(1.96).values,
                            window=10)

    # You typically want your plot to be ~1.33x wider than tall.
    # Common sizes: (10, 7.5) and (12, 9)
    plt.figure(figsize=(12, 9))

    # Remove the plot frame lines. They are unnecessary chartjunk.
    ax = plt.subplot(111)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

    # Ensure that the axis ticks only show up on the bottom and left of the plot.
    # Ticks on the right and top of the plot are generally unnecessary chartjunk.
    ax.get_xaxis().tick_bottom()
    ax.get_yaxis().tick_left()

    # Limit the range of the plot to only where the data is.
    # Avoid unnecessary whitespace.
    plt.ylim(63, 85)

    # Make sure your axis ticks are large enough to be easily read.
    # You don't want your viewers squinting to read your plot.
    plt.xticks(range(1850, 2011, 20), fontsize=14)
    plt.yticks(range(65, 86, 5), fontsize=14)

    # Along the same vein, make sure your axis labels are large
    # enough to be easily read as well. Make them slightly larger
    # than your axis tick labels so they stand out.
    plt.ylabel("Ply per Game", fontsize=16)

    # Use matplotlib's fill_between() call to create error bars.
    # Use the dark blue "#3F5D7D" as a nice fill color.
    plt.fill_between(years, mean_PlyCount - sem_PlyCount,
                     mean_PlyCount + sem_PlyCount, color="#3F5D7D")

    # Plot the means as a white line in between the error bars.
    # White stands out best against the dark blue.
    plt.plot(years, mean_PlyCount, color="white", lw=2)

    plt.savefig(filename, bbox_inches="tight");

def plot_shaded_lines(my_xticks, y1, y2, error1, error2, ylab, xlab, filename):
    plt.figure(figsize=(6,6))
    from matplotlib import rcParams
    rcParams.update({'figure.autolayout': True})

    x = range(0, len(y1))
    plt.plot(x, y1, 'k-', color="blue",  label='men')
    plt.fill_between(x, y1-error1, y1+error1, facecolor='blue', alpha=.2)

    plt.plot(x, y2, 'k-', color="red",  label='women')
    plt.fill_between(x, y2-error2, y2+error2, facecolor='red', alpha=.2)

    #if isinstance(x, (int, long, float, complex)):
    #    plt.xlim(np.min(x), np.max(x))
    plt.gcf().subplots_adjust(bottom=0.3)

    plt.xticks(x, my_xticks)
    plt.xticks(rotation=70, fontsize=14)
    plt.yticks(fontsize=14)
    #plt.setp(ax.get_xticklabels(), rotation='vertical', fontsize=14)
    plt.ylabel(ylab, fontsize=14)
    plt.xlabel(xlab, fontsize=14)
    plt.legend()

    plt.savefig(filename)

def plot_percentage(women_percentage, men_percentage, colors, xaxis_values, filename):
    """
    plots the percentage of people who show up in N language editions
    :param men:
    :param women:
    :param colors:
    :return:
    """
    fig = plt.figure(figsize=(6,6))
    if len(colors) == 0:
        plt.gca().set_color_cycle(['pink', 'blue', 'yellow', 'red', 'black'])
    else:
        plt.gca().set_color_cycle(colors)
    #print xaxis_values
    #print women_percentage
    #print len(xaxis_values)
    #print len(women_percentage)
    plt.plot(xaxis_values, women_percentage[:len(xaxis_values)], linewidth=1)
    plt.plot(xaxis_values, men_percentage[:len(xaxis_values)], linewidth=1)
    plt.xlabel("Log Num Editions")
    plt.ylabel("Log Proportion")
    plt.xscale("log")
    plt.yscale("log")
    #plt.ysca:len(xaxis_values)lotle("log")
    plt.savefig(filename)
    plt.close()

def plot_cdf(_itemFreqList, labels, colors, path, xlabel, plotCCDF, xlog, ylog):
    """
    :param _itemFreqList: list of 2D array which consist of items and frequency of that item
            labels: list of names for the lists
            colors: colors of the lists
            path: where to save the plot
            xlabel: name of xaxis; y-axis is always ccdf
            xlog: boolean indicating if x axis is log or not
    :return:
    """
    fig = plt.figure(figsize=(6,6))
    if len(colors) == 0:
        plt.gca().set_color_cycle(['pink', 'blue', 'yellow', 'red', 'black'])
    else:
        plt.gca().set_color_cycle(colors)

    for _itemFreq in _itemFreqList:
        _itemFreq_sorted = np.array(sorted(_itemFreq,key=lambda x: x[0]))

        numberOfFrequency_sorted = _itemFreq_sorted[:, 1]
        numberOfReviews_sorted = _itemFreq_sorted[:,0]

        cdf = numberOfFrequency_sorted.cumsum()/ float(numberOfFrequency_sorted.sum())

        if plotCCDF==1:
            ccdf= 1-cdf
            sign = "\geq"
        else:
            ccdf = cdf
            sign = "<"

        plt.plot(numberOfReviews_sorted, ccdf, linewidth=1)

    #plt.title("ccdf plot")
    if xlog==1:
        plt.xscale('log')
        plt.xlabel('$log \/ '+xlabel+' \/ x $')
    else:
        plt.xlabel('$'+xlabel+' \/ x $')
    if ylog == 1:
        plt.yscale('log')
        plt.ylabel("$ log \/ P(X "+sign+" x) $")
    else:
        plt.ylabel("$ P(X "+sign+" x)$")

    plt.legend(labels, loc='upper right')
    plt.grid(False)
    plt.savefig(path)

#def plot_rank_size(_itemFreqList, labels, colors, path, xlabel, ylabel, xlog, ylog):
    #TODO

def plot_facet_dist(data, facet_col, data_col, filename):
    fig = plt.figure(figsize=(6,6))
    g = sns.FacetGrid(data, col=facet_col, size=4, aspect=1.5)
    g.map(sns.distplot, data_col, kde=False)
    g.set(yscale='log')
    g.add_legend()
    plt.grid(False)
    plt.savefig(filename)
    plt.close()

def plot_numedition_size(_itemFreqList, labels, colors, path, xlabel, ylabel, xlog, ylog):

    fig = plt.figure(figsize=(6,6))

    if len(colors) == 0:
        plt.gca().set_color_cycle(['pink', 'blue', 'yellow', 'red', 'black'])
    else:
        plt.gca().set_color_cycle(colors)

    numberOfFrequency_sorted = {}
    for _itemFreq in _itemFreqList:
        _itemFreq_sorted = np.array(sorted(_itemFreq,key=lambda x: x[0]))

        numberOfFrequency_sorted = _itemFreq_sorted[:, 1]
        numberOfReviews_sorted = _itemFreq_sorted[:,0]

        plt.plot(numberOfReviews_sorted, numberOfFrequency_sorted, linewidth=1)

    #plt.title("ccdf plot")
    if xlog==1:
        plt.xscale('log')
        plt.xlabel('$log \/ '+xlabel+'  $')
    else:
        plt.xlabel(xlabel)
    if ylog == 1:
        plt.yscale('log')
        plt.ylabel("$ log "+ylabel)
    else:
        plt.ylabel(ylabel)

    plt.legend(labels, loc='upper right')
    plt.grid(False)
    plt.savefig(path, bbox_inches='tight')
    plt.close()

def plotratio(ratio, ratiofit, fakeratio, styles, labels, title, path, xlabel, ylabel, xlog, ylog):
    xval = list(ratio.keys())
    yval = list(ratio.values())
    fit_xval = list(ratiofit[:, 0])
    fit_yval = list(ratiofit[:, 1])

    fxval = list(fakeratio.keys())
    fyval = list(fakeratio.values())

    fig = plt.figure(figsize=(6,6))
    ax = fig.add_subplot(111)
    ax.set_autoscale_on(False)

    plt.title(title, loc='left', color='r')
    #plt.text(0.5, 0.5, title,
    # horizontalalignment='center',
    # verticalalignment='center',
    # transform = ax.transAxes)
    plotline(xval, yval, styles[0], path, xlabel, ylabel, xlog, ylog)
    plotline(fit_xval, fit_yval, styles[1], path, xlabel, ylabel, xlog, ylog)
    plotline(fxval, fyval, styles[2], path, xlabel, ylabel, xlog, ylog)
    ax.set_ylim([-1, np.max(yval)+3.5])
    ax.set_xlim([0, np.max(xval)])
    ax.legend(numpoints=3, scatterpoints=3)
    ax.legend(labels, loc='upper right')
    #plt.grid(False)
    plt.savefig(path, bbox_inches='tight')

def plotline(xvals, yvals, style, path, xlabel, ylabel, xlog, ylog):

    plt.plot(xvals, yvals, style, linewidth=2, markersize=5)

    #plt.title("ccdf plot")
    if xlog==1:
        plt.xscale('log')
        plt.xlabel('log \/ '+xlabel+'  ')
    else:
        plt.xlabel(xlabel)
    if ylog == 1:
        plt.yscale('log')
        plt.ylabel("log "+ylabel)
    else:
        plt.ylabel(ylabel)

def rank_size_plot(data, attr, label, filename):
    plt.figure(figsize=(6,6))

    plt.subplot(111)

    men = data[data.gender == 'male']
    rank_m = men[attr].sort(ascending=False, inplace=False)
    women = data[data.gender == 'female']
    rank_w = women.edition_count.sort(ascending=False, inplace=False)

    plt.plot(np.arange(rank_m.shape[0], dtype=np.float) / rank_m.shape[0] + 0.001, rank_m, label='Men', linestyle='none', marker='.', alpha=0.5)
    plt.plot(np.arange(rank_w.shape[0], dtype=np.float) / rank_w.shape[0] + 0.001, rank_w, label='Women', linestyle='none', marker='.', alpha=0.5)
    plt.xlabel('Normalized Rank')
    plt.ylabel(label)
    plt.yscale('log')
    plt.xscale('log')
    plt.legend()

    plt.xlim([0.001, 1.001])
    plt.savefig(filename, bbox_inches='tight')

def write_mannwithneyu(U, p, arr1, arr2, logfile):
    gsize1 = len(arr1)
    gsize2 = len(arr2)
    prob = U/float(gsize1*gsize2)
    if(p < 0.05):
        logfile.write("\n U=%s, p=%s, S1=%s, S2=%s, median1(mean1)=%s(%s), median2(mean2)=%s(%s), U/S1*S2 = %s  ******" % (U, p, gsize1, gsize2, np.median(arr1), np.mean(arr1), np.median(arr2), np.mean(arr2), prob))
        logfile.write("\n U statistic divide by the product of the two sample sizes --> probability that a score randomly drawn from population A will be greater than a score randomly drawn from population B is %s"% str(U/float(gsize1*gsize2)))
    else:
        logfile.write("\n  U=%s, p=%s, S1=%s, S2=%s, median1(mean1)=%s(%s), median2(mean2)=%s(%s), U/S1*S2 = %s " %  (U, p, gsize1, gsize2, np.median(arr1), np.mean(arr1), np.median(arr2), np.mean(arr2), prob))
        logfile.write("\n probability that a score randomly drawn from population A will be greater than a score randomly drawn from population B is %s"% str(U/float(gsize1*gsize2)))

#if __name__ == '__main__':
#    #test cdf function
#    from scipy.stats import itemfreq
#    item_frequency_female = itemfreq(np.array([2, 2, 2, 2, 2, 3, 3, 5, 6, 8]))
#    item_frequency_male = itemfreq(np.array([1, 1,1 , 1, 1, 1, 1, 1, 4, 4, 5]))
#    print item_frequency_female
#    print item_frequency_male
#    plot_cdf([item_frequency_female, item_frequency_male], ['female', 'male'], ['pink', 'blue'], "./test_ccdf.png", "num", 1, 0, 0)

max_stars_repo_name: clauwag/WikipediaGenderInequality

max_issues_repo_name: clauwag/WikipediaGenderInequality

max_forks_repo_name: clauwag/WikipediaGenderInequality


text:
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np

class crimemap:
    def __init__(self, map_dir):
        self.shapefile0 = gpd.read_file(map_dir)
        self.shapefile0.plot(figsize=(10, 10))

    def mapshow(self, x_min=-73.59, x_max=-73.55, y_min=45.49, y_max=45.53, grid_size=0.002):
        self.extent = [x_min, x_max, y_min, y_max]
        plt.xlim(x_min, x_max)
        plt.ylim(y_min, y_max)
        self.x_sticks = np.arange(x_min, x_max + grid_size/2, grid_size)
        self.y_sticks = np.arange(y_min, y_max + grid_size/2, grid_size)
        plt.xticks(self.x_sticks)
        plt.yticks(self.y_sticks)
        plt.grid()
        plt.show()

    def crime_calc(self, threshold=0.5):
        series_x = self.shapefile0.geometry.x
        series_y = self.shapefile0.geometry.y
        crime_hist, _, _ = np.histogram2d(series_y, series_x, bins=[self.y_sticks, self.x_sticks])
        crime_array = crime_hist.flatten()
        threshold_value = np.percentile(crime_array, threshold*100)
        crime_hist_cp = np.copy(crime_hist)
        crime_hist_cp[crime_hist_cp <= threshold_value] = 0
        crime_hist_cp[crime_hist_cp > threshold_value] = 1
        mean_ = np.mean(crime_hist)
        std_dev_ = np.std(crime_hist)
        print("Mean: {0:.4f}".format(mean_))
        print("Standard deviation: {0:.4f}".format(std_dev_))
        print("Threshold value: {0:.4f}".format(threshold_value))
        plt.imshow(crime_hist_cp, origin='lower', extent=self.extent, interpolation='nearest', aspect='auto')
        plt.show()

ext: py

max_stars_repo_name: TJ-Jimmy/Ai

max_issues_repo_name: TJ-Jimmy/Ai

max_forks_repo_name: TJ-Jimmy/Ai

reason: import numpy


text:
# Compared pedagogical method and neural network method by comparing data acquired from
# 'logicalRuleEvaluation.py' and 'generate_testing_dataSet.py (seen in enst cloud service)'
# %%
import torch
import numpy as np

import random
import os

# device=torch.device("cuda" if torch.cuda.is_available() else 'cpu')
device=torch.device('cpu')

# %%
# create a data storing path

num_epochs=100
steps=288

random.seed(0)

dirName='11082021/model_save'
if not os.path.exists(dirName):
    os.makedirs(dirName)

other_data_name='11082021/section6_seen'
if not os.path.exists(other_data_name):
    os.makedirs(other_data_name)

# if os.path.exists("./../v3/"+other_data_name+"/pedagogicalRuleEvaluation_v2_48_with_seen.log"):
#     f = open("./../v3/"+other_data_name+"/pedagogicalRuleEvaluation_v2_48_with_seen.log", "r+")
# else:
#     f = open("./../v3/"+other_data_name+"/pedagogicalRuleEvaluation_v2_48_with_seen.log", "w")

# import logging

# logging.basicConfig(level=logging.INFO, format='%(message)s')
# logger = logging.getLogger()
# logger.addHandler(logging.FileHandler("./../v3/"+other_data_name+"/pedagogicalRuleEvaluation_v2_48_with_seen.log", 'a'))
# print = logger.info

# %%
# load data for neural network
# totalEpochRewardsNN=np.load('./../v3/section6/not_seen_states'+'/totalEpochRewards_2_v3_with_not_seen_100.npy',allow_pickle=True)
# replayMemoryNN=np.load('./../v3/section6/not_seen_states'+'/replayMemory_2_v3_with_not_seen_100.npy',allow_pickle=True)

# # %%
# # load data for the rule system
# totalEpochRewardsRule=np.load('./../v3/section6/not_seen_states'+'/totalEpochRewards_2_v3_with_not_seen_RULES_100.npy',allow_pickle=True)
# replayMemoryRule=np.load('./../v3/section6/not_seen_states'+'/replayMemory_2_v3_with_not_seen_RULES_100.npy',allow_pickle=True)
# notSeenStatesRule=np.load('./../v3/section6/not_seen_states'+'/states_with_not_seen_RULES_100.npy',allow_pickle=True)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 10 epoch

# totalEpochRewardsNN=np.load('./../v3/section6/not_seen_states'+'/totalEpochRewards_2_v3_with_not_seen_100.npy',allow_pickle=True)
# replayMemoryNN=np.load('./../v3/section6/not_seen_states'+'/replayMemory_2_v3_with_not_seen_100.npy',allow_pickle=True)

# %%
# load data for the rule system
# totalEpochRewardsRule=np.load('./../v3/section6/not_seen_states'+'/totalEpochRewards_2_v3_with_not_seen_not_seen_RULES_100.npy',allow_pickle=True)
# replayMemoryRule=np.load('./../v3/section6/not_seen_states'+'/replayMemory_2_v3_with_not_seen_not_seen_RULES_100.npy',allow_pickle=True)
# notSeenStatesRule=np.load('./../v3/section6/not_seen_states'+'/states_with_not_seen_not_seen_RULES_100.npy',allow_pickle=True)

# %%
# %%
# %%
# load data for the rule system
# totalEpochRewardsRule=np.load('./../v3/section6/not_seen_decomp/totalEpochRewards_2_v3_with_not_seen_RULES_100.npy',allow_pickle=True)
# replayMemoryRule=np.load('./../v3/section6/not_seen_decomp/replayMemory_2_v3_with_not_seen_RULES_100.npy',allow_pickle=True)
# notSeenStatesRule=np.load('./../v3/section6/not_seen_decomp/states_with_not_seen_RULES_100.npy',allow_pickle=True)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# totalEpochRewardsNN=np.load('./../v3/'+other_data_name+'/totalEpochRewards_light_lstm_seen.npy',allow_pickle=True)
# replayMemoryNN=np.load('./../v3/'+other_data_name+'/replayMemory_light_lstm_seen.npy',allow_pickle=True)

# totalEpochRewardsRule=np.load('./../v3/'+other_data_name+'/totalEpochRewards_light_service_structure_lstm_RULES_100_seen.npy',allow_pickle=True)
# replayMemoryRule=np.load('./../v3/'+other_data_name+'/replayMemory_light_service_structure_lstm_RULES_100_seen.npy',allow_pickle=True)
# notSeenStatesRule=np.load('./../v3/'+other_data_name+'/notseen_states_light_service_structure_lstm_RULES_100_seen.npy',allow_pickle=True)

totalEpochRewardsNN=np.load('./../v3/'+other_data_name+'/totalEpochRewards_light_lstm_unseen.npy',allow_pickle=True)
replayMemoryNN=np.load('./../v3/'+other_data_name+'/replayMemory_light_lstm_unseen.npy',allow_pickle=True)

totalEpochRewardsRule=np.load('./../v3/'+other_data_name+'/totalEpochRewards_light_service_structure_lstm_RULES_100_unseen.npy',allow_pickle=True)
replayMemoryRule=np.load('./../v3/'+other_data_name+'/replayMemory_light_service_structure_lstm_RULES_100_unseen.npy',allow_pickle=True)
notSeenStatesRule=np.load('./../v3/'+other_data_name+'/notseen_states_light_service_structure_lstm_RULES_100_unseen.npy',allow_pickle=True)

# calculate epoch rewards

def judgePos(notSeenStatesRule,i,epoch):
    if i<=notSeenStatesRule.shape[1]-2:
        if notSeenStatesRule[0,i+1]>=(epoch+1)*steps:
            return True
    elif i==notSeenStatesRule.shape[1]-1:

        return True

    return False

seenEpochRewards=np.full([1,num_epochs],288)*4
idx_seen=np.full([1,num_epochs],None)

if notSeenStatesRule[0,0]!=None:
    totalEpochRewardsNN=np.full([1,num_epochs],0)

num=0
num_delete=0
for epoch in range(num_epochs):
    print(f'epoch: {epoch}')
    rewards=0
    restIdx=np.array([[i for i in range(epoch*steps,(epoch+1)*steps)]]).reshape(1,-1)

    idx_seen[0,epoch]=restIdx
    toDelete=np.full([1,1],None)

    if notSeenStatesRule[0,0]!=None:

        for i in range(notSeenStatesRule.shape[1]):

            if notSeenStatesRule[0,i]>=(epoch*steps) and notSeenStatesRule[0,i]<((epoch+1)*steps):

                if toDelete[0,0]==None:
                    toDelete[0,0]=notSeenStatesRule[0,i]-epoch*steps
                    toDelete=toDelete.astype(np.int)
                else:
                    toDeleteInt=np.array([[notSeenStatesRule[0,i]-epoch*steps]]).astype(np.int)
                    toDelete=np.concatenate((toDelete,toDeleteInt),axis=1)

            print(f'i: {notSeenStatesRule[0,i]}')
            result=judgePos(notSeenStatesRule,i,epoch)

            if result:
                if toDelete[0,0]!=None:
                    idx_seen[0,epoch]=np.delete(idx_seen[0,epoch],toDelete,axis=1)
                    num_not_seen=toDelete.shape[1]
                    seenEpochRewards[0,epoch]-=num_not_seen*4
                    num_delete+=toDelete.shape[1]
                    total=restIdx.shape[1]+toDelete.shape[1]

                    if total!=288:
                        print('================================================================')

                num=num+restIdx.shape[1]

                for idx in idx_seen[0,epoch][0]:

                    totalEpochRewardsNN[0,epoch]+=replayMemoryNN[idx,2].item()

                break

print(f'idx_not_seen num: {num_delete}')
print(f'idx_seen num: {num}')
if notSeenStatesRule[0,0]!=None:
    replayMemoryNN=np.delete(replayMemoryNN,notSeenStatesRule,axis=0)
print(f'replayMemoryNN shape: {replayMemoryNN.shape}')
print(f'replayMemoryRule shape: ',{replayMemoryRule.shape})

idx_seen_2=np.array([[i for i in range(0,100*288)]]).reshape(1,-1)
if notSeenStatesRule[0,0]!=None:
    idx_seen_2=np.delete(idx_seen_2,notSeenStatesRule,axis=1)
print(f'idx_seen_2 shape: {idx_seen_2.shape}')
print(f'notSeenStatesRule shape: {notSeenStatesRule.shape}')
print(f'notSeenStatesRule shape: {notSeenStatesRule}')

# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")

# print("================================================================================================")
# print("================================================================================================")
# print("================================================================================================")
print("================================================================================================")
print("================================================================================================")

print("================================================================================================")
num=0
for i in range(replayMemoryNN.shape[0]):
    # print('current step: ',i)
    if (replayMemoryNN[i,0]==replayMemoryRule[i,0]).all():
        if (replayMemoryNN[i,1][0]==replayMemoryRule[i,1][0]) and (replayMemoryNN[i,1][1]==replayMemoryRule[i,1][1]):
            num+=1

print(num)
print(f'replayMemoryNN shape: {replayMemoryNN.shape}')

# np.save('./../v3/section6/not_seen_states'+f'/totalEpochRewards_2_v3_with_not_seen_NN_100.npy',totalEpochRewardsNN)
# np.save('./../v3/section6/not_seen_states'+f'/replayMemory_2_v3_with_not_seen_NN_100.npy',replayMemoryNN)
# np.save('./../v3/section6/not_seen_states'+f'/seenEpochRewards_100.npy',seenEpochRewards)
# np.save('./../v3/section6/not_seen_states'+f'/idx_seen_100.npy',idx_seen)
# np.save('./../v3/section6/not_seen_states'+f'/idx_seen_100_2.npy',idx_seen_2)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# 10 epochs

# np.save('./../v3/'+other_data_name+f'/totalEpochRewards_light_service_structure_lstm_NN_100_seen.npy',totalEpochRewardsNN)
# np.save('./../v3/'+other_data_name+f'/replayMemory_light_service_structure_lstm_NN_100_seen.npy',replayMemoryNN)
# np.save('./../v3/'+other_data_name+f'/seenEpochRewards_light_service_structure_lstm_100_seen.npy',seenEpochRewards)
# np.save('./../v3/'+other_data_name+f'/idx_light_service_structure_lstm_100_seen.npy',idx_seen)
# np.save('./../v3/'+other_data_name+f'/idx_light_service_structure_lstm_100_2_seen.npy',idx_seen_2)

np.save('./../v3/'+other_data_name+f'/totalEpochRewards_light_service_structure_lstm_NN_100_unseen.npy',totalEpochRewardsNN)
np.save('./../v3/'+other_data_name+f'/replayMemory_light_service_structure_lstm_NN_100_unseen.npy',replayMemoryNN)
np.save('./../v3/'+other_data_name+f'/seenEpochRewards_light_service_structure_lstm_100_unseen.npy',seenEpochRewards)
np.save('./../v3/'+other_data_name+f'/idx_light_service_structure_lstm_100_unseen.npy',idx_seen)
np.save('./../v3/'+other_data_name+f'/idx_light_service_structure_lstm_100_2_unseen.npy',idx_seen_2)


text:
import torch
import torchvision.models as torch_models
import flaxvision.models as flax_models
from flax import nn
import numpy as np
import jax.numpy as jnp
from jax import random
from jax.config import config
config.enable_omnistaging()

import unittest, os
import logging

RNG = random.PRNGKey(0)

MODELS_LIST = [
    'vgg11', 'vgg11_bn', 'vgg13', 'vgg13_bn', 'vgg16', 'vgg16_bn', 'vgg19', 'vgg19_bn', 'resnet18', 'resnet34',
    'resnet50', 'resnet101', 'resnet152', 'resnext50_32x4d', 'resnext101_32x8d', 'wide_resnet50_2', 'wide_resnet101_2',
    'inception_v3', 'densenet121', 'densenet161', 'densenet169', 'densenet201', 'fcn_resnet50', 'fcn_resnet101',
    'deeplabv3_resnet50', 'deeplabv3_resnet101'
]

class TestPretrained(unittest.TestCase):

  def test_outputs(self):
    log = logging.getLogger(__name__)

    flax_input = jnp.ones((1, 224, 224, 3))
    torch_input = torch.ones([1, 3, 224, 224])

    flax_inception_input = jnp.ones((1, 299, 299, 3))
    torch_inception_input = torch.ones([1, 3, 299, 299])

    for key in MODELS_LIST:
      log.info(f'testing {key}')
      torch_model, (flax_model, flax_params) = self._get_model(key)
      torch_model.eval()

      if key == 'inception_v3':
        with torch.no_grad():
          torch_out = torch_model(torch_inception_input).detach().numpy()
        flax_out = flax_model.apply(flax_params, flax_inception_input, mutable=False)
        self.assertLess(
            np.mean(np.abs(flax_out[0] - torch_out)), 0.0001, "PyTorch and Flax models' outputs don't match.")
      else:
        with torch.no_grad():
          if 'fcn' in key or 'deeplab' in key:
            torch_out = torch_model(torch_input)['out'].detach().numpy()
          else:
            torch_out = torch_model(torch_input).detach().numpy()
        flax_out = flax_model.apply(flax_params, flax_input, mutable=False)
        self.assertLess(np.mean(np.abs(flax_out - torch_out)), 0.0001, "PyTorch and Flax models' outputs don't match.")

      del torch_model, flax_model, flax_params

  def _get_model(self, key):
    if key == 'vgg11':
      return (torch_models.vgg11(True), flax_models.vgg11(RNG))
    if key == 'vgg11_bn':
      return (torch_models.vgg11_bn(True), flax_models.vgg11_bn(RNG))
    if key == 'vgg13':
      return (torch_models.vgg13(True), flax_models.vgg13(RNG))
    if key == 'vgg13_bn':
      return (torch_models.vgg13_bn(True), flax_models.vgg13_bn(RNG))
    if key == 'vgg16':
      return (torch_models.vgg16(True), flax_models.vgg16(RNG))
    if key == 'vgg16_bn':
      return (torch_models.vgg16_bn(True), flax_models.vgg16_bn(RNG))
    if key == 'vgg19':
      return (torch_models.vgg19(True), flax_models.vgg19(RNG))
    if key == 'vgg19_bn':
      return (torch_models.vgg19_bn(True), flax_models.vgg19_bn(RNG))
    if key == 'resnet18':
      return (torch_models.resnet18(True), flax_models.resnet18(RNG))
    if key == 'resnet34':
      return (torch_models.resnet34(True), flax_models.resnet34(RNG))
    if key == 'resnet50':
      return (torch_models.resnet50(True), flax_models.resnet50(RNG))
    if key == 'resnet101':
      return (torch_models.resnet101(True), flax_models.resnet101(RNG))
    if key == 'resnet152':
      return (torch_models.resnet152(True), flax_models.resnet152(RNG))
    if key == 'resnext50_32x4d':
      return (torch_models.resnext50_32x4d(True), flax_models.resnext50_32x4d(RNG))
    if key == 'resnext101_32x8d':
      return (torch_models.resnext101_32x8d(True), flax_models.resnext101_32x8d(RNG))
    if key == 'wide_resnet50_2':
      return (torch_models.wide_resnet50_2(True), flax_models.wide_resnet50_2(RNG))
    if key == 'wide_resnet101_2':
      return (torch_models.wide_resnet101_2(True), flax_models.wide_resnet101_2(RNG))
    if key == 'inception_v3':
      return (torch_models.inception_v3(True), flax_models.inception_v3(RNG))
    if key == 'densenet121':
      return (torch_models.densenet121(True), flax_models.densenet121(RNG))
    if key == 'densenet161':
      return (torch_models.densenet161(True), flax_models.densenet161(RNG))
    if key == 'densenet169':
      return (torch_models.densenet169(True), flax_models.densenet169(RNG))
    if key == 'densenet201':
      return (torch_models.densenet201(True), flax_models.densenet201(RNG))
    if key == 'fcn_resnet50':
      return (torch_models.segmentation.fcn_resnet50(True), flax_models.fcn_resnet50(RNG))
    if key == 'fcn_resnet101':
      return (torch_models.segmentation.fcn_resnet101(True), flax_models.fcn_resnet101(RNG))
    if key == 'deeplabv3_resnet50':
      return (torch_models.segmentation.deeplabv3_resnet50(True), flax_models.deeplabv3_resnet50(RNG))
    if key == 'deeplabv3_resnet101':
      return (torch_models.segmentation.deeplabv3_resnet101(True), flax_models.deeplabv3_resnet101(RNG))

if __name__ == '__main__':
  logging.basicConfig(level=logging.INFO)
  unittest.main()

max_stars_repo_name: rolandgvc/flaxvision

max_issues_repo_name: rolandgvc/flaxvision

max_forks_repo_name: rolandgvc/flaxvision

reason: import numpy,import jax,from jax


text:
import gdxpds as gd
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

### This file reads optimal load scheduling from Julia ouput and generates aggregated
### consumption of alll categories, as in as in section 7.2.3 - step 2, following project report.
### Result is saved in csv and then converted to inc from xls2gdx ibrary in excel VBA.

# Sceanario and iteration
scenario='Expected'      ### choose between Expected, RapidScenario and SlowScenario
adoption=0.25            ### choose
iter='5it'               ### choose #it

# EV (categ,RRR,YYY,EH)
evnumbers=pd.read_csv('JuliatoBALMOREL/evsharesfin_'+scenario+'.csv',header=[0,1,2],index_col=0)
evnumbers.iloc[0:12]=evnumbers.iloc[0:12]*0.8 #houses only 80% at home
evnumbers.iloc[12:24]=evnumbers.iloc[12:24]*0.1 #apartments only 10% at home
categories=evnumbers.index
coefficients=pd.read_csv('JuliatoBALMOREL/evcoefficients.csv',skiprows=1,names=categories)/100 #[8*30]
years=['2020','2025','2030','2040']
lista=list(evnumbers.index.array)
lista.append('TOT')
EV=pd.DataFrame(index=range(8736),columns=range(1000))
EV.columns=pd.MultiIndex.from_product([lista,['DK1','DK2'],years,['BEV','BEVTOT','TOT','PHEV','PHEVTOT']])
TOTTOT=pd.DataFrame()
for DK in ['DK1','DK2']:
    for year in years:
        chargingBEV = pd.read_csv('JuliatoBALMOREL/Juliaout/FlexEV/FlexEV_' + DK + '_y' + year+ '.csv', header=None)
        chargingPHEV = pd.read_csv('JuliatoBALMOREL/Juliaout/FlexEV/FlexPHEV_' + DK + '_y' + year + '.csv', header=None)
        for category in evnumbers.index:
            weights=[evnumbers.loc[category][(DK,'BEV',year)],evnumbers.loc[category][(DK,'PHEV',year)]]
            EV[category,DK,year,'BEV']=sum(chargingBEV.iloc[:, i] * coefficients[category][i] for i in np.arange(20))
            EV[category,DK,year,'PHEV']=sum(chargingPHEV.iloc[:, i] * coefficients[category][i+20] for i in np.arange(10))
            EV[category, DK, year, 'TOT']=(EV[category,DK,year,'BEV'] * weights[0] +EV[category,DK,year,'PHEV'] * weights[1])/sum(weights)
            EV[category, DK, year, 'BEVTOT']=EV[category,DK,year,'BEV'] * evnumbers.loc[category][(DK,'BEV',year)]
            EV[category, DK, year, 'PHEVTOT']=EV[category,DK,year,'PHEV'] * evnumbers.loc[category][(DK,'PHEV',year)]
        EV['TOT', DK, year, 'BEVTOTT'] = EV[EV.columns[(EV.columns.get_level_values(1)==DK)&(EV.columns.get_level_values(2)==year) & (EV.columns.get_level_values(3)=='BEVTOT')]].sum(axis=1)
        EV['TOT', DK, year, 'PHEVTOTT'] = EV[EV.columns[(EV.columns.get_level_values(1)==DK)&(EV.columns.get_level_values(2)==year) & (EV.columns.get_level_values(3)=='PHEVTOT')]].sum(axis=1)
        TOTTOT[DK,year,'EV']=EV['TOT', DK, year, 'BEVTOTT']+EV['TOT', DK, year, 'PHEVTOTT']
# heavy
# EV.to_csv('JuliatoBALMOREL/EV'+scenario+str(adoption)+iter+'.csv')

# HP (categ,RRR,YYY,EH)
ehnumbers=pd.read_csv('JuliatoBALMOREL/ShareHP_categoriesfin_'+scenario+'.csv',header=[0,1,2],index_col=0)
EH=pd.DataFrame(index=range(8736),columns=range(800))
EH.columns=pd.MultiIndex.from_product([lista,['DK1','DK2'],years,['EH','AWTOT','AATOT','AASHTOT']]) # ('AA','AW','AASH')
for DK in ['DK1','DK2']:
    for year in years:
        AWprofile = pd.read_csv('JuliatoBALMOREL/Juliaout/FlexEH/FlexEHAW_' + DK + '_y' + year + '.csv', names=categories)
        AAprofile = pd.read_csv('JuliatoBALMOREL/Juliaout/FlexEH/FlexEHAA_' + DK + '_y' + year + '.csv', names=categories)
        AASHprofile = pd.read_csv('JuliatoBALMOREL/Juliaout/FlexEH/FlexEHSumHouse_' + DK + '_y' + year + '.csv',names=categories)
        for category in categories:
            weights=[ehnumbers.loc[category][DK, 'AW', year],ehnumbers.loc[category][DK, 'AA', year],ehnumbers.loc[category][DK, 'AASH', year]]
            #EH[category, DK, year, 'AA'] = (AWprofile[category] * weights[0] + AAprofile[category] * weights[1] +AASHprofile[category] * weights[2]) / sum(weights)  # change with share per each category
            #EH[category, DK, year, 'AW'] = (AWprofile[category] * weights[0] + AAprofile[category] * weights[1] +AASHprofile[category] * weights[2]) / sum(weights)  # change with share per each category
            #EH[category, DK, year, 'AASH'] = (AWprofile[category] * weights[0] + AAprofile[category] * weights[1] +AASHprofile[category] * weights[2]) / sum(weights)  # change with share per each category
            EH[category, DK, year, 'EH'] = (AWprofile[category] * weights[0] + AAprofile[category] * weights[1] +AASHprofile[category] * weights[2]) / sum(weights)  # change with share per each category
            EH[category, DK, year, 'AWTOT'] = AWprofile[category] * ehnumbers.loc[category][DK, 'AW', year]
            EH[category, DK, year, 'AATOT'] = AAprofile[category] * ehnumbers.loc[category][(DK, 'AA', year)]
            EH[category, DK, year, 'AASHTOT'] = AASHprofile[category] * ehnumbers.loc[category][(DK, 'AASH', year)]
        EH['TOT', DK, year, 'AWTOTT'] = EH[EH.columns[(EH.columns.get_level_values(1)==DK)&(EH.columns.get_level_values(2)==year) & (EH.columns.get_level_values(3)=='AWTOT')]].sum(axis=1)
        EH['TOT', DK, year, 'AATOTT'] = EH[EH.columns[(EH.columns.get_level_values(1)==DK)&(EH.columns.get_level_values(2)==year) & (EH.columns.get_level_values(3)=='AATOT')]].sum(axis=1)
        EH['TOT', DK, year, 'AASHTOTT'] = EH[EH.columns[(EH.columns.get_level_values(1)==DK)&(EH.columns.get_level_values(2)==year) & (EH.columns.get_level_values(3)=='AASHTOT')]].sum(axis=1)
        TOTTOT[DK,year,'EH']=EH['TOT', DK, year, 'AWTOTT']+EH['TOT', DK, year, 'AATOTT']+EH['TOT', DK, year, 'AASHTOTT']
# heavy
# EH.to_csv('JuliatoBALMOREL/EH'+scenario+str(adoption)+iter+'.csv')

# Yearly normalization
RESEVHP=pd.DataFrame()
ssn=[0,1680,3360,5040,8736]

RESEVHPDK1=pd.DataFrame()
RESEVHPDK2=pd.DataFrame()
for i in years:
    RESEVHPDK1[i]=pd.concat([TOTTOT['DK1',i,'EV'][ssn[p]:ssn[p+1]]-TOTTOT['DK1',i,'EV'][ssn[p]:ssn[p+1]].mean() for p in range(4)])+\
                    pd.concat([TOTTOT['DK1',i,'EH'][ssn[p]:ssn[p+1]]-TOTTOT['DK1',i,'EH'][ssn[p]:ssn[p+1]].mean() for p in range(4)])
    RESEVHPDK2[i]=pd.concat([TOTTOT['DK2',i,'EV'][ssn[p]:ssn[p+1]]-TOTTOT['DK2',i,'EV'][ssn[p]:ssn[p+1]].mean() for p in range(4)])+\
                    pd.concat([TOTTOT['DK2',i,'EH'][ssn[p]:ssn[p+1]]-TOTTOT['DK2',i,'EH'][ssn[p]:ssn[p+1]].mean() for p in range(4)])

RESEVHP['DK1']=pd.concat([RESEVHPDK1[i] for i in years])
RESEVHP['DK2']=pd.concat([RESEVHPDK2[i] for i in years])

# Adoption and to MW
# adoption= 0.40                 ### choose
RESEVHPOK=RESEVHP*adoption/1000

if scenario=='RapidScenario':
    RESEVHPOK['DK1']=RESEVHPOK['DK1']+TSRapDK1
    RESEVHPOK['DK2']=RESEVHPOK['DK2']+TSRapDK2
if scenario=='SlowScenario':
    RESEVHPOK['DK1']=RESEVHPOK['DK1']+TSSlowDK1
    RESEVHPOK['DK2']=RESEVHPOK['DK2']+TSSlowDK2

# Export
RESEVHPOK.to_csv('BALMOREL/RESEVHP files/xlsx/DE_RESEVHP'+scenario+str(adoption)+iter+'.csv') # normalized already but non adopted
RESEVHPOK.describe()

max_stars_repo_path: BALMOREL add-on/ResidEVHPaggr_step2.py

max_issues_repo_path: BALMOREL add-on/ResidEVHPaggr_step2.py

max_forks_repo_path: BALMOREL add-on/ResidEVHPaggr_step2.py


text:
# _*_ coding: utf-8 _*_

"""
Order number prediction of car parts based on time series method.

Author: StrongXGP (xgp1227@gmailcom)
Date:   2018/09/29
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, arma_order_select_ic

def is_stationary(ts, maxlag=15, autolag=None, regression='ct'):
    """Checks if the time series `ts` is stationary."""
    adf_result = adfuller(ts, maxlag=maxlag, autolag=autolag, regression=regression)
    return adf_result[0] < adf_result[4]['%5']

def d_param(ts, max_lag=12):
    """Finds out the differential order."""
    if is_stationary(ts):
        return 0
    else:
        for i in range(1, max_lag + 1):
            if is_stationary(ts.diff(i).dropna()):
                return i

def ARMA_params(ts, max_ar=4, max_ma=2, ic='aic'):
    """Finds out the AR order `p` and the ma order `q`."""
    if ic == 'aic':
        return arma_order_select_ic(ts.dropna(), max_ar=max_ar, max_ma=max_ma, ic=ic).aic_min_order
    elif ic == 'bic':
        return arma_order_select_ic(ts.dropna(), max_ar=max_ar, max_ma=max_ma, ic=ic).bic_min_order
    elif ic == 'hqic':
        return arma_order_select_ic(ts.dropna(), max_ar=max_ar, max_ma=max_ma, ic=ic).hqic_min_order

def main():
    """Process procedure."""

    # Load data
    data_path = "../data/"
    df = pd.read_csv()

max_stars_repo_name: Genpeng/auto-parts-time-series-forecast

max_issues_repo_name: Genpeng/auto-parts-time-series-forecast

max_forks_repo_name: Genpeng/auto-parts-time-series-forecast

reason: import numpy,from statsmodels


text:
import pandas as pd
import numpy as np
import os
import glob
import csv

#import shutil

def cleancsv(source):
    try:
        data = pd.read_csv(source, parse_dates=True)
        print('ok')
    except (FileNotFoundError, IOError):
        print('Wrong file or file path.')
        return
    if data.empty:
        print('data empty')
        return

    # shutil.copyfile(source, destination + source)
    # print('ok')

    data = data.drop_duplicates(subset='Date', keep='first')
    data['Date'] = pd.to_datetime(data['Date'])
    data.set_index('Date', drop=False, inplace=True)
    idx = pd.date_range(data.index.min(), data.index.max())
    indexed_data = data.reindex(index=idx, fill_value=np.nan)
    indexed_data = indexed_data.replace('0', np.nan)
    indexed_data = indexed_data.fillna(method='ffill')
    indexed_data = indexed_data.drop('Date', 1)
    print(indexed_data)
    indexed_data.to_csv(source, index_label='Date')
    print('job done')
    # #filename = destination + '/' + source
    # with open(destination, 'w') as f:
    #     indexed_data.to_csv(f, index_label='Date')

def calcopening(source):
    try:
        data = pd.read_csv(source, index_col=0, parse_dates=True)
    except (FileNotFoundError, IOError):
        print('Wrong file or file path.')
        return
    if data.empty:
        return

    data['Opening Price'] = data['Closing Price'].shift(1)
    # The Opening Price must be adjusted so that it is smaller than Maximum Price
    # and larger than Minimum Price
    data['Maximum Price'] = data[['Opening Price', 'Maximum Price', 'Minimum Price', 'Closing Price']].max(axis=1)
    data['Minimum Price'] = data[['Opening Price', 'Maximum Price', 'Minimum Price', 'Closing Price']].min(axis=1)
    data.set_value(data.index[0], 'Opening Price', data.get_value(data.index[0], 'Closing Price'))
    data.to_csv(source, index=True)
    print('complete')

def cleanall(source):

    os.chdir(source)
    for file in glob.glob("*.csv"):
        filename = os.path.basename(file)
        print('Cleaning ' + filename + '...\n')
        cleancsv(filename)

def applyfunc(func, source, *args, **kwargs):

   # os.chdir(source)
    for file in glob.glob("*.csv"):
        filename = os.path.basename(file)
        func(file, *args, **kwargs)

if __name__ == "__main__":
    cleanall('./data/')
    applyfunc(calcopening, './data/')

max_stars_repo_name: pudasainishushant/NEPSE-prediction

max_issues_repo_name: pudasainishushant/NEPSE-prediction

max_forks_repo_name: pudasainishushant/NEPSE-prediction


text:
# Created by Xingyu Lin, 11/12/2017

# coding: utf-8

import numpy as np
import matplotlib as mpl

mpl.use('Agg')
import matplotlib.pyplot as plt
from numpy import linalg as LA
import os
from os import path as osp
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib

def f(x, y):
  return np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)

def main():
  optimizers = ['L2L', 'Adam', 'Momentum', 'SGD', 'NAG', 'RMSProp']

  problem_path = './problems/quadratic.npz'
  npzfile = np.load(problem_path)
  problems_w, problems_b = npzfile['arr_0'], npzfile['arr_1']
  prob_num = len(problems_w)
  x = {}
  obj = {}

  for optimizer in optimizers:
    x[optimizer] = np.load(osp.join('./results', optimizer + '.npy'))

  for prob_idx in range(prob_num):
    plt.figure(figsize=(10, 6))
    for optimizer in optimizers:
      obj[optimizer] = list(
        map(lambda x: LA.norm(problems_w[prob_idx].dot(x) - np.transpose(problems_b[prob_idx])) ** 2, x[optimizer][prob_idx]))
      plt.plot(obj[optimizer], label=optimizer)
    plt.legend(loc='upper right')
    plt.xlabel('number of iterations')
    plt.ylabel('objective value')
    plt.savefig('./figs/loss_prob_{}.png'.format(prob_idx))

    # X = {}  # dictionary to store data
    # obj = {}  # dictionary to store obj value
    # W = 0.5 * np.eye(2)
    # Y = 0.5 * np.ones(2)
    # for name in names:
    #     X[name] = np.loadtxt(name + '.txt')
    #     obj[name] = list(map(lambda x: LA.norm(W.dot(x) - Y) ** 2, X[name]))
    #     plt.plot(obj[name], label=name)
    # plt.legend(loc='upper right')
    # plt.xlabel('number of iterations')
    # plt.ylabel('objective value')
    #

    # # plot level set and trajectory for SGD
    # W = 0.5 * np.eye(2)
    # Y = 0.5 * np.ones(2)
  optimizer = 'L2L'
  for prob_idx in range(prob_num):
    W = problems_w[prob_idx]
    Y = problems_b[prob_idx]
    A = W.transpose().dot(W)
    b = 2 * W.transpose().dot(Y)
    delta = 0.001
    minx = np.min(x[optimizer][prob_idx][:, 0])
    miny = np.min(x[optimizer][prob_idx][:, 1])
    maxx = np.max(x[optimizer][prob_idx][:, 0])
    maxy = np.max(x[optimizer][prob_idx][:, 1])
    min_plot = min(minx, miny)
    max_plot = max(maxx, maxy)
    t_min = min_plot - 0.3 * (max_plot - min_plot)
    t_max = max_plot + 0.3 * (max_plot - min_plot)
    x1 = np.arange(t_min, t_max, delta)
    x2 = np.arange(t_min, t_max, delta)
    X1, X2 = np.meshgrid(x1, x2)
    F = A[0, 0] * X1 ** 2 + A[1, 1] * X2 ** 2 + (A[0, 1] + A[1, 0]) * X1 * X2 - b[0] * X1 - b[1] * X2 + LA.norm(
      Y) ** 2
    plt.figure()
    plt.contourf(X1, X2, F, 100, cmap='RdGy')
    plt.colorbar()
    plt.axes().set_aspect('equal')
    plt.axis([t_min, t_max, t_min, t_max])

    plt.scatter(x[optimizer][prob_idx][:, 0], x[optimizer][prob_idx][:, 1], s=30, edgecolors='g', facecolors='none',
                marker='o')
    plt.plot(x[optimizer][prob_idx][:, 0], x[optimizer][prob_idx][:, 1], color='g')
    plt.savefig('./figs/contour_{}.png'.format(prob_idx))

if __name__ == "__main__":
  main()

reason: import numpy,from numpy


text:
#!/usr/bin/python

import rospy
import numpy as np
import math
import tf

from ackermann_msgs.msg import AckermannDriveStamped
from geometry_msgs.msg import PoseArray, Pose, PoseStamped
from std_msgs.msg import Bool

class GoToPointNode:
    targetPose = None
    drive = False
    reqed = False

    def __init__(self):
        # Default Settings
        self.k = 0.8
        self.topic_goal_in = "/move_base_simple/goal"
        self.topic_drive_out = "/vesc/ackermann_cmd_mux/input/nav"
        self.topic_plan_req = "goto_pt/plan_req"
        self.base_frame = "base_link"
        self.map_frame = "odom"

        # Initial Settings
        self.k_steer = 0.1
        self.kp = 0.6
        self.ki = 0
        self.kd = 0
        self.max_steering_angle = 0.3
        self.max_speed =.5
        self.lastDistance = 0
        self.lastTheta = 0
        self.distanceI = 0
        self.thetaI = 0
        self.threshold = 1
        self.rate=30

        # Param Settings
        self.k_steer = rospy.get_param('~k_steer', self.k_steer)
        self.kp = rospy.get_param('~k_dist', self.kp)
        self.topic_goal_in = rospy.get_param('~topic_goal_in', self.topic_goal_in)
        self.topic_plan_req = rospy.get_param('~topic_plan_req', self.topic_plan_req)
        self.topic_drive_out = rospy.get_param('~topic_drive_out', self.topic_drive_out)
        self.base_frame = rospy.get_param('~base_frame', self.base_frame)
        self.map_frame = rospy.get_param('~map_frame', self.map_frame)
        self.max_speed = rospy.get_param('~max_vel', self.max_speed)
        self.max_steering_angle = rospy.get_param('~max_theta', self.max_steering_angle)
        self.threshold = rospy.get_param('~threshold', self.threshold)
        self.rate = rospy.get_param('~rate', self.rate)

        # Pubs and Subs
        self.drive_pub = rospy.Publisher(self.topic_drive_out, AckermannDriveStamped, queue_size=1)
        self.req_pub = rospy.Publisher(self.topic_plan_req, Bool, queue_size=1)
        rospy.Subscriber(self.topic_goal_in, PoseStamped, self.new_dest_callback)

        self.listener = tf.TransformListener(True, rospy.Duration(10.0))

    def doDrive(self):
        # Get point and transform base_frame
        if(self.drive):
            try:
                self.targetPose.header.stamp = self.listener.getLatestCommonTime(self.base_frame,self.targetPose.header.frame_id)
                dest = self.listener.transformPose(self.base_frame, self.targetPose)
            except:
                print "TF Error!"
                return

            x = dest.pose.position.x
            y = dest.pose.position.y

            # Computer dist and theta
            distance = math.pow(math.pow(x,2) + math.pow(y,2),  0.5)
            theta = math.atan2(y,x)
            #rospy.loginfo("Distance: " + str(distance))

            if(distance < self.threshold * 2 and self.reqed == False):
                self.reqed = True
                self.req_pub.publish(True)

            if(distance < self.threshold):
                #rospy.loginfo("Goal Reached!")
                self.drive = False
                self.reqed = False
                self.req_pub.publish(True)
                return;

            # Prepare to Publish Drive Msg
            msg = AckermannDriveStamped()
            msg.header.stamp = rospy.Time.now()
            msg.header.frame_id = "base_link"

            # calc speed and theta
            # tested for points ahead/left, ahead/right, ahead, behind, and current loc

            self.thetaI=self.thetaI+theta
            tp= self.k_steer * theta
            ti= self.ki*(self.thetaI)
            td= self.kd* (theta-self.lastTheta)
            msg.drive.steering_angle= max(min(self.max_steering_angle,theta), -1*self.max_steering_angle)

            self.distanceI= self.distanceI+distance
            dp= self.kp * distance
            di= self.kd*(distance - self.lastDistance)
            dd= self.ki*(self.distanceI)
            self.lastDistance=distance
            msg.drive.speed= max(0.5,min(self.max_speed, dp+di+dd - 2*min(abs(theta), 0.3)))

            self.lastTheta=theta
            self.drive_pub.publish(msg)

    def new_dest_callback(self,data):

        self.targetPose = data;
        self.drive = True
        self.reqed = False

if __name__=="__main__":
    rospy.init_node("goto_pt")

    GotoPt = GoToPointNode()
    rate = rospy.Rate(GotoPt.rate)
    while not rospy.is_shutdown():
        GotoPt.doDrive();
        rate.sleep();

max_stars_repo_name: aybakana/ROS-Ogren

max_issues_repo_name: aybakana/ROS-Ogren

max_forks_repo_name: aybakana/ROS-Ogren


text:
import logging
import os

import pandas as pd
import numpy as np
import matplotlib

matplotlib.use('Agg')

from matplotlib import pyplot as plt
import matplotlib.colors as colors

from funfolding import binning, model, solution
from funfolding.visualization import visualize_classic_binning
from funfolding.visualization import visualize_tree_binning
from funfolding.visualization import visualize_model
from funfolding.visualization.visualize_llh import plot_llh_slice

import corner
from scipy import linalg

def generate_acceptance_correction(vec_f_truth,
                                   binning,
                                   logged_truth):
    e_min = 200
    e_max = 50000
    gamma = -2.7
    n_showers = 12000000
    if logged_truth:
        binning = np.power(10., binning)
    normalization = (gamma + 1) / (e_max ** (gamma + 1) - e_min ** (gamma + 1))
    corsika_cdf = lambda E: normalization * E ** (gamma + 1) / (gamma + 1)
    vec_acceptance = np.zeros_like(vec_f_truth, dtype=float)
    for i, vec_i_detected in enumerate(vec_f_truth):
        p_bin_i = corsika_cdf(binning[i + 1]) - corsika_cdf(binning[i])
        vec_acceptance[i] = p_bin_i * n_showers / vec_i_detected
    return vec_acceptance

if __name__ == '__main__':
    logging.captureWarnings(True)
    logging.basicConfig(
        format='%(processName)-10s %(name)s %(levelname)-8s %(message)s',
        level=logging.INFO)

    random_seed = 1340
    n_walkers = 100
    n_steps_used = 2000
    n_samples_test = 5000
    min_samples_leaf = 20
    tau = 1.
    binning_E = np.linspace(2.4, 4.2, 10)

    random_state = np.random.RandomState(random_seed)
    if not os.path.isfile('fact_simulations.hdf'):
        from get_fact_simulations import download
        logging.info('Downloading FACT simulations!')
        download()
    df = pd.read_hdf('fact_simulations.hdf', 'gamma_simulation')

    binned_E = np.digitize(df.loc[:, 'log10(energy)'],
                           binning_E)

    idx = np.arange(len(df))
    random_state.shuffle(idx)

    test_slice = slice(0, n_samples_test)
    binning_slice = slice(n_samples_test, n_samples_test + 10 * n_samples_test)
    A_slice = slice(n_samples_test + 10 * n_samples_test, None)

    idx_test = np.sort(idx[test_slice])
    idx_binning = np.sort(idx[binning_slice])
    idx_A = np.sort(idx[A_slice])

    binning_E = np.linspace(2.4, 4.2, 10)
    binned_E = np.digitize(df.loc[:, 'log10(energy)'], binning_E)
    binned_E_test = binned_E[idx_test]
    binned_E_binning = binned_E[idx_binning]
    binned_E_A = binned_E[idx_A]

    obs_array = df.get(['log10(ConcCore)', 'log10(E_RF)']).values

    X_A = obs_array[idx_A]
    X_binning = obs_array[idx_binning]
    X_test = obs_array[idx_test]

    classic_binning = binning.ClassicBinning(
        bins=[15, 25],
        random_state=random_state)

    classic_binning.fit(X_A)

    fig, ax = plt.subplots()
    visualize_classic_binning.plot_binning(ax,
                 classic_binning,
                 X_A,
                 log_c=False,
                 cmap='viridis')
    fig.savefig('05_fact_example_original_binning.png')

    tree_binning_2d = binning.TreeBinningSklearn(
        regression=False,
        min_samples_leaf=int(min_samples_leaf * 10.),
        random_state=random_state)

    tree_binning_2d.fit(X_binning,
                     binned_E_binning,
                     uniform=False)

    fig, ax = plt.subplots(figsize=(6, 6))
    visualize_tree_binning.plot_binning(ax,
                                        tree_binning_2d,
                                        limits=[-0.7,
                                                -0.2,
                                                2.7,
                                                4.2],
                                        X=X_A,
                                        linecolor='k',
                                        linewidth='0.3',
                                        log_c=False,
                                        as_hexbins=True,
                                        hex_kwargs={'rasterized': True,
                                                    'gridsize': 50},
                                        cmap='viridis')
    ax.set_ylabel('log10(Energy Estimator [GeV])')
    ax.set_xlabel('log10(Concentration [a.u.])')
    fig.savefig('05_fact_example_original_tree_binning.png', dpi=300)

    closest = classic_binning.merge(X_binning,
                                    min_samples=int(min_samples_leaf * 10.),
                                    mode='closest')
    fig, ax = plt.subplots()
    visualize_classic_binning.plot_binning(ax,
                 closest,
                 X_A,
                 log_c=False,
                 cmap='viridis')
    fig.savefig('05_fact_example_original_binning_closest.png')

    unmerged_model = model.LinearModel()
    binned_g_A = classic_binning.digitize(X_A)
    unmerged_model.initialize(digitized_obs=binned_g_A,
                              digitized_truth=binned_E_A)

    binned_g_A = closest.digitize(X_A)
    merged_model = model.LinearModel()
    merged_model.initialize(digitized_obs=binned_g_A,
                            digitized_truth=binned_E_A)

    single_obs_model = model.LinearModel()
    max_e = np.max(X_A[:, 1]) + 1e-3
    min_e = np.min(X_A[:, 1]) - 1e-3
    binning_E_obs = np.linspace(min_e, max_e, 11)
    binned_g_A = np.digitize(X_A[:, 1], binning_E_obs)
    single_obs_model.initialize(digitized_obs=binned_g_A,
                                digitized_truth=binned_E_A)

    single_obs_model_more_bins = model.LinearModel()
    max_e = np.max(X_A[:, 1]) + 1e-3
    min_e = np.min(X_A[:, 1]) - 1e-3
    binning_E_obs = np.linspace(min_e, max_e, closest.n_bins + 1)
    binned_g_A = np.digitize(X_A[:, 1], binning_E_obs)
    single_obs_model_more_bins.initialize(digitized_obs=binned_g_A,
                                          digitized_truth=binned_E_A)
    fig, ax = plt.subplots(figsize=(2, 6))
    visualize_model.plot_A(ax, merged_model)
    fig.savefig('05_A_single_obs_model.png', dpi=300)

    binned_g_A = tree_binning_2d.digitize(X_A)
    tree_2d_model = model.LinearModel()
    tree_2d_model.initialize(digitized_obs=binned_g_A,
                            digitized_truth=binned_E_A)
    fig, ax = plt.subplots(figsize=(2, 6))
    visualize_model.plot_A(ax, tree_2d_model)
    fig.savefig('05_A_tree_model.png', dpi=300)

    tree_obs = ["log10(E_RF)",
                "log10(Size)",
                "log10(ConcCore)",
                "log10(numPixelInShower)",
                "log10(Length)",
                "Width",
                "M3Trans",
                "M3Long",
                "m3l",
                "m3t",
                "Concentration_onePixel",
                "Concentration_twoPixel",
                "Leakage",
                "Leakage2",
                "concCOG",
                "numIslands",
                "phChargeShower_mean",
                "phChargeShower_variance",
                "phChargeShower_max"]

    obs_array = df.get(tree_obs).values

    X_tree_test = obs_array[idx_test]
    X_tree_binning = obs_array[idx_binning]
    X_tree_A = obs_array[idx_A]

    tree_binning_uniform = binning.TreeBinningSklearn(
        regression=False,
        min_samples_leaf=int(min_samples_leaf * 10.),
        random_state=random_state)

    tree_binning_uniform.fit(X_tree_binning,
                             binned_E_binning,
                             uniform=True)

    binned_g_A = tree_binning_uniform.digitize(X_tree_A)

    tree_model_uniform = model.LinearModel()
    tree_model_uniform.initialize(digitized_obs=binned_g_A,
                                  digitized_truth=binned_E_A)
    fig, ax = plt.subplots(figsize=(2, 6))
    visualize_model.plot_A(ax, tree_model_uniform)
    fig.savefig('05_A_tree_model_full_uniform.png', dpi=300)

    tree_binning = binning.TreeBinningSklearn(
        regression=False,
        min_samples_leaf=int(min_samples_leaf * 10.),
        random_state=random_state)

    tree_binning.fit(X_tree_binning,
                     binned_E_binning,
                     uniform=False)

    binned_g_A = tree_binning.digitize(X_tree_A)

    tree_model = model.LinearModel()
    tree_model.initialize(digitized_obs=binned_g_A,
                          digitized_truth=binned_E_A)
    visualize_model.plot_A(ax, tree_model)
    fig.savefig('05_A_tree_model_full.png', dpi=300)

    fig, ax = plt.subplots()
    svd_values = unmerged_model.evaluate_condition()
    bin_edges = np.linspace(0, len(svd_values), len(svd_values) + 1)
    bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2.
    ax.hist(bin_centers,
            bins=bin_edges,
            weights=svd_values,
            histtype='step',
            label='2 Observables (Unmerged; {} Bins)'.format(
                classic_binning.n_bins))

    svd_values = merged_model.evaluate_condition()
    ax.hist(bin_centers,
            bins=bin_edges,
            weights=svd_values,
            histtype='step',
            label='2 Observables (Merged; {} Bins)'.format(closest.n_bins))
#
#    svd_values = single_obs_model.evaluate_condition()
#    ax.hist(bin_centers,
#            bins=bin_edges,
#            weights=svd_values,
#            histtype='step',
#            label='Single Observable ({} Bins)'.format(closest.n_bins))
#
#    svd_values = single_obs_model_more_bins.evaluate_condition()
#    ax.hist(bin_centers,
#            bins=bin_edges,
#            weights=svd_values,
#            histtype='step',
#            label='Single Observable ({} Bins)'.format(closest.n_bins))

    svd_values = tree_2d_model.evaluate_condition()
    ax.hist(bin_centers,
            bins=bin_edges,
            weights=svd_values,
            histtype='step',
            label='Tree Based 2D ({} Bins)'.format(tree_binning.n_bins))

    svd_values = tree_model.evaluate_condition()
    ax.hist(bin_centers,
            bins=bin_edges,
            weights=svd_values,
            histtype='step',
            label='Tree Based ({} Bins)'.format(tree_binning.n_bins))
#
    svd_values = tree_model_uniform.evaluate_condition()
    ax.hist(bin_centers,
            bins=bin_edges,
            ooweights=svd_values,
            histtype='step',
            label='Tree Based ({} Bins; Uniform)'.format(tree_binning.n_bins))

    plt.legend(loc='lower left')
    ax.set_yscale("log", nonposy='clip')
    plt.savefig('05_condition.png')

    exit()

    binned_g_test = tree_binning.digitize(X_tree_test)
    vec_g, vec_f = tree_model.generate_vectors(binned_g_test,
                                               binned_E_test)

    vec_f_truth = np.array(np.bincount(binned_E), dtype=float)[1:]
    vec_acceptance = generate_acceptance_correction(vec_f_truth,
                                                    binning_E,
                                                    True)
    vec_f_truth /= np.sum(vec_f_truth)
    vec_f_truth *= n_samples_test

    llh = solution.StandardLLH(tau=tau,
                               log_f=True,
                               vec_acceptance=vec_acceptance,
                               C='thikonov')
    llh.initialize(vec_g=vec_g,
                   model=tree_model)

    sol_gd = solution.LLHSolutionGradientDescent(n_steps=500,
                                                 gamma=0.01)
    sol_gd.initialize(llh=llh, model=tree_model)
    sol_gd.set_x0_and_bounds()
    x, llh_values, gradient, hessian = sol_gd.fit()
    idx_best = np.argmax(llh_values)
    vec_f_str = ', '.join('{0:.2f}'.format(a)
                          for a in x[idx_best])
    logging.info('Best Fit (Gradient):\t{}\t(LLH: {})'.format(
        vec_f_str,
        llh_values[idx_best]))

    sol_mini = solution.LLHSolutionMinimizer()
    sol_mini.initialize(llh=llh, model=tree_model)
    sol_mini.set_x0_and_bounds(x0=x[idx_best])
    best_fit, mini_cov = sol_mini.fit(constrain_N=False)

    vec_f_str = ', '.join('{0:.2f}'.format(a)
                          for a in best_fit.x)
    logging.info('Best Fit (Minimizer):\t{}\t(LLH: {})'.format(
        vec_f_str,
        best_fit.fun))

    sol_mcmc = solution.LLHSolutionMCMC(n_burn_steps=100,
                                        n_used_steps=n_steps_used,
                                        n_walkers=n_walkers,
                                        n_threads=1,
                                        random_state=random_state)
    sol_mcmc.initialize(llh=llh, model=tree_model)
    sol_mcmc.set_x0_and_bounds(x0=best_fit.x)
    vec_f_est_mcmc, sigma_vec_f, sample, probs = sol_mcmc.fit()

    vec_f_str = ', '.join('{0:.2f}'.format(a)
                          for a in vec_f_est_mcmc)
    logging.info('Best Fit (MCMC):\t{}\t(LLH: {})'.format(
        vec_f_str,
        max(probs)))

    # sol_mcmc.n_threads = 9
    # logging.info('Calculating Eff sample size:')
    # n_eff = sol_mcmc.calc_effective_sample_size(sample, n_threads=9)
    # n_eff_str = ', '.join(str(n) for n in n_eff)
    # logging.info('per Walker:\t{} ({} Walker with {} steps)'.format(
    #     n_eff_str,
    #     n_walkers,
    #     n_steps_used))

    def create_llh_slice(llh, best_fit, selected_bin=None):
        if selected_bin is None:
            selected_bin = np.argmax(best_fit)
        points = np.linspace(0.9 * best_fit[selected_bin],
                             1.1 * best_fit[selected_bin],
                             31)
        llh_values = np.zeros_like(points)
        gradient_values = np.zeros_like(points)
        hessian_values = np.zeros_like(points)

        fig, [ax_grad, ax_hess] = plt.subplots(2, 1, figsize=(24, 18))
        diff = np.diff(points)[0] / 1.5
        for i, p_i in enumerate(points):
            best_fit[selected_bin] = p_i
            llh_values[i] = llh.evaluate_llh(best_fit)
            gradient_values[i] = llh.evaluate_gradient(best_fit)[selected_bin]
            hessian_values[i] = llh.evaluate_hessian(best_fit)[selected_bin,
                                                               selected_bin]
            lower_x = p_i - diff
            upper_x = p_i + diff

            grad_lower_y = llh_values[i] - (diff * gradient_values[i])
            grad_upper_y = llh_values[i] + (diff * gradient_values[i])

            hess_lower_y = gradient_values[i] - (diff * hessian_values[i])
            hess_upper_y = gradient_values[i] + (diff * hessian_values[i])

            if gradient_values[i] < 0:
                direction = -1.
            else:
                direction = 1.

            ax_hess.plot([lower_x, upper_x],
                         [hess_lower_y, hess_upper_y],
                         'k-')

        dy = gradient_values * diff
        dx = np.ones_like(points) * diff
        dx[gradient_values < 0] *= -1.
        dy[gradient_values < 0] *= -1.

        ax_grad.quiver(points, llh_values, dx, dy, angles='xy', scale_units='xy', scale=1.)
        ax_grad.plot(best_fit[selected_bin], llh_values[selected_bin], 'ro')
        ax_hess.plot(points, gradient_values, 'o')
        fig.savefig('05_llh_scan.png')
        plt.close(fig)
        return selected_bin

    logging.info('Creating plot of a LLH slice')
    fig = plot_llh_slice(llh, best_fit.x)
    fig.savefig('05_llh_slice.png')

    logging.info('Creating corner plot')
    corner_fig = corner.corner(sample,
                               truths=vec_f_est_mcmc,
                               truth_color='r',
                               rasterized=True)
    corner_fig.savefig('05_corner_fact.png')

    logging.info('Creating best fit plots')
    fig, ax = plt.subplots(1, 1, figsize=(6, 4))
    bin_mids = (binning_E[1:] + binning_E[:-1]) / 2.
    bin_width = (binning_E[1:] - binning_E[:-1]) / 2.
    plt.hist(bin_mids, bins=binning_E, weights=vec_f_truth, histtype='step')
    ax.errorbar(bin_mids,
                best_fit.x,
                yerr=np.sqrt(np.diag(np.absolute(mini_cov))),
                xerr=bin_width,
                ls="",
                color="k",
                label="Unfolding (Minimizer)")
    ax.set_yscale("log", nonposy='clip')
    ax.set_ylim([2e1, 2e3])
    fig.savefig('05_unfolding_minimizer.png')
    plt.close(fig)

    fig, ax = plt.subplots(1, 1, figsize=(6, 4))
    plt.hist(bin_mids, bins=binning_E, weights=vec_f_truth, histtype='step')
    ax.errorbar(bin_mids,
                vec_f_est_mcmc,
                yerr=[vec_f_est_mcmc - sigma_vec_f[0, :],
                      sigma_vec_f[1, :] - vec_f_est_mcmc],
                xerr=bin_width,
                ls="",
                color="r",
                label="Unfolding (MCMC)")
    ax.set_yscale("log", nonposy='clip')
    ax.set_ylim([2e1, 2e3])
    fig.savefig('05_unfolding_mcmc.png')
    plt.close(fig)

    logging.info('Creating LLH histogram')
    fig, ax = plt.subplots(1, 1, figsize=(6, 4.5))
    ax.hist(2*(np.max(probs) - probs),
            bins=50,
            weights=np.ones_like(probs) * 1./len(probs),
            histtype='step', lw=2)
    ax.set_yscale("log", nonposy='clip')
    ax.set_xlabel(r'$-2\cdot\ln\left(\frac{\mathdefault{LLH}}{\mathdefault{LLH}_{\mathdefault{Best Fit}}}\right)$')
    ax.set_ylabel(r'$\frac{\mathdefault{Bin}_i}{\sum_i \mathdefault{Bin}_i}$')
    plt.tight_layout()
    plt.savefig('05_hist_probs.png')

max_stars_repo_name: mbrner/funfolding

max_issues_repo_name: tudo-astroparticlephysics/funfolding

max_forks_repo_name: tudo-astroparticlephysics/funfolding

reason: import numpy,from scipy


text:
import argparse
from collections import namedtuple
import glob
import numpy as np
import os
import re
import struct
import json
import pandas as pd
from collections import OrderedDict
from tqdm.auto import tqdm

class Tags(object):
    """
    Represents tag name to id mapping from a header file.
    Preprocessor macros are held in the structure:
        tag_names: Dict<str, int>
    We also reverse index for faster subsequent parsing:
        tag_ids: Dict<int, str>
    """
    def __init__(self, fn: str):
        define = re.compile(r'^\s*#define\s+LOGGING_TAG_(\w+)\s+\(?(\w*?)\)?\s*$')
        self.tag_names = dict()
        with open(fn, 'r') as f:
            for line in f:
                match = define.match(line)
                if match is None:
                    continue
                try:
                    name, value = match.groups()
                    int_value = int(value, 0)
                    self.tag_names[name.lower()] = int_value
                except:
                    raise SyntaxError(line)
        self.tag_ids = {v: k for k, v in self.tag_names.items()}

    def __str__(self):
        return f'{{tag_names: {self.tag_names}, tag_ids: {self.tag_ids}}}'

class Log(object):
    """
    Represents a single log file.
        log: List<Tuple<int, int>>
    Also creates a human-readable version with tag names instead of IDs.
    This could hit performance issues for large log files.
    Preferably use IDs instead.
        log_h: List<Tuple<string, int>>
    """
    def __init__(self, fn: str, tags: Tags, verbose=0):
        self.fn = fn
        self.tags = tags

        if verbose > 0:
            print(f'Parsing {self.fn}...', end='')
        self.parse_log()
        if verbose > 0:
            print(' done')
        # print('log', self.log)

        self.make_human_readable()
        self.check_version()
        self.extract_constants()

    def make_human_readable(self):
        ''' Substitute tags for tag IDs to create a human-readable log '''
        assert all([tag_id in self.tags.tag_ids for (tag_id, value) in self.log])
        self.log_h = [(self.tags.tag_ids[tag_id], value) for (tag_id, value) in self.log]
        # for tag, value in self.log_h:
        #     print(f'{tag}: {value}')

    def check_version(self):
        """
        Assert logged version matches version from header file
        """
        try:
            self.version = next(value for (tag, value) in self.log_h if tag == 'version')
        except StopIteration:
            raise ValueError('no version tag found')
        if self.version is not self.tags.tag_names['version']:
            raise ValueError(f"expected version {self.tags.tag_names['version']}, got {self.version}")

    def extract_constants(self):
        self.clocks_per_sec = next(value for (tag, value) in self.log_h if tag == 'clocks_per_sec')

    def parse_log(self):
        with open(self.fn, 'rb') as f:
            d = f.read()
        self.log = list()
        offset = 0
        while (offset < len(d)):
            fmt = 'II'
            tag, value = struct.unpack_from(fmt, d, offset)
            offset += struct.calcsize(fmt)
            self.log.append((tag, value))

    def get_values(self, name: str):
        full_name = f'val_{name}'
        val_tag_id = self.tags.tag_names[full_name]
        values = [value for (tag_id, value) in self.log if tag_id == val_tag_id]
        return values

    def get_wall_clock_durations(self, name: str):
        ''' Automatically appends the wc_ prefix and _begin or _end suffixes'''
        begin_full_name = f'wc_{name}_begin'
        end_full_name = f'wc_{name}_end'

        begin_tag_id = self.tags.tag_names[begin_full_name]
        end_tag_id = self.tags.tag_names[end_full_name]

        begins = [value for (tag_id, value) in self.log if tag_id == begin_tag_id]
        ends = [value for (tag_id, value) in self.log if tag_id == end_tag_id]

        assert len(begins) == len(ends)
        # Some durations may wrap around, but check whether all fit in half-range
        diffs = [(ends[i] - begins[i]) % (1 << 32) for i in range(len(begins))]
        assert all([d < (1 << 31) for d in diffs])

        durations = [float(d) / 1e6 * 1e3 for d in diffs]
        return durations

    def get_wall_clock_tuples(self, name: str):
        begin_full_name = f'wc_{name}_begin'
        end_full_name = f'wc_{name}_end'

        begin_tag_id = self.tags.tag_names[begin_full_name]
        end_tag_id = self.tags.tag_names[end_full_name]

        begins = [value for (tag_id, value) in self.log if tag_id == begin_tag_id]
        ends = [value for (tag_id, value) in self.log if tag_id == end_tag_id]

        assert len(begins) == len(ends)
        return begins, ends

    def get_epoch_comp_comm_dataframe(self, migration_period: int):
        begin_tag_id = self.tags.tag_names['wc_epoch_begin']
        end_tag_id = self.tags.tag_names['wc_epoch_end']

        begins = [value for (tag_id, value) in self.log if tag_id == begin_tag_id]
        ends = [value for (tag_id, value) in self.log if tag_id == end_tag_id]

        assert len(begins) == len(ends)
        periods = range(0, len(begins)-migration_period, migration_period)
        compute_times = list()
        communication_times = list()
        for i in periods:
            compute_times.append((ends[i+migration_period-1] - begins[i]) / 1e3)
            communication_times.append((begins[i+migration_period] - ends[i+migration_period-1]) / 1e3)

        return list(zip(range(len(periods)), compute_times, communication_times))

    def get_cpu_clock_durations(self, name: str):
        ''' Automatically appends the cc_ prefix and _begin or _end suffixes'''
        begin_full_name = f'cc_{name}_begin'
        end_full_name = f'cc_{name}_end'

        begin_tag_id = self.tags.tag_names[begin_full_name]
        end_tag_id = self.tags.tag_names[end_full_name]

        begins = [value for (tag_id, value) in self.log if tag_id == begin_tag_id]
        ends = [value for (tag_id, value) in self.log if tag_id == end_tag_id]

        assert len(begins) == len(ends)
        # Some durations may wrap around, but check whether all fit in half-range
        diffs = [(ends[i] - begins[i]) % (1 << 32) for i in range(len(begins))]
        assert all([d < (1 << 31) for d in diffs])

        durations = [float(d) / self.clocks_per_sec * 1e3 for d in diffs]
        return durations

class Epochs(object):
    """
    Represents epochs across one or more logs.
    """
    def __init__(self, log: Log, tags: Tags):
        self.log = log
        self.tags = tags
        self.process_epochs_in_log()

    # class Epoch(object):
    #     self.begin = None
    #     self.end = None
    #     self.best_fitness = None

    # Old depreviated verison
    def process_epochs_in_log_legacy(self) -> None:
        Epoch = namedtuple('Epoch', ['begin', 'end', 'fitness'])

        # get relevant tag IDs for fast lookups
        begin_tag_id = self.tags.tag_names['wc_epoch_begin']
        end_tag_id = self.tags.tag_names['wc_epoch_end']
        fitness_tag_id = self.tags.tag_names['val_best_fitness']
        tag_ids = (begin_tag_id, end_tag_id, fitness_tag_id)

        it = iter(self.log.log)
        self.epochs: [Epoch] = []
        try:
            while True:
                begin = next(value for (tag_id, value) in it if tag_id == begin_tag_id)
                fitness = next(value for (tag_id, value) in it if tag_id == fitness_tag_id)
                end = next(value for (tag_id, value) in it if tag_id == end_tag_id)
                self.epochs.append(Epoch(
                    begin=float(begin)/1e3,
                    end=float(end)/1e3,
                    fitness=fitness
                ))
        except StopIteration:
            pass

    def process_epochs_in_log(self) -> None:
        Epoch = namedtuple('Epoch', ['end', 'fitness'])

        # get relevant tag IDs for fast lookups
        end_tag_id = self.tags.tag_names['wc_epoch_end']
        fitness_tag_id = self.tags.tag_names['val_best_fitness']

        it = iter(self.log.log)
        self.epochs: [Epoch] = []
        try:
            while True:
                fitness = next(value for (tag_id, value) in it if tag_id == fitness_tag_id)
                end = next(value for (tag_id, value) in it if tag_id == end_tag_id)
                self.epochs.append(Epoch(
                    end=float(end)/1e3,
                    fitness=fitness
                ))
        except StopIteration:
            pass

    def get_fitness_vs_time(self) -> ([int], [float]):
        return (
            [e.fitness for e in self.epochs],
            [e.end for e in self.epochs]
        )

    def get_fitness_vs_time_dataframe(self):
        return ((e.fitness, e.end, epoch) for epoch, e in enumerate(self.epochs))

def last_log():
    log_fns = glob.glob(os.path.join('logs', '*_tags.bin'))
    if len(log_fns) == 0:
        raise FileNotFoundError('No log files found')
    #print(log_fns)
    last_log_fn = max(log_fns)
    #print(last_log_fn)
    return last_log_fn

def logs_in_dir(path):
    log_fns = glob.glob(os.path.join(path, '*_tags.bin'))
    return log_fns

# Given a log_dir (generated with a leonhard run) and a name, saves a dataframe to name.gz
# That dataframe contains the epochs, wall clock times, fitness, rep, rank and all the variable parameters
def generate_fitness_wc_dataframe(log_dir, name, tag_loc="tags.hpp"):
    from tqdm import tqdm

    assert isinstance(name, str), "name must be a string"
    assert os.path.isdir(log_dir), "log_dir must be a directory"
    assert os.path.isfile(tag_loc), f"Could not find tag_loc {tag_loc}"

    # Extract tags
    tags = Tags("tags.hpp")

    all_names = os.listdir(log_dir)

    # Validate JSON
    json_file = list(filter(lambda x: ".json" in x, all_names))
    assert len(json_file) > 0, "Found no JSON file in the directory"
    assert len(json_file) <= 1, "Found multiple JSON files in the directory"
    json_file = json_file[0]

    # Extract repetitions
    with open(os.path.join(log_dir, json_file)) as file:
        json_file = json.load(file, object_pairs_hook=OrderedDict)
        repetitions = json_file["repetitions"]

    # Find unique runs (without repetitions)
    all_names = list(filter(lambda x: os.path.isdir(os.path.join(log_dir, x)), all_names))
    unique_names = list(set(map(lambda x: "_".join(x.split("_")[:-1]), all_names)))

    df = None
    for run_name in tqdm(unique_names):
        params = run_name.split("_")

        # Get parameter names
        param_names = list()
        for key, val in json_file["variable_params"].items():
            # Handle tuples
            if val["type"] == "tuple":
                param_names.extend(val["names"])
            else:
                param_names.append(key)
        param_names = list(map(lambda x: x.replace("-", ""), param_names))

        for repetition in tqdm(range(repetitions)):
            folder_name = run_name + "_" + str(repetition)
            folder_contents = os.listdir(os.path.join(log_dir, folder_name))
            folder_contents = list(filter(lambda x: ".bin" in x, folder_contents))
            try:
                for filename in folder_contents:
                    log = Log(os.path.join(log_dir, folder_name, filename), tags)
                    rank = int(filename.split("_")[-2])
                    epochs = Epochs(log, tags)
                    if df is None:
                        df = pd.DataFrame(epochs.get_fitness_vs_time_dataframe(), columns=["fitness", "wall clock time", "epoch"])
                        df["rank"] = rank
                        df["rep"] = repetition
                        for param, param_name in zip(params, param_names):
                            df[param_name] = param
                    else:
                        df2 = pd.DataFrame(epochs.get_fitness_vs_time_dataframe(), columns=["fitness", "wall clock time", "epoch"])
                        df2["rank"] = rank
                        df2["rep"] = repetition
                        for param, param_name in zip(params, param_names):
                            df2[param_name] = param
                        df = df.append(df2, ignore_index=True)
            except ValueError as e:
                print(f"Run {folder_name} probably crashed... ignoring")

    # Figure out correct file ending
    if name.endswith(".gz"):
        file_name = name
    else:
        file_name = name + ".gz"

    # Save to disk
    df.to_csv(file_name, compression="gzip", index=False)
    return df

def generate_periods_dataframe(log_dir, name, tag_loc="tags.hpp"):
    from itertools import count

    assert isinstance(name, str), "name must be a string"
    assert os.path.isdir(log_dir), "log_dir must be a directory"
    assert os.path.isfile(tag_loc), f"Could not find tag_loc {tag_loc}"

    # Extract tags
    tags = Tags("tags.hpp")

    all_names = os.listdir(log_dir)

    # Validate JSON
    json_file = list(filter(lambda x: ".json" in x, all_names))
    assert len(json_file) > 0, "Found no JSON file in the directory"
    assert len(json_file) <= 1, "Found multiple JSON files in the directory"
    json_file = json_file[0]

    # Get parameters from JSON
    with open(os.path.join(log_dir, json_file)) as file:
        json_file = json.load(file, object_pairs_hook=OrderedDict)

        # Extract repetitions
        repetitions = json_file["repetitions"]

    # Find unique runs (without repetitions)
    all_names = list(filter(lambda x: os.path.isdir(os.path.join(log_dir, x)), all_names))
    unique_names = list(set(map(lambda x: "_".join(x.split("_")[:-1]), all_names)))

    df = None
    for run_name in tqdm(unique_names):
        params = run_name.split("_")

        # Get parameter names
        param_names = list()
        for key, val in json_file["variable_params"].items():
            # Handle tuples
            if val["type"] == "tuple":
                param_names.extend(val["names"])
            else:
                param_names.append(key)
        param_names = list(map(lambda x: x.replace("-", ""), param_names))

        for repetition in tqdm(range(repetitions)):
            folder_name = run_name + "_" + str(repetition)
            folder_contents = os.listdir(os.path.join(log_dir, folder_name))
            folder_contents = list(filter(lambda x: ".bin" in x, folder_contents))
            try:
                for filename in folder_contents:
                    log = Log(os.path.join(log_dir, folder_name, filename), tags)
                    rank = int(filename.split("_")[-2])

                    val_tag_id = tags.tag_names["val_computation"]
                    comp_values = [value for (tag_id, value) in log.log if tag_id == val_tag_id]
                    val_tag_id = tags.tag_names["val_communication"]
                    comm_values = [value for (tag_id, value) in log.log if tag_id == val_tag_id]

                    if df is None:
                        df = pd.DataFrame(zip(count(), comp_values, comm_values), columns=["period", "computation time", "communication time"])
                        df["rank"] = rank
                        df["rep"] = repetition
                        for param, param_name in zip(params, param_names):
                            df[param_name] = param
                    else:
                        df2 = pd.DataFrame(zip(count(), comp_values, comm_values), columns=["period", "computation time", "communication time"])
                        df2["rank"] = rank
                        df2["rep"] = repetition
                        for param, param_name in zip(params, param_names):
                            df2[param_name] = param
                        df = df.append(df2, ignore_index=True)
            except ValueError as e:
                print(f"Run {folder_name} probably crashed... ignoring")

    comp_df = df.drop(columns="communication time")
    comp_df = comp_df.rename(columns={"computation time" : "time"})
    comp_df["type"] = "computation"
    comm_df = df.drop(columns="computation time")
    comm_df = comm_df.rename(columns={"communication time" : "time"})
    comm_df["type"] = "communication"
    new_df = comp_df.append(comm_df, ignore_index=True)

    # Figure out correct file ending
    if name.endswith(".gz"):
        file_name = name + "_periods"
    else:
        file_name = name + "_periods.gz"

    # Save to disk
    new_df.to_csv(file_name, compression="gzip", index=False)

    return new_df

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='Process GA logs')
    parser.add_argument('--tags', help='path to C header file with log defines')
    parser.add_argument('--log', help='path to log file')
    parser.add_argument('--dir', help='path to directory with log files')
    args = parser.parse_args()

    # Parse header file with tag definitions
    tags_fn = os.path.join('logging', 'tags.hpp') if not args.tags else args.tags
    tags = Tags(tags_fn)
    # print(tags)

    # parse binary log
    log_fns = (
        [args.log] if args.log else
        logs_in_dir(args.dir) if args.dir else
        [last_log()]
    )
    logs = [Log(log_fn, tags) for log_fn in log_fns]

    epochss = [Epochs(log, tags) for log in logs]
    from matplotlib import pyplot as plt
    for epochs in epochss:
        fitness, time = epochs.get_fitness_vs_time()
        plt.plot([t/1e3 for t in time], fitness)
    plt.xlabel('time [s]')
    plt.ylabel('distance')
    plt.show()

    for log in logs:
        print(log.fn)
        print('    total:',
            f"wall clock {log.get_wall_clock_durations('logging')[0]:.3f}ms",
            f"CPU clock {log.get_cpu_clock_durations('logging')[0]:.3f}ms"
        )
        def print_stats(name):
            ri_arr = log.get_wall_clock_durations(name)
            print(f'    {name}:',
                f"length {len(ri_arr)},",
                f"mean {np.mean(ri_arr):.3f}ms,",
                f"std {np.std(ri_arr):.3f}ms"
        )
        print_stats('epoch')
        print_stats('rank_individuals')
        print_stats('breed_population')
        print_stats('mutate_population')