text:
# -*- coding: utf-8 -*-
"""
Created on Tue Mar 16 11:13:37 2021

@author: dv516
"""

import numpy as np
import pickle

import pyro
pyro.enable_validation(True)  # can help with debugging
pyro.set_rng_seed(1)

from algorithms.PyBobyqa_wrapped.Wrapper_for_pybobyqa import PyBobyqaWrapper
from algorithms.Bayesian_opt_Pyro.utilities_full import BayesOpt
from algorithms.nesterov_random.nesterov_random import nesterov_random
from algorithms.simplex.simplex_method import simplex_method
from algorithms.CUATRO.CUATRO import CUATRO
from algorithms.Finite_differences.Finite_differences import finite_Diff_Newton
from algorithms.Finite_differences.Finite_differences import Adam_optimizer
from algorithms.Finite_differences.Finite_differences import BFGS_optimizer
from algorithms.SQSnobfit_wrapped.Wrapper_for_SQSnobfit import SQSnobFitWrapper
from algorithms.DIRECT_wrapped.Wrapper_for_Direct import DIRECTWrapper

from test_functions import quadratic_constrained

import matplotlib.pyplot as plt

def trust_fig(oracle, bounds):
    N = 200
    lim = 2
    x = np.linspace(-lim, lim, N)
    y = np.linspace(-lim, lim, N)
    X,Y = np.meshgrid(x, y)
    Z = oracle.sample_obj(X,Y)
    constr = oracle.sample_constr(X,Y)

    levels_list = np.logspace(-2, 1.5, 10)

    fig = plt.figure(figsize = (6,4))
    ax = fig.add_subplot()

    ax.contour(X,Y,Z*constr, levels = levels_list)
    ax.plot([bounds[0,0], bounds[0, 1]], [bounds[1,0], bounds[1, 0]], c = 'k')
    ax.plot([bounds[0,0], bounds[0, 1]], [bounds[1,1], bounds[1, 1]], c = 'k')
    ax.plot([bounds[0,0], bounds[0, 0]], [bounds[1,0], bounds[1, 1]], c = 'k')
    ax.plot([bounds[0,1], bounds[0, 1]], [bounds[1,0], bounds[1, 1]], c = 'k')

    return ax, fig

def average_from_list(solutions_list):
    N = len(solutions_list)
    f_best_all = np.zeros((N, 100))
    for i in range(N):
        f_best = np.array(solutions_list[i]['f_best_so_far'])
        x_ind = np.array(solutions_list[i]['samples_at_iteration'])
        for j in range(100):
            ind = np.where(x_ind <= j+1)
            if len(ind[0]) == 0:
                f_best_all[i, j] = f_best[0]
            else:
                f_best_all[i, j] = f_best[ind][-1]
    f_median = np.median(f_best_all, axis = 0)
    # f_av = np.average(f_best_all, axis = 0)
    # f_std = np.std(f_best_all, axis = 0)
    f_min = np.min(f_best_all, axis = 0)
    f_max = np.max(f_best_all, axis = 0)
    return f_best_all, f_median, f_min, f_max

class RB:
    def __init__(self, objective, ineq = []):
        self.obj = objective ; self.ieq = ineq
    def sample_obj(self, x, y):
        return self.obj(x, y)
    def sample_constr(self, x, y):
        if self.ieq == []:
            if (type(x) == float) or (type(x) == int):
                return 1
            else:
                return np.ones(len(x))
        elif (type(x) == float) or (type(x) == int):
            temporary = [int(g(x, y)) for g in self.ieq]
            return np.product(np.array(temporary))
        else:
            temporary = [g(x, y).astype(int) for g in self.ieq]
            return np.product(np.array(temporary), axis = 0)

def fix_starting_points(complete_list, x0, init_out, only_starting_point = False):
    if only_starting_point:
        for i in range(len(complete_list)):
            dict_out = complete_list[i]
            f_arr = dict_out['f_best_so_far']
            N_eval = len(f_arr)
            g_arr = dict_out['g_best_so_far']
            dict_out['x_best_so_far'][0] = np.array(x0)
            dict_out['f_best_so_far'][0] = init_out[0]
            dict_out['g_best_so_far'][0] = np.array(init_out[1])
            complete_list[i] = dict_out
    else:
        for i in range(len(complete_list)):
            dict_out = complete_list[i]
            f_arr = dict_out['f_best_so_far']
            N_eval = len(f_arr)
            g_arr = dict_out['g_best_so_far']
            dict_out['x_best_so_far'][0] = np.array(x0)
            dict_out['f_best_so_far'][0] = init_out[0]
            dict_out['g_best_so_far'][0] = np.array(init_out[1])

            for j in range(1, N_eval):
                if (g_arr[j] > 1e-3).any() or (init_out[0] < f_arr[j]):
                    dict_out['x_best_so_far'][j] = np.array(x0)
                    dict_out['f_best_so_far'][j] = init_out[0]
                    dict_out['g_best_so_far'][j] = np.array(init_out[1])
            complete_list[i] = dict_out

    return complete_list

def Problem_quadratic(x):
    f = quadratic_constrained.quadratic_f
    g = quadratic_constrained.quadratic_g

    return f(x), [g(x)]

bounds = np.array([[-1.5,1.5],[-1.5,1.5]])
# x0 = np.array([1, 1])
x0 = np.array([-0.5, 1.5])

max_f_eval = 100
max_it = 100

initial_output = Problem_quadratic(x0)

quadraticTight_pybobyqa = PyBobyqaWrapper().solve(Problem_quadratic, x0, bounds=bounds.T, \
                                      maxfun= max_f_eval, constraints=2, \
                                      seek_global_minimum= True, \
                                      objfun_has_noise=False)
# print(x0)

N = 10
quadraticTight_Nest_list = []
for i in range(N):
    rnd_seed = i
    quadratic_Nest = nesterov_random(Problem_quadratic, x0, bounds, max_iter = 100, \
                          constraints = 1, rnd_seed = i, alpha = 1e-3, mu = 1e-2, \
                           max_f_eval = max_f_eval)
    quadraticTight_Nest_list.append(quadratic_Nest)
print('10 Nesterov iterations completed')
# print(x0)

N = 10
quadraticTight_simplex_list = []
for i in range(N):
    rnd_seed = i
    quadratic_simplex = simplex_method(Problem_quadratic, x0, bounds, max_iter = 100, \
                            constraints = 1, rnd_seed = i, max_f_eval = max_f_eval)
    quadraticTight_simplex_list.append(quadratic_simplex)
print('10 simplex iterations completed')

quadraticTight_FiniteDiff = finite_Diff_Newton(Problem_quadratic, x0, bounds = bounds, \
                                   con_weight = 100)

quadraticTight_BFGS = BFGS_optimizer(Problem_quadratic, x0, bounds = bounds, \
                         con_weight = 100)

quadraticTight_Adam = Adam_optimizer(Problem_quadratic, x0, method = 'forward', \
                                      bounds = bounds, alpha = 0.4, \
                                      beta1 = 0.2, beta2  = 0.1, \
                                      max_f_eval = 100, con_weight = 100)

N_min_s = 15
init_radius = 2
method = 'Discrimination'
N = 10
quadraticTight_CUATRO_global_list = []
for i in range(N):
    rnd_seed = i
    quadratic_CUATRO_global = CUATRO(Problem_quadratic, x0, init_radius, bounds = bounds, \
                          N_min_samples = N_min_s, tolerance = 1e-10,\
                          beta_red = 0.9, rnd = rnd_seed, method = 'global', \
                          constr_handling = method)
    quadraticTight_CUATRO_global_list.append(quadratic_CUATRO_global)
print('10 CUATRO global iterations completed')

N_min_s = 6
init_radius = 0.5
method = 'Fitting'
N = 10
quadraticTight_CUATRO_local_list = []
for i in range(N):
    rnd_seed = i
    quadratic_CUATRO_local = CUATRO(Problem_quadratic, x0, init_radius, bounds = bounds, \
                          N_min_samples = N_min_s, tolerance = 1e-10,\
                          beta_red = 0.9, rnd = rnd_seed, method = 'local', \
                          constr_handling = method)
    quadraticTight_CUATRO_local_list.append(quadratic_CUATRO_local)
print('10 CUATRO local iterations completed')

N = 10
quadraticTight_SQSnobFit_list = []
for i in range(N):
    quadratic_SQSnobFit = SQSnobFitWrapper().solve(Problem_quadratic, x0, bounds, \
                                   maxfun = max_f_eval, constraints=2)
    quadraticTight_SQSnobFit_list.append(quadratic_SQSnobFit)
print('10 SnobFit iterations completed')

N = 10
quadraticTight_DIRECT_list = []
quadratic_DIRECT_f = lambda x, grad: Problem_quadratic(x)
for i in range(N):
    quadratic_DIRECT =  DIRECTWrapper().solve(quadratic_DIRECT_f, x0, bounds, \
                                   maxfun = max_f_eval, constraints=2)
    quadraticTight_DIRECT_list.append(quadratic_DIRECT)
print('10 DIRECT iterations completed')

# with open('BayesQuadratic_list.pickle', 'rb') as handle:
#     quadratic_Bayes_list = pickle.load(handle)

# quadratic_Bayes_list = fix_starting_points(quadratic_Bayes_list, x0, initial_output)
quadraticTight_DIRECT_list = fix_starting_points(quadraticTight_DIRECT_list, x0, initial_output)
quadraticTight_simplex_list = fix_starting_points(quadraticTight_simplex_list, x0, initial_output)
quadraticTight_pybobyqa['x_best_so_far'][0] = np.array(x0)
quadraticTight_pybobyqa['f_best_so_far'][0] = initial_output[0]
quadraticTight_pybobyqa['g_best_so_far'][0] = np.array(initial_output[1])

x_best_pyBbyqa = np.array(quadraticTight_pybobyqa['x_best_so_far'])
f_best_pyBbyqa = np.array(quadraticTight_pybobyqa['f_best_so_far'])
# x_ind_pyBbyqa = np.array(RB_pybobyqa['samples_at_iteration'])
# nbr_feval_pyBbyqa = len(RB_pybobyqa['f_store'])

x_best_finDiff = np.array(quadraticTight_FiniteDiff['x_best_so_far'])
f_best_finDiff = np.array(quadraticTight_FiniteDiff['f_best_so_far'])
x_ind_findDiff = np.array(quadraticTight_FiniteDiff['samples_at_iteration'])

x_best_BFGS = np.array(quadraticTight_BFGS['x_best_so_far'])
f_best_BFGS = np.array(quadraticTight_BFGS['f_best_so_far'])
x_ind_BFGS = np.array(quadraticTight_BFGS['samples_at_iteration'])

x_best_Adam = np.array(quadraticTight_Adam['x_best_so_far'])
f_best_Adam = np.array(quadraticTight_Adam['f_best_so_far'])
x_ind_Adam = np.array(quadraticTight_Adam['samples_at_iteration'])

plt.rcParams["font.family"] = "Times New Roman"
ft = int(15)
font = {'size': ft}
plt.rc('font', **font)
params = {'legend.fontsize': 12.5,
              'legend.handlelength': 2}
plt.rcParams.update(params)

fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax1.step(np.arange(len(f_best_pyBbyqa)), f_best_pyBbyqa, where = 'post', \
         label = 'PyBobyqa')
ax1.step(x_ind_findDiff, f_best_finDiff, where = 'post', \
         label = 'Newton Fin. Diff.')
ax1.step(x_ind_BFGS, f_best_BFGS, where = 'post', \
         label = 'BFGS')
ax1.step(x_ind_Adam, f_best_Adam, where = 'post', \
         label = 'Adam')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax1.legend()
ax1.set_yscale('log')
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_Deterministic_Convergence_plot.svg', format = "svg")

quadratic_f = lambda x, y: x**2 + 10*y**2 + x*y
quadratic_g = lambda x, y: 1 - x - y <= 0

oracle = RB(quadratic_f, ineq = [quadratic_g])

ax2, fig2 = trust_fig(oracle, bounds)
ax2.plot(x_best_pyBbyqa[:,0], x_best_pyBbyqa[:,1], '--x', \
         label = 'PyBobyqa')
ax2.plot(x_best_finDiff[:,0], x_best_finDiff[:,1], '--x', \
         label = 'Newton Fin. Diff.')
ax2.plot(x_best_BFGS[:,0], x_best_BFGS[:,1], '--x', \
         label = 'BFGS')
ax2.plot(x_best_Adam[:,0], x_best_Adam[:,1], '--x', \
         label = 'Adam')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.legend()
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_Deterministic_2Dspace_plot.svg', format = "svg")

fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(len(quadraticTight_CUATRO_global_list)):
    x_best = np.array(quadraticTight_CUATRO_global_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_CUATRO_global_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_CUATRO_global_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'CUATRO_g'+str(i))
    # ax1.plot(x_ind, f_best, label = 'CUATRO_g'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'CUATRO_g'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.legend()
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_CUATROg_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_CUATROg_2Dspace_plot.svg', format = "svg")

fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(len(quadraticTight_CUATRO_local_list)):
    x_best = np.array(quadraticTight_CUATRO_local_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_CUATRO_local_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_CUATRO_local_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'CUATRO_l'+str(i))
    # ax1.plot(x_ind, f_best, label = 'CUATRO_l'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'CUATRO_l'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.legend()
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_CUATROl_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_CUATROl_2Dspace_plot.svg', format = "svg")

# fig1 = plt.figure()
# ax1 = fig1.add_subplot()
# ax2, fig2 = trust_fig(oracle, bounds)
# for i in range(len(quadraticTight_Bayes_list)):
#     x_best = np.array(quadraticTight_Bayes_list[i]['x_best_so_far'])
#     f_best = np.array(quadraticTight_Bayes_list[i]['f_best_so_far'])
#     nbr_feval = len(quadraticTight_Bayes_list[i]['f_store'])
#     ax1.step(np.arange(len(f_best)), f_best, where = 'post', \
#           label = 'BO'+str(i)+'; #f_eval: ' + str(nbr_feval))
#     ax2.plot(x_best[:,0], x_best[:,1], '--', \
#           label = 'BO'+str(i)+'; #f_eval: ' + str(nbr_feval))
# ax1.legend()
# ax1.set_yscale('log')
# ax1.set_xlabel('Nbr. of function evaluations')
# ax1.set_ylabel('Best function evaluation')
# ax2.set_xlabel('$x_1$')
# ax2.set_ylabel('$x_2$')
# ax2.legend()
# ax2.set_xlim(bounds[0])
# ax2.set_ylim(bounds[1])
# fig1.savefig('QuadraticTight_plots/Quadratic_BO_Convergence_plot.svg', format = "svg")
# fig2.savefig('QuadraticTight_plots/Quadratic_BO_2Dspace_plot.svg', format = "svg")

fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(len(quadraticTight_simplex_list)):
    x_best = np.array(quadraticTight_simplex_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_simplex_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_simplex_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'Simplex'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'Simplex'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
ax2.legend()
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_simplex_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_simplex_2Dspace_plot.svg', format = "svg")

## Change to x_best_So_far
fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(N):
    x_best = np.array(quadraticTight_Nest_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_Nest_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_Nest_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'Nest.'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'Nest.'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
ax2.legend()
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_Nesterov_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_Nesterov_2Dspace_plot.svg', format = "svg")

## Change to x_best_So_far
fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(N):
    x_best = np.array(quadraticTight_SQSnobFit_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_SQSnobFit_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_SQSnobFit_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'Snobfit'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'Snobfit'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
ax2.legend()
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_SQSnobFit_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_SQSnobFit_2Dspace_plot.svg', format = "svg")

## Change to x_best_So_far
fig1 = plt.figure()
ax1 = fig1.add_subplot()
ax2, fig2 = trust_fig(oracle, bounds)
for i in range(N):
    x_best = np.array(quadraticTight_DIRECT_list[i]['x_best_so_far'])
    f_best = np.array(quadraticTight_DIRECT_list[i]['f_best_so_far'])
    x_ind = np.array(quadraticTight_DIRECT_list[i]['samples_at_iteration'])
    ax1.step(x_ind, f_best, where = 'post', label = 'DIRECT'+str(i))
    ax2.plot(x_best[:,0], x_best[:,1], '--', label = 'DIRECT'+str(i))
ax1.legend()
ax1.set_yscale('log')
ax1.set_xlabel('Nbr. of function evaluations')
ax1.set_ylabel('Best function evaluation')
ax2.set_xlabel('$x_1$')
ax2.set_ylabel('$x_2$')
ax2.set_xlim(bounds[0])
ax2.set_ylim(bounds[1])
ax2.legend()
fig1.savefig('Quadratic_plots_Tight/QuadraticTight_DIRECT_Convergence_plot.svg', format = "svg")
fig2.savefig('Quadratic_plots_Tight/QuadraticTight_DIRECT_2Dspace_plot.svg', format = "svg")

sol_Cg = average_from_list(quadraticTight_CUATRO_global_list)
test_CUATROg, test_av_CUATROg, test_min_CUATROg, test_max_CUATROg = sol_Cg
sol_Cl = average_from_list(quadraticTight_CUATRO_local_list)
test_CUATROl, test_av_CUATROl, test_min_CUATROl, test_max_CUATROl = sol_Cl
sol_Nest = average_from_list(quadraticTight_Nest_list)
test_Nest, test_av_Nest, test_min_Nest, test_max_Nest = sol_Nest
sol_Splx = average_from_list(quadraticTight_simplex_list)
test_Splx, test_av_Splx, test_min_Splx, test_max_Splx = sol_Splx
sol_SQSF = average_from_list(quadraticTight_SQSnobFit_list)
test_SQSF, test_av_SQSF, test_min_SQSF, test_max_SQSF = sol_SQSF
sol_DIR = average_from_list(quadraticTight_DIRECT_list)
test_DIR, test_av_DIR, test_min_DIR, test_max_DIR = sol_DIR
# sol_BO = average_from_list(quadraticTight_Bayes_list)
# test_BO, test_av_BO, test_min_BO, test_max_BO = sol_BO

fig = plt.figure()
ax = fig.add_subplot()
ax.step(np.arange(1, 101), test_av_CUATROg, where = 'post', label = 'CUATRO_g', c = 'b')
ax.fill_between(np.arange(1, 101), test_min_CUATROg, \
                test_max_CUATROg, color = 'b', alpha = .5, step = 'post')
ax.step(np.arange(1, 101), test_av_CUATROl, where = 'post', label = 'CUATRO_l', c = 'c')
ax.fill_between(np.arange(1, 101), test_min_CUATROl, \
                test_max_CUATROl, color = 'c', alpha = .5, step = 'post')
ax.step(np.arange(1, 101), test_av_SQSF, where = 'post', label = 'Snobfit', c = 'orange')
ax.fill_between(np.arange(1, 101), test_min_SQSF, \
                test_max_SQSF, color = 'orange', alpha = .5, step = 'post')
ax.step(np.arange(len(f_best_pyBbyqa)), f_best_pyBbyqa, where = 'post', \
          label = 'Py-BOBYQA', c = 'green')
# ax.step(np.arange(1, 101), test_av_BO, where = 'post', \
#           label = 'BO', c = 'r')
# ax.fill_between(np.arange(1, 101), test_min_BO, \
#                 test_max_BO, color = 'r', alpha = .5, step = 'post')

ax.legend()
ax.set_xlabel('Number of function evaluations')
ax.set_ylabel('Best function evaluation')
ax.set_yscale('log')
ax.set_xlim([1, 100])
ax.set_ylim([0.8, 25])
fig.savefig('Quadratic_publication_plots/QuadraticTight_Model.svg', format = "svg")

fig = plt.figure()
ax = fig.add_subplot()
ax.step(np.arange(1, 101), test_av_Nest, where = 'post', label = 'Nesterov', c = 'brown')
ax.fill_between(np.arange(1, 101), test_min_Nest, \
                test_max_Nest, color = 'brown', alpha = .5, step = 'post')
ax.step(np.arange(1, 101), test_av_Splx, where = 'post', label = 'Simplex', c = 'green')
ax.fill_between(np.arange(1, 101), test_min_Splx, \
                test_max_Splx, color = 'green', alpha = .5, step = 'post')
ax.step(x_ind_findDiff, f_best_finDiff, where = 'post', \
          label = 'Newton', c = 'black')
ax.step(x_ind_BFGS, f_best_BFGS, where = 'post', \
          label = 'BFGS', c = 'orange')
ax.step(x_ind_Adam, f_best_Adam, where = 'post', \
          label = 'Adam', c = 'blue')
ax.step(np.arange(1, 101), test_av_DIR, where = 'post', label = 'DIRECT', c = 'violet')
ax.fill_between(np.arange(1, 101), test_min_DIR, \
                test_max_DIR, color = 'violet', alpha = .5, step = 'post')

ax.legend()
ax.set_xlabel('Number of function evaluations')
ax.set_ylabel('Best function evaluation')
ax.set_yscale('log')
ax.set_xlim([1, 100])
ax.set_ylim([0.8, 25])
ax.legend(loc = 'upper right')
fig.savefig('Quadratic_publication_plots/QuadraticTight_Others.svg', format = "svg")


text:
from __future__ import division
import numpy as np
import pandas
import json as js
import tarfile as tf
import os
from DSTK.GAM.utils.shape_function import ShapeFunction
from DSTK.utils.function_helpers import sigmoid
import abc

class BaseGAM(object):
    __metaclass__ = abc.ABCMeta

    def __init__(self):
        self.shapes = dict()
        self.is_fit = False
        self.initialized = False
        self.n_features = None
        self.feature_names = None

    def _get_index_for_feature(self, feature_name):
        return self.feature_names.index(feature_name)

    def logit_score(self, vec):
        return np.sum([func.get_value(vec[self._get_index_for_feature(feat)]) for feat, func in self.shapes.iteritems()])

    def _score_single_record(self, vec):
        return [sigmoid(-2 * np.sum([func.get_value(vec[self._get_index_for_feature(feat)]) for feat, func in self.shapes.iteritems()])), \
                sigmoid( 2 * np.sum([func.get_value(vec[self._get_index_for_feature(feat)]) for feat, func in self.shapes.iteritems()]))]

    def _get_feature_value_pair_single_record(self, vec):
        return sorted([(feat, func.get_value(vec[self._get_index_for_feature(feat)])) for feat, func in self.shapes.iteritems()], key=lambda x: x[1])

    def feature_value_pairs(self, X):
        if isinstance(X, pandas.core.frame.DataFrame):
            data = X.as_matrix()
        else:
            data = X

        if data.ndim == 1:
            return self._get_feature_value_pair_single_record(data)
        if data.ndim == 2:
            return [self._get_feature_value_pair_single_record(vec) for vec in data]

    def score(self, X):
        if isinstance(X, pandas.core.frame.DataFrame):
            data = X.as_matrix()
        else:
            data = X

        if data.ndim == 1:
            return self._score_single_record(data)
        if data.ndim == 2:
            return np.asarray([self._score_single_record(vec) for vec in data])

    @staticmethod
    def _check_input(data, labels):
        if not np.isfinite(data).ravel().any():
            raise ValueError("Encountered invalid value in the training data")
        if not np.isfinite(labels).any():
            raise ValueError("Encountered invalid value in the target data")

        assert len(data) == len(labels), "Data and Targets have different lentgth."

    @abc.abstractmethod
    def train(self, data, targets, **kwargs):
        raise NotImplementedError

    @abc.abstractmethod
    def _get_metadata_dict(self):
        raise NotImplementedError

    def serialize(self, model_name, file_path=None):

        if not self.is_fit:
            raise RuntimeError("Cannot serialize untrained GAM.")

        if file_path is None:
            file_path = os.getcwd()

        folder_path = '{}/{}'.format(file_path, model_name)
        if not os.path.exists(folder_path):
            os.makedirs(folder_path)

        for shape in self.shapes.itervalues():
            shape.serialize('{}/{}/shapes'.format(file_path, model_name), meta_data={'model_name': model_name})

        meta_data_dct = self._get_metadata_dict()

        dct = {
            'model_name': model_name,
            'num_features': self.n_features,
            'feature_names': self.feature_names,
            'metadata': meta_data_dct,
            'shape_data': './shapes'
            }

        with open('{}/{}/{}.json'.format(file_path, model_name, model_name), 'w') as fp:
            js.dump(dct, fp, sort_keys=True, indent=2, separators=(',', ': '))

        self._make_tarfile('{}/{}.tar.gz'.format(file_path, model_name), '{}/{}'.format(file_path, model_name))

    @staticmethod
    def _make_tarfile(output_filename, source_dir):
        with tf.open(output_filename, "w:gz") as tar:
            tar.add(source_dir, arcname=os.path.basename(source_dir))

def load_from_tar(file_name):

    model_name = file_name.split('/')[-1].replace('.tar.gz', '')

    gam = DeserializedGAM()

    with tf.open(file_name, "r:gz") as tar:
        f = tar.extractfile('{}/{}.json'.format(model_name, model_name))
        content = js.loads(f.read())
        gam._n_features = content['num_features']
        gam.feature_names = content['feature_names']

        for member in tar.getmembers():
            if member.isfile() and (member.name != '{}/{}.json'.format(model_name, model_name)):
                f = tar.extractfile(member.path)
                content = js.loads(f.read())
                gam.shapes.update({content['feature_name']: ShapeFunction(content['splits'],
                                                                          content['values'],
                                                                          content['feature_name'])})

    assert set(gam.shapes.keys()) == set(gam.feature_names), 'feature names and shape names do not match up'
    return gam

class DeserializedGAM(BaseGAM):
    def _get_metadata_dict(self):
        pass

    def train(self, data, targets, **kwargs):
        pass

    def __init__(self, **kwargs):
        super(DeserializedGAM, self).__init__()

ext: py

max_stars_repo_name: jotterbach/dstk

max_issues_repo_name: jotterbach/dstk

max_forks_repo_name: jotterbach/dstk

reason: import numpy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np

from configs.constants import Constants

class ShapeGenerator:

    TEST_RANDOM_VALUES = [[0, 0.5], [0.3, 0.9], [0, 0.8], [0.5, 0.6], [0.5, 0.9],
                          [0.8, 0], [0.5, 0.1], [1, 0.4], [0.9, 0.2], [1, 0]]

    def __init__(self):
        basic_parameters = Constants.get_basic_parameters()
        physics_parameters = basic_parameters["physics"]
        self.min_permittivity = physics_parameters["min_permittivity"]
        self.max_permittivity = physics_parameters["max_permittivity"]

    def generate_shapes(self, no_of_shapes, test, image_i):
        pass

    def get_test_shape_parameters(self, max_side, image_i, shape_i):
        center_range = [-1 + max_side + 0.05, 1 - max_side - 0.05]
        center_x = center_range[0] + (center_range[1] - center_range[0]) * self.TEST_RANDOM_VALUES[image_i - 1][shape_i]
        center_y = center_range[0] + (center_range[1] - center_range[0]) * self.TEST_RANDOM_VALUES[image_i - 1][shape_i]
        relative_permittivity = self.min_permittivity + (
                self.max_permittivity - self.min_permittivity) * self.TEST_RANDOM_VALUES[image_i - 1][shape_i]

        return center_x, center_y, relative_permittivity

    def get_shape_parameters(self, max_side):
        center_range = [-1 + max_side + 0.05, 1 - max_side - 0.05]
        center_x = center_range[0] + (center_range[1] - center_range[0]) * np.random.uniform()
        center_y = center_range[0] + (center_range[1] - center_range[0]) * np.random.uniform()
        relative_permittivity = self.min_permittivity + (
                    self.max_permittivity - self.min_permittivity) * np.random.uniform()

        return center_x, center_y, relative_permittivity

    def get_shape_name(self):
        pass

max_stars_repo_name: estanislaoledesma/genper

max_issues_repo_name: estanislaoledesma/genper

max_forks_repo_name: estanislaoledesma/genper


text:
import time

import numpy as np
import torch.nn as nn
import torch

class Batch:
    """
    Object for holding a batch of data with mask during training.
    """

    def __init__(self, src, trg=None, pad=0):
        self.src = src
        self.src_mask = (src != pad).unsqueeze(-2)
        if trg is not None:
            # in the decoder, we need predict next word base on the what we have generate.
            # trg remove the last word and trg_y remove the first word, so when decoder generate trg[i],
            # the next word is expected as trg_y[i]
            self.trg = trg[:, :-1]
            self.trg_y = trg[:, 1:]
            self.trg_mask = self.make_std_mask(self.trg, pad)
            self.n_tokens = (self.trg_y != pad).data.sum()

    @staticmethod
    def make_std_mask(tgt, pad):
        """
        Create a mask to hide padding and future words.
        tgt的形状是（batch_size, token_length）而在后续的处理中，经过词嵌入以后形状会变为（batch_size, token_length, d_model）
        在计算multi-head attention的情况下，会变成四维向量：（batch_size, h, token_length, d_model）
        由于不同的head的mask是相同的，所以只需要在第一个维度添加一个占位维度即可。
        所以在生成tgt_mask时，需要在（batch_size, token_length）的中间加一个占位维度。
        """
        tgt_mask = (tgt != pad).unsqueeze(-2)
        res = subsequent_mask(tgt.size(-1))
        tgt_mask = tgt_mask & res
        return tgt_mask

def subsequent_mask(size):
    """
    Mask out subsequent positions.
    """
    # make the num matrix to a bool matrix
    return np.triu(np.ones((1, size, size), dtype=np.int8), k=1) == 0

global max_src_in_batch, max_tgt_in_batch

def batch_size_fn(new, count):
    """
    Keep augmenting batch and calculate total number of tokens + padding.
    """
    global max_src_in_batch, max_tgt_in_batch
    if count == 1:
        max_src_in_batch = 0
        max_tgt_in_batch = 0
    max_src_in_batch = max(max_src_in_batch,  len(new.src))
    max_tgt_in_batch = max(max_tgt_in_batch,  len(new.trg) + 2)
    src_elements = count * max_src_in_batch
    tgt_elements = count * max_tgt_in_batch
    return max(src_elements, tgt_elements)

class NoamOpt:
    """
    Optimizer wrapper that implements rate.
    """
    def __init__(self, model_size, factor, warmup, optimizer):
        self.optimizer = optimizer
        self._step = 0
        self.warmup = warmup
        self.factor = factor
        self.model_size = model_size
        self._rate = 0

    def step(self):
        """
        Update parameters and rate
        """
        self._step += 1
        rate = self.rate()
        for p in self.optimizer.param_groups:
            p['lr'] = rate
        self._rate = rate
        self.optimizer.step()

    def rate(self, step=None):
        """
        Implement `l_rate` above"
        """
        if step is None:
            step = self._step
        return self.factor * (self.model_size ** (-0.5) * min(step ** (-0.5), step * self.warmup ** (-1.5)))

def get_std_opt(model):
    return NoamOpt(model_size=model.src_embed[0].d_model,
                   factor=2,
                   warmup=4000,
                   optimizer=torch.optim.Adam(model.parameters(), lr=0, betas=(0.9, 0.98), eps=1e-9))

def data_gen(V, batch, batches):
    """
    Generate random data for a src-tgt copy task.
    """
    for i in range(batches):
        data = np.random.randint(1, V, size=(batch, 10), dtype=np.int64)
        data[:, 0] = 1
        src = torch.from_numpy(data)
        tgt = torch.from_numpy(data)
        yield Batch(src, tgt, 0)

class SimpleLossCompute:
    """
    A simple loss compute and train function.
    """

    def __init__(self, generator, criterion, opt=None):
        self.generator = generator
        self.criterion = criterion
        self.opt = opt

    def __call__(self, x, y, norm):
        x = self.generator(x)

        loss = self.criterion(x.contiguous().view(-1, x.size(-1)), y.contiguous().view(-1)) / norm
        loss.backward()
        if self.opt is not None:
            self.opt.step()
            self.opt.optimizer.zero_grad()
        return loss.data * norm

class LabelSmoothing(nn.Module):
    """
    Implement label smoothing.
    """

    def __init__(self, size, padding_idx, smoothing=0.0):
        super(LabelSmoothing, self).__init__()
        self.criterion = nn.KLDivLoss(size_average=False)
        self.padding_idx = padding_idx
        self.confidence = 1.0 - smoothing
        self.smoothing = smoothing
        self.size = size
        self.true_dist = None

    def forward(self, x, target):
        assert x.size(1) == self.size
        true_dist = x.data.clone()
        true_dist.fill_(self.smoothing / (self.size - 2))
        true_dist.scatter_(1, target.data.unsqueeze(1), self.confidence)
        true_dist[:, self.padding_idx] = 0
        mask = torch.nonzero(target.data == self.padding_idx)
        if mask.dim() > 0:
            true_dist.index_fill_(0, mask.squeeze(), 0.0)
        self.true_dist = true_dist
        return self.criterion(x, true_dist)

def run_epoch(data_iter, model, loss_compute):
    """
    Standard Training and Logging Function
    """
    start = time.time()
    total_tokens = 0
    total_loss = 0
    tokens = 0
    for i, batch in enumerate(data_iter):
        out = model.forward(batch.src, batch.trg, batch.src_mask, batch.trg_mask)
        loss = loss_compute(out, batch.trg_y, batch.n_tokens)
        total_loss += loss
        total_tokens += batch.n_tokens
        tokens += batch.n_tokens
        if i % 50 == 1:
            elapsed = time.time() - start
            print("Epoch Step: %d Loss: %f Tokens per Sec: %f" % (i, loss / batch.n_tokens, tokens / elapsed))
            start = time.time()
            tokens = 0
    return total_loss / total_tokens

def greedy_decode(model, src, src_mask, max_len, start_symbol):
    memory = model.encode(src, src_mask)
    ys = torch.ones(1, 1).fill_(start_symbol).type_as(src.data)
    for i in range(max_len-1):
        out = model.decode(memory, src_mask, ys, subsequent_mask(ys.size(1)))
        prob = model.generator(out[:, -1])
        _, next_word = torch.max(prob, dim=1)
        next_word = next_word.data[0]
        ys = torch.cat([ys,
                        torch.ones(1, 1).type_as(src.data).fill_(next_word)], dim=1)
    return ys

max_stars_repo_name: MartrixG/transformer

max_issues_repo_name: MartrixG/transformer

max_forks_repo_name: MartrixG/transformer


text:
'''
This son of a gun tests the restructuring of TF tensors for an SOM project.

bmu_index = tf.argmin(tf.sqrt(tf.reduce_sum(
tf.pow(tf.sub(self._weightage_vects, tf.pack(
    [self._vect_input for i in range(m*n)])), 2), 1)),
                  0)

RESULTS in above: 2 = the exp for tf.pow,
                  1 = sum across the RGB channels,
                  0 = min index across the 1D (0th dimension) vector

Hey also, these are just notes, so none of this stuff runs together or anything like that. I'm just using pieces of this script in various situations.

'''

import numpy as np
import tensorflow as tf

#### 2D output / 1D input example #####################################################
'''
ns = 4
dim = 3

m = 8
n = 9

temp = np.zeros((ns,dim))
for i in range(ns):
	for j in range(dim):
		temp[i,j] = i*10 + j

weights = tf.random_normal([m*n, dim])
inp = tf.constant(temp[0], dtype=tf.float32)

bmuIdx = tf.argmin(tf.sqrt(tf.reduce_sum(tf.pow(tf.sub(weights, tf.pack([inp for i in range(m*n)])), 2), 1)), 0)
'''
### 3D output / 1D input example #####################################################
ns = 4
dim = 3

l = 7
m = 8
n = 9

temp = np.zeros((ns,dim))
for i in range(ns):
	for j in range(dim):
		temp[i,j] = i*10 + j

#print('original = \n', temp, '\n')

weights = np.zeros((l*m*n, dim))
for i in range(l*m*n):
	for j in range(dim):
		weights[i,j] = (i - (l*m*n/2))**2 + j

#weights = tf.random_normal([l*m*n, dim])
inp = tf.constant(temp[0])

bmuIdx = tf.argmin(tf.sqrt(tf.reduce_sum(tf.pow(tf.sub(weights, tf.pack([inp for i in range(l*m*n)])), 2), 1)), 0)

### neuron locations class ######################################################
#3D output
def neuron_locations(l, m, n):
	for k in range(l):
		for i in range(m):
			for j in range(n):
				yield np.array([k, i, j])

#2D output
'''
def neuron_locations(m, n):
	for i in range(m):
		for j in range(n):
			yield np.array([i, j])
'''
### slice stuff #################################################################
#3D output
locs = tf.constant(np.array(list(neuron_locations(l, m, n))))

slce = tf.pad(tf.reshape(bmuIdx, [1]), np.array([[0,1]]))

bmu_loc = tf.reshape(tf.slice(locs, slce, tf.constant(np.array([1, 3]))), [3])

#2D output
'''
locs = tf.constant(np.array(list(neuron_locations(m, n))))

slce = tf.pad(tf.reshape(bmuIdx, [1]), np.array([[0,1]]))

bmu_loc = tf.reshape(tf.slice(locs, slce, tf.constant(np.array([1, 2]))), [2])
'''

sess = tf.Session()
print('bmu index = ', sess.run(bmuIdx))
print('slice = ', sess.run(slce))
#print(sess.run(locs))
print('bmu location = ', sess.run(bmu_loc))

max_stars_repo_name: gganssle/organize-yaself

max_issues_repo_name: gganssle/organize-yaself

max_forks_repo_name: gganssle/organize-yaself


text:
"""

获取通用数据的工具库

"""

import pandas as pd
import numpy as np
from sklearn.datasets import load_iris, load_boston, load_breast_cancer, load_diabetes

import tushare as ts

from mllib.utils import dateutils

import os.path

CURRENT_PATH = os.path.abspath(os.path.dirname(__file__))

# 链家网数据
def get_lianjia_data():
    path = os.path.join(CURRENT_PATH, "../data/lj_data.csv")
    return pd.read_csv(path, encoding='gb18030')

# 正太分布数据
def get_normalize_data(size = 1000):
    x = np.random.normal(size)
    y = np.random.normal(size)

    return (x, y)

# 均匀分布
def get_linspace_data(a, b, num=5):
    return np.linspace(a, b, num)

# flower
def get_iris_flower_data():
    return load_iris()

def get_boston_house_price():
    return load_boston()

# 乳谢癌数据
def get_breast_cancer_disease_data():
    return load_breast_cancer()

# 糖尿病数据
def get_diabetes_disease_data():
    return load_diabetes()

# 小费
def get_tips_data():
    path = os.path.join(CURRENT_PATH, "../data/tips.csv")
    return pd.read_csv(path)

# 坦坦尼克号死亡情况数据
def get_titanic_data():
    path = os.path.join(CURRENT_PATH, "../data/titanic.csv")
    return pd.read_csv(path)

# 银行用户信用数据
def get_bank_credit_data():
    path = os.path.join(CURRENT_PATH, "../data/credit-data.csv")
    return pd.read_csv(path)

# 航班用户数据
def get_flights_passengers_data():
    path = os.path.join(CURRENT_PATH, "../data/flights.csv")
    return pd.read_csv(path)

# Google stock data
def get_stocks_google():
    path = os.path.join(CURRENT_PATH, "../data/GOOG.csv")
    return pd.read_csv(path)

# apple stock data
def get_stocks_nvidia():
    path = os.path.join(CURRENT_PATH, "../data/NVDA.csv")
    return pd.read_csv(path)

# msft stock data
def get_stocks_msft():
    path = os.path.join(CURRENT_PATH, "../data/MSFT.csv")
    return pd.read_csv(path)

# apple stock data
def get_stocks_apple():
    path = os.path.join(CURRENT_PATH, "../data/AAPL.csv")
    return pd.read_csv(path)

# 80条最新财经新闻
def get_latest_finance_news():
    return ts.get_latest_news()  # 默认获取最近80条新闻数据，只提供新闻类型、链接和标题

# 实时电影票房数据
def get_realtime_movies_boxoffice():
    return ts.realtime_boxoffice()

# 昨日票房
def get_yesterday_movies_boxoffice():
    return ts.day_boxoffice()

# 上月票房
def get_last_month_movies_boxoffice():
    return ts.month_boxoffice()

# 股票 行业分类
def get_finance_industry_classify():
    return ts.get_industry_classified()

# 股票 概念分类
def get_finance_concept_classify():
    return ts.get_concept_classified()

# 股票龙虎榜数据
def get_finance_top_list():
    lst =  ts.top_list(dateutils.get_today_str())
    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(1))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(2))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(3))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(4))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(5))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(6))

    if len(lst) == 0:
        lst = ts.top_list(dateutils.get_backdays_str(7))

    return lst

# 营业部上榜数据
def get_finance_broker_top_list():
    return ts.broker_tops()

# yahoo 收盘价s
# see here: http://www.learndatasci.com/python-finance-part-yahoo-finance-api-pandas-matplotlib/
def get_yahoo_finance_close_data(froms='2016-01-01',tos='2017-08-28', tickers=['AAPL', 'MSFT', 'SPY']):
    from pandas_datareader import data

    # Define the instruments to download. We would like to see Apple, Microsoft and the S&P500 index.
    #tickers = ['AAPL', 'MSFT', 'SPY']

    # Define which online source one should use
    data_source = 'google'

    # We would like all available data from 01/01/2000 until 12/31/2016.
    start_date = froms
    end_date = tos

    # User pandas_reader.data.DataReader to load the desired data. As simple as that.
    panel_data = data.DataReader(tickers, data_source, start_date, end_date)

    # Getting just the adjusted closing prices. This will return a Pandas DataFrame
    # The index in this DataFrame is the major index of the panel_data.
    close = panel_data.ix['Close']

    # Getting all weekdays between 01/01/2000 and 12/31/2016
    all_weekdays = pd.date_range(start=start_date, end=end_date, freq='B')

    # How do we align the existing prices in adj_close with our new set of dates?
    # All we need to do is reindex close using all_weekdays as the new index
    close = close.reindex(all_weekdays)

    return close

if __name__ == '__main__':
    df = get_lianjia_data()
    print(df)

max_stars_repo_name: foxlog/mllib

max_issues_repo_name: foxlog/mllib

max_forks_repo_name: foxlog/mllib


text:
# Nathaniel Gordon, CS4100 Final Project
import json
import numpy as np
import collections
from scipy import optimize
import random
from multiprocessing import Pool

# Evaluation function
def score(comb, db, pool):
    selected = []
    for c in comb:
        selected.append(int(c))

    # Return 0 if duplicate cards
    if len(selected) != len(set(selected)):
        return 0

    # Count how many cards break the curve
    overCurve = 0
    curveTally = np.zeros(7)
    idealCurve = [45, 2, 7, 6, 4, 3, 1]
    for selind in selected:
        cmc = db[pool[selind]]['convertedManaCost']
        if cmc < 6:
            if curveTally[cmc] == idealCurve[cmc]:
                overCurve = overCurve + 1
            else:
                curveTally[cmc] = curveTally[cmc] + 1
        else:
            if curveTally[6] == idealCurve[6]:
                overCurve = overCurve + 1
            else:
                curveTally[6] = curveTally[6] + 1
    curveScore = (len(comb) - overCurve) / len(comb)

    # Sum the total value of all cards in the deck
    cardScore = 0
    for selind in selected:
        cardScore = cardScore + db[pool[selind]]['value']

    return -1*(curveScore*cardScore)

# Similar to score(), but returns more feedback
def getDeckStats(deck, db, pool):
    # Count how many cards break the curve
    overCurve = 0
    curveTally = np.zeros(7)
    idealCurve = [45, 2, 7, 6, 4, 3, 1]
    for card in deck:
        cmc = db[card]['convertedManaCost']
        if cmc < 6:
            if curveTally[cmc] == idealCurve[cmc]:
                overCurve = overCurve + 1
            else:
                curveTally[cmc] = curveTally[cmc] + 1
        else:
            if curveTally[6] == idealCurve[6]:
                overCurve = overCurve + 1
            else:
                curveTally[6] = curveTally[6] + 1
    curveScore = float(len(deck) - overCurve) / float(len(deck))

    # Sum the total value of all cards in the deck
    cardScore = 0
    for card in deck:
        cardScore = cardScore + db[card]['value']

    return (cardScore, curveScore)

# Uniqueness constraint for optimizers
def uniqueConstraint(x):
    y=[]
    for xi in x:
        y.append(int(xi))
    if len(y) != len(set(y)):
        return 0
    return 1

# Perform an experimental trial
def runExp(key):
    key = str(key)
    print('Running trial', key)

    pool = pools[key]
    expDict = {'Initial': {}, 'Annealing': {}, 'Genetic': {}, 'Random': {}}

    expDict['Initial']['Deck']=[]
    for card in pool:
        expDict['Initial']['Deck'].append(card)

    expDict['Initial']['Score']=0

    # Remove basic lands
    tbr = []
    for card in pool:
        if 'Basic Land' in db[card]['type']:
            tbr.append(card)
    for card in tbr:
        pool.remove(card)

    # Remove underrepresented colors
    colorMembers = []
    for card in pool:
        for cardColor in db[card]['colorIdentity']:
            colorMembers.append(cardColor)

    colorRef = ['R', 'W', 'G', 'B', 'U']
    colorCounts = [colorMembers.count(c) for c in colorRef]
    removeColors = []
    for i in range(3):
        minColorInd = np.argmin(colorCounts)
        removeColors.append(colorRef[minColorInd])
        del colorRef[minColorInd]
        del colorCounts[minColorInd]

    tbr = []
    for card in pool:
        if len(list(set(removeColors).intersection(db[card]['colorIdentity']))) > 0:
            tbr.append(card)
    for card in tbr:
        pool.remove(card)

    varbound=np.array([[0,len(pool)-1]]*23)

    nlc = optimize.NonlinearConstraint(uniqueConstraint, .5, 1.5)

    annealRes = optimize.dual_annealing(score, args=(db, pool), bounds=varbound, maxiter=10000)
    if annealRes.fun == 0:
        print("Annealing Score:", annealRes.fun)
        return None

    genRes = optimize.differential_evolution(score, varbound, args=(db, pool), constraints=(nlc))
    randResSel = random.sample(range(len(pool)), 23)

    print("Annealing Score:", annealRes.fun)
    print("Genetic Score:", genRes.fun)

    if annealRes.fun < 0 and genRes.fun < 0:

        expDict['Annealing']['Deck']=[]
        for x in annealRes.x:
            expDict['Annealing']['Deck'].append(db[pool[int(x)]]['name'])

        expDict['Annealing']['Score']=annealRes.fun

        expDict['Annealing']['Card Score'], expDict['Annealing']['Curve Score'] = getDeckStats(expDict['Annealing']['Deck'], db, pool)

        expDict['Genetic']['Deck']=[]
        for x in genRes.x:
            expDict['Genetic']['Deck'].append(db[pool[int(x)]]['name'])

        expDict['Genetic']['Score']=genRes.fun

        expDict['Genetic']['Card Score'], expDict['Genetic']['Curve Score'] = getDeckStats(expDict['Genetic']['Deck'], db, pool)

        expDict['Random']['Deck']=[]
        for x in randResSel:
            expDict['Random']['Deck'].append(db[pool[int(x)]]['name'])

        expDict['Random']['Card Score'], expDict['Random']['Curve Score'] = getDeckStats(expDict['Random']['Deck'], db, pool)

        expDict['Random']['Score']= - expDict['Random']['Card Score']*expDict['Random']['Curve Score']

        return expDict

if __name__ == "__main__":

    # Load card database
    with open("data/m19cards.json", "r") as read_file:
        db = json.load(read_file)

    # Load trial pools
    with open("data/trialPools.json", "r") as read_file:
        pools = json.load(read_file)

    # Run experiment trials with pool
    trials = 2000
    resDict = {}

    mp = Pool()
    results = mp.map(runExp, range(trials))

    # Convert results to serializable JSON dict
    for i in range(len(results)):
        if results[i]:
            resDict[str(i)] = results[i]

    # Serializing json
    json_object = json.dumps(resDict, indent = 4)

    # Writing to sample.json
    with open("algOut.json", "w") as outfile:
        outfile.write(json_object)

reason: import numpy,from scipy


text:
import random
import math
import numpy as np
from GraphEngine import GraphicEngine

class FourierTransform:
    def __init__(self, input_signal_vector):
        self.number_of_data = len(input_signal_vector)
        self.number_of_transform = int(self.number_of_data / 2)

        self.input_signal_vector = input_signal_vector

        self.component_real = np.zeros((self.number_of_transform, 1), dtype=float)
        self.component_imaginary = np.zeros((self.number_of_transform, 1), dtype=float)

        self.component_phases = np.zeros((self.number_of_transform, 1), dtype=float)
        self.component_magnitude = np.zeros((self.number_of_transform, 1), dtype=float)

        self.synthesis_signal = np.zeros((self.number_of_data, 1), dtype=float)

        self.graph_engine = GraphicEngine()

    def fourier_transformation(self):
        for k in range(self.number_of_transform):
            self.component_real[k] = sum([self.input_signal_vector[i] * math.cos(2 * math.pi * k * i / self.number_of_data) for i in range(self.number_of_data)])
            self.component_imaginary[k] = sum([-self.input_signal_vector[i] * math.sin(2 * math.pi * k * i / self.number_of_data) for i in range(self.number_of_data)])

    def fourier_calculate_phase_and_mag(self):
        # self.component_phases = math.atan(self.component_imaginary/self.component_real)
        self.component_magnitude = [math.sqrt(self.component_real[i] ** 2 + self.component_imaginary[i] ** 2) for i in range(self.number_of_transform)]

    def fourier_synthesis(self):
        for i in range(self.number_of_data):
            self.synthesis_signal[i] = sum(
                [self.component_real[k] / self.number_of_data * math.cos(2 * math.pi * k * i / self.number_of_data) - self.component_imaginary[k] / self.number_of_data * math.sin(
                    2 * math.pi * k * i / self.number_of_data) for k in
                 range(self.number_of_transform)])

    def plot_signal(self):
        self.graph_engine.plot_signal([self.input_signal_vector, self.component_real, self.component_imaginary, self.synthesisize_signal], [-1, 1])

if __name__ == "__main__":
    input_vector_signal = np.zeros((100, 1), dtype=float)
    for i in range(100):
        input_vector_signal[i] = random.randrange(-1, 1)

    fourier_transform = FourierTransform(input_vector_signal)
    fourier_transform.fourier_transformation()
    fourier_transform.fourier_calculate_phase_and_mag()
    fourier_transform.fourier_synthesis()
    fourier_transform.plot_signal()


text:
import numpy as np
import unittest

import chainer
from chainer import testing
from chainer.testing import attr

from chainercv.links import SSD300
from chainercv.links import SSD512

@testing.parameterize(*testing.product({
    'insize': [300, 512],
    'n_fg_class': [1, 5, 20],
}))
class TestSSDVGG16(unittest.TestCase):

    def setUp(self):
        if self.insize == 300:
            self.link = SSD300(n_fg_class=self.n_fg_class)
            self.n_bbox = 8732
        elif self.insize == 512:
            self.link = SSD512(n_fg_class=self.n_fg_class)
            self.n_bbox = 24564

    def _check_call(self):
        x = self.link.xp.array(
            np.random.uniform(-1, 1, size=(1, 3, self.insize, self.insize)),
            dtype=np.float32)

        loc, conf = self.link(x)

        self.assertIsInstance(loc, chainer.Variable)
        self.assertIsInstance(loc.array, self.link.xp.ndarray)
        self.assertEqual(loc.shape, (1, self.n_bbox, 4))
        self.assertIsInstance(conf, chainer.Variable)
        self.assertIsInstance(conf.array, self.link.xp.ndarray)
        self.assertEqual(conf.shape, (1, self.n_bbox, self.n_fg_class + 1))

    @attr.slow
    def test_call_cpu(self):
        self._check_call()

    @attr.gpu
    @attr.slow
    def test_call_gpu(self):
        self.link.to_gpu()
        self._check_call()

@testing.parameterize(*testing.product({
    'model': [SSD300, SSD512],
    'n_fg_class': [None, 10, 20],
    'pretrained_model': ['voc0712', 'imagenet'],
}))
class TestSSDVGG16Pretrained(unittest.TestCase):

    @attr.slow
    def test_pretrained(self):
        kwargs = {
            'n_fg_class': self.n_fg_class,
            'pretrained_model': self.pretrained_model,
        }

        if self.pretrained_model == 'voc0712':
            valid = self.n_fg_class in {None, 20}
        elif self.pretrained_model == 'imagenet':
            valid = self.n_fg_class is not None

        if valid:
            self.model(**kwargs)
        else:
            with self.assertRaises(ValueError):
                self.model(**kwargs)

testing.run_module(__name__, __file__)

max_stars_repo_name: souravsingh/chainercv

max_issues_repo_name: souravsingh/chainercv

max_forks_repo_name: souravsingh/chainercv


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py:percent
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.13.7
#   kernelspec:
#     display_name: Python [conda env:bandit_38]
#     language: python
#     name: conda-env-bandit_38-py
# ---

# %% language="javascript"
# IPython.notebook.kernel.restart()

# %%
import datetime
import numpy as np
import pandas as pd
import re
import socket
import sys

from io import StringIO

from urllib.request import urlopen, Request
from urllib.error import HTTPError

# %%
workdir = '/Users/pnorton/tmp/streamflow_work'
stnfile = f'{workdir}/something'

base_url = 'http://waterservices.usgs.gov/nwis'

t1 = re.compile('^#.*$\n?', re.MULTILINE)   # remove comment lines
t2 = re.compile('^5s.*$\n?', re.MULTILINE)  # remove field length lines

# %%
# Set timeout in seconds - if not set defaults to infinite time for response
timeout = 30
socket.setdefaulttimeout(timeout)

# %%
# class WaterServices:
#     def __init__(self):
#         pass

#     def get():
#         # Get something from water services
#         response = urlopen(f'{BASE_NWIS_URL}/dv/{url_final}')
#         encoding = response.info().get_param('charset', failobj='utf8')
#         streamgage_obs_page = response.read().decode(encoding)

# %%
def nwis_load_daily_statistics(src_dir):
    col_names = ['agency_cd', 'site_no', 'parameter_cd', 'ts_id', 'loc_web_ds', 'month_nu', 'day_nu',
                 'begin_yr', 'end_yr', 'count_nu', 'mean_va']
    col_types = [np.str_, np.str_, np.str_, np.int, np.str_, np.int, np.int, np.int, np.int, np.int, np.float]
    cols = dict(zip(col_names, col_types))

    # Start with an empty dataframe
    nwis_daily = pd.DataFrame(columns=col_names)

    for region in range(1, 19):
        # region = '01'
        sys.stdout.write(f'\rRegion: {region:02}')
        sys.stdout.flush()
        # print(f'Region {region+1:02}')

        # Read the rdb file into a dataframe
        df = pd.read_csv(f'{src_dir}/conus_daily_HUC_{region:02}_obs.tab', sep='\t', dtype=cols)

        nwis_daily = nwis_daily.append(df, ignore_index=True)
    print('')
    return nwis_daily

# %%
nwis_daily = nwis_load_daily_statistics('/Users/pnorton/Projects/National_Hydrology_Model/datasets/streamflow')

# %%
nwis_cnt = nwis_daily.groupby(['site_no'])['loc_web_ds'].nunique()
# nwis_cnt = nwis_cnt[nwis_cnt < 2]

# %%
nwis_cnt['07311630']

# %%
nwis_daily[nwis_daily['site_no'] == '07311630']

# %%
nwis_daily

# %%

max_stars_repo_path: 10_NWIS/Retrieve daily streamflow.py

max_stars_repo_name: pnorton-usgs/notebooks

max_issues_repo_path: 10_NWIS/Retrieve daily streamflow.py

max_issues_repo_name: pnorton-usgs/notebooks

max_forks_repo_path: 10_NWIS/Retrieve daily streamflow.py

max_forks_repo_name: pnorton-usgs/notebooks


text:
import sys
import numpy as np
from model import *
import pickle

from utils import get_gen_data

if __name__ == "__main__":

    np.random.seed(2)

    try:
        # get an index as an argument
        idx = int(sys.argv[1])
    except:
        idx = 0

    n = 1000
    f_res = {'m1': [], 'm2': [], 'm3': [], 'y': []}
    obs_cov_res = {'obs': [], 'int0': [], 'int1': [], 'int2': []}
    int_cov_res = {'m1': [], 'm2': [], 'm3': []}
    for i in range(n):
        cov = np.array([[1, .7, .3, .4], [.7, 1, .5, .3], [.3, .5, 1, .1], [.4, .3, .1, 1]])
        f = {'m2': lambda m1: 2 * m1,
             'm3': lambda m1, m2: 1 * m1 - 1 * m2,
             'y': lambda m1, m2, m3: 1.5 * m1 - 2 * m2 - 1.5 * m3}
        gen_data = get_gen_data(cov, f)

        m1, m2, m3, y = gen_data(n)
        x1 = np.random.randn(n)
        x2 = np.random.randn(n) * 2
        x3 = np.random.randn(n) * 2
        do_m1, m2_dom1, m3_dom1, y_dom1 = gen_data(n, x1=x1)
        m1_dom2, do_m2, m3_dom2, y_dom2 = gen_data(n, x2=x2)
        m1_dom3, m2_dom3, do_m3, y_dom3 = gen_data(n, x3=x3)

        g = ig.Graph(directed=True)
        g.add_vertices([0, 1, 2, 3])
        g.add_edges([(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)])
        model = JointInterventionModel('test', g)
        obs_data = torch.cat([m1, m2, m3, y], dim=1)
        int_data = {1: torch.cat([m1_dom2, do_m2, m3_dom2, y_dom2], dim=1),
                    0: torch.cat([do_m1, m2_dom1, m3_dom1, y_dom1], dim=1),
                    2: torch.cat([m1_dom3, m2_dom3, do_m3, y_dom3], dim=1)}
        model.train(obs_data, int_data, outer_n_iter=10, inner_n_iter=1000, lr=0.01)

        for f, k in zip(model.F, f_res.keys()):
            f_res[k] += [f.squeeze().detach().numpy()]
        for cov, k in zip(model.int_covs.values(), int_cov_res.keys()):
            int_cov_res[k] += [cov.squeeze().detach().numpy()]

        obs_cov_res['obs'] += [model.obs_cov.squeeze().detach().numpy()]

    for k in f_res.keys():
        f_res[k] = np.stack(f_res[k])
    for k in int_cov_res.keys():
        int_cov_res[k] = np.stack(int_cov_res[k])
    obs_cov_res['obs'] = np.stack(obs_cov_res['obs'])

    # save the results
    with open("../_synthetic_results/params_res{}.pickle".format(idx), 'wb') as handle:
        pickle.dump(f_res, handle, protocol=pickle.HIGHEST_PROTOCOL)
        pickle.dump(int_cov_res, handle, protocol=pickle.HIGHEST_PROTOCOL)
        pickle.dump(obs_cov_res, handle, protocol=pickle.HIGHEST_PROTOCOL)