text:
import numpy as np
import pandas as pd
import argparse
from copy import deepcopy
import torch
from torch.optim import Adam
import torch.nn.functional as F
import torch.nn as nn
from torch.autograd import Variable

import torch.optim.lr_scheduler as Scheduler
from torch.utils.tensorboard import SummaryWriter

from evaluator import Evaluator
from memory import EpisodicMemory
from PER_memory import PrioritizedReplayBuffer
from agent import Agent
from util import *
# from util import prGreen, prYellow, prRed, to_tensor, soft_update,
import matplotlib.pylab as plt

DEMO_flag = 1

class RDPG(object):
    def __init__(self, demo_env, env, args):

        ### Create Environment Classes ###
        self.env = env
        self.demo_env = demo_env

        ##### Create Replay Buffer #####
        self.is_PER_replay = args.is_PER_replay
        if self.is_PER_replay:

            ### Porioritized experience replay (PER) buffer ###
            self.memory = PrioritizedReplayBuffer(args.PER_size, args.seed, alpha=args.p_alpha,
                                                  beta_init=1.0, beta_inc_n=100, max_t=args.exp_traj_len)

        else:
            ### Original replay buffer ###
            self.memory = EpisodicMemory(capacity=args.rmsize,
                                         max_train_traj_len=args.exp_traj_len,  # 整段episode被分段存取
                                         window_length=args.window_length)

        ### Evaluator is used for test trained agent ###
        self.evaluate = Evaluator(args.test_episodes,  #我們test那年有243個交易日
                                  max_episode_length=args.max_episode_length, #每天最長是240分鐘，最多交易240次
                                  args=args)

        ##### Model Setting #####
        self.rnn_mode = args.rnn_mode
        self.seq_len = args.seq_len
        self.hidden_rnn = args.hidden_rnn
        self.num_layer = args.num_rnn_layer
        self.agent = Agent(args)
        if torch.cuda.is_available() :
            self.agent.cuda()

        ### Hyper-parameters
        self.batch_size = args.bsize
        self.exp_traj_len = args.exp_traj_len
        self.max_episode_length = args.max_episode_length
        self.tau = args.tau
        self.discount = args.discount
        self.warmup = args.warmup
        self.a_update_freq = args.a_update_freq

        ##### Behavior Cloning Setting #####
        self.is_BClone = args.is_BClone
        self.is_Qfilt = args.is_Qfilt
        self.use_Qfilt = args.use_Qfilt
        if self.is_BClone:
            self.lambda_Policy = args.lambda_Policy
            self.lambda_BC = 1-self.lambda_Policy
        else:
            self.lambda_Policy = 1
            self.lambda_BC = 1-self.lambda_Policy
        # self.lambda_BC = args.lambda_BC
        self.BC_loss_func = nn.MSELoss(reduce=False)
        # self.BC_loss_func = nn.BCELoss(reduce=False)

        ##### PER demonstration Setting #####
        self.lambda_balance = args.lambda_balance
        self.small_const = args.small_const
        self.priority_const = args.priority_const

        self.is_demo_warmup = args.is_demo_warmup
        if self.is_demo_warmup:
            demo_protect_size = (self.max_episode_length /self.exp_traj_len) *self.warmup
            self.memory.set_protect_size(int(demo_protect_size))

        self.is_pretrain = args.is_pretrain
        self.pretrain_itrs = args.Pretrain_itrs

        ### Optimizer and LR_scheduler ###
        beta1 = args.beta1
        beta2 = args.beta2
        self.rnn_optim = Adam(self.agent.rnn.parameters(), lr=args.r_rate, betas=(beta1, beta2))
        self.rnn_scheduler = Scheduler.StepLR(self.rnn_optim, step_size=args.sch_step_size, gamma=args.sch_gamma)
        self.critic_optim  = Adam(self.agent.critic.parameters(), lr=args.c_rate, betas=(beta1, beta2))
        self.critic_scheduler = Scheduler.StepLR(self.critic_optim, step_size=args.sch_step_size, gamma=args.sch_gamma)
        self.actor_optim  = Adam(self.agent.actor.parameters(), lr=args.a_rate, betas=(beta1, beta2))
        self.actor_scheduler = Scheduler.StepLR(self.actor_optim, step_size=args.sch_step_size, gamma=args.sch_gamma)

        ### initialized values
        self.demoN_ratio = 0
        self.priority = 0
        self.actor_loss = 0
        self.BC_loss = 0
        self.BC_loss_Qf = 0
        self.tot_policy_loss = 0
        self.critic_loss = 0

        ### other setting ###
        if args.seed > 0:
            self.seed(args.seed)
        self.writer = SummaryWriter(args.logdir)
        self.is_training = True
        self.save_threshold = args.save_threshold
        self.date = args.date

    def train(self, num_episodes, checkpoint_path, debug):

        epi_idx = None #training時因為episode是random選date，並非照time_order，所以在此設成none沒關係。

        self.agent.is_training = True
        step = episode_steps = trajectory_steps = 0
        episode_reward = 0.
        state0 = None
        ewma_reward = 0
        episode = 1
        train_epi_reward = []
        train_ewma_reward = []
        train_actor_loss = []
        train_bc_loss = []
        train_bcQf_loss = []
        train_totPolicy_loss = []
        train_critic_loss = []
        demoN_ratio_batch = []

        while episode <= num_episodes:
            episode_steps = 1
            while episode_steps <= self.max_episode_length:

                if self.is_demo_warmup:
                #################### warmup adopt expert policy (Dual Thrust Strategy) ####################
                    if episode <= self.warmup:  # Note: warmup generate demonstrations, so here use demo_env
                        if state0 is None:
                            self.agent.reset()
                            state0 = deepcopy(self.demo_env.reset(epi_idx))  #training時的env reset是random選date，所以此時epi_idx=None沒關係。
                            state0 = state0.values

                        ### Note the following action will be from demonstration, not random.
                        action, epsilon = self.agent.random_action() #其實改成demo後，此action到demo_env也會變成demo的action，而非random
                        action_bc, next_state, reward, done, infos = self.demo_env.step(np.argmax(action))

                    else:  # normal training, without demonstration, so her use training env
                        if state0 is None:
                            self.agent.reset()
                            state0 = deepcopy(self.env.reset(epi_idx))  #training時的env reset是random選date，所以此時epi_idx=None沒關係。
                            state0 = state0.values

                        state0_cuda = to_tensor(np.array([state0])).cuda()
                        # state0_cuda = to_tensor(state0).cuda()
                        action, epsilon = self.agent.select_action(state0_cuda)
                        action_bc, next_state, reward, done, infos = self.env.step(np.argmax(action))
                ######################################################################

                else:
                ############## original random warmup experiences ##############
                    #### reset if it is the start of episode
                    if state0 is None:
                        self.agent.reset()
                        state0 = deepcopy(self.env.reset(epi_idx))  #training時的env reset是random選date，所以此時epi_idx=None沒關係。
                        state0 = state0.values
                    state0_cuda = to_tensor(np.array([state0])).cuda()

                    if episode <= self.warmup:
                        action, epsilon = self.agent.random_action()
                    else: # 正式training，新append experience為agent data
                        action, epsilon = self.agent.select_action(state0_cuda)

                    act = action
                    action_bc, next_state, reward, done, infos = self.env.step(np.argmax(act))
                ###############################################################

                next_state = next_state.values
                next_state = deepcopy(next_state)

                ###### agent observe and update policy #####

                #if episode <= self.warmup or episode-self.warmup > self.pretrain:
                if self.is_PER_replay == False:
                    self.memory.append(action_bc, state0, action, reward, done)
                else:
                    self.memory.add((torch.from_numpy(state0).float(),
                                           torch.from_numpy(action).float(),
                                           torch.from_numpy(action_bc).float(),
                                           torch.tensor([reward]).float(),
                                           torch.from_numpy(next_state).float(),
                                           torch.tensor([0.95]),
                                           int(episode<=self.warmup)),
                                           trajectory_steps)

                # update
                step += 1
                episode_steps += 1
                trajectory_steps += 1
                episode_reward += reward
                state0 = deepcopy(next_state)

                ##### 此exp_traj_len，目前設定成設定每16steps update一次 ##### (但每次episode還是跑到done才會結束)
                if trajectory_steps >= self.exp_traj_len:
                    ### 以下設定是為了讓hidden_state繼續往下一個step傳遞 ###
                    self.agent.reset_rnn_hidden(done=False)  #注意done對應model.py裡的reset_lstm_hidden_state()
                    trajectory_steps = 0
                    if episode > self.warmup:
                        self.update_policy(done)

                if done: # end of episod
                    # reset
                    state0 = None
                    ewma_reward = 0.05 * episode_reward + 0.95 * ewma_reward
                    (lr_Critic, lr_Actor, lr_RNN) = (self.critic_scheduler.get_last_lr()[0], self.actor_scheduler.get_last_lr()[0], self.rnn_scheduler.get_last_lr()[0])
                    if debug: prYellow('[Step:{}, Episode:{}, Len:{}] [lr:{}, epsl:{}] [Reward={:.3f}  Ewma_Reward={:.3f}]'.format(step, episode, episode_steps, lr_Critic, np.round(epsilon,3), episode_reward, ewma_reward))
                    if self.is_BClone:
                        print('A_loss={:.3f}\tBC_loss={:.3f}\tTotP_loss={:.3f}\tC_loss={:.3f}'.format(self.actor_loss, self.BC_loss, self.tot_policy_loss, self.critic_loss))
                    else:
                        print('TotP_loss={:.3f}\tC_loss={:.3f}'.format(self.tot_policy_loss, self.critic_loss))

                    ##### log data for Tensorboard #####
                    self.writer.add_scalar('Train/Episode Reward', episode_reward, episode)
                    self.writer.add_scalar('Train/EWMA Reward', ewma_reward, episode)
                    self.writer.add_scalar('Train/Actor Loss', self.actor_loss, episode)
                    self.writer.add_scalar('Train/BC Loss', self.BC_loss, episode)
                    self.writer.add_scalar('Train/Qf Loss', self.BC_loss_Qf, episode)
                    self.writer.add_scalar('Train/tot Policy Loss', self.tot_policy_loss, episode)
                    self.writer.add_scalar('Train/Critic Loss', self.critic_loss, episode)
                    self.writer.add_scalar('Train/Demon ratio in batch', self.demoN_ratio, episode)

                    ##### append training info for plot #####
                    train_epi_reward.append(episode_reward)
                    train_ewma_reward.append(ewma_reward)
                    train_actor_loss.append(self.actor_loss)
                    train_bc_loss.append(self.BC_loss)
                    train_bcQf_loss.append(self.BC_loss_Qf)
                    train_totPolicy_loss.append(self.tot_policy_loss)
                    train_critic_loss.append(self.critic_loss)
                    demoN_ratio_batch.append(self.demoN_ratio)

                    if episode == self.warmup:
                        ite=1
                        if self.is_pretrain:
                            while(ite<=self.pretrain_itrs):
                                self.update_policy(done)
                                if ite%(self.pretrain_itrs/10) == 0:
                                    print('demoN_ratio=',self.demoN_ratio)
                                    print('Pretrain',ite,':','TotP_loss={:.3f}\tC_loss={:.3f}\tDemo_ratio={:.3f}'.format(self.tot_policy_loss, self.critic_loss, self.demoN_ratio))
                                ite+=1

                    episode_reward = 0.
                    episode += 1
                    self.agent.reset_rnn_hidden(done=False)  #注意done對應model.py裡的reset_hidden_state()
                    break

            ##### Save models #####
            if (episode-1) >= 150 or (episode-1) % 100 == 0 or ewma_reward > self.save_threshold:
                self.agent.save_model(checkpoint_path, (episode-1), ewma_reward)

            ##### Plot Training Curves #####
            if (episode-1) % 100 == 0:
                if self.is_BClone:
                    self.train_plot_bc(episode-1, train_epi_reward, train_ewma_reward,
                                       train_totPolicy_loss, train_critic_loss,
                                       train_actor_loss, train_bc_loss, train_bcQf_loss)
                    self.train_demoN_ratio(episode-1, demoN_ratio_batch)
                else:
                    self.train_plot(episode-1, train_epi_reward, train_ewma_reward,
                                    train_totPolicy_loss, train_critic_loss)
                    self.train_demoN_ratio(episode-1, demoN_ratio_batch)

            ##### Apply Q-filter to BC loss #####
            if (episode-1) >= (self.warmup+self.use_Qfilt):
                self.is_Qfilt=True

#            if step >= args.warmup and episode > args.bsize:
#                # Update weights
#                agent.update_policy()

    def update_policy(self, done):
        ### Sample batch of trajectories
        t_len = 0
        if self.is_PER_replay == False:
            experiences = self.memory.sample(self.batch_size)
            if len(experiences) == 0: # not enough samples
                return
            t_len = len(experiences)
        else:
            (state0s, actions, action_bcs, rewards, state1s, batch_gammas, batch_flagss), \
                                           weights, idxes = self.memory.sample(self.batch_size)
            t_len = len(state0s)

        actor_loss_total = 0  #actor loss
        BC_loss_total = 0  #BC loss
        BC_loss_Qf_total = 0  #BC loss after Q-filter
        policy_loss_total = 0  #policy loss
        value_loss_total = 0  #critic loss
        demo_cnt = []
        for t in range(t_len): # iterate over episodes
            if self.is_PER_replay == False and t == t_len-1:
                break
            a_target_cx = Variable(torch.zeros(self.num_layer, self.batch_size, self.hidden_rnn)).type(FLOAT).cuda()
            a_target_hx = Variable(torch.zeros(self.num_layer, self.batch_size, self.hidden_rnn)).type(FLOAT).cuda()

            a_cx = Variable(torch.zeros(self.num_layer, self.batch_size, self.hidden_rnn)).type(FLOAT).cuda()
            a_hx = Variable(torch.zeros(self.num_layer, self.batch_size, self.hidden_rnn)).type(FLOAT).cuda()

            if self.is_PER_replay == False:
                action_bc = np.stack((trajectory.action_bc for trajectory in experiences[t]))

                state0 = np.stack((trajectory.state0 for trajectory in experiences[t]))
                action = np.stack((trajectory.action for trajectory in experiences[t]))
                action = to_tensor(action)
                reward = np.expand_dims(np.stack((trajectory.reward for trajectory in experiences[t])), axis=1)
                reward = to_tensor(reward)
                state1 = np.stack((trajectory.state0 for trajectory in experiences[t+1]))

                state0_cuda = to_tensor(state0).cuda()
                state1_cuda = to_tensor(state1).cuda()
            else:
                state0 = state0s[t]
                action = actions[t]
                reward = rewards[t]
                state1 = state1s[t]
                batch_flags = batch_flagss[t]
                action_bc = action_bcs[t]
                state0_cuda = state0.cuda()
                state1_cuda = state1.cuda()

                ##### calculate demonstration ratio in a batch #####
                d_flags = torch.from_numpy(batch_flags)
                demo_select = d_flags == DEMO_flag
                N_act = demo_select.sum().item()
                demo_cnt.append(N_act/self.batch_size)

            ######################## critic loss calculation ########################
            # with torch.no_grad():
            if self.rnn_mode == 'lstm':
                xh0, _ = self.agent.rnn(state0_cuda, (a_hx, a_cx))
                current_q = self.agent.critic([xh0, action.cuda()])

                with torch.no_grad():
                    xh1, _ = self.agent.rnn_target(state0_cuda, (a_hx, a_cx))
                    target_action = self.agent.actor_target(xh1)
                    target_action = target_action.detach()
                    next_q_value = self.agent.critic_target([xh1, target_action])

            elif self.rnn_mode == 'gru':
                xh0, _ = self.agent.rnn(state0_cuda, a_hx)
                current_q = self.agent.critic([xh0, action.cuda()])

                with torch.no_grad():
                    xh1, _ = self.agent.rnn_target(state1_cuda, a_hx)
                    target_action = self.agent.actor_target(xh1)
                    target_action = target_action.detach()
                    next_q_value = self.agent.critic_target([xh1, target_action])

            target_q = reward + (1-done) * self.discount * next_q_value.cpu()

            value_loss = 0
            if self.is_PER_replay == False:
                value_loss = F.smooth_l1_loss(current_q, target_q.cuda())
            else:
                value_loss = (F.smooth_l1_loss(current_q, target_q.cuda()) * torch.tensor(weights).cuda()).mean()
            value_loss /= t_len # divide by experience length
            value_loss_total += value_loss
            ####### update Critic per step #######
            self.agent.rnn.zero_grad()
            self.agent.actor.zero_grad()
            self.agent.critic.zero_grad()
            value_loss.backward()
            self.critic_optim.step()
            self.rnn_optim.step()

            ########################## Actor loss calculation & Update Actor ##########################
            ################## Preliminary for loss calculation ##################:
            if t % self.a_update_freq ==0: # update Actor per 3-steps
                if self.rnn_mode == 'lstm':
                    xh_b0, _ = self.agent.rnn(state0_cuda, (a_hx, a_cx))
                    behavior_action = self.agent.actor(xh_b0)
                    ### Estimate actor action ###
                    q_action = self.agent.critic([xh_b0, action.cuda()])  #這邊的critic是behavior_critic，有別於target_critic

                    ### Calculate Actor loss based on Q-value ###
                    actor_loss = -self.agent.critic([xh_b0, behavior_action])

                    ### calculate q_actor_loss for priority
                    q_actor_loss = self.agent.critic([xh_b0, behavior_action])

                    ##### Behavior Cloning Loss #####
                    if self.is_BClone:
                        ### Estimate prophetic action ###
                        q_action_bc = self.agent.critic([xh_b0, action_bc.cuda()])

                        ### Q_filter & BC_loss ###
                        BC_loss = self.BC_loss_func(behavior_action, action_bc.cuda())
                        BC_loss = torch.sum(BC_loss,dim=1).unsqueeze(1)

                        Q_filter = torch.gt(q_action_bc, q_action)
                        BC_loss_Qf = BC_loss * (Q_filter.detach())
                        if self.is_Qfilt:
                            ### modified Policy loss ###
                            policy_loss = (self.lambda_Policy*actor_loss) + (self.lambda_BC*BC_loss_Qf)
                        else:
                            ### modified Policy loss ###
                            policy_loss = (self.lambda_Policy*actor_loss) + (self.lambda_BC*BC_loss)

                    else:  ### Original Policy loss ###
                        policy_loss = actor_loss

                elif self.rnn_mode == 'gru':
                    xh_b0, _ = self.agent.rnn(state0_cuda, a_hx)
                    behavior_action = self.agent.actor(xh_b0)
                    ### Estimate actor action ###
                    q_action = self.agent.critic([xh_b0, action.cuda()])  #這邊的critic是behavior_critic，有別於target_critic

                    ### Calculate Actor loss based on Q-value ###
                    behavior_action = self.agent.actor(xh_b0)
                    actor_loss = -self.agent.critic([xh_b0, behavior_action])

                    ### calculate q_actor_loss for priority
                    q_actor_loss = self.agent.critic([xh_b0, behavior_action])

                    ##### Behavior Cloning Loss #####
                    if self.is_BClone:
                        ### Estimate prophetic action ###
                        q_action_bc = self.agent.critic([xh_b0, action_bc.cuda()])

                        ### Q_filter & BC_loss ###
                        BC_loss = self.BC_loss_func(behavior_action, action_bc.cuda())
                        BC_loss = torch.sum(BC_loss,dim=1).unsqueeze(1)

                        Q_filter = torch.gt(q_action_bc, q_action)
                        BC_loss_Qf = BC_loss * (Q_filter.detach())
                        if self.is_Qfilt:
                            ### modified Policy loss ###
                            policy_loss = (self.lambda_Policy*actor_loss) + (self.lambda_BC*BC_loss_Qf)
                        else:
                            ### modified Policy loss ###
                            policy_loss = (self.lambda_Policy*actor_loss) + (self.lambda_BC*BC_loss)

                    else:  ### Original Policy loss ###
                        policy_loss = actor_loss

                ################## Actor loss calculation ##################
                if self.is_BClone:
                    BC_loss /= t_len
                    BC_loss_total +=  BC_loss.mean()  #BC loss
                    BC_loss_Qf  /= t_len
                    BC_loss_Qf_total += BC_loss_Qf.mean()
                    actor_loss /= t_len
                    actor_loss_total += actor_loss.mean()   #actor loss
                else:
                    BC_loss_total = torch.zeros(1)
                    BC_loss_Qf_total = torch.zeros(1)
                    actor_loss_total = torch.zeros(1)

                policy_loss /= t_len # divide by experience length
                policy_loss_total += policy_loss.mean()

                ####### Update Actor ###########
                self.agent.rnn.zero_grad()
                self.agent.actor.zero_grad()
                self.agent.critic.zero_grad()
                policy_loss = policy_loss.mean()
                policy_loss.backward()
                self.actor_optim.step()
                self.rnn_optim.step()

        ##### Learning rate Scheduling #####
        self.rnn_scheduler.step()
        self.critic_scheduler.step()
        self.actor_scheduler.step()

        ###Update priority###
        if self.is_PER_replay == True:
            TDerror_square = (target_q.cpu() - current_q.cpu()).pow(2)
            loss2actor_square = q_actor_loss.cpu().pow(2)

            self.priority = (TDerror_square + self.lambda_balance*loss2actor_square).detach().numpy().ravel() + self.small_const
            self.priority[batch_flags == DEMO_flag] += self.priority_const

            self.memory.update_priorities(idxes, self.priority)
            self.demoN_ratio = np.sum(demo_cnt) / t_len

        ########### Record all losses ############
        self.actor_loss = actor_loss_total.item()
        self.BC_loss = BC_loss_total.item()
        self.BC_loss_Qf = BC_loss_Qf_total.item()
        self.tot_policy_loss = policy_loss_total.item()
        self.critic_loss = value_loss_total.item()
        ##### 以下的batch期望值是整條segmt_traj，但因為有照time_order後，會不會就不符i.i.d.了？ #####
        ##### update once after experience (segmt_traj) ##### (雖有問題，但也值得try看看)
#        policy_loss_total /= self.batch_size # divide by number of trajectories
#        value_loss_total /= self.batch_size # divide by number of trajectories
#
#        self.agent.critic.zero_grad()
#        value_loss_total.backward()
#        self.critic_optim.step()
#
#        self.agent.actor.zero_grad()
#        policy_loss_total.backward()
#        self.actor_optim.step()

        ##### Target_Net update #####
        soft_update(self.agent.rnn_target, self.agent.rnn, self.tau)
        soft_update(self.agent.actor_target, self.agent.actor, self.tau)
        soft_update(self.agent.critic_target, self.agent.critic, self.tau)

    def test(self, model_path, model_fn, description, lackM=True, debug=False):

        if self.agent.load_weights(model_path, model_fn) == False:
            prRed("model path not found")
            return

        self.agent.is_training = False
        self.agent.eval()

        with torch.no_grad():
            test_mean_reward = self.evaluate(self.env, self.agent, description, lackM=lackM, debug=debug, save=False)
            if debug: prYellow('[Evaluate]: mean_reward:{}'.format(test_mean_reward))

    def seed(self,s):
        torch.manual_seed(s)
        if USE_CUDA:
            torch.cuda.manual_seed(s)

    def train_plot_bc(self, episode, train_epi_reward, train_ewma_reward,
                      train_totPolicy_loss, train_critic_loss,
                      train_actor_loss, train_bc_loss, train_bcQf_loss):
        font_size = 16
        plt.figure(num=1, figsize=(12, 6))

        plt.subplot(321)
        plt.title('Episode Reward', fontsize=font_size)
        plt.plot(train_epi_reward)

        plt.subplot(322)
        plt.title('EWMA Reward', fontsize=font_size)
        plt.plot(train_ewma_reward)

        plt.subplot(323)
        plt.title('total Policy Loss', fontsize=font_size)
        plt.plot(train_totPolicy_loss)

        plt.subplot(324)
        plt.title('Critic Loss', fontsize=font_size)
        plt.plot(train_critic_loss)

        plt.subplot(325)
        plt.title('Actor Loss', fontsize=font_size)
        plt.plot(train_actor_loss)

        plt.subplot(326)
        plt.title('BC Loss', fontsize=font_size)
        bc, = plt.plot(train_bc_loss, label='BC')
        Qf, = plt.plot(train_bcQf_loss, label='Qf')
        plt.legend(handles=[bc, Qf], loc='upper center', fontsize=10)
        # plt.legend([bc, Qf], ["BC", "Qf"], loc='upper left') #, facecolor='blue')
        # plt.legend(bbox_to_anchor=(1.05, 1.0, 0.3, 0.2), loc='upper left')

        plt.tight_layout()
        train_his_fn = 'lamBC_' +str(np.round(self.lambda_BC,3)) +'_' +self.rnn_mode +'_' +str(self.date)
        plt.savefig('results/TrainCurve_epi' +str(episode) +'_' +train_his_fn +'.jpg')
        plt.close()

        ##### Save Training History to csv file #####
        # epi_r = np.array(train_epi_reward)
        # ewma_r = np.array(train_ewma_reward)
        # totP_loss = np.array(train_totPolicy_loss)
        # c_loss = np.array(train_critic_loss)
        # a_loss = np.array(train_actor_loss)
        # bc_loss = np.array(train_bc_loss)
        # qf_loss = np.array(train_bcQf_loss)

        dic = {'Episode Reward':train_epi_reward,
               'EWMA Reward':train_ewma_reward,
               'total Policy Loss':train_totPolicy_loss,
               'Critic Loss':train_critic_loss,
               'Actor Loss':train_actor_loss,
               'BC Loss':train_bc_loss,
               'BC_Qf Loss':train_bcQf_loss
              }
        train_history = pd.DataFrame(dic)
        train_history.to_csv('results/TrainHis_' +train_his_fn +'.csv',  index=False)

    def train_plot(self, episode, train_epi_reward, train_ewma_reward, train_totPolicy_loss, train_critic_loss):
        font_size = 16
        plt.figure(num=1, figsize=(12, 7))

        ax1 = plt.subplot(221)
        ax1.set_title('Episode Reward', fontsize=font_size)
        plt.plot(train_epi_reward)

        plt.subplot(222)
        plt.title('EWMA Reward', fontsize=font_size)
        plt.plot(train_ewma_reward)

        plt.subplot(223)
        plt.title('total Policy Loss', fontsize=font_size)
        plt.plot(train_totPolicy_loss)

        plt.subplot(224)
        plt.title('Critic Loss', fontsize=font_size)
        plt.plot(train_critic_loss)

        plt.tight_layout()
        train_his_fn = '_' +self.rnn_mode +'_' +str(self.date)
        plt.savefig('results/TrainCurve_epi' +str(episode) +train_his_fn +'.jpg')
        plt.close()

        ##### Save Training History to csv file #####
        dic = {'Episode Reward':train_epi_reward,
               'EWMA Reward':train_ewma_reward,
               'total Policy Loss':train_totPolicy_loss,
               'Critic Loss':train_critic_loss,
              }
        train_history = pd.DataFrame(dic)
        train_history.to_csv('results/TrainHis_' +train_his_fn +'.csv',  index=False)

    def train_demoN_ratio(self, episode, demoN_ratio_batch):
        plt.xlabel('# of episode')
        plt.ylabel('batch_demoN_ratio')
        plt.plot(demoN_ratio_batch)
        fig_fn = 'batch_demoN_ratio_ep' +'.jpg'
        plt.savefig(fig_fn)
        plt.close()


text:
# Copyright 2017 Division of Medical Image Computing, German Cancer Research Center (DKFZ)
#
# 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.

from __future__ import print_function

from typing import List, Union

from future import standard_library
import threading

standard_library.install_aliases()
from builtins import range
from multiprocessing import Process
from multiprocessing import Queue
from queue import Queue as thrQueue
import numpy as np
import sys
import logging
from multiprocessing import Event
from time import sleep, time
from threadpoolctl import threadpool_limits

try:
    import torch
except ImportError:
    torch = None

def producer(queue, data_loader, transform, thread_id, seed, abort_event, wait_time: float = 0.02):
    np.random.seed(seed)
    data_loader.set_thread_id(thread_id)
    item = None

    try:
        while True:
            # check if abort event was set
            if not abort_event.is_set():
                # print("worker %d event not set" % thread_id)
                if item is None:
                    try:
                        item = next(data_loader)
                        if transform is not None:
                            item = transform(**item)
                    except StopIteration:
                        item = "end"

                if not queue.full():
                    queue.put(item)
                    item = None
                else:
                    sleep(wait_time)
            else:
                # print("worder %d event is now set, exiting" % thread_id)
                return
    except KeyboardInterrupt:
        abort_event.set()
        return
    except Exception as e:
        print("Exception in background worker %d:\n" % thread_id, e)
        abort_event.set()
        return

def results_loop(in_queues: List[Queue], out_queue: thrQueue, abort_event: Event, pin_memory: bool,
                 gpu: Union[int, None] = None, wait_time: float = 0.02):
    do_pin_memory = torch is not None and pin_memory and gpu is not None and torch.cuda.is_available()

    if do_pin_memory:
        print('using pin_memory on device', gpu)
        torch.cuda.set_device(gpu)

    item = None
    queue_ctr = 0
    end_ctr = 0

    while True:
        # if abort_event is set we need to clean up. This is where it hangs sometimes so it makes sense to drain all
        # the incoming queues and ignore all the errors occuring during this process.
        try:
            if abort_event.is_set():
                # print('abort event is set')
                return

            # if we don't have an item we need to fetch it first. If the queue we want to get it from it empty, try
            # again later
            if item is None:
                current_queue = in_queues[queue_ctr % len(in_queues)]
                if not current_queue.empty():
                    # get the item
                    item = current_queue.get()
                    # if we do pin memory, do it now, otherwise skip this
                    if do_pin_memory:
                        if isinstance(item, dict):
                            for k in item.keys():
                                if isinstance(item[k], torch.Tensor):
                                    item[k] = item[k].pin_memory()
                    queue_ctr += 1

                    if isinstance(item, str) and item == 'end':
                        end_ctr += 1
                    if end_ctr == len(in_queues):
                        end_ctr = 0
                        queue_ctr = 0

                else:
                    sleep(wait_time)
                    continue

            # we only arrive here if item is not None. Now put item in to the out_queue
            if not out_queue.full():
                out_queue.put(item)
                item = None
            else:
                sleep(wait_time)
                continue
        except KeyboardInterrupt:
            abort_event.set()
            raise KeyboardInterrupt

class MultiThreadedAugmenter(object):
    """ Makes your pipeline multi threaded. Yeah!
    If seeded we guarantee that batches are retunred in the same order and with the same augmentation every time this
    is run. This is realized internally by using une queue per worker and querying the queues one ofter the other.
    Args:
        data_loader (generator or DataLoaderBase instance): Your data loader. Must have a .next() function and return
        a dict that complies with our data structure
        transform (Transform instance): Any of our transformations. If you want to use multiple transformations then
        use our Compose transform! Can be None (in that case no transform will be applied)
        num_processes (int): number of processes
        num_cached_per_queue (int): number of batches cached per process (each process has its own
        multiprocessing.Queue). We found 2 to be ideal.
        seeds (list of int): one seed for each worker. Must have len(num_processes).
        If None then seeds = range(num_processes)
        pin_memory (bool): set to True if all torch tensors in data_dict are to be pinned. Pytorch only.
        timeout (int): How long do we wait for the background workers to do stuff? If timeout seconds have passed and
        self.__get_next_item still has not gotten an item from the workers we will perform a check whether all
        background workers are still alive. If all are alive we wait, if not we set the abort flag.
        wait_time (float): set this to be lower than the time you need per iteration. Don't set this to 0,
        that will come with a performance penalty. Default is 0.02 which will be fine for 50 iterations/s
    """

    def __init__(self, data_loader, transform, num_processes, num_cached_per_queue=2, seeds=None, pin_memory=False,
                 timeout=10, wait_time=0.02):
        self.timeout = timeout
        self.pin_memory = pin_memory
        self.transform = transform
        if seeds is not None:
            assert len(seeds) == num_processes
        else:
            seeds = [None] * num_processes
        self.seeds = seeds
        self.generator = data_loader
        self.num_processes = num_processes
        self.num_cached_per_queue = num_cached_per_queue
        self._queues = []
        self._processes = []
        self._end_ctr = 0
        self._queue_ctr = 0
        self.pin_memory_thread = None
        self.pin_memory_queue = None
        self.abort_event = Event()
        self.wait_time = wait_time

    def __iter__(self):
        return self

    def next(self):
        return self.__next__()

    def __get_next_item(self):
        item = None

        while item is None:
            if self.abort_event.is_set():
                self._finish()
                raise RuntimeError("MultiThreadedAugmenter.abort_event was set, something went wrong. Maybe one of "
                                   "your workers crashed")

            if not self.pin_memory_queue.empty():
                item = self.pin_memory_queue.get()
            else:
                sleep(self.wait_time)

        return item

    def __next__(self):
        if len(self._queues) == 0 or len(self._processes) != self.num_processes or not all([i.is_alive() for i
                                                                                in self._processes]):
            self._start()
        try:
            item = self.__get_next_item()

            while isinstance(item, str) and (item == "end"):
                self._end_ctr += 1
                if self._end_ctr == self.num_processes:
                    self._end_ctr = 0
                    self._queue_ctr = 0
                    logging.debug("MultiThreadedGenerator: finished data generation")
                    raise StopIteration

                item = self.__get_next_item()

            return item

        except KeyboardInterrupt:
            logging.error("MultiThreadedGenerator: caught exception: {}".format(sys.exc_info()))
            self.abort_event.set()
            self._finish()
            raise KeyboardInterrupt

    def _start(self):
        if len(self._processes) == 0:
            self.abort_event.clear()

            logging.debug("starting workers")
            self._queue_ctr = 0
            self._end_ctr = 0

            if hasattr(self.generator, 'was_initialized'):
                self.generator.was_initialized = False

            with threadpool_limits(limits=1, user_api="blas"):
                for i in range(self.num_processes):
                    self._queues.append(Queue(self.num_cached_per_queue))
                    self._processes.append(Process(target=producer, args=(
                        self._queues[i], self.generator, self.transform, i, self.seeds[i], self.abort_event)))
                    self._processes[-1].daemon = True
                    self._processes[-1].start()

            if torch is not None and torch.cuda.is_available():
                gpu = torch.cuda.current_device()
            else:
                gpu = None

            # more caching = more performance. But don't cache too much or your RAM will hate you
            self.pin_memory_queue = thrQueue(max(3, self.num_cached_per_queue * self.num_processes // 2))

            self.pin_memory_thread = threading.Thread(target=results_loop, args=(
                self._queues, self.pin_memory_queue, self.abort_event, self.pin_memory, gpu, self.wait_time))

            self.pin_memory_thread.daemon = True
            self.pin_memory_thread.start()
        else:
            logging.debug("MultiThreadedGenerator Warning: start() has been called but workers are already running")

    def _finish(self, timeout=10):
        self.abort_event.set()

        start = time()
        while self.pin_memory_thread is not None and self.pin_memory_thread.is_alive() and start + timeout > time():
            sleep(0.2)

        if len(self._processes) != 0:
            logging.debug("MultiThreadedGenerator: shutting down workers...")
            [i.terminate() for i in self._processes]

            for i, p in enumerate(self._processes):
                self._queues[i].close()
                self._queues[i].join_thread()

            self._queues = []
            self._processes = []
            self._queue = None
            self._end_ctr = 0
            self._queue_ctr = 0

            del self.pin_memory_queue

    def restart(self):
        self._finish()
        self._start()

    def __del__(self):
        logging.debug("MultiThreadedGenerator: destructor was called")
        self._finish()

ext: py

max_stars_repo_name: zhongzisha/batchgenerators

max_issues_repo_name: zhongzisha/batchgenerators

max_forks_repo_name: zhongzisha/batchgenerators

reason: import numpy


text:
"""Functions for analyzing output results dictionary"""
import itertools

import numpy as np

from gradient_loop import make_hash

def calculate_new_quantity(in_keys: list, out_name: str, f, results: dict, hp_file):
    """Calculate new quantity for each data point.

    Calculates the function ``f`` applied to the input ``keys`` and creates a new entry in the
    dictionary with name ``q_name``.
    """
    for i, values in enumerate(itertools.product(*hp_file.hyperparameters.values())):
        hp = {**dict(zip(hp_file.hyperparameters.keys(), values)), **hp_file.log_parameters}
        name = make_hash(hp)

        try:
            r = results[name]
        except KeyError:
            continue

        in_q = [r[k] for k in in_keys]
        out_q = f(*in_q)
        results[name][out_name] = out_q

def avg_quantities(keys: list, results: dict, hp_file):
    """Average quantities specified in ``keys`` over the range of seed values."""
    hps = hp_file.hyperparameters.copy()  # Load hyperparameters
    seed = hps["seed"]  # Record seed list
    del hps["seed"]  # Remove seed from hyperparameters

    results_avg = {}  # Create new dictionary to hold averaged quantities

    for i, values in enumerate(itertools.product(*hps.values())):
        hps_no_seed = dict(zip(hps.keys(), values))
        hps_no_seed_plus_log = {
            **hps_no_seed,
            **hp_file.log_parameters,
        }  # Add back in log parameters needed for hash

        q = []

        for s in seed:
            hp_with_seed = {**{"seed": s}, **hps_no_seed_plus_log}  # Add back in seed within loop
            name = make_hash(hp_with_seed)  # Find corresponding hash
            try:
                q.append([results[name][k] for k in keys])
            except KeyError:
                continue

        hp_not_looped = [
            k for k, v in hps.items() if len(v) == 1
        ]  # Find hyperparameters not looped over
        for h in hp_not_looped:  # Remove hyperparameters from dictionary that are not looped over
            del hps_no_seed[h]

        q = np.array(q)
        q_mean = np.mean(q, axis=0)

        results_avg[make_hash(hps_no_seed)] = {
            **hps_no_seed,
            **{k: q_mean[i] for i, k in enumerate(keys)},
        }

    return results_avg

def access_quantities(keys: list, results: dict, hp_file, average: bool = False):
    """Access quantities specified in ``keys`` over the range of seed values."""
    hps = hp_file.hyperparameters.copy()  # Load hyperparameters
    seed = hps["seed"]  # Record seed list
    del hps["seed"]  # Remove seed from hyperparameters

    results_avg = {}  # Create new dictionary to hold averaged quantities

    for i, values in enumerate(itertools.product(*hps.values())):
        hps_no_seed = dict(zip(hps.keys(), values))
        hps_no_seed_plus_log = {
            **hps_no_seed,
            **hp_file.log_parameters,
        }  # Add back in log parameters needed for hash

        q = []

        for s in seed:
            hp_with_seed = {**{"seed": s}, **hps_no_seed_plus_log}  # Add back in seed within loop
            name = make_hash(hp_with_seed)  # Find corresponding hash
            try:
                q.append([results[name][k] for k in keys])
            except KeyError:
                continue

        hp_not_looped = [
            k for k, v in hps.items() if len(v) == 1
        ]  # Find hyperparameters not looped over
        for h in hp_not_looped:  # Remove hyperparameters from dictionary that are not looped over
            del hps_no_seed[h]

        q = np.array(q)
        if average:
            q = np.mean(q, axis=0)

        results_avg[make_hash(hps_no_seed)] = {
            **hps_no_seed,
            **{k: q[:, i] for i, k in enumerate(keys)},
        }

    return results_avg

def calculate_slice(pinned: dict, results: dict):
    """Extracts all data points from ``results`` with pinned values specified by ``pinned``."""
    slice_d = {}

    for k, v in results.items():
        v_c = v.copy()
        contained = True

        for k_p, v_p in pinned.items():
            if v_c[k_p] != v_p:
                contained = False
            else:
                del v_c[k_p]

        if contained:
            slice_d.update({k: v_c})

    return slice_d

def make_numpy(results_slice: dict, x: str, y: str):
    """Make a slice into two numpy arrays suitable for plotting, with the x and y axes specified
    as key strings."""
    slice_list = sorted([(v[x], v[y]) for k, v in results_slice.items()])
    return np.array([s_l[0] for s_l in slice_list]), np.array([s_l[1] for s_l in slice_list])

def make_numpy2d(results_slice: dict, x: str, y: str, z: str, z_elem: int = None):
    """Make a slice into three numpy arrays suitable for 2D plotting, with the x, y, and z axes
    specified as key strings."""
    slice_list = sorted([(v[x], v[y], v[z]) for k, v in results_slice.items()])
    if z_elem:
        return (
            np.array([s_l[0] for s_l in slice_list]),
            np.array([s_l[1] for s_l in slice_list]),
            np.array([s_l[2][z_elem] for s_l in slice_list]),
        )

    return (
        np.array([s_l[0] for s_l in slice_list]),
        np.array([s_l[1] for s_l in slice_list]),
        np.array([s_l[2] for s_l in slice_list]),
    )


text:
import numpy as np
import pytest

from psyrun.pspace import AmbiguousOperationError, dict_concat, missing, Param

@pytest.mark.parametrize('args,result', [
    (({0: 0}, {0: 1}, {0: 2}), {0: [0, 1, 2]}),
    (({0: 2},), {0: [2]}),
    (({},), {}),
    (({0: 0, 1: 1}, {0: 0, 1: 1}), {0: [0, 0], 1: [1, 1]}),
    ((), {}),
    (({0: 0}, {1: 1}), {0: [0, None], 1: [None, 1]})
])
def test_dict_concat(args, result):
    assert dict_concat(args) == result

class TestParam(object):
    def test_empty(self):
        assert Param().build() == {}

    def test_single(self):
        space = Param(a=[1, 2]).build()
        assert sorted(space['a']) == [1, 2]

    def test_multiple(self):
        space = Param(a=[1, 2], b=[3, 4]).build()
        assert sorted(space['a']) == [1, 2]
        assert sorted(space['b']) == [3, 4]

    def test_multiple_unequal_length_raises_error(self):
        with pytest.raises(ValueError):
            Param(a=[1, 2], b=[3])

    def test_only_scalars(self):
        space = Param(a=1, b=2).build()
        assert space['a'] == [1]
        assert space['b'] == [2]

    def test_treats_string_not_as_sequence(self):
        assert Param(a='item').build()['a'] == ['item']
        assert Param(a=['item']).build()['a'] == ['item']

    def test_scalars_are_broadcasted(self):
        space = Param(a=[1, 2], b=3).build()
        assert sorted(space['a']) == [1, 2]
        assert sorted(space['b']) == [3, 3]

    def test_knows_about_its_keys(self):
        assert sorted(Param(a=[0], b=[1]).keys()) == ['a', 'b']

    def test_length(self):
        assert len(Param(a=[1, 2, 3], b=3)) == 3

    def test_keeps_dtype(self):
        space = Param(a=[1], b=[2.]).build()
        assert np.issubdtype(type(space['a'][0]), np.signedinteger)
        assert np.issubdtype(type(space['b'][0]), np.floating)

    def test_from_dict(self):
        d = {'a': [1, 2], 'b.x': 3, 'b.y': [4, 5]}
        space = Param(**d).build()
        assert space['a'] == [1, 2]
        assert space['b.x'] == [3, 3]
        assert space['b.y'] == [4, 5]

    def test_hierarchical(self):
        d1 = {'root.a': [1, 2], 'b': [3, 4]}
        d2 = {'root.b': [5, 6]}
        space = (Param(**d1) * Param(**d2)).build()
        assert space['root.a'] == [1, 1, 2, 2]
        assert space['b'] == [3, 3, 4, 4]
        assert space['root.b'] == [5, 6, 5, 6]

    def test_str(self):
        p = Param(**{'a': [1, 2], 'b.b': [3, 4]})
        assert str(p) == "Param(a=[1, 2], b.b=[3, 4])"

    def test_rep(self):
        p = Param(**{'a': [1, 2], 'b.b': [3, 4]})
        assert repr(p) == "Param(**{'a': [1, 2], 'b.b': [3, 4]})"

class TestProduct(object):
    def test_product(self):
        space = (Param(a=[1, 2, 3]) * Param(b=[4, 5])).build()
        assert space['a'] == [1, 1, 2, 2, 3, 3]
        assert space['b'] == [4, 5, 4, 5, 4, 5]

    def test_overlapping_product(self):
        with pytest.raises(AmbiguousOperationError):
            _ = Param(a=[1, 2]) * Param(a=[2, 3], b=[4, 5])

    def test_product_with_empty_right(self):
        space = (Param(a=[1, 2, 3]) * Param()).build()
        assert sorted(space['a']) == [1, 2, 3]

    def test_product_with_empty_left(self):
        space = (Param() * Param(a=[1, 2, 3])).build()
        assert sorted(space['a']) == [1, 2, 3]

    def test_product_with_zero_elements(self):
        space = (Param(a=[1, 2, 3]) * Param(b=[])).build()
        assert space['a'] == []
        assert space['b'] == []

    def test_length(self):
        assert len(Param(a=[1, 2, 3]) * Param(b=[1, 2])) == 6
        assert len(Param(a=[1, 2]) * Param()) == 2
        assert len(Param() * Param(a=[1, 2])) == 2
        assert len(Param(a=[1, 2]) * Param(b=[])) == 0
        assert len(Param(a=[]) * Param(b=[1, 2])) == 0

    def test_multidimensional(self):
        space = (Param(a=[[1, 1], [2, 2]]) * Param(b=[1, 2])).build()
        assert space['a'] == [[1, 1], [1, 1], [2, 2], [2, 2]]
        assert space['b'] == [1, 2, 1, 2]

class TestSum(object):
    def test_sum(self):
        space = (Param(a=[1]) + Param(a=[2])).build()
        assert sorted(space['a']) == [1, 2]

    def test_sum_distinct_params(self):
        space = (Param(a=[1]) + Param(b=[2])).build()
        assert space['a'][0] == 1
        assert space['b'][1] == 2
        assert np.isnan(space['a'][1])
        assert np.isnan(space['b'][0])

    def test_sum_with_empty(self):
        space = (Param(a=[1]) + Param()).build()
        assert sorted(space['a']) == [1]

    def test_sum_with_zero_elements(self):
        space = (Param(a=[1]) + Param(b=[])).build()
        assert space['a'] == [1]
        assert np.isnan(space['b'][0])

    def test_length(self):
        assert len(Param(a=[1]) + Param(a=[2, 3])) == 3

    def test_multidimensional(self):
        space = (Param(a=[[1, 1], [2, 2]]) + Param(a=[[3, 3]])).build()
        assert space['a'] == [[1, 1], [2, 2], [3, 3]]

        space = (Param(a=[[1, 1], [2, 2]]) + Param(a=[3, 4])).build()
        assert space['a'] == [[1, 1], [2, 2], 3, 4]

class TestDifference(object):
    def test_difference(self):
        space = (Param(a=[1, 2, 3]) - Param(a=[2])).build()
        assert sorted(space['a']) == [1, 3]

    def test_multiple_builds(self):
        space = (Param(a=[1, 2, 3]) - Param(a=[2]))
        assert sorted(space.build()['a']) == [1, 3]
        assert sorted(space.build()['a']) == [1, 3]

    def test_param_missing_in_subtrahend(self):
        space = (Param(a=[1, 2, 2], b=[4, 5, 6]) - Param(a=[2])).build()
        assert sorted(space['a']) == [1]
        assert sorted(space['b']) == [4]

    def test_param_missign_in_minuend(self):
        with pytest.raises(AmbiguousOperationError):
            _ = Param(a=[1, 2, 3]) - Param(a=[2], b=[3])

    def test_empty_subtrahend(self):
        space = (Param(a=[1, 2, 3]) - Param()).build()
        assert sorted(space['a']) == [1, 2, 3]

    def test_length(self):
        assert len(Param(a=[1, 2, 2], b=[4, 5, 6]) - Param(a=[2])) == 1

    def test_multidimensional(self):
        space = (Param(a=[[1, 1], [2, 2]]) - Param(a=[[2, 2]])).build()
        assert space['a'] == [[1, 1]]

class TestMissing(object):
    def test_missing(self):
        space = missing(Param(a=[1, 2, 3]), Param(a=[2])).build()
        assert sorted(space['a']) == [1, 3]

    def test_param_missing_in_subtrahend(self):
        with pytest.raises(AmbiguousOperationError):
            missing(Param(a=[1, 2, 2], b=[4, 5, 6]), Param(a=[2])).build()

    def test_param_missign_in_minuend(self):
        space = missing(Param(a=[1, 2, 3]), Param(a=[2], b=[3])).build()
        assert sorted(space['a']) == [1, 3]

    def test_empty_subtrahend(self):
        space = missing(Param(a=[1, 2, 3]), Param()).build()
        assert sorted(space['a']) == [1, 2, 3]

    def test_length(self):
        assert len(missing(Param(a=[1, 2, 3]), Param(a=[2]))) == 2

    def test_multidimensional(self):
        space = missing(Param(a=[[1, 1], [2, 2]]), Param(a=[[1, 1]])).build()
        assert space['a'] == [[2, 2]]

max_stars_repo_name: jgosmann/psyrun

max_issues_repo_name: jgosmann/psyrun

max_forks_repo_name: jgosmann/psyrun


text:
import numpy as np
from numpy import array, asarray, pi, sin, cos
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay

nm = 2
ni = 6

radius = 1.
angle = pi / 3.

center = array([0,0])
centers = [center]

points = []

def create_points(X, pts):
    for m in range(1, nm):
        for i in range(0, ni):
            P = m * radius * array([cos(i*angle)-X[0], sin(i*angle)-X[1]])
            pts.append(P)
    return pts

centers = create_points(center, centers)
for X in centers:
    points = create_points(X, points)

points = array(points)
tri = Delaunay(points)

barycenters = []
for T in tri.simplices:
    pts = points[T]
    A = array(pts[0,:]) ; B = array(pts[1,:]) ; C = array(pts[2,:])
    b = (A+B+C) / 3.
    barycenters.append(b)

barycenters = array(barycenters)

tri_b = Delaunay(barycenters)

plt.triplot(points[:,0], points[:,1], tri.simplices.copy())
plt.plot(points[:,0], points[:,1], 'or')
plt.plot(barycenters[:,0], barycenters[:,1], 'ob')

xmin = points[:,0].min()  ;   xmax = points[:,0].max()
ymin = points[:,1].min()  ;   ymax = points[:,1].max()

plt.xlim([xmin-0.2, xmax+0.2])
plt.ylim([ymin-0.2, ymax+0.2])
plt.show()

max_stars_repo_name: ratnania/caid

max_issues_repo_name: pyccel/caid

max_forks_repo_name: pyccel/caid

reason: import numpy,from numpy,from scipy


text:
import logging
import pickle
from time import perf_counter

import numpy as np
import pandas as pd

logger = logging.getLogger(__name__)

def read_table_csv(table_obj, csv_seperator=','):
    """
    Reads csv from path, renames columns and drops unnecessary columns
    """
    df_rows = pd.read_csv(table_obj.csv_file_location, header=None, escapechar='\\', encoding='utf-8', quotechar='"',
                          sep=csv_seperator)
    df_rows.columns = [table_obj.table_name + '.' + attr for attr in table_obj.attributes]

    for attribute in table_obj.irrelevant_attributes:
        df_rows = df_rows.drop(table_obj.table_name + '.' + attribute, axis=1)

    return df_rows.apply(pd.to_numeric, errors="ignore")
    # return df_rows.convert_objects()

def find_relationships(schema_graph, table, incoming=True):
    relationships = []

    for relationship_obj in schema_graph.relationships:

        if relationship_obj.end == table and incoming:
            relationships.append(relationship_obj)
        if relationship_obj.start == table and not incoming:
            relationships.append(relationship_obj)

    return relationships

def prepare_single_table(schema_graph, table, path, max_distinct_vals=10000, csv_seperator=',',
                         max_table_data=20000000):
    """
    Reads table csv. Adds multiplier fields, missing value imputation, dict for categorical data. Adds null tuple tables.

    :param schema_graph:
    :param table:
    :return:
    """
    table_meta_data = dict()
    table_obj = schema_graph.table_dictionary[table]
    table_data = read_table_csv(table_obj, csv_seperator=csv_seperator)
    table_sample_rate = table_obj.sample_rate

    relevant_attributes = [x for x in table_obj.attributes if x not in table_obj.irrelevant_attributes]

    table_meta_data['hdf_path'] = path
    table_meta_data['incoming_relationship_means'] = {}

    # manage functional dependencies
    logger.info(f"Managing functional dependencies for table {table}")
    table_meta_data['fd_dict'] = dict()
    cols_to_be_dropped = []
    for attribute_wo_table in table_obj.attributes:
        attribute = table + '.' + attribute_wo_table
        fd_children = table_obj.children_fd_attributes(attribute)
        if len(fd_children) > 0:
            for child in fd_children:
                logger.info(f"Managing functional dependencies for {child}->{attribute}")
                distinct_tuples = table_data.drop_duplicates([attribute, child])[[attribute, child]].values
                reverse_dict = {}
                for attribute_value, child_value in distinct_tuples:
                    if reverse_dict.get(attribute_value) is None:
                        reverse_dict[attribute_value] = []
                    reverse_dict[attribute_value].append(child_value)
                if table_meta_data['fd_dict'].get(attribute) is None:
                    table_meta_data['fd_dict'][attribute] = dict()
                table_meta_data['fd_dict'][attribute][child] = reverse_dict
            # remove from dataframe and relevant attributes
            cols_to_be_dropped.append(attribute)
            relevant_attributes.remove(attribute_wo_table)
    table_data.drop(columns=cols_to_be_dropped, inplace=True)

    # add multiplier fields
    logger.info("Preparing multipliers for table {}".format(table))
    incoming_relationships = find_relationships(schema_graph, table, incoming=True)

    for relationship_obj in incoming_relationships:
        logger.info("Preparing multiplier {} for table {}".format(relationship_obj.identifier, table))

        neighbor_table = relationship_obj.start
        neighbor_table_obj = schema_graph.table_dictionary[neighbor_table]
        neighbor_sample_rate = neighbor_table_obj.sample_rate

        left_attribute = table + '.' + relationship_obj.end_attr
        right_attribute = neighbor_table + '.' + relationship_obj.start_attr

        neighbor_table_data = read_table_csv(neighbor_table_obj, csv_seperator=csv_seperator).set_index(right_attribute,
                                                                                drop=False)
        table_data = table_data.set_index(left_attribute, drop=False)

        assert len(table_obj.primary_key) == 1, \
            "Currently, only single primary keys are supported for table with incoming edges"
        table_primary_key = table + '.' + table_obj.primary_key[0]
        assert table_primary_key == left_attribute, "Currently, only references to primary key are supported"

        # fix for new pandas version
        table_data.index.name = None
        neighbor_table_data.index.name = None
        muls = table_data.join(neighbor_table_data, how='left')[[table_primary_key, right_attribute]] \
            .groupby([table_primary_key]).count()

        mu_nn_col_name = relationship_obj.end + '.' + relationship_obj.multiplier_attribute_name_nn
        mu_col_name = relationship_obj.end + '.' + relationship_obj.multiplier_attribute_name

        muls.columns = [mu_col_name]
        # if we just have a sample of the neighbor table we assume larger multipliers
        muls[mu_col_name] = muls[mu_col_name] * 1 / neighbor_sample_rate
        muls[mu_nn_col_name] = muls[mu_col_name].replace(to_replace=0, value=1)

        table_data = table_data.join(muls)

        relevant_attributes.append(relationship_obj.multiplier_attribute_name)
        relevant_attributes.append(relationship_obj.multiplier_attribute_name_nn)

        table_meta_data['incoming_relationship_means'][relationship_obj.identifier] = table_data[mu_nn_col_name].mean()

    # save if there are entities without FK reference (e.g. orders without customers)
    outgoing_relationships = find_relationships(schema_graph, table, incoming=False)
    for relationship_obj in outgoing_relationships:
        fk_attribute_name = table + '.' + relationship_obj.start_attr

        table_meta_data[relationship_obj.identifier] = {
            'fk_attribute_name': fk_attribute_name,
            'length': table_data[fk_attribute_name].isna().sum() * 1 / table_sample_rate,
            'path': None
        }

    # null value imputation and categorical value replacement
    logger.info("Preparing categorical values and null values for table {}".format(table))
    table_meta_data['categorical_columns_dict'] = {}
    table_meta_data['null_values_column'] = []
    del_cat_attributes = []

    for rel_attribute in relevant_attributes:

        attribute = table + '.' + rel_attribute

        # categorical value
        if table_data.dtypes[attribute] == object:

            logger.debug("\t\tPreparing categorical values for column {}".format(rel_attribute))

            distinct_vals = table_data[attribute].unique()

            if len(distinct_vals) > max_distinct_vals:
                del_cat_attributes.append(rel_attribute)
                logger.info("Ignoring column {} for table {} because "
                            "there are too many categorical values".format(rel_attribute, table))
            # all values nan does not provide any information
            elif not table_data[attribute].notna().any():
                del_cat_attributes.append(rel_attribute)
                logger.info("Ignoring column {} for table {} because all values are nan".format(rel_attribute, table))
            else:
                if not table_data[attribute].isna().any():
                    val_dict = dict(zip(distinct_vals, range(1, len(distinct_vals) + 1)))
                    val_dict[np.nan] = 0
                else:
                    val_dict = dict(zip(distinct_vals, range(1, len(distinct_vals) + 1)))
                    val_dict[np.nan] = 0
                table_meta_data['categorical_columns_dict'][attribute] = val_dict

                table_data[attribute] = table_data[attribute].map(val_dict.get)
                # because we are paranoid
                table_data[attribute] = table_data[attribute].fillna(0)
                # apparently slow
                # table_data[attribute] = table_data[attribute].replace(val_dict)
                table_meta_data['null_values_column'].append(val_dict[np.nan])

        # numerical value
        else:

            logger.debug("\t\tPreparing numerical values for column {}".format(rel_attribute))

            # all nan values
            if not table_data[attribute].notna().any():
                del_cat_attributes.append(rel_attribute)
                logger.info("Ignoring column {} for table {} because all values are nan".format(rel_attribute, table))
            else:
                contains_nan = table_data[attribute].isna().any()

                # not the best solution but works
                unique_null_val = table_data[attribute].mean() + 0.0001
                assert not (table_data[attribute] == unique_null_val).any()

                table_data[attribute] = table_data[attribute].fillna(unique_null_val)
                table_meta_data['null_values_column'].append(unique_null_val)

                if contains_nan:
                    assert (table_data[attribute] == unique_null_val).any(), "Null value cannot be found"

    # remove categorical columns with too many entries from relevant tables and dataframe
    relevant_attributes = [x for x in relevant_attributes if x not in del_cat_attributes]
    logger.info("Relevant attributes for table {} are {}".format(table, relevant_attributes))
    logger.info("NULL values for table {} are {}".format(table, table_meta_data['null_values_column']))
    del_cat_attributes = [table + '.' + rel_attribute for rel_attribute in del_cat_attributes]
    table_data = table_data.drop(columns=del_cat_attributes)

    assert len(relevant_attributes) == len(table_meta_data['null_values_column']), \
        "Length of NULL values does not match"
    table_meta_data['relevant_attributes'] = relevant_attributes
    table_meta_data['relevant_attributes_full'] = [table + '.' + attr for attr in relevant_attributes]
    table_meta_data['length'] = len(table_data) * 1 / table_sample_rate

    assert not table_data.isna().any().any(), "Still contains null values"

    # save modified table
    if len(table_data) < max_table_data:
        table_data.to_hdf(path, key='df', format='table')
    else:
        table_data.sample(max_table_data).to_hdf(path, key='df', format='table')

    # add table parts without join partners
    logger.info("Adding table parts without join partners for table {}".format(table))
    for relationship_obj in incoming_relationships:
        logger.info("Adding table parts without join partners "
                    "for table {} and relationship {}".format(table, relationship_obj.identifier))

        neighbor_table = relationship_obj.start
        neighbor_table_obj = schema_graph.table_dictionary[neighbor_table]
        neighbor_primary_key = neighbor_table + '.' + neighbor_table_obj.primary_key[0]

        left_attribute = table + '.' + relationship_obj.end_attr
        right_attribute = neighbor_table + '.' + relationship_obj.start_attr

        table_data = table_data.set_index(left_attribute, drop=False)
        neighbor_table_data = read_table_csv(neighbor_table_obj, csv_seperator=csv_seperator).set_index(right_attribute,
                                                                                drop=False)
        null_tuples = table_data.join(neighbor_table_data, how='left')
        null_tuples = null_tuples.loc[null_tuples[neighbor_primary_key].isna(),
                                      [table + '.' + attr for attr in relevant_attributes]]
        if len(null_tuples) > 0 and neighbor_table_obj.sample_rate < 1:
            logger.warning(f"For {relationship_obj.identifier} {len(null_tuples)} tuples without a join partner were "
                           f"found. This is potentially due to the sampling rate of {neighbor_table_obj.sample_rate}.")

        if len(null_tuples) > 0:
            null_tuple_path = path + relationship_obj.start + relationship_obj.start_attr + '.hdf'
            table_meta_data[relationship_obj.identifier] = {
                'length': len(null_tuples) * 1 / table_sample_rate,
                'path': null_tuple_path
            }
            null_tuples.to_hdf(null_tuple_path, key='df', format='table')

    return table_meta_data

def prepare_all_tables(schema_graph, path, csv_seperator=',', max_table_data=20000000):
    prep_start_t = perf_counter()
    meta_data = {}
    for table_obj in schema_graph.tables:
        table = table_obj.table_name
        logger.info("Preparing hdf file for table {}".format(table))
        meta_data[table] = prepare_single_table(schema_graph, table, path + '/' + table + '.hdf',
                                                csv_seperator=csv_seperator, max_table_data=max_table_data)

    with open(path + '/meta_data.pkl', 'wb') as f:
        pickle.dump(meta_data, f, pickle.HIGHEST_PROTOCOL)
    prep_end_t = perf_counter()

    with open(path + '/build_time_hdf.txt', "w") as text_file:
        text_file.write(str(round(prep_end_t-prep_start_t)))

    return meta_data

max_stars_repo_name: LumingSun/deepdb-public

max_issues_repo_name: LumingSun/deepdb-public

max_forks_repo_name: jschj/deepdb-public


text:
"""
Actor-Critic
-------------
It uses TD-error as the Advantage.

Actor Critic History
----------------------
A3C > DDPG > AC

Advantage
----------
AC converge faster than Policy Gradient.

Disadvantage (IMPORTANT)
------------------------
The Policy is oscillated (difficult to converge), DDPG can solve
this problem using advantage of DQN.

Reference
----------
paper: https://papers.nips.cc/paper/1786-actor-critic-algorithms.pdf
View more on MorvanZhou's tutorial page: https://morvanzhou.github.io/tutorials/

Environment
------------
CartPole-v0: https://gym.openai.com/envs/CartPole-v0

A pole is attached by an un-actuated joint to a cart, which moves along a
frictionless track. The system is controlled by applying a force of +1 or -1
to the cart. The pendulum starts upright, and the goal is to prevent it from
falling over.

A reward of +1 is provided for every timestep that the pole remains upright.
The episode ends when the pole is more than 15 degrees from vertical, or the
cart moves more than 2.4 units from the center.

Prerequisites
--------------
tensorflow >=2.0.0a0
tensorlayer >=2.0.0

"""
import argparse
import time

import gym
import numpy as np
import tensorflow as tf

import tensorlayer as tl
from common.buffer import *
from common.networks import *
from common.utils import *
from tensorlayer.models import Model

tl.logging.set_verbosity(tl.logging.DEBUG)

###############################  Actor-Critic  ####################################

class Actor(object):

    def __init__(self, n_features, n_actions, lr, hidden_dim, hidden_layer):

        self.model = DeterministicPolicyNetwork(n_features, n_actions, hidden_dim, hidden_layer).model()
        self.model.train()
        self.optimizer = tf.optimizers.Adam(lr)

    def learn(self, s, a, td):
        with tf.GradientTape() as tape:
            _logits = self.model(np.array([s]))
            ## cross-entropy loss weighted by td-error (advantage),
            # the cross-entropy mearsures the difference of two probability distributions: the predicted logits and sampled action distribution,
            # then weighted by the td-error: small difference of real and predict actions for large td-error (advantage); and vice versa.
            _exp_v = tl.rein.cross_entropy_reward_loss(logits=_logits, actions=[a], rewards=td[0])
        grad = tape.gradient(_exp_v, self.model.trainable_weights)
        self.optimizer.apply_gradients(zip(grad, self.model.trainable_weights))
        return _exp_v

    def choose_action(self, s):
        _logits = self.model(np.array([s]))
        _probs = tf.nn.softmax(_logits).numpy()
        return tl.rein.choice_action_by_probs(_probs.ravel())  # sample according to probability distribution

    def choose_action_greedy(self, s):
        _logits = self.model(np.array([s]))  # logits: probability distribution of actions
        _probs = tf.nn.softmax(_logits).numpy()
        return np.argmax(_probs.ravel())

    def save_ckpt(self):  # save trained weights
        save_model(self.model, 'model_actor', 'AC')

    def load_ckpt(self):  # load trained weights
        load_model(self.model, 'model_actor', 'AC')

class Critic(object):

    def __init__(self, n_features, gamma, lr, hidden_dim, hidden_layer):
        self.GAMMA = gamma

        self.model = ValueNetwork(n_features, hidden_dim, hidden_layer).model()  # from common.networks
        self.model.train()

        self.optimizer = tf.optimizers.Adam(lr)

    def learn(self, s, r, s_):
        v_ = self.model(np.array([s_]))
        with tf.GradientTape() as tape:
            v = self.model(np.array([s]))
            ## TD_error = r + lambd * V(newS) - V(S)
            td_error = r + self.GAMMA * v_ - v
            loss = tf.square(td_error)
        grad = tape.gradient(loss, self.model.trainable_weights)
        self.optimizer.apply_gradients(zip(grad, self.model.trainable_weights))

        return td_error

    def save_ckpt(self):  # save trained weights
        save_model(self.model, 'model_critic', 'AC')

    def load_ckpt(self):  # load trained weights
        load_model(self.model, 'model_critic', 'AC')

def learn(
        env_id, train_episodes, test_episodes=1000, max_steps=1000, gamma=0.9, actor_lr=1e-3, critic_lr=1e-2,
        actor_hidden_dim=30, actor_hidden_layer=1, critic_hidden_dim=30, critic_hidden_layer=1, seed=2,
        save_interval=100, mode='train', render=False
):
    '''
    parameters
    -----------
    env: learning environment
    train_episodes:  total number of episodes for training
    test_episodes:  total number of episodes for testing
    max_steps:  maximum number of steps for one episode
    number_workers: manually set number of workers
    gamma: reward discount factor
    actor_lr: learning rate for actor
    critic_lr: learning rate for critic
    save_interval: timesteps for saving the weights and plotting the results
    mode: train or test

    '''

    env = make_env(env_id)
    env.seed(seed)  # reproducible
    np.random.seed(seed)
    tf.random.set_seed(seed)  # reproducible
    # env = env.unwrapped
    N_F = env.observation_space.shape[0]
    # N_A = env.action_space.shape[0]
    N_A = env.action_space.n

    print("observation dimension: %d" % N_F)  # 4
    print("observation high: %s" % env.observation_space.high)  # [ 2.4 , inf , 0.41887902 , inf]
    print("observation low : %s" % env.observation_space.low)  # [-2.4 , -inf , -0.41887902 , -inf]
    print("num of actions: %d" % N_A)  # 2 : left or right

    actor = Actor(n_features=N_F, n_actions=N_A, lr=actor_lr, hidden_dim=actor_hidden_dim,\
    hidden_layer = actor_hidden_layer)
    # we need a good teacher, so the teacher should learn faster than the actor
    critic = Critic(n_features=N_F, gamma=gamma, lr=critic_lr, hidden_dim=critic_hidden_dim,\
    hidden_layer = critic_hidden_layer)

    if mode == 'train':
        t0 = time.time()
        rewards = []
        for i_episode in range(train_episodes):
            # episode_time = time.time()
            s = env.reset().astype(np.float32)
            t = 0  # number of step in this episode
            all_r = []  # rewards of all steps

            while True:

                if render: env.render()

                a = actor.choose_action(s)

                s_new, r, done, info = env.step(a)
                s_new = s_new.astype(np.float32)

                if done: r = -20

                all_r.append(r)

                td_error = critic.learn(
                    s, r, s_new
                )  # learn Value-function : gradient = grad[r + lambda * V(s_new) - V(s)]
                try:
                    actor.learn(s, a, td_error)  # learn Policy : true_gradient = grad[logPi(s, a) * td_error]
                except KeyboardInterrupt:  # if Ctrl+C at running actor.learn(), then save model, or exit if not at actor.learn()
                    actor.save_ckpt()
                    critic.save_ckpt()
                    # logging

                s = s_new
                t += 1

                if done or t >= max_steps:
                    ep_rs_sum = sum(all_r)

                    if 'running_reward' not in globals():
                        running_reward = ep_rs_sum
                    else:
                        running_reward = running_reward * 0.95 + ep_rs_sum * 0.05

                    rewards.append(running_reward)
                    # start rending if running_reward greater than a threshold
                    # if running_reward > DISPLAY_REWARD_THRESHOLD: RENDER = True
                    # print("Episode: %d reward: %f running_reward %f took: %.5f" % \
                    #     (i_episode, ep_rs_sum, running_reward, time.time() - episode_time))
                    print('Episode: {}/{}  | Episode Reward: {:.4f}  | Running Time: {:.4f}'\
                    .format(i_episode, train_episodes, ep_rs_sum, time.time()-t0 ))

                    # Early Stopping for quick check
                    if t >= max_steps:
                        print("Early Stopping")
                        s = env.reset().astype(np.float32)
                        rall = 0
                        while True:
                            env.render()
                            # a = actor.choose_action(s)
                            a = actor.choose_action_greedy(s)  # Hao Dong: it is important for this task
                            s_new, r, done, info = env.step(a)
                            s_new = np.concatenate((s_new[0:N_F], s[N_F:]), axis=0).astype(np.float32)
                            rall += r
                            s = s_new
                            if done:
                                # print("reward", rall)
                                s = env.reset().astype(np.float32)
                                rall = 0
                    break

            if i_episode % save_interval == 0:
                actor.save_ckpt()
                critic.save_ckpt()
                plot(rewards, Algorithm_name='AC', Env_name=env_id)
        actor.save_ckpt()
        critic.save_ckpt()

    if mode == 'test':
        actor.load_ckpt()
        critic.load_ckpt()
        t0 = time.time()

        for i_episode in range(test_episodes):
            episode_time = time.time()
            s = env.reset().astype(np.float32)
            t = 0  # number of step in this episode
            all_r = []  # rewards of all steps
            while True:
                if render: env.render()
                a = actor.choose_action(s)
                s_new, r, done, info = env.step(a)
                s_new = s_new.astype(np.float32)
                if done: r = -20

                all_r.append(r)
                s = s_new
                t += 1

                if done or t >= max_steps:
                    ep_rs_sum = sum(all_r)

                    if 'running_reward' not in globals():
                        running_reward = ep_rs_sum
                    else:
                        running_reward = running_reward * 0.95 + ep_rs_sum * 0.05
                    # start rending if running_reward greater than a threshold
                    # if running_reward > DISPLAY_REWARD_THRESHOLD: RENDER = True
                    # print("Episode: %d reward: %f running_reward %f took: %.5f" % \
                    #     (i_episode, ep_rs_sum, running_reward, time.time() - episode_time))
                    print('Episode: {}/{}  | Episode Reward: {:.4f}  | Running Time: {:.4f}'\
                    .format(i_episode, test_episodes, ep_rs_sum, time.time()-t0 ))

                    # Early Stopping for quick check
                    if t >= max_steps:
                        print("Early Stopping")
                        s = env.reset().astype(np.float32)
                        rall = 0
                        while True:
                            env.render()
                            # a = actor.choose_action(s)
                            a = actor.choose_action_greedy(s)  # Hao Dong: it is important for this task
                            s_new, r, done, info = env.step(a)
                            s_new = np.concatenate((s_new[0:N_F], s[N_F:]), axis=0).astype(np.float32)
                            rall += r
                            s = s_new
                            if done:
                                # print("reward", rall)
                                s = env.reset().astype(np.float32)
                                rall = 0
                    break

max_stars_repo_name: Helilysyt/tensorlayer

max_issues_repo_name: Helilysyt/tensorlayer

max_forks_repo_name: Helilysyt/tensorlayer


text:
#Konstantin Ntokas
#
#
# 2019-11-17
# -----------------------------------------------------------------------------
# This function calcultes the input data for the reliability diagramm
#
# input:
#           calculation:    mxn matrix; m = number of simulations
#                                       n = number of member in ensemble
#           observation:    mx1 vector; m = number of records
#
# output:
#          nomi:            nominal probability of the intervals
#          Meff:            mean effective probability of the intervals
#          Mlen:            mean effective length of the intervals
# -----------------------------------------------------------------------------
import numpy as np

def Reliability(calculation, observation):
    # transform input into numpy array
    calculation = np.array(calculation, dtype='float64')
    observation = np.array(observation, dtype='float64')
    dim1 = calculation.shape
    if len(dim1) == 1:
        calculation = calculation.reshape((1,dim1[0]))
    dim2 = observation.shape
    if len(dim2) == 0:
        observation = observation.reshape((1,1))
    elif len(dim2) == 1:
        observation = observation.reshape((dim2[0],1))

    # Nominal probability of the intervals
    bin_start = np.arange(0.05, 0.5, 0.05)
    bin_end = np.arange(0.95, 0.5, -0.05)
    bins = np.concatenate((bin_start, bin_end))
    nb = len(bin_start)
    nomi = bin_end - bin_start

    # initialisation
    L = np.size(calculation, axis=0)
    length = np.zeros((L,nb))
    eff = np.zeros((L,nb))

    for i in range(L):
        if ~np.isnan(observation[i]):
            # get quantile for each bin and the median
            q = np.quantile(calculation[i,:], bins)
            qmed = np.median(calculation[i,:])

            # Compute lenghts of intervals (for Luc Perreault)
            length[i,:] = q[nb:] - q[:nb]

            # Locate observation in the ensemble
            if observation[i] <= qmed:
                eff[i,:] = q[:nb] <= observation[i]
            else:
                eff[i,:] = q[nb:] >= observation[i]
        else:
            length[i, :] = np.nan
            eff[i, :] = np.nan

    # Compute averages
    Meff = np.nanmean(eff, axis=0)
    Mlen = np.nanmean(length, axis=0)  # (For Luc)

    return nomi, Meff, Mlen


text:
from Stream import Stream, TimeAndValue
from Operators import stream_agent, tf, stream_func, ef
from examples_element_wrapper import print_stream
import numpy as np
from random import randint
from Stream import _no_value
"""Simple example of the window function wf with
window size and step size of 3. Output is the mean
of the clipped window with a_min=9, a_max=11. These
values are used in this example so that you can
check the results easily.

"""
class tv(object):
    def __init__(self, time, value):
        self.time = time
        self.value = value

def print_stream_tv(stream):
    def print_element(v, count):
        print '{0}[{1}] time = {2}, value = {3}'.format(stream.name, count, v.time, v.value)
        return (count+1)

    ef(inputs=stream, outputs=None, func=print_element, state=0)

x = Stream('x')
y = Stream('y')
z = Stream('z')
print_stream(x)
print_stream(y)
print_stream(z)

def h(list_of_timed_windows, threshold):
    return sum([sum([w[1] for w in timed_window])
             for timed_window in list_of_timed_windows])

tf([x,y], z, h, 10, 10, threshold=20)

x.extend([(2, 10), (6, 11), (12, 12), (14, 13), (19, 12), (22, 13), (60, 23)])
y.extend([(3, 11), (4, 15), (14, 12), (16, 13), (21, 12), (22, 13), (50, 23)])

max_stars_repo_name: AnomalyInc/StreamPy

max_issues_repo_name: StreamPy/StreamPy

max_forks_repo_name: StreamPy/StreamPy


text:
import numpy
from Util import Util

class YuwanaSeborg:
    def __init__(self, data):
        self.data = data
        self.estimative = numpy.zeros(len(self.data))

        k_c = 1.0
        setpoint = 1.0

        [y_p1, y_m1, y_p2] = Util.findCriticalPoints(self.data, 3)

        print("y_p1 = " + str(y_p1))
        print("y_m1 = " + str(y_m1))
        print("y_p2 = " + str(y_p2))

        # Finding y when steady-state
        y_ss = (y_p2 * y_p1 - (y_m1**2))/(y_p1 + y_p2 - 2*y_m1)

        # Computing zeta process which is important to find other process parameters
        aux1 = numpy.log((y_ss - y_m1)/(y_p1 - y_ss))

        aux2 = numpy.log((y_p2 - y_ss)/(y_p1 - y_ss))

        zeta_1 = -aux1 / numpy.sqrt(numpy.pi**2 + aux1**2)
        zeta_2 = -aux2 / numpy.sqrt(4 * (numpy.pi**2) + aux2**2)

        zeta_m = (zeta_1 + zeta_2)/2

        # Determining overshootDuration
        deltaT = self.__getOvershootDuration(data, y_ss)

        print("deltaT = " + str(deltaT))

        k_m = y_ss / (k_c * (setpoint - y_ss))
        k_f = k_c * k_m

        # Finding process parameters
        aux1 = numpy.sqrt((1 - zeta_m**2) * (k_f + 1))
        aux2 = zeta_m * numpy.sqrt(k_f + 1) + numpy.sqrt((zeta_m**2) * (k_f + 1) + k_f)

        theta_m = 2 * (deltaT / numpy.pi) * (aux1 / aux2)

        tau_m = (deltaT / numpy.pi) * aux1 * aux2

        print("k_m = " + str(k_m))
        print("tau_m = " + str(tau_m))

        # Estimating transfer function parameters

        self.k = k_f / (k_f + 1)

        self.tau = numpy.sqrt((theta_m * tau_m) / (2 * (k_f + 1)))

        self.zeta = (tau_m + 0.5 * theta_m * (1 - k_f)) / numpy.sqrt(2 * theta_m * tau_m * (k_f + 1))

        self.theta = theta_m

        self.__estimate()

    def __estimate(self):
        if self.zeta < 1:
            beta = numpy.sqrt(1 - self.zeta ** 2 )
            phi = numpy.arccos(self.zeta)
            gama = self.theta / (2 * self.tau)

            for i in range(len(self.data)):
                t = self.data[i, 0]
                alpha = (beta * t / self.tau) + phi

                self.estimative[i] = self.k * (1 \
                    + numpy.exp(-self.zeta * t / self.tau) \
                    * (gama * numpy.cos(alpha) \
                    - (1 + self.zeta * gama) * numpy.sin(alpha) / beta))

    def __getOvershootDuration(self, data, y_ss):
        t1 = data[0, 0]
        t2 = data[-1, 0]
        flag = True

        for value in data:
            if flag:
                if value[1] >= y_ss:
                    t1 = value[0]
                    flag = False
            else:
                if value[1] <= y_ss:
                    t2 = value[0]
                    break

        return (t2 - t1)

    def showResults(self):
        print("The transfer function")
        print("G(s) = " + str(self.k) + " * (1 - 0.5 * " + str(self.theta)+ " * s)")
        print("---------------------------")
        print(str(self.tau) + "**2 * s**2 + 2 * "+ str(self.zeta) +" * " + str(self.tau) + " * s + 1")

max_stars_repo_name: glauberrleite/system-identifier

max_issues_repo_name: glauberrleite/system-identifier

max_forks_repo_name: glauberrleite/system-identifier


text:
# Copyright (C) 2020. Huawei Technologies Co., Ltd. All rights reserved. THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

import numpy as np
import torch

# normalize to [0, 1]
class Normalize():
    def __init__(self, worker=None):
        pass

    def _is_numpy_image(self, img):
        return isinstance(img, np.ndarray) and (img.ndim in {2, 3})

    def __call__(self, input_dict):
        input = input_dict['rgb'].astype(np.float32)
        max_value = np.iinfo(input_dict['rgb'].dtype).max

        for i in range(input.shape[0]):
            pic = input[i, ...]
            if not self._is_numpy_image(pic):
                raise TypeError('pic should be ndarray. Got {}'.format(type(pic)))

            # handle numpy array
            if pic.ndim == 2:#shape属性长度，灰度图，处理为三维[h,w,c],同opencv
                pic = pic[:, :, None]

            # normalize
            pic = pic.astype(np.float32) / max_value

            input[i, ...] = pic

        input = input.transpose((0, 3, 1, 2))#[b,c,h,w]
        input_dict['rgb'] = input

        return input_dict


text:
import krpc
import numpy as np
#import time
import Global

conn = Global.conn
vessel = conn.space_center.active_vessel

def norm(x):
    return np.linalg.norm(x)

def r2d(x):
    return x*180.0/np.pi

def d2r(x):
    return x*np.pi/180.0

#function to calculate distance between 2 (lat,lon) coords
def db2(lat1, lon1, lat2, lon2):
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = np.sin(0.5*dlat)**2.0 + np.cos(lat1)*np.cos(lat2)*(np.sin(0.5*dlon)**2.0)
    c = 2.0*np.arcsin(min(1.0, np.sqrt(a)))
    d = vessel.orbit.body.equatorial_radius*c
    return d

#function to calculate (lat, long) from state vectors
def latlon(x):
    lon = np.arctan2(x[1],x[0])
    p = (x[0]**2.0 + x[1]** 2.0)**0.5
    lat = np.arctan2(x[2],p)
    lat = r2d(lat)
    lon = r2d(lon)
    return np.array((lat, lon))

def Toggle(init_x, init_v):
    bmu = vessel.orbit.body.gravitational_parameter
    br = vessel.orbit.body.equatorial_radius

    r = init_x
    v = init_v

    r = [r[0], r[2], r[1]]

    v = [v[0], v[2], v[1]]

    #calculate orbital elements from state vectors
    h = np.cross(r, v)
    e = np.subtract(np.divide(np.cross(v, h), bmu), np.divide(r, norm(r)))
    ec = norm(e)

    n = np.transpose([-h[1], h[0], 0])

    if np.dot(r, v)>=0:
        nu = np.arccos(np.dot(e, r)/(ec*norm(r)))
    else:
        nu = 2*np.pi - np.arccos(np.dot(e, r)/(ec*norm(r)))

    i = np.arccos(h[2]/norm(h))
    if n[1]>=0:
        RAAN = np.arccos(n[0]/norm(n))
    else:
        RAAN = 2*np.pi - np.arccos(n[0]/norm(n))

    if e[2]>=0:
        w = np.arccos(np.dot(n, e)/(ec*norm(n)))
    else:
        w = 2*np.pi - np.arccos(np.dot(n, e)/(ec*norm(n)))

    a = (2.0/norm(r) - (norm(v)**2.0)/bmu)**-1

    elem = [a, ec, i, w, RAAN, nu]

    #calculate new true anomaly when distance to centre of earth is 6371km on the orbit
    elem[5] = -np.arccos(((1.0-elem[1]**2.0)*(elem[0]/br)-1.0)/elem[1])

    #convert back to state vectors
    X = np.cos(elem[4])*np.cos(elem[3]+elem[5]) - np.sin(elem[4])*np.sin(elem[3]+elem[5])*np.cos(elem[2])
    Y = np.sin(elem[4])*np.cos(elem[3]+elem[5]) + np.cos(elem[4])*np.sin(elem[3]+elem[5])*np.cos(elem[2])
    Z = np.sin(elem[2])*np.sin(elem[3]+elem[5])

    X *= br
    Y *= br
    Z *= br

    R = [X, Y, Z]

    return R


text:
import math

import numpy as np

from pymoo.model.crossover import Crossover
from pymoo.operators.crossover.util import crossover_mask
from pymoo.rand import random

class HalfUniformCrossover(Crossover):

    def __init__(self, prob=0.5):
        super().__init__(2, 2)
        self.prob = prob

    def _do(self, problem, pop, parents, **kwargs):

        # get the X of parents and count the matings
        X = pop.get("X")[parents.T]
        _, n_matings, n_var = X.shape

        # the mask do to the crossover
        M = np.full((n_matings, n_var), False)

        not_equal = X[0] != X[1]

        # create for each individual the crossover range
        for i in range(n_matings):
            I = np.where(not_equal[i])[0]

            n = math.ceil(len(I) / 2)
            if n > 0:
                _I = I[random.perm(len(I))[:n]]
                M[i, _I] = True

        _X = crossover_mask(X, M)
        return pop.new("X", _X)

max_stars_repo_name: temaurer/pymoo

max_issues_repo_name: temaurer/pymoo

max_forks_repo_name: temaurer/pymoo