text:
# PPO-LSTM
import gym
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.distributions import Categorical

from envs.acrobot_simulator import AcrobotSimulator
from envs.acrobot_simulator_po import AcrobotSimulator_po

from torchdiffeq import odeint as odeint
import time
import numpy as np

device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

# Hyperparameters
learning_rate = 0.0005
gamma = 0.98
lmbda = 0.95
eps_clip = 0.1
K_epoch = 2
T_horizon = 20
ode_method = "explicit_adams"
atol = rtol = 1e-8

class ODEFunc(nn.Module):
    def __init__(self, input_dim, ode_dim=20):
        super(ODEFunc, self).__init__()
        self.l1 = nn.Linear(input_dim, ode_dim)
        self.l2 = nn.Linear(ode_dim, ode_dim)
        self.l3 = nn.Linear(ode_dim, input_dim)
        self.relu = nn.ReLU()

    def forward(self, t, x):
        """
        Perform one step in solving ODE.
        """
        x = self.l1(x)
        x = self.relu(x)
        x = self.l2(x)
        x = self.relu(x)
        x = self.l3(x)
        return x

class PPO(nn.Module):
    def __init__(self, action_dim=3, state_dim=4, latent_dim=32):
        super(PPO, self).__init__()
        self.data = []
        self.action_dim = action_dim

        self.fc1 = nn.Linear(state_dim, 64)
        self.lstm = nn.LSTM(64, latent_dim)
        self.fc_pi = nn.Linear(latent_dim, action_dim)
        self.fc_v = nn.Linear(latent_dim, 1)
        self.optimizer = optim.Adam(self.parameters(), lr=learning_rate)

        self.fc_xi = nn.Linear(latent_dim * 2 + action_dim, latent_dim * 2)
        self.ode_func = ODEFunc(latent_dim * 2, ode_dim=20)

    def xi(self, hidden_tild, a):
        a = torch.tensor(a).to(device)
        a = F.one_hot(a, num_classes=self.action_dim)
        input = torch.cat((hidden_tild[0].flatten(), hidden_tild[1].flatten(), a), dim=-1) # cell state, hidden_state
        hidden_prime = self.fc_xi(input)
        return hidden_prime

    def pi(self, obs, hidden):
        obs = F.relu(self.fc1(obs))
        obs = obs.view(-1, 1, 64)
        obs, hidden_tild = self.lstm(obs, hidden)
        obs = self.fc_pi(obs)
        prob = F.softmax(obs, dim=2)
        return prob, hidden_tild

    def v(self, obs, hidden):
        obs = F.relu(self.fc1(obs))
        obs = obs.view(-1, 1, 64)
        obs, lstm_hidden = self.lstm(obs, hidden)
        v = self.fc_v(obs)
        return v

    def put_data(self, transition):
        self.data.append(transition)

    def make_batch(self):
        s_lst, a_lst, r_lst, s_prime_lst, prob_a_lst, h_in_lst, h_out_lst, done_lst = [], [], [], [], [], [], [], []
        for transition in self.data:
            s, a, r, s_prime, prob_a, h_in, h_out, done = transition

            s_lst.append(s)
            a_lst.append([a])
            r_lst.append([r])
            s_prime_lst.append(s_prime)
            prob_a_lst.append([prob_a])
            h_in_lst.append(h_in)
            h_out_lst.append(h_out)
            done_mask = 0 if done else 1
            done_lst.append([done_mask])

        s, a, r, s_prime, done_mask, prob_a = torch.tensor(s_lst, dtype=torch.float), torch.tensor(a_lst), \
                                              torch.tensor(r_lst), torch.tensor(s_prime_lst, dtype=torch.float), \
                                              torch.tensor(done_lst, dtype=torch.float), torch.tensor(prob_a_lst)
        self.data = []
        return s, a, r, s_prime, done_mask, prob_a, h_in_lst[0], h_out_lst[0]

    def train_net(self):
        s, a, r, s_prime, done_mask, prob_a, (h1_in, h2_in), (h1_out, h2_out) = self.make_batch()
        s = torch.FloatTensor(s).to(device)
        r = r.float().to(device)
        a = a.to(device)
        s_prime = torch.FloatTensor(s_prime).to(device)
        done_mask = torch.FloatTensor(done_mask).to(device)
        prob_a = torch.FloatTensor(prob_a).to(device)

        first_hidden  = (h1_in.detach(), h2_in.detach())
        second_hidden = (h1_out.detach(), h2_out.detach())

        for i in range(K_epoch):
            v_prime = self.v(s_prime, second_hidden).squeeze(1)
            td_target = r + gamma * v_prime * done_mask
            v_s = self.v(s, first_hidden).squeeze(1)
            delta = td_target - v_s
            delta = delta.detach().cpu().numpy()

            advantage_lst = []
            advantage = 0.0
            for item in delta[::-1]:
                advantage = gamma * lmbda * advantage + item[0]
                advantage_lst.append([advantage])
            advantage_lst.reverse()
            advantage = torch.tensor(advantage_lst, dtype=torch.float).to(device)

            pi, _ = self.pi(s, first_hidden)
            pi_a = pi.squeeze(1).gather(1, a)
            ratio = torch.exp(torch.log(pi_a) - torch.log(prob_a))  # a/b == log(exp(a)-exp(b))

            # import pdb; pdb.set_trace()
            surr1 = ratio * advantage
            surr2 = torch.clamp(ratio, 1 - eps_clip, 1 + eps_clip) * advantage
            loss = -torch.min(surr1, surr2) + F.smooth_l1_loss(v_s, td_target.detach())

            self.optimizer.zero_grad()
            loss.mean().backward(retain_graph=True)
            self.optimizer.step()

def main():
    exp_name = "lstmppo_acrobot_po_continuous"
    env = AcrobotSimulator_po(continuous_time=True)
    model = PPO(action_dim=3, state_dim=2).to(device)

    score = 0.0
    print_interval = 10
    results = []

    for n_epi in range(10000):
        z = (torch.zeros([1, 1, 32], dtype=torch.float).to(device), torch.zeros([1, 1, 32], dtype=torch.float).to(device))
        s = env.reset()
        done = False
        t = 0
        while not done:
            for t in range(T_horizon):
                h_in = z
                prob, z_tild = model.pi(torch.from_numpy(s).float().to(device), h_in)
                prob = prob.view(-1)
                m = Categorical(prob)
                a = m.sample().item()
                s_prime, r, done, info = env.step(a)
                prev_t = t
                t = prev_t + info["dt"]
                z_prime = model.xi(z_tild, a)
                time = torch.FloatTensor([prev_t, t]).to(device)
                z = odeint(model.ode_func, z_prime, time, rtol=rtol, atol=atol, method=ode_method)[-1]
                z = torch.split(z.reshape(1, 1, -1), split_size_or_sections=32, dim=-1)
                # import pdb;pdb.set_trace()
                model.put_data((s, a, r / 100.0, s_prime, prob[a].item(), h_in, z_tild, done))
                s = s_prime

                score += r
                if done:
                    break

            model.train_net()

        if n_epi % 10 == 0 and n_epi != 0:
            print("# of episode :{}, avg score : {:.1f}".format(n_epi, score / print_interval))
            results.append([n_epi, score / print_interval])
            np.save(exp_name, np.array(results))
            score = 0.0

    env.close()

if __name__ == '__main__':
    main()


text:
import numpy as np
import os
import sys
sys.path.append(os.path.abspath(r'../lib'))
import NumCpp  # noqa E402

np.random.seed(666)

####################################################################################
def test_num2str():
    """Tests the NumCpp Utils"""
    value = np.random.randint(1, 100, [1, ], dtype=np.int8).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.int16).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.int32).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.int64).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.uint8).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.uint16).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.uint32).item()
    assert NumCpp.num2str(value) == str(value)

    value = np.random.randint(1, 100, [1, ], dtype=np.uint64).item()
    assert NumCpp.num2str(value) == str(value)

####################################################################################
def test_sqr():
    value = np.random.randint(1, 12, [1, ], dtype=np.int8).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.int16).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.int32).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.int64).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 15, [1, ], dtype=np.uint8).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.uint16).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.uint32).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ], dtype=np.uint64).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ]).astype(np.double).item()
    assert NumCpp.sqr(value) == value ** 2

    value = np.random.randint(1, 100, [1, ]).astype(np.float32).item()
    assert NumCpp.sqr(value) == value ** 2

####################################################################################
def test_cube():
    value = np.random.randint(1, 6, [1, ], dtype=np.int8).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 32, [1, ], dtype=np.int16).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ], dtype=np.int32).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ], dtype=np.int64).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 7, [1, ], dtype=np.uint8).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 41, [1, ], dtype=np.uint16).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ], dtype=np.uint32).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ], dtype=np.uint64).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ]).astype(np.double).item()
    assert NumCpp.cube(value) == value ** 3

    value = np.random.randint(1, 100, [1, ]).astype(np.float32).item()
    assert NumCpp.cube(value) == value ** 3

####################################################################################
def test_power():
    value = np.random.randint(1, 4, [1, ], dtype=np.int8).item()
    power = np.random.randint(1, 4, dtype=np.uint8).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int16).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int32).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int64).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint8).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint16).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint32).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint64).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ]).astype(np.double).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ]).astype(np.float32).item()
    assert NumCpp.power(value, power) == np.power(value, power)

    real = np.random.rand(1).astype(np.double).item()
    imag = np.random.rand(1).astype(np.double).item()
    value = np.complex128(complex(real, imag))
    assert np.round(NumCpp.power(value, power), 5) == np.round(np.power(value, power), 5)

    real = np.random.rand(1).astype(np.double).item()
    imag = np.random.rand(1).astype(np.double).item()
    value = np.complex64(complex(real, imag))
    assert np.round(NumCpp.power(value, power), 5) == np.round(np.power(value, power), 5)

####################################################################################
def test_powerf():
    value = np.random.randint(1, 5, [1, ], dtype=np.int8).item()
    power = np.random.rand(1).item() * 5
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int16).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int32).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.int64).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint8).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint16).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint32).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.randint(1, 5, [1, ], dtype=np.uint64).item()
    assert NumCpp.powerf(value, power) == np.power(value, power)

    value = np.random.rand(1).astype(np.double).item()
    assert np.round(NumCpp.powerf(value, power), 5) == \
        np.round(np.power(value, power), 5)

    value = np.random.rand(1).astype(np.float32).item()
    assert np.round(NumCpp.powerf(value, np.float32(power)), 5) == \
        np.round(np.power(value, np.float32(power)), 5)

    realPower = np.random.rand(1).astype(np.double).item()
    imagPower = np.random.rand(1).astype(np.double).item()
    power = np.complex128(complex(realPower, imagPower))
    real = np.random.rand(1).astype(np.double).item()
    imag = np.random.rand(1).astype(np.double).item()
    value = np.complex128(complex(real, imag))
    assert np.round(NumCpp.powerf_complex(value, power), 5) == np.round(np.float_power(value, power), 5)

    realPower = np.random.rand(1).astype(np.double).item()
    imagPower = np.random.rand(1).astype(np.double).item()
    power = np.complex64(complex(realPower, imagPower))
    real = np.random.rand(1).astype(np.double).item()
    imag = np.random.rand(1).astype(np.double).item()
    value = np.complex64(complex(real, imag))
    assert np.round(NumCpp.powerf_complex(value, power), 5) == np.round(np.float_power(value, power), 5)

ext: py

max_stars_repo_name: leannejdong/NumCpp

max_issues_repo_name: leannejdong/NumCpp

max_forks_repo_name: leannejdong/NumCpp

reason: import numpy


text:
from pb2 import common_pb2
import matplotlib.cm as cm
import matplotlib
import os
import sys
from bh_tsne import bhtsne
matplotlib.use('Agg')
from matplotlib import pyplot as plt
matplotlib.rcParams.update({'font.size': 17})
matplotlib.rcParams.update({'pdf.fonttype': 42})
matplotlib.rcParams.update({'ps.fonttype': 42})
import numpy as np

def make_N_colors(cmap_name, N):
  cmap = cm.get_cmap(cmap_name, N)
  return cmap(np.arange(N))

def plot2d(x, y, label, fname):
  plt.clf()
  color=make_N_colors('gist_rainbow', 10)
  marker=["o","x","d","s","+","v",">","p",'*','.']
  count=[0]*10
  for i,lb in enumerate(label):
    plt.scatter(x[i],y[i], marker=marker[lb], s=40, c='w', edgecolor=color[lb])
    count[lb]+=1
  #print count
  #plt.show()
  plt.savefig(fname)
  plt.clf()

def plot_all_feature(infile, outfolder):
  fd = open(infile, 'rb')
  bps = common_pb2.BlobProtos()
  bps.ParseFromString(fd.read())
  labels = None
  outprefix = os.path.join(outfolder, os.path.splitext(os.path.split(infile)[1])[0])
  for (name, blob) in zip(bps.name, bps.blob):
    if 'label' in name:
      assert blob.shape[0] == len(blob.data)
      labels = np.asarray(blob.data, dtype = np.int)
  for (name, blob) in zip(bps.name, bps.blob):
    if len(blob.shape) >1 :
      s = (blob.shape[0], len(blob.data) / blob.shape[0])
      fea = np.asarray(blob.data, dtype = np.float32).reshape(s)
      X = []
      Y = []
      if fea.shape[1] == 2:
        plot2d(fea[:,0], fea[:,1], labels, outprefix + '-' + name + '.png')
      else:
        for point in bhtsne.bh_tsne(fea):
          X.append(point[0])
          Y.append(point[1])
        plot2d(X, Y, labels, outprefix + '-' + name + '.png')

if __name__ == "__main__":
  if len(sys.argv) != 3:
    print 'Usage: python plot.py <path> <out_folder>'
    print 'the program use <path.dat> as input to generate a picture at <out_folder/>'
    sys.exit()

  plot_all_feature(sys.argv[1], sys.argv[2])

max_stars_repo_name: chonho/singa-web

max_issues_repo_name: chonho/singa-web

max_forks_repo_name: chonho/singa-web


text:
import matplotlib.pyplot as plt
import math as m
import sympy as sp

x = sp.Symbol('x')

"""Now that we just need to make the differentiation to work"""

# x - distance along beam (m)
# y - deflection (m)
# theta(x) - slope(m/m)
# E - modulus of elasticity (Pa = N/m**2)
# I - Moment of inertia (m**4)
# M(x) - moment (Nm)
# V(x) - shear (N)
# w(x) - distributed load (N/m)

E = 200
I = 0.0003
w_0 = 2.5
y_0 = 0
y_L = 0
x = 1
dx = 0.125
L = 3

def theta_x(L):
    return w_0/(120*E*I*L)*(-5*x**4+6*L**2*x**2 - L**4)

def M_x():
    return w_0/(120*L)*(-20*x*3+6*L**2*2*x)

def V_x():
    return w_0/(120*L)*(-60*x*2+6*L**2*2)

def w_x():
    return -(w_0/(120*L)*(-120*x))

V_x_values = []
theta_values = []
M_x_values = []
w_x_values = []

while L>=0:
    if L == 0:
        break
    w_x_values.append(w_x())
    M_x_values.append(M_x())
    theta_values.append(theta_x(L))
    V_x_values.append(V_x())
    L +=-dx
    print(L)

print("Theta values")
print(theta_values)
print("Moment values")
print(M_x_values)
print("Shear Values")
print(V_x_values)
print("Distributed Load")
print(w_x_values)
# dy/dx = theta(x)
# dtheta/dx = M(x)/(EI)
# dM/dx = v(x)
# dV/dx = -w(x)

# plt.plot(x,slope)
# plt.show()

reason: import sympy


text:
import numpy as np
import pytest
from spector import indices, vector

def test_basic():
    ind = indices()
    assert str(ind) == 'indices([])'
    assert len(ind) == 0
    assert 0 not in ind
    assert list(ind) == []
    assert indices(np.array([False]))
    assert indices(vector([0]))
    with pytest.raises(TypeError):
        indices(np.array([0.0]))

    assert ind.add(0)
    assert not ind.add(0)
    assert str(ind) == 'indices([0])'
    assert len(ind) == 1
    assert 0 in ind
    assert list(ind) == [0]

    assert ind.discard(0)
    assert not ind.discard(0)
    ind.update(range(3))
    assert set(ind) == set(indices(np.array(ind))) == {0, 1, 2}
    np.array(ind).sum() == 3
    ind.update(indices([3]))
    assert 3 in ind
    ind.clear()
    assert not ind

def test_cmp():
    ind = indices([0, 1])
    assert ind == ind == indices([0, 1]) == indices([1, 0])
    assert ind != {} and not (ind == {})
    assert not (ind != ind)
    assert ind != indices([0, 2])
    assert not (ind <= indices([0]))
    assert ind <= ind
    assert not ind < ind
    assert ind < indices([0, 1, 2])
    with pytest.raises(TypeError):
        ind < None
    with pytest.raises(TypeError):
        ind <= None
    assert not ind.isdisjoint(indices([0])) and ind.isdisjoint(indices([2]))
    assert ind @ indices([0]) == 1
    assert ind @ indices(range(3)) == 2

def test_sets():
    x, y = indices([0, 1]), indices([1, 2])
    assert x | y == indices([0, 1, 2])
    assert x & y == indices([1])
    assert x - y == indices([0])
    assert x ^ y == indices([0, 2])
    assert x.union() == x and x.union() is not x
    assert x.union([1], [2]) == indices([0, 1, 2])
    assert x.intersection() == x and x.intersection() is not x
    assert x.intersection(y, [1]) == indices([1])
    assert x.difference() == x and x.difference() is not x
    assert x.difference(y, [3]) == indices([0])
    assert x.dot() == 2
    assert x.dot(list(y)) == x.dot(y) == x.dot(y, y) == 1

    z = x
    z ^= y
    assert z is x and z == indices([0, 2])
    z |= y
    assert z is x and z == indices([0, 1, 2])
    z -= y
    z -= y  # smaller vector is iterated
    assert z is x and z == indices([0])
    z &= y
    assert z is x and z == indices()
    assert not y & z

def test_dense():
    arr = np.zeros(4, bool)
    arr[0] = arr[2] = True
    ind = indices.fromdense(arr)
    assert set(ind) == {0, 2}
    assert np.array_equal(ind.todense(), arr[:3])
    assert np.array_equal(ind.todense(4), arr)

max_stars_repo_name: coady/spector

max_issues_repo_name: coady/spector

max_forks_repo_name: coady/spector


text:
import numpy as np
import networkx as nx

with open("day22-input.txt") as f:
    depth = int(f.readline().split(" ")[1])
    x, y = (int(n) for n in f.readline().split(" ")[1].split(','))

xdim, ydim = int(x*3), int(y*3)

cave = np.empty((xdim, ydim), dtype=int)

xx, yy = np.meshgrid(np.arange(xdim), np.arange(ydim), indexing="ij")

cave[xx == 0] = (yy[xx == 0] * 48271 + depth) % 20183
cave[yy == 0] = (xx[yy == 0] * 16807 + depth) % 20183

for i in range(1, xdim):
    for j in range(1, ydim):
        cave[i, j] = ((cave[i-1, j] % 20183) * (cave[i, j-1] % 20183) + depth) % 20183

cave[x, y] = depth

cave_type = cave % 3

print(f"Answer part  I: {cave_type[:x+1, :y+1].sum()}")

grid = nx.Graph()

common_states = {
    2: [0, 1],
    3: [1],
    4: [1, 2],
    5: [0],
    6: [2],
    8: [0, 2]
    }

for (r, c), t in np.ndenumerate(cave_type):
    layers = [(r, c, s) for s in common_states[2 * 2**t]]
    grid.add_edge(*layers, weight=7)
    if r < (xdim - 1):
        t2 = cave_type[r+1, c]
        for s in common_states[2**t + 2**t2]:
            grid.add_edge((r, c, s), (r+1, c, s), weight=1)
    if c < (ydim - 1):
        t2 = cave_type[r, c+1]
        for s in common_states[2**t + 2**t2]:
            grid.add_edge((r, c, s), (r, c+1, s), weight=1)

print(f"Answer part II: {nx.shortest_path_length(grid, (0, 0, 0), (x, y, 0), weight='weight')}")

reason: import numpy,import networkx


text:
import logging
from itertools import product

import numpy as np
import pymc3 as pm
import theano
import theano.tensor as tt
from pymc3.variational.callbacks import CheckParametersConvergence

import csrank.numpy_util as npu
import csrank.theano_util as ttu
from csrank.learner import Learner
from csrank.util import print_dictionary
from .discrete_choice import DiscreteObjectChooser
from .likelihoods import likelihood_dict, LogLikelihood, create_weight_dictionary, fit_pymc3_model

class MixedLogitModel(DiscreteObjectChooser, Learner):
    def __init__(self, n_object_features, n_mixtures=4, loss_function='', regularization='l2', **kwargs):
        """
            Create an instance of the Mixed Logit model for learning the discrete choice function. In this model we
            assume weights of this model to be random due to which this model can learn different variations in choices
            amongst the individuals. The utility score for each object in query set :math:`Q` is defined as
            :math:`U_r(x) = w_r \\cdot x`, where :math:`w_r` is the k-th sample weight vector from the underlying distribution
            The probability of choosing an object :math:`x_i` is defined by taking softmax over the
            utility scores of the objects:

            .. math::

                P(x_i \\lvert Q) = \\frac{1}{R} \\sum_{r=1}^R \\frac{exp(U_r(x_i))}{\\sum_{x_j \\in Q} exp(U_r(x_j))}

            The discrete choice for the given query set :math:`Q` is defined as:

            .. math::

                dc(Q) := \\operatorname{argmax}_{x_i \\in Q }  \\; P(x_i \\lvert Q)

            Parameters
            ----------
            n_object_features : int
                Number of features of the object space
            n_mixtures: int (range : [2, inf])
                The number of logit models (:math:`R`) which are used to estimate the choice probability
            loss_function : string , {‘categorical_crossentropy’, ‘binary_crossentropy’, ’categorical_hinge’}
                Loss function to be used for the discrete choice decision from the query set
            regularization : string, {‘l1’, ‘l2’}, string
               Regularizer function (L1 or L2) applied to the `kernel` weights matrix
            random_state : int or object
                Numpy random state
            **kwargs
                Keyword arguments for the algorithms

            References
            ----------
                [1] Kenneth E Train. „Discrete choice methods with simulation“. In: Cambridge university press, 2009. Chap Mixed Logit, pp. 153–172.

                [2] Kenneth Train. Qualitative choice analysis. Cambridge, MA: MIT Press, 1986

                [3] Daniel McFadden and Kenneth Train. „Mixed MNL models for discrete response“. In: Journal of applied Econometrics 15.5 (2000), pp. 447–470
        """
        self.logger = logging.getLogger(MixedLogitModel.__name__)
        self.n_object_features = n_object_features
        self.loss_function = likelihood_dict.get(loss_function, None)
        if regularization in ['l1', 'l2']:
            self.regularization = regularization
        else:
            self.regularization = 'l2'
        self._config = None
        self.n_mixtures = n_mixtures
        self.model = None
        self.trace = None
        self.trace_vi = None
        self.Xt = None
        self.Yt = None
        self.p = None

    @property
    def model_configuration(self):
        """
            Constructs the dictionary containing the priors for the weight vectors for the model according to the
            regularization function. The parameters are:
                * **weights** : Distribution of the weigh vectors to evaluates the utility of the objects

            For ``l1`` regularization the priors are:

            .. math::

                \\text{mu}_w \\sim \\text{Normal}(\\text{mu}=0, \\text{sd}=5.0) \\\\
                \\text{b}_w \\sim \\text{HalfCauchy}(\\beta=1.0) \\\\
                \\text{weights} \\sim \\text{Laplace}(\\text{mu}=\\text{mu}_w, \\text{b}=\\text{b}_w)

            For ``l2`` regularization the priors are:

            .. math::

                \\text{mu}_w \\sim \\text{Normal}(\\text{mu}=0, \\text{sd}=5.0) \\\\
                \\text{sd}_w \\sim \\text{HalfCauchy}(\\beta=1.0) \\\\
                \\text{weights} \\sim \\text{Normal}(\\text{mu}=\\text{mu}_w, \\text{sd}=\\text{sd}_w)
        """
        if self._config is None:
            if self.regularization == 'l2':
                weight = pm.Normal
                prior = 'sd'
            elif self.regularization == 'l1':
                weight = pm.Laplace
                prior = 'b'
            self._config = {
                'weights': [weight, {'mu': (pm.Normal, {'mu': 0, 'sd': 5}), prior: (pm.HalfCauchy, {'beta': 1})}]}
            self.logger.info('Creating model with config {}'.format(print_dictionary(self._config)))
        return self._config

    def construct_model(self, X, Y):
        """
            Constructs the mixed logit model by applying priors on weight vectors **weights** as per
            :meth:`model_configuration`. The probability of choosing the object :math:`x_i` from the query set
            :math:`Q = \\{x_1, \\ldots ,x_n\\}` assuming we draw :math:`R` samples of the weight vectors is:

            .. math::

                P(x_i \\lvert Q) = \\frac{1}{R} \\sum_{r=1}^R \\frac{exp(U_r(x_i))}{\\sum_{x_j \\in Q} exp(U_r(x_j))}

            Parameters
            ----------
            X : numpy array
                (n_instances, n_objects, n_features)
                Feature vectors of the objects
            Y : numpy array
                (n_instances, n_objects)
                Preferences in the form of discrete choices for given objects

            Returns
            -------
             model : pymc3 Model :class:`pm.Model`
        """
        with pm.Model() as self.model:
            self.Xt = theano.shared(X)
            self.Yt = theano.shared(Y)
            shapes = {'weights': (self.n_object_features, self.n_mixtures)}
            weights_dict = create_weight_dictionary(self.model_configuration, shapes)
            utility = tt.dot(self.Xt, weights_dict['weights'])
            self.p = tt.mean(ttu.softmax(utility, axis=1), axis=2)
            yl = LogLikelihood('yl', loss_func=self.loss_function, p=self.p, observed=self.Yt)
        self.logger.info("Model construction completed")

    def fit(self, X, Y, sampler='variational', tune=500, draws=500,
            vi_params={"n": 20000, "method": "advi", "callbacks": [CheckParametersConvergence()]}, **kwargs):
        """
            Fit a mixed logit model on the provided set of queries X and choices Y of those objects. The provided
            queries and corresponding preferences are of a fixed size (numpy arrays). For learning this network
            the categorical cross entropy loss function for each object :math:`x_i \\in Q` is defined as:

            .. math::

                C_{i} =  -y(i)\\log(P_i) \\enspace,

            where :math:`y` is ground-truth discrete choice vector of the objects in the given query set :math:`Q`.
            The value :math:`y(i) = 1` if object :math:`x_i` is chosen else :math:`y(i) = 0`.

            Parameters
            ----------
            X : numpy array (n_instances, n_objects, n_features)
                Feature vectors of the objects
            Y : numpy array (n_instances, n_objects)
                Choices for given objects in the query
            sampler : {‘variational’, ‘metropolis’, ‘nuts’}, string
                The sampler used to estimate the posterior mean and mass matrix from the trace

                    * **variational** : Run inference methods to estimate posterior mean and diagonal mass matrix
                    * **metropolis** : Use the MAP as starting point and Metropolis-Hastings sampler
                    * **nuts** : Use the No-U-Turn sampler
            vi_params : dict
                The parameters for the **variational** inference method
            draws : int
                The number of samples to draw. Defaults to 500. The number of tuned samples are discarded by default
            tune : int
                Number of iterations to tune, defaults to 500. Ignored when using 'SMC'. Samplers adjust
                the step sizes, scalings or similar during tuning. Tuning samples will be drawn in addition
                to the number specified in the `draws` argument, and will be discarded unless
                `discard_tuned_samples` is set to False.
            **kwargs :
                Keyword arguments for the fit function of :meth:`pymc3.fit`or :meth:`pymc3.sample`
        """
        self.construct_model(X, Y)
        fit_pymc3_model(self, sampler, draws, tune, vi_params, **kwargs)

    def _predict_scores_fixed(self, X, **kwargs):
        summary = dict(pm.summary(self.trace)['mean'])
        weights = np.zeros((self.n_object_features, self.n_mixtures))
        for i, k in product(range(self.n_object_features), range(self.n_mixtures)):
            weights[i][k] = summary['weights[{},{}]'.format(i, k)]
        utility = np.dot(X, weights)
        p = np.mean(npu.softmax(utility, axis=1), axis=2)
        return p

    def predict(self, X, **kwargs):
        return super().predict(X, **kwargs)

    def predict_scores(self, X, **kwargs):
        return super().predict_scores(X, **kwargs)

    def predict_for_scores(self, scores, **kwargs):
        return DiscreteObjectChooser.predict_for_scores(self, scores, **kwargs)

    def set_tunable_parameters(self, n_mixtures=4, loss_function='', regularization="l1", **point):
        """
            Set tunable parameters of the Mixed Logit model to the values provided.

            Parameters
            ----------
            n_mixtures: int (range : [2, inf])
                The number of logit models (:math:`R`) which are used to estimate the choice probability
            loss_function : string , {‘categorical_crossentropy’, ‘binary_crossentropy’, ’categorical_hinge’}
                Loss function to be used for the discrete choice decision from the query set
            regularization : string, {‘l1’, ‘l2’}, string
               Regularizer function (L1 or L2) applied to the `kernel` weights matrix
            point: dict
                Dictionary containing parameter values which are not tuned for the network
        """
        if loss_function in likelihood_dict.keys():
            self.loss_function = likelihood_dict.get(loss_function, None)
        self.n_mixtures = n_mixtures
        self.regularization = regularization
        self.model = None
        self.trace = None
        self.trace_vi = None
        self.Xt = None
        self.Yt = None
        self.p = None
        self._config = None
        if len(point) > 0:
            self.logger.warning('This ranking algorithm does not support'
                                ' tunable parameters called: {}'.format(print_dictionary(point)))

max_stars_repo_name: helegraf/cs-ranking

max_issues_repo_name: helegraf/cs-ranking

max_forks_repo_name: helegraf/cs-ranking

reason: import numpy,import theano,import pymc3,from pymc3


text:
import numpy as np

import sys
sys.path.append('./GMRES_API')
import GMRES
sys.path.append('./RestartAlgorithm_API')
import RestartAlgorithm

from matplotlib import pyplot as plt

def GMRES_test( A_mat, b_mat, x_mat, iterations_between_restarts,restarting_iteration_ending_threshold ):
    print("\n")

    # The restatrt algorithm of GMRES
    #=========================================================================================================
    restarted_GMRES = RestartAlgorithm.RestartAlgorithm()
    restarted_GMRES.kernel_algorithm_register( GMRES.GMRES_API( A_mat, b_mat, iterations_between_restarts ) )
    restarted_GMRES.restart_initial_input( x_mat )
    restarted_GMRES.maximum_restarting_iteration_register( 500 )
    restarted_GMRES.restarting_iteration_ending_threshold_register( restarting_iteration_ending_threshold )
    # x_final, r_trend = restarted_GMRES.run_restart()
    return  restarted_GMRES.run_restart()
    #=========================================================================================================

def main():

    A_mat = np.array( [[4.00, 1.00, 1.00, 1.11],
                       [1.50, 3.00, 1.00, 1.73],
                       [1.50, 2.00, 3.00, 1.17],
                       [0.30, 0.50, 3.00, 3.00]] )

    b_mat = np.array( [3.0, 2.0, 1.0, 1.5] )

    x_mat = np.array( [1.0, 1.0, 1.0, 1.0] )
    print("x  =", x_mat)

    # The phenomena described in
    # M. Embree
    # The tortoise and the hare restart GMRES
    # SIAM Review, 45 (2) (2003), pp. 256-266
    #==========================================================================
    # GMRES with restart, 1 iterations in each restart ( GMRES(1) )
    x_final_1, r_trend_1 = GMRES_test( A_mat, b_mat, x_mat, 1, 1.0e-14 )
    # Converged, the fastest at first

    # GMRES with restart, 2 iterations in each restart ( GMRES(2) )
    x_final_2, r_trend_2 = GMRES_test( A_mat, b_mat, x_mat, 2, 1.0e-14 )
    # Converged, but not the fastest at first, faster later.

    # GMRES with restart, 3 iterations in each restart ( GMRES(3) )
    #x_final_3, r_trend_3 = GMRES_test( A_mat, b_mat, x_mat, 3, 1.0e-10 )
    # This will explode if set the threshold too small.

    # GMRES with restart, 4 iterations in each restart ( GMRES(4) )
    #GMRES_test_itr2 = GMRES.GMRES_API( A_mat, b_mat, 4) # This will explode
    #x_final_4, r_trend_4 = GMRES_test( A_mat, b_mat, x_mat, 4, 1.0e-10 )
    #==========================================================================

    # Draw the residual trend by the sequence of restarts
    #====================================================
    plt.title("restarted_GMRES_residual_trend")
    plt.xlabel("restart")
    plt.ylabel("residual")
    min_restart_shown =0    #0   #18
    max_restart_shown =50   #500 #23
    plt.plot(r_trend_1[min_restart_shown:max_restart_shown])
    plt.plot(r_trend_2[min_restart_shown:max_restart_shown])
    #plt.plot(r_trend_3[min_restart_shown:max_restart_shown])
    #plt.plot(r_trend_4[min_restart_shown:max_restart_shown])
    plt.show()
    #====================================================

    xx = np.matmul( np.linalg.inv(A_mat), b_mat )
    print("ANS : xx =", xx)

if __name__ == '__main__':
    main()


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

import numpy as np
from datetime import date

''' RINEX NAWIGACYJNY'''

def readrnxnav(file):
    m=1
    nav=np.zeros((2000,37))
    inav=np.zeros((2000))
    n=-1
    with open(file, "r") as f:
        for s in f:
            answer = s.find('END OF HEADER') # skip header
            if answer != -1:
                break
        for s in f:
            s = s.replace('D', 'E')
            if m==1:
                prn=int(s2n(s,1,2))
                a = np.empty((1,6))
                a[:] = np.NaN
                a[0,0:6]=np.array(s2e(s,4,23))
            else:
                a=np.append(a,s2n(s,4,19))
            for x in range (3):
                p=23+x*19
                a=np.append(a,s2n(s,p,19))
            if m<8:
                m+=1
            else:
                n+=1
                nav[n,:]=a
                inav[n]=prn
                m=1
        nav=nav[0:n+1,:]
        inav=inav[0:n+1]
        inav = inav.astype(int)
    f.close()
    return nav, inav

def readrnxobs(file, time_start, time_end, GNSS):
    with open(file, "r") as f:
        for s in f:
            label = s[59:]
            if label.find('SYS / # / OBS TYPES') == 1:
                if s[0] == GNSS:
                    p = 7
                    types_header = []
                    for i in range(int(s[4:4+2])):
                        if p > 58:
                            p = 7
                            s = next(f)
                        types_header.append(s[p:p+3])
                        p += 4

            elif label.find('END OF HEADER') == 1:
                break
            types_of_obs = ['C1C']
        ind = np.zeros((len(types_header)))
        for n in range(len(types_of_obs)):
            i=(types_header.index(types_of_obs[n])) if types_of_obs[n] in types_header else -1#np.empty((0))
            if i>-1:
                ind[i]=n+1

        obs = np.zeros((150000, len(types_of_obs)))*np.nan
        iobs = np.zeros((150000, 3))
        n=0
        for s in f:
            label = s[0]
            if label == '>':
                epoch = s2e(s,2,29)
                y = epoch[0]
                # tt = (date.toordinal(date(epoch[0],epoch[1],epoch[2]))+366-t0)*86400+np.dot((epoch[3:6]), ([3600,60,1])) + 6*86400
                tt = date2tow(epoch)[1] - date2tow(epoch)[2] * 86400
                if tt > (date2tow(time_end)[1] - date2tow(epoch)[2] * 86400):
                    break
                else:
                    flag = int(round(tt))>=(date2tow(time_start)[1] - date2tow(time_start)[2] * 86400)
                if flag:
                    number_of_all_sats = int(s[33:33+2])
                    iobs[n+np.arange(0,number_of_all_sats),1] = tt
                    iobs[n+np.arange(0,number_of_all_sats),2] = date2tow(epoch)[1]
                    for sat in range(number_of_all_sats):
                        s = next(f)
                        p = 3
                        if s[0] == GNSS:
                            for i in range(len(types_header)):
                                if ind[i] != 0:
                                    obs[n+sat, int(ind[i] - 1)] = s2n(s,p,15)
                                    iobs[n+sat,0] = s2n(s,1,2)
                                p+=16
                    n += number_of_all_sats
        obs = obs[0:n, :]
        iobs = iobs[0:n,:]
        obs = np.delete(obs,iobs[:,0]==0, axis=0)
        iobs = np.delete(iobs,iobs[:,0]==0, axis=0)
        f.close()
        iobs = iobs.astype(int)
    return obs, iobs

def s2e(s,p,n):
    epoch = [int(s[p:p+4]), int(s[p+5:p+5+2]), int(s[p+8:p+8+2]), int(s[p+11:p+11+2]), int(s[p+14:p+14+2]), float(s[p+17:n])]
    return epoch

def date2tow(data):
    dday=date.toordinal(date(data[0],data[1],data[2])) - (date.toordinal(date(1980,1,6)))
    week = dday//7
    dow = dday%7
    tow = dow * 86400 + data[3] * 3600 + data[4] * 60 + data[5]
    return week, tow, dow

def s2n(s,p,n):
    a = s[p:p+n]
    if (not (a and not a.isspace())):
        a = np.nan
    else:
        a = float(a)
    return a

max_stars_repo_path: Single Point Positioning/readrnx.py

max_issues_repo_path: Single Point Positioning/readrnx.py

max_forks_repo_path: Single Point Positioning/readrnx.py


text:
import h5py
import numpy as np

# Read hdf5 file
hdf = h5py.File('config/model.h5','r')
# Get keys of the hdf5 file
data = list(hdf.keys())
# Read data in each group
for key in data:
    # Get the values
    data = list(hdf[key].values())
    weights = list(data[0]['kernel:0'])
    # Convert to numpy array
    weights = np.array(weights)
    # Print shape of array to see mapping of the layers
    # Format: (input, output) for each layer
    print(weights.shape)

max_stars_repo_name: nikhilraghava/diabetes-risk-predictor

max_issues_repo_name: nikhilraghava/diabetes-risk-predictor

max_forks_repo_name: nikhilraghava/diabetes-risk-predictor


text:
import numpy as np

def LUdec(matrix):
    '''
    :param matrix: matriz a descomponer en l y u
    :return L, U
    '''

    n = matrix.shape[0]
    U = matrix
    L = np.eye(n)

    for k in range(n-1):
        for l in range(k+1, n):
            m = U[l, k] / U[k, k];
            U[l, k] = 0;
            U[l, k + 1:] = U[l, k+1:n] - m * U[k, k+1:]
            L[l, k] = m;

    return L, U

L, U = LUdec(np.array([[3, 5], [6, 7]]))
print(L, U)

max_stars_repo_name: mzahnd/algoritmos-mn

max_issues_repo_name: mzahnd/algoritmos-mn

max_forks_repo_name: mzahnd/algoritmos-mn


text:
import datetime
import numpy as np

from hyperopt import hp, fmin, tpe
import os
import sys

from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression

sys.path.insert(0, os.getcwd())

import workdir.classes.config
from qml.cv import QCV
from qml.models import QXgb, QXgb2
from workdir.classes.models import qm

if __name__ == "__main__":
    CV_SCORE_TO_STOP = 0.556
    DATAS = [219, 266, 269]

    EVALS_ROUNDS = 4000

    rounds = EVALS_ROUNDS

    cv = QCV(qm)
    counter = 0
    def fn(params):
        global counter

        counter +=1

        params['C'] = params['C'] / 10
        params['max_iter'] = int(1.3 ** params['max_iter'])

        model_id = qm.add_by_params(
            LogisticRegression(
                n_jobs=-1,
                max_iter=params['max_iter'],
                solver=params['solver'],
                C=params['C']
            ),
            'hyperopt log_regr',
            predict_fn='predict_proba'
        )
        res = cv.cross_val(model_id, params['data_id'], seed=1000, early_stop_cv=lambda x: x>CV_SCORE_TO_STOP)
        res = np.float64(res)
        res_arr = [res]
        # if res < CV_SCORE_TO_STOP:
        #     for i in range(7):
        #         res = cv.cross_val(model_id, data_id, seed=1001 + i, force=True)
        #         res = np.float64(res)
        #         res_arr.append(res)

        print(params['data_id'], model_id, "{}/{}".format(counter, rounds), res_arr, datetime.datetime.now(),  params)
        return np.mean(res_arr)
    space = {
        'max_iter': hp.quniform('max_iter', 18, 28, 1),
        'solver': hp.choice('solver', ['newton-cg', 'lbfgs', 'liblinear', 'sag']),
        'C': hp.quniform('C', 5, 40,1),
    }

    counter = 0
    space['data_id'] = hp.choice('data_id', DATAS)
    rounds = EVALS_ROUNDS
    fmin(fn, space, algo=tpe.suggest, max_evals=rounds)


text:
import numpy as np

from moviepy.editor import *

w,h = moviesize = (720,380)

duration = 1

def f(t,size, a = np.pi/3, thickness = 20):
    w,h = size
    v = thickness* np.array([np.cos(a),np.sin(a)])[::-1]
    center = [int(t*w/duration),h/2]
    return biGradientScreen(size,center,v,0.6,0.0)

logo = ImageClip("../../videos/logo_descr.png").\
         resize(width=w/2).\
         set_mask(mask)

screen = logo.on_color(moviesize, color = (0,0,0), pos='center')

shade = ColorClip(moviesize,color=(0,0,0))
mask_frame = lambda t : f(t,moviesize,duration)
shade.mask = VideoClip(ismask=True, get_frame = mask_frame)

cc = CompositeVideoClip([im.set_pos(2*["center"]),shade],
                         size = moviesize)

cc.subclip(0,duration).write_videofile("moviepy_logo.avi",fps=24)

max_stars_repo_name: sruon/moviepy

max_issues_repo_name: sruon/moviepy

max_forks_repo_name: sruon/moviepy


text:
"""Functions that generate artificial test cases for machine learning.

Copyright by Nikolaus Ruf
Released under the MIT license - see LICENSE file for details
"""

import numpy as np

def create_checkerboard(size=8, samples_per_square=100, random_state=None):  # pragma: no cover
    """Create samples for the 'checkerboard' problem.

    :param size: integer greater or equal to 2; number of rows and columns making up the checkerboard
    :param samples_per_square: positive integer; average number of samples per square
    :param random_state: instance of numpy.random.RandomState, integer, or None; a random state is used for sampling,
        while any other argument is used as seed to create a random state
    :return: 2D numpy float array and 1d numpy int array; the first value is the (random) feature matrix; the second
        value is the binary target, which is a deterministic function of the features; the first dimension of each value
        has length equal to samples_per_square * size ** 2
    """
    if size < 2:
        raise ValueError("Parameter size must be at least 2.")
    if samples_per_square <= 0:
        raise ValueError("Parameter samples_per_square must positive.")
    if not isinstance(random_state, np.random.RandomState):
        random_state = np.random.RandomState(random_state)
    num_samples = samples_per_square * size ** 2
    features = random_state.uniform(low=0.0, high=1.0, size=(num_samples, 2))
    target = np.int64(np.floor(size * features[:, 0]) + np.floor(size * features[:, 1])) % 2
    return features, target

def create_continuous_xor(size=6, samples_per_orthant=100, random_state=None):  # pragma: no cover
    """Create samples for the 'continuous XOR' problem.

    :param size: integer greater or equal to 1; number of features
    :param samples_per_orthant: positive integer; average number of features per orthant
    :param random_state: instance of numpy.random.RandomState, integer, or None; a random state is used for sampling,
        while any other argument is used as seed to create a random state
    :return: 2D numpy float array and 1d numpy int array; the first value is the (random) feature matrix; the second
        value is the binary target, which is a deterministic function of the features; the first dimension of each value
        has length equal to samples_per_orthant * 2 ** size
    """
    if size < 1:
        raise ValueError("Parameter size must be at least 1.")
    if samples_per_orthant <= 0:
        raise ValueError("Parameter samples_per_orthant must positive.")
    if not isinstance(random_state, np.random.RandomState):
        random_state = np.random.RandomState(random_state)
    num_samples = samples_per_orthant * 2 ** size
    features = random_state.uniform(low=-1.0, high=1.0, size=(num_samples, size))
    target = np.int64(np.prod(features, axis=1) >= 0)
    return features, target


text:
import os
from itertools import islice

import nltk
import numpy as np
import pandas as pd
import streamlit as st
import yake
from nltk import FreqDist, StanfordTagger, ngrams, word_tokenize
from sklearn.feature_extraction.text import TfidfVectorizer

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def get_tfidf_features(sentences, sort=True, filter_pos_tags=[]):
    new_text = sentences
    if filter_pos_tags != []:
        new_text = ""
        pos_tags = get_stanford_tags(sentences)
        for sentence in pos_tags:
            new_text += (
                " ".join(
                    [
                        word
                        for word, word_class in sentence
                        if word_class in filter_pos_tags
                    ]
                )
                + "\n"
            )

    lines = new_text.split("\n")
    vectorizer = TfidfVectorizer()
    X = vectorizer.fit_transform(lines).toarray()
    X = np.mean(X, axis=0)
    X = X.reshape(-1)
    result_dict = dict(zip(list(vectorizer.get_feature_names()), list(X)))
    if sort:
        result_dict_tuple = sorted(
            result_dict.items(), key=lambda x: x[1], reverse=True
        )
        result_dict = {k: v for k, v in result_dict_tuple}
    return result_dict

from .misc import length_removed_keywords, remove_keywords_util

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def get_stanford_tags(sentences):
    pos_tags = []
    for line in sentences.split("\n"):
        pos_tagged = nltk.pos_tag(word_tokenize(line))
        pos_tags.append(pos_tagged)
    return pos_tags

def find_ngrams_for_sentences(sentences, n=1):
    sentence_list = sentences.split("\n")
    ngrams_list = []
    for sentence in sentence_list:
        ngrams_list += list(ngrams(sentence.split(), n))
    return ngrams_list

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def find_freq(text, n=1, sort=False, filter_pos_tags=[]):
    new_text = text
    if filter_pos_tags != []:
        new_text = ""
        pos_tags = get_stanford_tags(text)
        for sentence in pos_tags:
            new_text += (
                " ".join(
                    [
                        word
                        for word, word_class in sentence
                        if word_class in filter_pos_tags
                    ]
                )
                + "\n"
            )

    ngrams_lst = find_ngrams_for_sentences(new_text, n)
    gram_count_mapping = dict(FreqDist(ngrams_lst))
    gram_count_mapping = {" ".join(k): v for k, v in gram_count_mapping.items()}
    if sort:
        sorted_gram_count_tuple = sorted(
            gram_count_mapping.items(), key=lambda x: x[1], reverse=True
        )
        gram_count_mapping = {k: v for k, v in sorted_gram_count_tuple}
    return gram_count_mapping

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def find_norm_freq(text, n=1, sort=False, filter_pos_tags=[]):
    gram_count_mapping = find_freq(
        text=text, n=n, sort=sort, filter_pos_tags=filter_pos_tags
    )
    norm_factor = sum(list(gram_count_mapping.values()))
    gram_count_mapping = {k: v / norm_factor for k, v in gram_count_mapping.items()}
    return gram_count_mapping

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def find_productivity(words, text, n=2, normalize=False, filter_pos_tags=[]):
    if normalize:
        fdist = find_norm_freq(
            text=text, n=n, sort=True, filter_pos_tags=filter_pos_tags
        )
    else:
        fdist = find_freq(text=text, n=n, sort=True, filter_pos_tags=filter_pos_tags)

    ngrams_lst = list(fdist.keys())
    prods = {}
    for word in words:
        ngrams_lst_having_word = []
        for ngram in ngrams_lst:
            if word in ngram:
                ngrams_lst_having_word.append(ngram)

        ngrams_lst_having_word = list(set(ngrams_lst_having_word))

        # find f_m_i for every ngram
        f_m_is = []
        for i in range(len(ngrams_lst_having_word)):
            f_m_i = fdist[ngrams_lst_having_word[i]]
            f_m_is.append(f_m_i)

        total_f_m_i = sum(f_m_is)
        p_m_is = [i / total_f_m_i for i in f_m_is]
        prod = -np.sum(np.multiply(np.log2(p_m_is), p_m_is))
        prods[word] = prod

    return prods

@st.cache(persist=eval(os.getenv("PERSISTENT")))
def freq_top_k(
    text,
    top_k=20,
    n=1,
    normalize=True,
    filter_pos_tags=[],
    tfidf=False,
    remove_keywords_path="./removed_keywords/removedphrases.txt",
):
    if tfidf:
        sorted_gram_count_mapping = get_tfidf_features(
            text, sort=True, filter_pos_tags=filter_pos_tags
        )
    else:
        if normalize:
            sorted_gram_count_mapping = find_norm_freq(
                text, n=n, sort=True, filter_pos_tags=filter_pos_tags
            )
        else:
            sorted_gram_count_mapping = find_freq(
                text, n=n, sort=True, filter_pos_tags=filter_pos_tags
            )

    if remove_keywords_path is not None and os.path.isfile(remove_keywords_path):
        sorted_gram_count_mapping = remove_keywords_util(
            remove_keywords_path, sorted_gram_count_mapping
        )

    if top_k < len(sorted_gram_count_mapping):
        sorted_gram_count_mapping = dict(
            islice(sorted_gram_count_mapping.items(), top_k)
        )

    return sorted_gram_count_mapping

def yake_keyword_extraction(
    text_file,
    top_k=20,
    language="en",
    max_ngram_size=2,
    window_size=2,
    deduplication_threshold=0.9,
    deduplication_algo="seqm",
    remove_keywords_path="./removed_keywords/removedphrases.txt",
):
    with open(text_file, "r") as f:
        text = f.read()
    if remove_keywords_path is not None and os.path.isfile(remove_keywords_path):
        len_rem_keywords = length_removed_keywords(remove_keywords_path)
    else:
        len_rem_keywords = 0
    custom_kw_extractor = yake.KeywordExtractor(
        lan=language,
        n=max_ngram_size,
        dedupLim=deduplication_threshold,
        dedupFunc=deduplication_algo,
        windowsSize=window_size,
        top=top_k + len_rem_keywords,
        features=None,
    )
    keywords = custom_kw_extractor.extract_keywords(text)
    if remove_keywords_path is not None and os.path.isfile(remove_keywords_path):
        keywords = remove_keywords_util(remove_keywords_path, dict(keywords))
    x = list(keywords.keys())
    y = [1 / (1e5 * keywords[key]) for key in keywords.keys()]

    df = pd.DataFrame()
    df["ngram"] = x
    df["score"] = y
    return df

max_stars_repo_name: gchhablani/DRIFT

max_issues_repo_name: gchhablani/DRIFT

max_forks_repo_name: gchhablani/DRIFT


text:
"""
Inspired by:

    1. RLKit: https://github.com/rail-berkeley/rlkit [X]
    2. StabelBaselines: https://github.com/hill-a/stable-baselines
    3. SpinningUp OpenAI: https://github.com/openai/spinningup
    4. CleanRL: https://github.com/vwxyzjn/cleanrl

"""

import os, subprocess, sys
import argparse
import importlib
import datetime
import random

import time
import wandb

import numpy as np
import torch as T
import torch.nn.functional as F
from torch.distributions.normal import Normal
nn = T.nn

from rl.algorithms.mfrl.mfrl import MFRL
# from rl.algorithms.mfrl.reinforce import Reinforce
from rl.control.policy import NPGPolicy, StochasticPolicy
from rl.value_functions.v_function import VFunction

class ActorCritic: # Done
    """
    Actor-Critic
        An entity contains both the actor (policy) that acts on the environment,
        and a critic (V-function) that evaluate that state given a policy.
    """
    def __init__(self,
                 obs_dim, act_dim,
                 act_up_lim, act_low_lim,
                 configs, seed, device
                 ) -> None:
        print('Initialize AC!')
        super(ActorCritic, self).__init__()
        self.obs_dim, self.act_dim = obs_dim, act_dim
        self.act_up_lim, self.act_low_lim = act_up_lim, act_low_lim
        self.configs, self.seed = configs, seed
        self._device_ = device

        self.actor, self.critic = None, None
        self._build()

    def _build(self):
        self.actor = self._set_actor()
        self.critic = self._set_critic()

    def _set_actor(self):
        net_configs = self.configs['actor']['network']
        # return PPOPolicy(
        #     self.obs_dim, self.act_dim,
        #     self.act_up_lim, self.act_low_lim,
        #     net_configs, self._device_, self.seed)
        return StochasticPolicy(
            self.obs_dim, self.act_dim,
            self.act_up_lim, self.act_low_lim,
            net_configs, self._device_, self.seed)

    def _set_critic(self):
        net_configs = self.configs['critic']['network']
        return VFunction(
            self.obs_dim, self.act_dim,
            net_configs, self._device_, self.seed)

    def get_v(self, o):
        return self.critic(o)

    def get_pi(self, o, a=None, reparameterize=False, deterministic=False, return_log_pi=True):
        action, log_pi, entropy = self.actor(o, a, reparameterize, deterministic, return_log_pi)
        return action, log_pi, entropy

    # def get_action(self, o, a=None):
    #     o = T.Tensor(o)
    #     if a: a = T.Tensor(a)
    #     with T.no_grad(): action, _, _ = self.actor(T.Tensor(o), a)
    #     return action.cpu().numpy()

    def get_action(self, o, a=None, reparameterize=False, deterministic=False, return_log_pi=False):
        o = T.Tensor(o)
        if a: a = T.Tensor(a)
        with T.no_grad(): a, _, _ = self.actor(o, a, reparameterize, deterministic, return_log_pi)
        return a.cpu()

    def get_action_np(self, o, a=None, reparameterize=False, deterministic=False, return_log_pi=False):
        return self.get_action(o, a, reparameterize, deterministic, return_log_pi).numpy()

    def get_pi_and_v(self, o, a=None, reparameterize=False, deterministic=False, return_log_pi=True):
        action, log_pi, entropy = self.actor(o, a, reparameterize, deterministic, return_log_pi)
        return action, log_pi, entropy, self.critic(o)

    def get_a_and_v_np(self, o, a=None, reparameterize=False, deterministic=False, return_log_pi=True):
        o = T.Tensor(o)
        if a: a = T.Tensor(a)
        with T.no_grad(): a, log_pi, entropy = self.actor(o, a, reparameterize, deterministic, return_log_pi)
        return a.cpu().numpy(), log_pi.cpu().numpy(), self.critic(o).cpu().numpy()

class NPG(MFRL):
    """
    Algorithm: Natural Policy Gradient (On-Policy, Model-Free)

        01. Initialize: Models parameters( Policy net πφ, Value net Vψ )
        02. Initialize: Trajectory buffer Dτ = {}
        03. Hyperparameters: Disc. factor γ, GAE λ, num traj's Nτ, rollout horizon H
        04. for n = 0, 1, 2, ..., N:
        05.     Collect set of traj's {τ^πk} by πk = π(φk) for e = 0, 1, 2, ..., E
        06.     Aggregate the traj's in traj buffer, Dτ = Dτ U {τ^πk}
        07.     Compute RTG R^_t, GAE A^_t based on Vk = V(θk)
        08.     for g = 0, 1, 2, ..., G do
        09.         Update πφ by maxz Jπ
                        φ = arg max_φ {(1/T|Dk|) sum sum min((π/πk), 1 +- eps) Aπk }
        10.         Fit Vθ by MSE(Jv)
                        ψ = arg min_ψ {(1/T|Dk|) sum sum (Vψ(st) - RTG)^2 }
        11.     end for
        12. end for
    """
    def __init__(self, exp_prefix, configs, seed, device, wb) -> None:
        super(NPG, self).__init__(exp_prefix, configs, seed, device)
        # print('init NPG Algorithm!')
        self.configs = configs
        self.seed = seed
        self._device_ = device
        self.WandB = wb
        self._build()

    def _build(self):
        super(NPG, self)._build()
        self._build_npg()

    def _build_npg(self):
        self._set_actor_critic()

    def _set_actor_critic(self):
        self.actor_critic = ActorCritic(
            self.obs_dim, self.act_dim,
            self.act_up_lim, self.act_low_lim,
            self.configs, self.seed, self._device_)

    def learn(self):
        N = self.configs['algorithm']['learning']['epochs']
        NT = self.configs['algorithm']['learning']['epoch_steps']
        Ni = self.configs['algorithm']['learning']['init_epochs']
        Nx = self.configs['algorithm']['learning']['expl_epochs']

        E = self.configs['algorithm']['learning']['env_steps']
        F = self.configs['algorithm']['learning']['train_AC_freq']
        G = self.configs['algorithm']['learning']['grad_AC_steps']

        batch_size = self.configs['data']['batch_size']
        mini_batch_size = self.configs['data']['mini_batch_size']

        global_step = 0
        start_time = time.time()
        o, d, Z, el, t = self.learn_env.reset(), 0, 0, 0, 0
        oldJs = [0, 0]
        JVList, JPiList = [0]*Ni, [0]*Ni
        logs = dict()
        lastEZ, lastES = 0, -2

        start_time_real = time.time()
        for n in range(1, N + 1):
            if self.configs['experiment']['print_logs']:
                print('=' * 80)
                if n > Nx:
                    print(f'\n[ Epoch {n}   Learning ]')
                elif n > Ni:
                    print(f'\n[ Epoch {n}   Exploration + Learning ]')
                else:
                    print(f'\n[ Epoch {n}   Inintial Exploration ]')

            nt = 0
            # self.buffer.reset()
            learn_start_real = time.time()
            while nt < NT:
                # Interaction steps (On-Policy)
                for e in range(1, E+1):
                    print('t: ', t, end='\r')
                    # global_step += 1 #* num_envs
                    o, d, Z, el, t = self.internact_op(n, o, d, Z, el, t)

                    if t % 1000 == 0:
                        # print(f"Training: global_step={global_step}, return={round(Z, 2)}, ep_length={el}")
                        EZ, ES, EL = self.evaluate_op()
                        print(f"Evaluation: global_step={t}, return={round(np.mean(EZ), 2)}, ep_length={np.mean(EL)}")
                        logs['evaluation/episodic_return_mean'] = np.mean(EZ)
                        logs['evaluation/episodic_length_mean'] = np.mean(EL)
                        if self.WandB:
                            wandb.log(logs)

                # with T.no_grad(): v = self.actor_critic.get_v(T.Tensor(o))
                # self.buffer.traj_tail(d, v)

                # Taking gradient steps after exploration
                if n > Ni and n % F == 0:
                    # print('updateAC')
                    with T.no_grad(): v = self.actor_critic.get_v(T.Tensor(o)).cpu()
                    self.buffer.traj_tail(d, v)
                    # Optimizing policy and value networks
                    # b_inds = np.arange(batch_size)
                    for g in range(1, G+1):
                        # NPG-P >>>>
                        for b in range(0, batch_size, mini_batch_size):
                            # print('ptr: ', self.buffer.ptr)
                            mini_batch = self.buffer.sample_batch(mini_batch_size, device=self._device_)
                            Jv, Jpi, stop_pi = self.trainAC(g, mini_batch, oldJs)
                            oldJs = [Jv, Jpi]
                            JVList.append(Jv.item())
                            JPiList.append(Jpi.item())
                        # NPG-P <<<<

                        # np.random.shuffle(b_inds)
                        # for start in range(0, batch_size, mini_batch_size):
                        #     end = start + mini_batch_size
                        #     mb_inds = b_inds[start:end]
                        #     mini_batch = self.buffer.sample_inds(mb_inds)
                        #     self.trainAC(mini_batch)
                    self.buffer.reset()
                nt += E

            # logs['time/training                  '] = time.time() - learn_start_real
            # logs['training/npg/Jv                '] = np.mean(JVList)
            # logs['training/npg/Jpi               '] = np.mean(JPiList)
            #
            # eval_start_real = time.time()
            # EZ, ES, EL = self.evaluate_op()
            #
            # logs['time/evaluation                '] = time.time() - eval_start_real
            #
            # if self.configs['environment']['type'] == 'mujoco-pddm-shadowhand':
            #     logs['evaluation/episodic_score_mean '] = np.mean(ES)
            #     logs['evaluation/episodic_score_std  '] = np.std(ES)
            # else:
            #     logs['evaluation/episodic_return_mean'] = np.mean(EZ)
            #     logs['evaluation/episodic_return_std '] = np.std(EZ)
            # logs['evaluation/episodic_length_mean'] = np.mean(EL)
            #
            # logs['time/total                     '] = time.time() - start_time_real
            #
            # if n > (N - 50):
            #     if self.configs['environment']['type'] == 'mujoco-pddm-shadowhand':
            #         if np.mean(ES) > lastES:
            #             print(f'[ Epoch {n}   Agent Saving ]                    ')
            #             env_name = self.configs['environment']['name']
            #             alg_name = self.configs['algorithm']['name']
            #             T.save(self.actor_critic.actor,
            #             f'./saved_agents/{env_name}-{alg_name}-seed:{self.seed}-epoch:{n}.pth.tar')
            #             lastES = np.mean(ES)
            #     else:
            #         if np.mean(EZ) > lastEZ:
            #             print(f'[ Epoch {n}   Agent Saving ]                    ')
            #             env_name = self.configs['environment']['name']
            #             alg_name = self.configs['algorithm']['name']
            #             T.save(self.actor_critic.actor,
            #             f'./saved_agents/{env_name}-{alg_name}-seed:{self.seed}-epoch:{n}.pth.tar')
            #             lastEZ = np.mean(EZ)
            #
            # # Printing logs
            # if self.configs['experiment']['print_logs']:
            #     for k, v in logs.items():
            #         print(f'{k}  {round(v, 2)}')
            #
            # # WandB
            # if self.WandB:
            #     wandb.log(logs)

        self.learn_env.close()
        self.eval_env.close()

    def trainAC(self, g, batch, oldJs):
        Jv = self.updateV(batch, oldJs[0])
        Jpi, stop_pi = self.updatePi(batch, oldJs[1])
        return Jv, Jpi, stop_pi

    def updateV(self, batch, Jv_old):
        """"
        Jv(θ) =
        """
        max_grad_norm = kl_targ = self.configs['critic']['network']['max_grad_norm']

        observations, _, returns, _, _, _ = batch.values()
        v = self.actor_critic.get_v(observations)

        Jv = 0.5 * ( (v - returns) ** 2 ).mean(axis=0)

        self.actor_critic.critic.optimizer.zero_grad()
        Jv.backward()
        nn.utils.clip_grad_norm_(self.actor_critic.critic.parameters(), max_grad_norm)
        self.actor_critic.critic.optimizer.step()

        return Jv

    def updatePi(self, batch, Jpi_old, stop_pi=False):
        """
        Jπ(φ) =
        """

        clip_eps = self.configs['actor']['clip_eps']
        kl_targ = self.configs['actor']['kl_targ']
        normz_step_size = self.configs['actor']['normz_step_size']

        entropy_coef = self.configs['actor']['entropy_coef']
        max_grad_norm = self.configs['actor']['network']['max_grad_norm']

        observations, actions, _, _, advantages, log_pis_old = batch.values()

        _, log_pi, entropy = self.actor_critic.get_pi(observations, actions)
        logratio = log_pi - log_pis_old
        ratio = logratio.exp()

        with T.no_grad():
            # calculate approx_kl http://joschu.net/blog/kl-approx.html
            approx_kl_old = (-logratio).mean()
            approx_kl = ((ratio - 1) - logratio).mean()

        # clipped_ratio = T.clamp(ratio, 1-clip_eps, 1+clip_eps)
        # Jpg = - ( T.min(ratio * advantages, clipped_ratio * advantages) ).mean(axis=0)
        # Jentropy = entropy_coef * entropy.mean()
        # Jpi = Jpg + Jentropy

        # if approx_kl_old <= kl_targ:
        clipped_ratio = T.clamp(ratio, 1-clip_eps, 1+clip_eps)
        Jpg = - ( T.min(ratio * advantages, clipped_ratio * advantages) ).mean(axis=0)
        Jentropy = entropy_coef * 0. #entropy.mean()
        Jpi = Jpg + Jentropy

        # Jvp =
        # Jnp =

        self.actor_critic.actor.optimizer.zero_grad()
        Jpi.backward()
        # nn.utils.clip_grad_norm_(self.actor_critic.actor.parameters(), max_grad_norm)
        self.actor_critic.actor.optimizer.step()

        return Jpi, stop_pi

def main(exp_prefix, config, seed, device, wb):

    print('Start an NPG experiment...')
    print('\n')

    configs = config.configurations

    if seed:
        random.seed(seed), np.random.seed(seed), T.manual_seed(seed)

    alg_name = configs['algorithm']['name']
    env_name = configs['environment']['name']
    env_type = configs['environment']['type']

    group_name = f"{env_name}-{alg_name}-X"
    exp_prefix = f"seed:{seed}"

    if wb:
        wandb.init(
            name=exp_prefix,
            group=group_name,
            # project='test',
            project='AMMI-RL-2022',
            config=configs
        )

    agent = NPG(exp_prefix, configs, seed, device, wb)

    agent.learn()

    print('\n')
    print('... End the NPG experiment')

if __name__ == "__main__":

    import argparse

    parser = argparse.ArgumentParser()

    parser.add_argument('-exp_prefix', type=str)
    parser.add_argument('-cfg', type=str)
    parser.add_argument('-seed', type=str)
    parser.add_argument('-device', type=str)
    parser.add_argument('-wb', type=str)

    args = parser.parse_args()

    exp_prefix = args.exp_prefix
    sys.path.append(f"{os.getcwd()}/configs")
    config = importlib.import_module(args.cfg)
    seed = int(args.seed)
    device = args.device
    wb = eval(args.wb)

    main(exp_prefix, config, seed, device, wb)

max_stars_repo_name: RamiSketcher/AMMI-RL

max_issues_repo_name: RamiSketcher/AMMI-RL

max_forks_repo_name: RamiSketcher/AMMI-RL


text:
"""
CRAYNN - yet another NN toolkit.
"""

from setuptools import setup, find_packages
from codecs import open
import os.path as osp
import numpy as np

here = osp.abspath(osp.dirname(__file__))

with open(osp.join(here, 'README.md'), encoding='utf-8') as f:
  long_description = f.read()

setup(
  name = 'craynn',

  version='0.9.1',

  description="""Yet another neural network toolkit.""",

  long_description = long_description,

  url='https://github.com/maxim-borisyak/craynn',

  author='Maxim Borisyak and contributors.',
  author_email='mborisyak at hse dot ru',

  maintainer = 'Maxim Borisyak',
  maintainer_email = 'mborisyak at hse dot ru',

  license='MIT',

  classifiers=[
    'Development Status :: 4 - Beta',

    'Intended Audience :: Science/Research',

    'Topic :: Scientific/Engineering :: Image Recognition',

    'License :: OSI Approved :: MIT License',

    'Programming Language :: Python :: 3',
  ],

  keywords='neural networks',

  packages=find_packages(exclude=['contrib', 'examples', 'docs', 'tests']),

  extras_require={
    'dev': ['check-manifest'],
    'test': ['nose>=1.3.0'],
  },

  install_requires=[
    'numpy',
    'scipy',
    'matplotlib',
    'theano',
    'lasagne',
    'pydotplus',
  ],

  include_package_data=True,
  ext_modules = [],
  include_dirs = [np.get_include()]
)

max_stars_repo_name: maxim-borisyak/craynn

max_issues_repo_name: maxim-borisyak/craynn

max_forks_repo_name: maxim-borisyak/craynn


text:
from sympy.vector import CoordSys3D
from symplyphysics import (
    symbols, Function, Eq, pretty, sin, pi
)
from symplyphysics.laws.nuclear.buckling import geometric_buckling_from_neutron_flux

# Description
## Neutron flux formula for the uniform spherical reactor. The spherical reactor is situated in spherical
## geometry at the origin of coordinates.

## Law: Ф(r) = C1 * (sin(Pi * r / R) / r)
## Where:
## C1 - neutron flux power constant.
## r - distance from the center of sphere.
## Pi - Pi constant.
## R - sphere radius.
## Ф(r) - neutron flux density.

neutron_flux_power_constant = symbols('C1', constant=True)
distance_from_center = symbols('distance_from_center')
sphere_radius = symbols('sphere_radius')
neutron_flux_function = symbols('neutron_flux_function', cls = Function)

# This constant is being used for geometric buckling calculation
# See: [geometric buckling for uniform sphere](geometric_buckling_for_uniform_sphere.py)
radial_constant = pi / sphere_radius

law = Eq(neutron_flux_function(distance_from_center),
    neutron_flux_power_constant * sin(radial_constant * distance_from_center) / distance_from_center)

# Check the solution by passing the known neutron flux to the geometric_buckling_from_neutron_flux.
# Neutron flux is a function of radius in the spherical coordinates.

# Boundary conditions:
# - vacuum boundary condition: Ф(R + d) = Ф(Re) = 0
# - finite flux condition: 0 <= Ф(r) < ∞
# - interface condition: the neutron flux and the normal component of the neutron current must be continuous
# - source condition: all neutrons flowing through the bounding area of the source must come from the neutron source
# - albedo boundary condition: Ф(Ralbedo) = 0

# define flux function in spherical coordinates as a function of sphere radius
spherical_coordinates = CoordSys3D('spherical_coordinates', transformation='spherical')
neutron_flux_function_spherical = law.subs(distance_from_center, spherical_coordinates.r)

solved = geometric_buckling_from_neutron_flux.apply_neutron_flux_function(neutron_flux_function_spherical.rhs)

# check with the derived law: Bg^2 = radial_constant**2
assert solved.rhs == radial_constant**2

def print():
    return pretty(law, use_unicode=False)

# There is no calculate() method. Neutron flux is usually being used internally to pass to other laws.

max_stars_repo_name: blackyblack/symplyphysics

max_issues_repo_name: blackyblack/symplyphysics

max_forks_repo_name: blackyblack/symplyphysics

reason: from sympy


text:
#!/usr/bin/env python

'produce a solution file for kaggle from category predictions files'

import os
import sys
import numpy as np
from glob import glob

input_dir = sys.argv[1]
output_file = sys.argv[2]

print "loading predictions..."

os.chdir( input_dir )
files = glob( '*.txt' )

all_p = np.empty(( 0, 2 ))
for f in files:
	print f
	p = np.loadtxt( f )
	all_p = np.vstack(( all_p, p ))

# for category 3, a missing example prediction from p_kaggle_no_holdout
all_p = np.vstack(( all_p, np.array(( -11.858253, 89655507 ))))

assert( all_p.shape[0] == 1351242 )

# sort by prediction and reverse
i = np.argsort( all_p[:,0] )[::-1]

# ids of most certain predictions
p_ids = all_p[i,1].astype ( int )

print "writing..."

np.savetxt( output_file, p_ids, fmt = '%d', header = 'Id', comments = '' )

max_stars_repo_name: zygmuntz/kaggle-avito

max_issues_repo_name: zygmuntz/kaggle-avito

max_forks_repo_name: zygmuntz/kaggle-avito


text:
#!/usr/bin/env python

import io, yaml, os
import matplotlib.pyplot as plt
import numpy as np
import warnings, dask, netCDF4
from datetime import datetime, date
import xarray as xr
from mom6_tools.DiagsCase import DiagsCase
from ncar_jobqueue import NCARCluster
from dask.distributed import Client
from mom6_tools import m6plot
from mom6_tools  import m6toolbox
from mom6_tools.MOM6grid import MOM6grid

def options():
  try: import argparse
  except: raise Exception('This version of python is not new enough. python 2.7 or newer is required.')
  parser = argparse.ArgumentParser(description='''Script for plotting poleward heat transport.''')
  parser.add_argument('diag_config_yml_path', type=str, help='''Full path to the yaml file  \
    describing the run and diagnostics to be performed.''')
  parser.add_argument('-v', '--variables', nargs='+', default=['T_ady_2d', 'T_diffy_2d', 'T_lbd_diffy_2d'],
                     help='''Variables to be processed (default=['T_ady_2d', 'T_diffy_2d', 'T_lbd_diffy_2d'])''')
  parser.add_argument('-sd','--start_date', type=str, default='',
                      help='''Start year to compute averages. Default is to use value set in diag_config_yml_path''')
  parser.add_argument('-ed','--end_date', type=str, default='',
                      help='''End year to compute averages. Default is to use value set in diag_config_yml_path''')
  parser.add_argument('-nw','--number_of_workers',  type=int, default=2,
                      help='''Number of workers to use (default=2).''')
  parser.add_argument('-debug',   help='''Add priting statements for debugging purposes''',
                      action="store_true")

  cmdLineArgs = parser.parse_args()
  return cmdLineArgs

def main(stream=False):
  # Get options
  args = options()
  nw = args.number_of_workers
  if not os.path.isdir('PNG/HT'):
    print('Creating a directory to place figures (PNG/HT)... \n')
    os.system('mkdir -p PNG/HT')
  if not os.path.isdir('ncfiles'):
    print('Creating a directory to place figures (ncfiles)... \n')
    os.system('mkdir ncfiles')

  # Read in the yaml file
  diag_config_yml = yaml.load(open(args.diag_config_yml_path,'r'), Loader=yaml.Loader)

  # Create the case instance
  dcase = DiagsCase(diag_config_yml['Case'])
  args.case_name = dcase.casename
  args.savefigs = True; args.outdir = 'PNG/HT'
  RUNDIR = dcase.get_value('RUNDIR')
  print('Run directory is:', RUNDIR)
  print('Casename is:', dcase.casename)
  print('Variables to be processed:', args.variables)
  print('Number of workers to be used:', nw)

  # set avg dates
  avg = diag_config_yml['Avg']
  if not args.start_date : args.start_date = avg['start_date']
  if not args.end_date : args.end_date = avg['end_date']

  # read grid info
  grd = MOM6grid(RUNDIR+'/'+dcase.casename+'.mom6.static.nc')
  depth = grd.depth_ocean
  # remote Nan's, otherwise genBasinMasks won't work
  depth[np.isnan(depth)] = 0.0
  basin_code = m6toolbox.genBasinMasks(grd.geolon, grd.geolat, depth)
  parallel, cluster, client = m6toolbox.request_workers(nw)
  print('Reading dataset...')
  startTime = datetime.now()
  variables = args.variables

  def preprocess(ds):
    ''' Compute montly averages and return the dataset with variables'''
    for var in variables:
      print('Processing {}'.format(var))
      if var not in ds.variables:
        print('WARNING: ds does not have variable {}. Creating dataarray with zeros'.format(var))
        jm, im = grd.geolat.shape
        tm = len(ds.time)
        da = xr.DataArray(np.zeros((tm, jm, im)), dims=['time','yq','xh'], \
             coords={'yq' : grd.yq, 'xh' : grd.xh, 'time' : ds.time}).rename(var)
        ds = xr.merge([ds, da])
    #return ds[variables].resample(time="1Y", closed='left', \
    #       keep_attrs=True).mean(dim='time', keep_attrs=True)
    return ds[variables]

  if parallel:
    ds = xr.open_mfdataset(RUNDIR+'/'+dcase.casename+'.mom6.hm_*.nc', \
         parallel=True, data_vars='minimal', chunks={'time': 12},\
         coords='minimal', compat='override', preprocess=preprocess)
  else:
    ds = xr.open_mfdataset(RUNDIR+'/'+dcase.casename+'.mom6.hm_*.nc', \
         data_vars='minimal', coords='minimal', compat='override', \
         preprocess=preprocess)
  print('Time elasped: ', datetime.now() - startTime)

  print('Selecting data between {} and {}...'.format(args.start_date, args.end_date))
  startTime = datetime.now()
  ds_sel = ds.sel(time=slice(args.start_date, args.end_date))
  print('Time elasped: ', datetime.now() - startTime)

  print('Computing yearly means...')
  startTime = datetime.now()
  ds_sel = ds_sel.resample(time="1Y", closed='left',keep_attrs=True).mean('time',keep_attrs=True)
  print('Time elasped: ', datetime.now() - startTime)

  print('Computing time mean...')
  startTime = datetime.now()
  ds_sel = ds_sel.mean('time').load()
  print('Time elasped: ', datetime.now() - startTime)

  if parallel:
    print('Releasing workers...')
    client.close(); cluster.close()

  varName = 'T_ady_2d'
  print('Saving netCDF files...')
  attrs = {'description': 'Time-mean poleward heat transport by components ', 'units': ds[varName].units,
       'start_date': args.start_date, 'end_date': args.end_date, 'casename': dcase.casename}
  m6toolbox.add_global_attrs(ds_sel,attrs)

  ds_sel.to_netcdf('ncfiles/'+dcase.casename+'_heat_transport.nc')
  # create a ndarray subclass
  class C(np.ndarray): pass

  if varName in ds.variables:
    tmp = np.ma.masked_invalid(ds_sel[varName].values)
    tmp = tmp[:].filled(0.)
    advective = tmp.view(C)
    advective.units = ds[varName].units
  else:
    raise Exception('Could not find "T_ady_2d" in file "%s"'%(args.infile+args.monthly))

  varName = 'T_diffy_2d'
  if varName in ds.variables:
    tmp = np.ma.masked_invalid(ds_sel[varName].values)
    tmp = tmp[:].filled(0.)
    diffusive = tmp.view(C)
    diffusive.units = ds[varName].units
  else:
    diffusive = None
    warnings.warn('Diffusive temperature term not found. This will result in an underestimation of the heat transport.')

  varName = 'T_lbd_diffy_2d'
  if varName in ds.variables:
    tmp = np.ma.masked_invalid(ds_sel[varName].values)
    tmp = tmp[:].filled(0.)
    lbd = tmp.view(C)
    #lbd.units = ds[varName].units
  else:
    lbd = None
    warnings.warn('Lateral boundary mixing term not found. This will result in an underestimation of the heat transport.')

  plt_heat_transport_model_vs_obs(advective, diffusive, lbd, basin_code, grd, args)
  return

def plt_heat_transport_model_vs_obs(advective, diffusive, lbd, basin_code, grd, args):
  """Plots model vs obs poleward heat transport for the global, Pacific and Atlantic basins"""
  # Load Observations
  fObs = netCDF4.Dataset('/glade/work/gmarques/cesm/datasets/Trenberth_and_Caron_Heat_Transport.nc')
  # POP JRA-55, 31 year (years 29-59)
  pop = xr.open_dataset('/glade/u/home/bryan/MOM6-modeloutputanalysis/MHT_mean.g210.GIAF_JRA.v13.gx1v7.01.nc')
  # Estimate based on the JRA-55 v1.3 forcing (Tsujino et al, 2019)
  # basin = 0 is Global; basin = 1 is Atlantic and basin = 2 is IndoPacific
  jra = xr.open_dataset('/glade/work/gmarques/cesm/datasets/Heat_transport/jra55fcst_v1_3_annual_1x1/nht_jra55do_v1_3.nc')
  #Trenberth and Caron
  yobs = fObs.variables['ylat'][:]
  NCEP = {}; NCEP['Global'] = fObs.variables['OTn']
  NCEP['Atlantic'] = fObs.variables['ATLn'][:]; NCEP['IndoPac'] = fObs.variables['INDPACn'][:]
  ECMWF = {}; ECMWF['Global'] = fObs.variables['OTe'][:]
  ECMWF['Atlantic'] = fObs.variables['ATLe'][:]; ECMWF['IndoPac'] = fObs.variables['INDPACe'][:]

  #G and W
  Global = {}
  Global['lat'] = np.array([-30., -19., 24., 47.])
  Global['trans'] = np.array([-0.6, -0.8, 1.8, 0.6])
  Global['err'] = np.array([0.3, 0.6, 0.3, 0.1])

  Atlantic = {}
  Atlantic['lat'] = np.array([-45., -30., -19., -11., -4.5, 7.5, 24., 47.])
  Atlantic['trans'] = np.array([0.66, 0.35, 0.77, 0.9, 1., 1.26, 1.27, 0.6])
  Atlantic['err'] = np.array([0.12, 0.15, 0.2, 0.4, 0.55, 0.31, 0.15, 0.09])

  IndoPac = {}
  IndoPac['lat'] = np.array([-30., -18., 24., 47.])
  IndoPac['trans'] = np.array([-0.9, -1.6, 0.52, 0.])
  IndoPac['err'] = np.array([0.3, 0.6, 0.2, 0.05,])

  GandW = {}
  GandW['Global'] = Global
  GandW['Atlantic'] = Atlantic
  GandW['IndoPac'] = IndoPac

  if args.case_name != '':  suptitle = args.case_name
  else: suptitle = ''

  # Global Heat Transport
  plt.figure(figsize=(12,10))
  HTplot = heatTrans(advective,diffusive,lbd)
  yy = grd.geolat_c[:,:].max(axis=-1)
  ave_title = ', averaged from {} to {}'.format(args.start_date, args.end_date)
  plotHeatTrans(yy,HTplot,title='Global Y-Direction Heat Transport [PW]'+ave_title)
  plt.plot(pop.lat_aux_grid.values,pop.MHT_global.values,'orange',linewidth=1,label='POP')
  jra_mean_global = jra.nht[:,0,:].mean('time').values
  jra_std_global = jra.nht[:,0,:].std('time').values
  plt.plot(jra.lat, jra_mean_global,'k', label='JRA-55 v1.3', color='#1B2ACC', lw=1)
  plt.fill_between(jra.lat, jra_mean_global-jra_std_global, jra_mean_global+jra_std_global,
    alpha=0.25, edgecolor='#1B2ACC', facecolor='#089FFF')
  plt.plot(yobs,NCEP['Global'],'k--',linewidth=0.5,label='NCEP')
  plt.plot(yobs,ECMWF['Global'],'k.',linewidth=0.5,label='ECMWF')
  plotGandW(GandW['Global']['lat'],GandW['Global']['trans'],GandW['Global']['err'])
  plt.xlabel(r'Latitude [$\degree$N]',fontsize=10)
  plt.suptitle(suptitle)
  plt.legend(loc=0,fontsize=10)
  annotateObs()
  if diffusive is None: annotatePlot('Warning: Diffusive component of transport is missing.')
  if lbd is None: annotatePlot('Warning: LBD component of transport is missing.')

  if args.savefigs:
    objOut = args.outdir+'/'+args.case_name+'_HeatTransport_global.png'
    plt.savefig(objOut); plt.close()
  else:
    plt.show()
  # Atlantic Heat Transport
  m = 0*basin_code; m[(basin_code==2) | (basin_code==4) | (basin_code==6) | (basin_code==7) | (basin_code==8)] = 1
  plt.figure(figsize=(12,10))
  HTplot = heatTrans(advective, diffusive, lbd, vmask=m*np.roll(m,-1,axis=-2))
  yy = grd.geolat_c[:,:].max(axis=-1)
  HTplot[yy<-34] = np.nan
  plotHeatTrans(yy,HTplot,title='Atlantic Y-Direction Heat Transport [PW]'+ave_title)
  plt.plot(pop.lat_aux_grid.values,pop.MHT_atl.values,'orange',linewidth=1,label='POP')
  jra_mean_atl = jra.nht[:,1,:].mean('time').values
  jra_std_atl = jra.nht[:,1,:].std('time').values
  plt.plot(jra.lat, jra_mean_atl,'k', label='JRA-55 v1.3', color='#1B2ACC', lw=1)
  plt.fill_between(jra.lat, jra_mean_atl-jra_std_atl, jra_mean_atl+jra_std_atl,
    alpha=0.25, edgecolor='#1B2ACC', facecolor='#089FFF')
  plt.plot(yobs,NCEP['Atlantic'],'k--',linewidth=0.5,label='NCEP')
  plt.plot(yobs,ECMWF['Atlantic'],'k.',linewidth=0.5,label='ECMWF')
  plotGandW(GandW['Atlantic']['lat'],GandW['Atlantic']['trans'],GandW['Atlantic']['err'])
  plt.xlabel(r'Latitude [$\degree$N]',fontsize=10)
  plt.suptitle(suptitle)
  plt.legend(loc=0,fontsize=10)
  annotateObs()
  if diffusive is None: annotatePlot('Warning: Diffusive component of transport is missing.')
  if lbd is None: annotatePlot('Warning: LBD component of transport is missing.')
  if args.savefigs:
    objOut = args.outdir+'/'+args.case_name+'_HeatTransport_Atlantic.png'
    plt.savefig(objOut); plt.close()
  else:
    plt.show()
  # Indo-Pacific Heat Transport
  m = 0*basin_code; m[(basin_code==3) | (basin_code==5)] = 1
  plt.figure(figsize=(12,10))
  HTplot = heatTrans(advective, diffusive, lbd, vmask=m*np.roll(m,-1,axis=-2))
  yy = grd.geolat_c[:,:].max(axis=-1)
  HTplot[yy<-34] = np.nan
  plotHeatTrans(yy,HTplot,title='Indo-Pacific Y-Direction Heat Transport [PW]'+ave_title)
  plt.plot(pop.lat_aux_grid.values,(pop.MHT_global-pop.MHT_atl).values,'orange',linewidth=1,label='POP')
  jra_mean_indo = jra.nht[:,2,:].mean('time').values
  jra_std_indo = jra.nht[:,2,:].std('time').values
  plt.plot(jra.lat, jra_mean_indo,'k', label='JRA-55 v1.3', color='#1B2ACC', lw=1)
  plt.fill_between(jra.lat, jra_mean_indo-jra_std_indo, jra_mean_indo+jra_std_indo,
    alpha=0.25, edgecolor='#1B2ACC', facecolor='#089FFF')
  plt.plot(yobs,NCEP['IndoPac'],'k--',linewidth=0.5,label='NCEP')
  plt.plot(yobs,ECMWF['IndoPac'],'k.',linewidth=0.5,label='ECMWF')
  plotGandW(GandW['IndoPac']['lat'],GandW['IndoPac']['trans'],GandW['IndoPac']['err'])
  plt.xlabel(r'Latitude [$\degree$N]',fontsize=10)
  annotateObs()
  if diffusive is None: annotatePlot('Warning: Diffusive component of transport is missing.')
  if lbd is None: annotatePlot('Warning: LBD component of transport is missing.')
  plt.suptitle(suptitle)
  plt.legend(loc=0,fontsize=10)
  if args.savefigs:
    objOut = args.outdir+'/'+args.case_name+'_HeatTransport_IndoPacific.png'
    plt.savefig(objOut); plt.close()
  else:
    plt.show()
  return

def heatTrans(advective, diffusive=None, lbd=None, vmask=None):
  """Converts vertically integrated temperature advection into heat transport"""
  HT = advective[:]
  if diffusive is not None:
    HT = HT + diffusive[:]
  if lbd is not None:
    HT = HT + lbd[:]
  if len(HT.shape) == 3:
    HT = HT.mean(axis=0)
  if advective.units == "Celsius meter3 second-1":
    rho0 = 1.035e3
    Cp = 3992.
    HT = HT * (rho0 * Cp)
    HT = HT * 1.e-15  # convert to PW
  elif advective.units == "W":
    HT = HT * 1.e-15
  else:
    print('Unknown units')
  if vmask is not None: HT = HT*vmask
  HT = HT.sum(axis=-1); HT = HT.squeeze() # sum in x-direction
  return HT

def plotHeatTrans(y, HT, title, xlim=(-80,90), ylim=(-2.5,3.0)):
  plt.plot(y, y*0., 'k', linewidth=0.5)
  plt.plot(y, HT, 'r', linewidth=1.5,label='Model')
  plt.xlim(xlim); plt.ylim(ylim)
  plt.title(title)
  plt.grid(True)

def annotatePlot(label):
  fig = plt.gcf()
  #fig.text(0.1,0.85,label)
  fig.text(0.535,0.12,label)

def annotateObs():
  fig = plt.gcf()
  fig.text(0.13,0.85,r"Trenberth, K. E. and J. M. Caron, 2001: Estimates of Meridional Atmosphere and Ocean Heat Transports. J.Climate, 14, 3433-3443.", fontsize=8)
  fig.text(0.13,0.825,r"Ganachaud, A. and C. Wunsch, 2000: Improved estimates of global ocean circulation, heat transport and mixing from hydrographic data.", fontsize=8)
  fig.text(0.13,0.8,r"Nature, 408, 453-457", fontsize=8)

def plotGandW(lat,trans,err):
  low = trans - err
  high = trans + err
  for n in range(0,len(low)):
    if n == 0:
      plt.plot([lat[n],lat[n]], [low[n],high[n]], 'c', linewidth=2.0, label='G&W')
    else:
      plt.plot([lat[n],lat[n]], [low[n],high[n]], 'c', linewidth=2.0)
  plt.scatter(lat,trans,marker='s',facecolor='cyan')

if __name__ == '__main__':
  main()


text:
import typing
import sys
import numpy as np
import numba as nb

@nb.njit
def pow(mod: int, x: int, n: int) -> int:
    y = 1
    while n:
        if n & 1: y = y * x % mod
        x = x * x % mod
        n >>= 1
    return y

@nb.njit((nb.i8, ), cache=True)
def solve(l: int) -> typing.NoReturn:
    mod = 10 ** 9 + 7
    m = 10 << 17
    f = np.zeros(m, np.int64)
    g = np.zeros(m, np.int64)
    h = np.zeros(m, np.int64)
    f[0] = 1
    g[0] = 2
    h[0] = 1
    for i in range(m - 1):
        f[i + 1] = (3 * f[i] + pow(mod, g[i], 3)) % mod
        x = (g[i] * g[i] - h[i] * h[i]) % mod
        g[i + 1] = (x + g[i]) * g[i] % mod
        h[i + 1] = x * h[i] % mod
    print(f[l])

def main() -> typing.NoReturn:
    l = int(input())
    solve(l)

main()

max_stars_repo_name: kagemeka/competitive-programming

max_issues_repo_name: kagemeka/competitive-programming

max_forks_repo_name: kagemeka/competitive-programming

reason: import numpy,import numba


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

"""
Earthquake Predictor Neural Network Experiment
Prediction Map
"""
import numpy as np
from PIL import Image

def load_world_map(scale):
    print('Loading world map...')
    if scale==6:
        earth_map = np.array(Image.open('map/1080.jpg'))
    if scale==4:
        earth_map = np.array(Image.open('map/720.jpg'))
    if scale==2:
        earth_map = np.array(Image.open('map/360.jpg'))
    return earth_map

def merge_maps(background_image, earthquake_matrix, scale, min_mag):
    print('Merging earthquakes with world map...')
    earth_map = background_image
    quake_map = Image.new('RGBA',(181*scale,91*scale),(0,0,0,0))
    border_map = Image.new('RGBA',(181*scale,91*scale),(0,0,0,0))
    quake_map = np.array(quake_map)
    border_map = np.array(border_map)
    for y in range(91):
        for x in range(181):
            if earthquake_matrix[90-y,x] >= (min_mag + 1) / 11 and earthquake_matrix[90-y,x] > 0:
                for i in range(scale):
                    for j in range(scale):
                        quake_map[y*scale+i,x*scale+j] = (0,0,0,earthquake_matrix[90-y,x])
                        if scale > 4:
                            if i == 0 or i == scale-1 or j == 0 or j == scale-1:
                                border_map[y*scale+i,x*scale+j] = (0,0,0,0)
                            elif i == 1 or i == scale-2 or j == 1 or j == scale-2:
                                border_map[y*scale+i,x*scale+j] = (0,0,0,255)
                        elif scale > 2:
                            if i == 0 or i == scale-1 or j == 0 or j == scale-1:
                                border_map[y*scale+i,x*scale+j] = (0,0,0,255)
    earth_map = Image.fromarray(earth_map)
    yellow_map = Image.new('RGB',(181*scale,91*scale),(255,255,0))
    red_map = Image.new('RGB',(181*scale,91*scale),(255,0,0))
    border_map = Image.fromarray(border_map)
    quake_map = Image.fromarray(quake_map)
    earth_map = Image.composite(red_map, earth_map, quake_map)
    earth_map = Image.composite(yellow_map, earth_map, border_map)
    earth_map = np.array(earth_map)
    return earth_map

def create_legend(image, scale):
    x1 = 5
    x2 = 45
    y1 = 81
    y2 = 84
    legend = Image.new('RGBA',(181*scale,91*scale),(0,0,0,0))
    border = Image.new('RGBA',(181*scale,91*scale),(0,0,0,0))
    legend = np.array(legend)
    border = np.array(border)
    for y in range(91):
        for x in range(181):
            if x >= x1 and x <= x2 and y >= y1 and y <= y2:
                for i in range(scale):
                    for j in range(scale):
                        legend[y*scale+i,x*scale+j] = (0,0,0, (float(x-x1)/float(x2-x1))*255)
                        if x == x1 and j == 0 or x == x2 and j == scale-1 or y == y1 and i == 0 or y == y2 and i == scale-1:
                            border[y*scale+i,x*scale+j] = (0,0,0,255)
                        elif j == (float(x-x1)/float(x2-x1))*(scale-1) and (i == 1 or i == scale-2):
                            border[y*scale+i,x*scale+j] = (0,0,0,255)
                        elif (float(x-x1) % (float(x2-x1) / 10)) == 0 and (y == y1 and i == 1 or y == y2 and i == scale-2):
                            border[y*scale+i,x*scale+j] = (0,0,0,255)
    red = Image.new('RGB',(181*scale,91*scale),(255,0,0))
    yellow = Image.new('RGB',(181*scale,91*scale),(255,255,0))
    legend = Image.fromarray(legend)
    border = Image.fromarray(border)
    image = Image.fromarray(image)
    image = Image.composite(red,image,legend)
    image = Image.composite(yellow,image,border)
    image = np.array(image)
    return image

def create_map(earthquake_img, scale=6, min_mag=-1.0,filename='prediction/prediction_map.png'):
    print('Creating map of earthquakes...')
    bg = load_world_map(scale)
    eq = np.array(Image.open(earthquake_img))
    earthquake_map = merge_maps(bg, eq, scale, min_mag)
    earthquake_map = create_legend(earthquake_map, scale)
    Image.fromarray(earthquake_map).save(filename)

max_stars_repo_name: cognoscitive/earthquakepredictor

max_issues_repo_name: cognoscitive/earthquakepredictor

max_forks_repo_name: cognoscitive/earthquakepredictor