text:
import numpy as np
import time
import pybullet as p
import pybullet_data
import os
import random
from gym import error, spaces, utils

from gym_pybullet_drones.envs.BaseAviary import DroneModel, Physics, BaseAviary

######################################################################################################################################################
#### Single tethered drone environment class for reinforcement learning applications (in this implementation, at some pose) ##########################
######################################################################################################################################################
class RLTetherAviary(BaseAviary):

    ####################################################################################################
    #### Initialize the environment ####################################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - drone_model (DroneModel)         desired drone type (associated to an .urdf file) ###########
    #### - num_drones (int)                 desired number of drones in the aviary #####################
    #### - visibility_radius (float)        used to compute the drones' adjacency matrix, in meters ####
    #### - initial_xyzs ((3,1) array)       initial XYZ position of the drones #########################
    #### - initial_rpys ((3,1) array)       initial orientations of the drones (radians) ###############
    #### - physics (Physics)                desired implementation of physics/dynamics #################
    #### - freq (int)                       the frequency (Hz) at which the physics engine advances ####
    #### - aggregate_phy_steps (int)        number of physics updates within one call of .step() #######
    #### - gui (bool)                       whether to use PyBullet's GUI ##############################
    #### - record (bool)                    whether to save a video of the simulation ##################
    #### - obstacles (bool)                 whether to add obstacles to the simulation #################
    ####################################################################################################
    def __init__(self, drone_model: DroneModel=DroneModel.CERLAB, num_drones: int=1,
                        neighbourhood_radius: float=np.inf, initial_xyzs=None, initial_rpys=None,
                        physics: Physics=Physics.PYB_TETHER, freq: int=200, aggregate_phy_steps: int=1,
                        gui=False, record=False, obstacles=False, user_debug_gui=True):
        if num_drones!=1: print("[ERROR] in RLTakeoffAviary.__init__(), RLTakeoffAviary only accepts num_drones=1" ); exit()
        super().__init__(drone_model=drone_model, neighbourhood_radius=neighbourhood_radius,
            initial_xyzs=initial_xyzs, initial_rpys=initial_rpys, physics=physics, freq=freq, aggregate_phy_steps=aggregate_phy_steps,
            gui=gui, record=record, obstacles=obstacles, user_debug_gui=user_debug_gui)

        self.N_ACTIONS = 5

        self.INIT_X_BOUND = 0.2
        self.INIT_Y_BOUND = 0.2
        self.INIT_Z_BOUND = 0.2
        self.INIT_R_BOUND = 0.1
        self.INIT_P_BOUND = 0.1
        self.INIT_YAW_BOUND = 0

        self.geoFenceMax_XY = 0.99
        self.geoFenceMax_Z  = 1.5
        self.geoFenceMin_Z  = 0

        self.penaltyPosition = 10
        self.penaltyAngle = 10
        self.penaltyVelocity = 50
        self.penaltyAngularVelocity = 50
        self.penaltyFlag = 1000
        self.penaltyTetherUsage = 0
        self.penaltyDiffTetherUsage = 5

        self.positionThreshold = 0.15
        self.angleThreshold = np.pi/18
        self.angularVelocityThreshold = 0.05
        self.velocityThreshold = 0.1

        self.rewardGoal = 20

        self.TETHER_MIN_LENGTH = 0.20
        self.MAX_TETHER_FORCE = 0.25*(4*(self.MAX_RPM**2 * self.KF))
        print("[INFO] Max tether force set to ", round(self.MAX_TETHER_FORCE,4), "N")

    def _housekeeping(self):
        #### Initialize/reset counters and zero-valued variables ###########################################
        self.RESET_TIME = time.time(); self.step_counter = 0; self.first_render_call = True
        self.X_AX = -1*np.ones(self.NUM_DRONES); self.Y_AX = -1*np.ones(self.NUM_DRONES); self.Z_AX = -1*np.ones(self.NUM_DRONES)
        if self.PHYSICS == Physics.PYB_TETHER:
            self.TETHER_AX = -1*np.ones(self.NUM_DRONES)
        self.GUI_INPUT_TEXT = -1*np.ones(self.NUM_DRONES); self.USE_GUI_RPM=False; self.last_input_switch = 0
        self.last_action = -1*np.ones((self.NUM_DRONES,self.N_ACTIONS))
        self.last_clipped_action = np.zeros((self.NUM_DRONES,self.N_ACTIONS)); self.gui_input = np.zeros(self.N_ACTIONS)
        self.no_pybullet_dyn_accs = np.zeros((self.NUM_DRONES,3))

        #### Initialize the drones kinematic information ##################################################
        unitVector = np.random.rand(3)-np.array([0.5,0.5,0])
        startingPoint = np.ones((self.NUM_DRONES,3))*np.array([0,0,0.5])
        self.pos = 0.5*np.random.rand(1)*unitVector/np.linalg.norm(unitVector)*np.ones((self.NUM_DRONES,3))+startingPoint

        self.quat = np.zeros((self.NUM_DRONES,4)); self.rpy = np.zeros((self.NUM_DRONES,3))
        self.vel = np.zeros((self.NUM_DRONES,3)); self.ang_v = np.zeros((self.NUM_DRONES,3))

        #### Set PyBullet's parameters #####################################################################
        p.setGravity(0, 0, -self.G, physicsClientId=self.CLIENT)
        p.setRealTimeSimulation(0, physicsClientId=self.CLIENT)
        p.setTimeStep(self.TIMESTEP, physicsClientId=self.CLIENT)
        p.setAdditionalSearchPath(pybullet_data.getDataPath(), physicsClientId=self.CLIENT)
        #### Load ground plane, drone and obstacles models #################################################
        self.PLANE_ID = p.loadURDF("plane.urdf", physicsClientId=self.CLIENT)
        self.DRONE_IDS = np.array([p.loadURDF(os.path.dirname(os.path.abspath(__file__))+"/../assets/"+self.URDF, self._setInitialXYZ(self.INIT_XYZS[i,:]), p.getQuaternionFromEuler(self._setInitialRPY(self.INIT_RPYS[i,:])), physicsClientId=self.CLIENT) for i in range(self.NUM_DRONES)])
        for i in range(self.NUM_DRONES):
            #### Show the frame of reference of the drone, note thet it can severly slow down the GUI ##########
            if self.GUI and self.USER_DEBUG: self._showDroneLocalAxes(i)
            #### Disable collisions between drones' and the ground plane, e.g., to start a drone at [0,0,0] ####
            p.setCollisionFilterPair(bodyUniqueIdA=self.PLANE_ID, bodyUniqueIdB=self.DRONE_IDS[i], linkIndexA=-1, linkIndexB=-1, enableCollision=0, physicsClientId=self.CLIENT)
        if self.OBSTACLES: self._addObstacles()

    ####################################################################################################
    #### Return the action space of the environment, a Box(4,) #########################################
    ####################################################################################################
    #### Additional component in action space for tether force magnitude ###############################
    ####################################################################################################
    def _actionSpace(self):
        #### Action vector ######## P0            P1            P2            P3            P4
        act_lower_bound = np.array([-1,           -1,           -1,           -1,           0])
        act_upper_bound = np.array([ 1,            1,            1,            1,           1])
        return spaces.Box( low=act_lower_bound, high=act_upper_bound, dtype=np.float32 )

    ####################################################################################################
    #### Return the observation space of the environment, a Box(20,) ###################################
    ####################################################################################################
    def _observationSpace(self):
        #### Observation vector ### X        Y        Z      R       P       Y       VX       VY       VZ       WR       WP       WY       P0            P1            P2            P3
        obs_lower_bound = np.array([-1,      -1,      0,     -1,     -1,     -1,     -1,      -1,      -1,      -1,      -1,      -1,      -1,           -1,           -1,           -1])
        obs_upper_bound = np.array([ 1,       1,      1,      1,      1,      1,      1,       1,       1,       1,       1,       1,       1,            1,            1,            1])
        return spaces.Box( low=obs_lower_bound, high=obs_upper_bound, dtype=np.float32 )

    ####################################################################################################
    #### Return the current observation of the environment #############################################
    ####################################################################################################
    #### Returns #######################################################################################
    #### - obs (20,) array                  for its content see _observationSpace() ####################
    ####################################################################################################
    def _computeObs(self):
        droneState = self._getDroneStateVector(0)
        droneState = np.hstack([droneState[0:3],droneState[7:20]])
        droneState.reshape(16,)
        return self._clipAndNormalizeState(droneState)

    ####################################################################################################
    #### Preprocess the action passed to step() ########################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - action ((N_ACTIONS,1) array)     unclipped RPMs commanded to the 4 motors of the drone ######
    ####################################################################################################
    #### Returns #######################################################################################
    #### - clipped_action ((N_ACTIONS,1) array)  clipped RPMs commanded to the 4 motors of the drone ###
    ####################################################################################################
    def _preprocessAction(self, action):
        rpm, mag = self._normalizedActionToRPM(action)
        return np.hstack((np.clip(np.array(rpm), 0, self.MAX_RPM), np.clip(np.array(mag), 0, self.MAX_TETHER_FORCE)))

    ####################################################################################################
    #### Compute the current reward value(s) ###########################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - obs (..)                         the return of _computeObs() ################################
    ####################################################################################################
    #### Returns #######################################################################################
    #### - reward (..)                      the reward(s) associated to the current obs/state ##########
    ####################################################################################################
    def _computeReward(self, obs):
        # x = obs[0]
        # y = obs[1]
        # z = obs[2]

        # dx = obs[3]
        # dy = obs[4]
        # dz = obs[5]

        # roll = obs[6]
        # pitch = obs[7]
        # yaw = obs[8]

        # droll = obs[9]
        # dpitch = obs[10]
        # dyaw = obs[11]

        xy = obs[0:2]
        rp = obs[3:5]
        v = obs[6:9]
        # angle = obs[3:6]
        omega = obs[9:12]
        # actions = obs[12:16]
        # tetherUsage = self.tether_force_mag
        deltaTetherUsage = np.abs(self.tether_force_mag - self.last_tether_force_mag)

        errorPosition = np.sum(np.square(xy))                                       # error in x and y
        errorAngle = np.sum(np.square(rp))                                          # error in r and p
        errorVelocity = np.sum(np.square(v))                                        # error in dx, dy, dz
        errorAngularVelocity = np.sum(np.square(omega[0:3]))                        # error in droll, dpitch, dyaw

        penaltyPosition = errorPosition*self.penaltyPosition                        # get cost from position        [error * penalty]
        penaltyAngle = errorAngle*self.penaltyAngle                                 # get cost from angle           [error^2 * penalty]
        penaltyVelocity = errorVelocity*self.penaltyVelocity                        # get cost from velovity        [error * penalty]
        penaltyAngularVelocity = errorAngularVelocity*self.penaltyAngularVelocity   # get cost from ang vel         [error * penalty]
        # penaltyTetherUsage = tetherUsage * self.penaltyTetherUsage                  # get cost from tether force    [usage * penalty]
        penaltyDiffTetherUsage = deltaTetherUsage * self.penaltyDiffTetherUsage     # get cost from dtether force   [change in usage * penalty]

        # check if out of geofence
        outOfGeoFence = np.linalg.norm(obs[0:2]) > self.geoFenceMax_XY or obs[2] > self.geoFenceMax_Z or obs[2] < self.geoFenceMin_Z

        crashed = True if obs[2]<self.COLLISION_H else False                        # check if drone is crashed
        penaltyFlag = self.penaltyFlag if outOfGeoFence or crashed else 0           # set penaltyFlag if geoFence or crashed

        rewardPosition = np.sqrt(errorPosition)<self.positionThreshold                          # determine if position is within threshold
        rewardVelocity = np.sqrt(errorVelocity)<self.velocityThreshold                          # determine if lin vel is within threshold
        rewardAngularVelocity = np.sqrt(errorAngularVelocity)<self.angularVelocityThreshold     # determine if ang vel is within threshold

        if all([rewardPosition, rewardVelocity, rewardAngularVelocity]):            # if all conditions are met, add goal reward, else don't
            rewardGoal = self.rewardGoal
        else:
            rewardGoal = 0

        #actions = np.abs(actions)
        #penaltyControl = np.sum(np.square(actions-self.HOVER_RPM))*self.penaltyControl
        #penaltyDiffControl = np.square(np.max(actions)-np.min(actions))*self.penaltyDiffControl

        # sum and return reward
        return rewardGoal - penaltyPosition - penaltyAngle - penaltyVelocity - penaltyAngularVelocity \
               - penaltyFlag + penaltyDiffTetherUsage# - penaltyTetherUsage - penaltyControl - penaltyDiffControl

    ####################################################################################################
    #### Compute the current done value(s) #############################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - obs (..)                         the return of _computeObs() ################################
    ####################################################################################################
    #### Returns #######################################################################################
    #### - done (..)                        the done value(s) associated to the current obs/state ######
    ####################################################################################################
    def _computeDone(self, obs):
        outOfGeoFence = np.linalg.norm(obs[0:2]) > self.geoFenceMax_XY or \
                        obs[2] > self.geoFenceMax_Z or \
                        obs[2] < self.geoFenceMin_Z
        outOfTime = True if (self.step_counter/self.SIM_FREQ > 10) else False
        crashed = True if obs[2]<self.COLLISION_H else False
        return outOfGeoFence or outOfTime or crashed

    ####################################################################################################
    #### Compute the current info dict(s) ##############################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - obs (..)                         the return of _computeObs() ################################
    ####################################################################################################
    #### Returns #######################################################################################
    #### - info (..)                        the info dict(s) associated to the current obs/state #######
    ####################################################################################################
    def _computeInfo(self, obs):
        return {"answer": 42} #### Calculated by the Deep Thought supercomputer in 7.5M years

    ####################################################################################################
    #### Normalize the 20 values in the simulation state to the [-1,1] range ###########################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - state ((20,1) array)             raw simulation state #######################################
    ####################################################################################################
    #### Returns #######################################################################################
    #### - normalized state ((20,1) array)  clipped and normalized simulation state ####################
    ####################################################################################################
    def _clipAndNormalizeState(self, state):
        clipped_pos = np.clip(state[0:3], -1, 1)
        clipped_rp = np.clip(state[3:5], -np.pi/3, np.pi/3)
        clipped_vel = np.clip(state[6:9], -5, 5)
        clipped_ang_vel_rp = np.clip(state[9:11], -10*np.pi, 10*np.pi)
        clipped_ang_vel_y = np.clip(state[11], -20*np.pi, 20*np.pi)
        # if self.GUI: self._clipAndNormalizeStateWarning(state, clipped_pos, clipped_rp, clipped_vel, clipped_ang_vel_rp, clipped_ang_vel_y)
        normalized_pos = clipped_pos
        normalized_rp = clipped_rp/(np.pi/3)
        normalized_y = state[5]/np.pi
        normalized_vel = clipped_vel/5
        normalized_ang_vel_rp = clipped_ang_vel_rp/(10*np.pi)
        normalized_ang_vel_y = clipped_ang_vel_y/(20*np.pi)
        return np.hstack([normalized_pos, normalized_rp, normalized_y, normalized_vel, normalized_ang_vel_rp, normalized_ang_vel_y, state[12:16] ]).reshape(16,)

    ####################################################################################################
    #### Print a warning if any of the 20 values in a state vector is out of the normalization range ###
    ####################################################################################################
    def _clipAndNormalizeStateWarning(self, state, clipped_pos, clipped_rp, clipped_vel, clipped_ang_vel_rp, clipped_ang_vel_y):
        if not(clipped_pos==np.array(state[0:3])).all(): print("[WARNING] it", self.step_counter, "in RLTakeoffAviary._clipAndNormalizeState(), out-of-bound position [{:.2f} {:.2f} {:.2f}], consider a more conservative implementation of RLTakeoffAviary._computeDone()".format(state[0], state[1], state[2]))
        if not(clipped_rp==np.array(state[3:5])).all(): print("[WARNING] it", self.step_counter, "in RLTakeoffAviary._clipAndNormalizeState(), out-of-bound roll/pitch [{:.2f} {:.2f}], consider a more conservative implementation of RLTakeoffAviary._computeDone()".format(state[3], state[4]))
        if not(clipped_vel==np.array(state[6:9])).all(): print("[WARNING] it", self.step_counter, "in RLTakeoffAviary._clipAndNormalizeState(), out-of-bound velocity [{:.2f} {:.2f} {:.2f}], consider a more conservative implementation of RLTakeoffAviary._computeDone()".format(state[6], state[7], state[8]))
        if not(clipped_ang_vel_rp==np.array(state[9:11])).all(): print("[WARNING] it", self.step_counter, "in RLTakeoffAviary._clipAndNormalizeState(), out-of-bound angular velocity [{:.2f} {:.2f} {:.2f}], consider a more conservative implementation of RLTakeoffAviary._computeDone()".format(state[9], state[10], state[11]))
        if not(clipped_ang_vel_y==np.array(state[11])): print("[WARNING] it", self.step_counter, "in RLTakeoffAviary._clipAndNormalizeState(), out-of-bound angular velocity [{:.2f} {:.2f} {:.2f}], consider a more conservative implementation of RLTakeoffAviary._computeDone()".format(state[9], state[10], state[11]))

    ####################################################################################################
    #### Denormalize the [-1,1] range to the [0, MAX RPM] & [0, MAX_TETHER_FORCE] range ################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - action ((N_ACTIONS,1) array)     normalized [-1,1] actions applied to the 4 motors ##########
    ####################################################################################################
    #### Returns #######################################################################################
    #### - rpm ((N_ACTIONS,1) array)        RPM values to apply to the 4 motors ########################
    #### - tether_force_mag                 force applied by tether ####################################
    ####################################################################################################
    def _normalizedActionToRPM(self, action):
        if np.any(np.abs(action[0:self.N_ACTIONS-1]))>1:
            print("\n[ERROR] it", self.step_counter, "in BaseAviary._normalizedActionToRPM(), out-of-bound action")
        rpm = np.where(action[0:self.N_ACTIONS-1] <= 0, (action[0:self.N_ACTIONS-1]+1)*self.HOVER_RPM, action[0:self.N_ACTIONS-1]*self.MAX_RPM) # Non-linear mapping: -1 -> 0, 0 -> HOVER_RPM, 1 -> MAX_RPM
        tether = action[self.N_ACTIONS-1]*self.MAX_TETHER_FORCE # Non-linear mapping: 0 -> 0, 1 -> MAX_TETHER_FORCE
        return rpm, tether

    ####################################################################################################
    #### Return the state vector of the nth drone ######################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - nth_drone (int)                  order position of the drone in list self.DRONE_IDS #########
    ####################################################################################################
    #### Returns #######################################################################################
    #### - state ((21,) array)              the state vector of the nth drone ##########################
    ####################################################################################################
    def _getDroneStateVector(self, nth_drone):
        state = np.hstack([self.pos[nth_drone,:], self.quat[nth_drone,:], self.rpy[nth_drone,:], self.vel[nth_drone,:], self.ang_v[nth_drone,:], self.last_action[nth_drone,:]])
        return state.reshape(21,)

    ####################################################################################################
    #### Set initial xyz for current episode ###########################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - initial_xyzs ((3,1) array)     initial x, y, and z position of drone in episode #############
    ####################################################################################################
    def _setInitialXYZ(self, initial_xyzs):
        return np.vstack([
                            random.uniform(-self.INIT_X_BOUND, self.INIT_X_BOUND),
                            random.uniform(-self.INIT_Y_BOUND, self.INIT_Y_BOUND),
                            random.uniform(-self.INIT_Z_BOUND, self.INIT_Z_BOUND)+1.0
                                                                ])

    ####################################################################################################
    #### Set initial rpy for current episode ###########################################################
    ####################################################################################################
    #### Arguments #####################################################################################
    #### - initial_rpys ((3,1) array)     initial r, p, and y orientations of drone in episode #########
    ####################################################################################################
    def _setInitialRPY(self, initial_rpys):
        return np.vstack([
                            random.uniform(-self.INIT_R_BOUND, self.INIT_R_BOUND),
                            random.uniform(-self.INIT_P_BOUND, self.INIT_P_BOUND),
                            random.uniform(-self.INIT_YAW_BOUND, self.INIT_YAW_BOUND)
                                                                ])


text:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import scipy

train_data = pd.read_csv('../input/train.csv')
mean_train_data = train_data.groupby('place_id').mean()
std_train_data = train_data.groupby('place_id').std()

acc_df = pd.concat([mean_train_data['accuracy'],std_train_data['x'],std_train_data['y']], axis=1)
acc_df.rename(columns={'accuracy':'mean_accuracy','x':'std_x','y':'std_y'}, inplace=True)

acc_df.fillna(0, inplace=True)

p = plt.hist(acc_df.mean_accuracy, bins=np.arange(min(acc_df.mean_accuracy), max(acc_df.mean_accuracy) + 1, 1))
plt.xlabel('Mean Accuracy')
plt.ylabel('Count')
plt.title('Counts of mean accuracy for places')

plt.show()

p = plt.hist(acc_df.std_x, bins=np.arange(min(acc_df.std_x), max(acc_df.std_x) + .01, .01))
plt.xlabel('Std(y)')
plt.ylabel('Count')
plt.title('Counts of std(x) for places')

plt.show()

p = plt.hist(acc_df.std_y, bins=np.arange(min(acc_df.std_y), max(acc_df.std_y) + .001, .001))
plt.xlabel('Std(y)')
plt.ylabel('Count')
plt.title('Counts of std(y) for places')

plt.show()

p = plt.scatter(acc_df.mean_accuracy,acc_df.std_x)
plt.xlabel('Mean Accuracy')
plt.ylabel('Std(x)')
plt.title('std(x) vs Mean accuracy for each place')

plt.show()

p = plt.scatter(acc_df.mean_accuracy, acc_df.std_y)
plt.xlabel('Mean Accuracy')
plt.ylabel('Std(y)')
plt.title('std(y) vs mean accuracy for each place')

plt.show()

p = plt.scatter(acc_df.mean_accuracy,acc_df.std_y)
l = plt.plot([i for i in range(0,1000)], [func(i, -0.01277189,  1.4337799 ,  2.71943575) for i in range(0,1000)], 'r', linewidth=1)
plt.xlabel('Mean Accuracy')
plt.ylabel('Std(y)')
plt.title('std(y) vs Mean accuracy for each place (exp decay)')

plt.show()
popt, pcov = scipy.optimize.curve_fit(func, acc_df.mean_accuracy, acc_df.std_y)
popt, pcov = scipy.optimize.curve_fit(func4, acc_df.mean_accuracy, acc_df.std_y)
p = plt.scatter(acc_df.mean_accuracy,acc_df.std_y)
l = plt.plot([i for i in range(0,1000)], [func(i, 221.24964964,   12.02131292,  804.60004683) for i in range(0,1000)], 'r', linewidth=1)
plt.xlabel('Mean Accuracy')
plt.ylabel('Std(y)')
plt.title('std(y) vs Mean accuracy for each place (1/x)')

plt.show()

ext: py

max_stars_repo_name: khushhallchandra/Deep-Learning

max_issues_repo_name: khushhallchandra/Deep-Learning

max_forks_repo_name: khushhallchandra/Deep-Learning

reason: import numpy,import scipy


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fftpack import fft, fft2, ifft
import scipy.io as sio
import copy

## general simulation parameters

fs = 1024
npnts = fs*5 # 5 seconds# parameters
freq = 4 # peak frequency

# centered time vector
timevec = np.arange(0,npnts)/fs
timevec = timevec - np.mean(timevec)

# for power spectrum
hz = np.linspace(0,fs/2,int(np.floor(npnts/2)+1))

## Mexican hat wavelet

# the wavelet
s = .4
MexicanWavelet = (2/(np.sqrt(3*s)*np.pi**.25)) * (1- (timevec**2)/(s**2) ) * np.exp( (-timevec**2)/(2*s**2) )

# amplitude spectrum
MexicanPow = np.abs(fft(MexicanWavelet)/npnts)

# time-domain plotting
plt.subplot(211)
plt.plot(timevec,MexicanWavelet,'k')
plt.xlabel('Time (sec.)')
plt.title('Mexican wavelet in time domain')

# frequency-domain plotting
plt.subplot(212)
plt.plot(hz,MexicanPow[:len(hz)],'k')
plt.xlim([0,freq*3])
plt.xlabel('Frequency (Hz)')
plt.title('Mexican wavelet in frequency domain')
plt.show()

reason: import numpy,import scipy,from scipy


text:
import numpy as np
from math import gcd
# def calculate(l,g,s,r,n1,n2,n3):
#     if g[2]:
#         if l > g[1]:
#             temp_g = [g[0],g[1],g[2]-1]
#             return calculate(l-g[1],temp_g,s,r,n1+1,n2,n3)
#         else:
#             m = n1*g[0] + n2*s[0] + n3*r[0]
#             l2,g2,s2,r2,n1_2,n2_2,n3_2 =
#     return

def treasures(info, limit):
    l = int(limit*1000)
    g = [info['golden coin']['price'], info['golden coin']['weight'], info['golden coin']['amount']]
    s = [info['silver coin']['price'], info['silver coin']['weight'], info['silver coin']['amount']]
    r = [info['ruby']['price'], info['ruby']['weight'], info['ruby']['amount']]

    d = l
    for i in (g,s,r):
        d = gcd(d,i[1])
    l = int(l/d)
    g[1] = int(g[1]/d)
    s[1] = int(s[1]/d)
    r[1] = int(r[1]/d)

    result = np.zeros((g[2]+s[2]+r[2],l+1,4),dtype=int)
    for i in range(l+1):
        if i >= g[1]:
            result[(0,i,0)] = g[0]
            result[(0,i,1)] = 1

    for i in range(1,g[2]):
        for j in range(1,l+1):
            if j < g[1]:
                result[(i,j)] = result[(i-1,j)]
            elif result[(i-1,j,0)] < result[(i-1,j-g[1],0)]+g[0]:
                result[(i,j)] = result[(i-1,j-g[1])]
                result[(i,j,0)] += g[0]
                result[(i,j,1)] += 1
            else:
                result[(i,j)] = result[(i-1,j)]

    for i in range(g[2],g[2]+s[2]):
        for j in range(1,l+1):
            if j < s[1]:
                result[(i,j)] = result[(i-1,j)]
            elif result[(i-1,j,0)] < result[(i-1,j-s[1],0)]+s[0]:
                result[(i,j)] = result[(i-1,j-s[1])]
                result[(i,j,0)] += s[0]
                result[(i,j,2)] += 1
            else:
                result[(i,j)] = result[(i-1,j)]

    for i in range(g[2]+s[2],g[2]+s[2]+r[2]):
        for j in range(1,l+1):
            if j < r[1]:
                result[(i,j)] = result[(i-1,j)]
            elif result[(i-1,j,0)] < result[(i-1,j-r[1],0)]+r[0]:
                result[(i,j)] = result[(i-1,j-r[1])]
                result[(i,j,0)] += r[0]
                result[(i,j,3)] += 1
            else:
                result[(i,j)] = result[(i-1,j)]

    n1 = result[(-1,-1,1)]
    n2 = result[(-1,-1,2)]
    n3 = result[(-1,-1,3)]

    result_list = []
    if n1:
        result_list.append('golden coin: '+str(n1))
    if n2:
        result_list.append('silver coin: '+str(n2))
    if n3:
        result_list.append('ruby: '+str(n3))

    return result_list

if __name__ == '__main__':
    print("Example:")
    print(treasures({'golden coin':
                        {'price': 100, 'weight': 50, 'amount': 200},
                     'silver coin':
                        {'price': 10, 'weight': 20, 'amount': 1000},
                     'ruby':
                        {'price': 1000, 'weight': 200, 'amount': 2}}, 5))

    #These "asserts" using only for self-checking and not necessary for auto-testing
    assert treasures({'golden coin':
                         {'price': 100, 'weight': 50, 'amount': 200},
                      'silver coin':
                         {'price': 10, 'weight': 20, 'amount': 1000},
                      'ruby':
                         {'price': 1000, 'weight': 200, 'amount': 2}}, 5) == [
                          'golden coin: 92', 'ruby: 2']
    assert treasures({'golden coin':
                         {'price': 100, 'weight': 50, 'amount': 100},
                      'silver coin':
                         {'price': 10, 'weight': 20, 'amount': 100},
                      'ruby':
                         {'price': 1000, 'weight': 200, 'amount': 1}}, 7.5) == [
                          'golden coin: 100', 'silver coin: 100', 'ruby: 1']
    print("Coding complete? Click 'Check' to earn cool rewards!")

reason: import numpy


text:
import numpy as np

def read_motion(key_frames):
    key_motion = open("key_data_collection/Talk_05_Key.bvh", "a")
    key_motion.write("Frames: "+str(len(key_frames))+"\n")
    key_motion.write("Frame Time: 0.016667\n")

    whole_motion = open("key_data_collection/NaturalTalking_005.bvh", "r")
    check = False
    for i, line in enumerate(whole_motion):
        if line == "MOTION\n":
            check = True
            key_frames += i
            key_frames += 2
        if check:
            if i in key_frames:
                key_motion.write(line)

if __name__ == "__main__":
    motion_001 = np.array([1, 48, 98, 187, 262, 290, 340, 690, 700, 710, 717, 786, 847, 1044, 1250, 1469, 1600, 1653, 1735, 1916, 2364, 2419, 2605, 2619, 2635, 2676, 2967, 3578, 3591, 3616, 3764, 4097, 4210, 5047, 5179, 5311, 5320, 5346, 5396, 5770, 5824, 6086, 6135, 6157, 6171, 6746, 7315, 7439, 7812, 8911, 8923, 9267, 9417, 9460, 9490, 9623, 9638, 9663, 9695, 9828, 9922, 9940, 9971, 10120, 10262, 10316, 10914, 10931, 10938, 10943, 10948, 10956, 11392, 11398, 11492, 11499, 11504, 11512, 11531, 11550, 11605, 11630, 11650, 11742, 12002])
    motion_002 = np.array([1, 24, 27, 35, 115, 200, 239, 334, 550, 594, 646, 740, 755, 1289, 1359, 1575, 1957, 2960, 3199, 3246, 6827, 7018, 9500, 9518, 9534, 9560, 9591, 9640, 9699, 11020, 11122, 11179])
    motion_004 = np.array([1, 335, 377, 848, 1069, 1560, 1979, 2164, 2229, 2366, 2721, 2895, 3640, 5183, 6331, 9180, 10485, 11119, 11791, 12820, 13413, 14640, 16579, 17726, 18557, 19428, 20259, 20575, 21327, 26471, 27935, 33118, 34780, 35136, 35532, 36205])
    motion_005 = np.array([225, 1102, 3868, 22138, 24072, 31562])

    key_frames = motion_005
    read_motion(key_frames)

max_stars_repo_name: JasperRice/NAO-G

max_issues_repo_name: JasperRice/NAO-G

max_forks_repo_name: JasperRice/NAO-G


text:
import os, time, sys, struct

import numpy as np

import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib.image import NonUniformImage
import matplotlib.ticker as ticker
from matplotlib import cm

import scipy.interpolate as interp
import scipy.integrate as integrate

import astropy.constants as const
import astropy.units as u
from astropy.cosmology import z_at_value
from astropy.cosmology import WMAP9 as cosmo

from fractions import Fraction

import StrainandNoise_v2 as SnN

current_path = os.getcwd()
splt_path = current_path.split("/")
top_path_idx = splt_path.index('DetectorDesignSensitivities')
top_directory = "/".join(splt_path[0:top_path_idx+1])
load_directory = top_directory + '/LoadFiles/InstrumentFiles/'

def getSNRMatrix(source,instrument,var_x,sampleRate_x,var_y,sampleRate_y):
    # # Setting Up SNR Calculation
    # Uses the variable given and the data range to sample the space either logrithmically or linearly based on the
    # selection of variables. Then it computes the SNR for each value.
    # Returns the variable ranges used to calculate the SNR for each matrix, then returns the SNRs with size of the sample_yXsample_x
    #

    #Get PhenomD waveform
    [phenomD_f,phenomD_h] = source.Get_PhenomD_Strain()
    source.Set_Instrument(instrument)

    #Get Samples for source variables (will return None if they arent var_x or var_y)
    [sample_x,sample_y,recalculate_strain] = Get_Samples(source,var_x,sampleRate_x,var_y,sampleRate_y)
    #Get instrument noise and frequency
    instrument.Get_ASD()
    instrument.Set_f_opt()

    #Check if either sample is not a source variable, if not use instrument samples
    #Very sloppy way of doing it...
    recalculate_noise = 'neither' #Don't need to update inst_var_dict
    if len(sample_x) == 0 and len(sample_y) == 0:
        recalculate_noise = 'both' #Need to update inst_var_dict on x-axis and y-axis
        [sample_x,sample_y,_] = Get_Samples(instrument,var_x,sampleRate_x,var_y,sampleRate_y)
    elif len(sample_x) == 0 and len(sample_y) != 0:
        recalculate_noise = 'x' #Need to update inst_var_dict on x-axis
        [sample_x,_,_] = Get_Samples(instrument,var_x,sampleRate_x,var_y,sampleRate_y)
    elif len(sample_x) != 0 and len(sample_y) == 0:
        recalculate_noise = 'y' #Need to update inst_var_dict on y-axis
        [_,sample_y,_] = Get_Samples(instrument,var_x,sampleRate_x,var_y,sampleRate_y)

    #Make sure samples have the correct astropy units
    [sample_x,sample_y] = Handle_Units(sample_x,var_x,sample_y,var_y)

    sampleSize_x = len(sample_x)
    sampleSize_y = len(sample_y)
    SNRMatrix = np.zeros((sampleSize_x,sampleSize_y))
    tmpSNRMatrix=np.zeros((sampleSize_x,sampleSize_y))

    for i in range(sampleSize_x):
        for j in range(sampleSize_y):

            if recalculate_noise == 'x':
                instrument.Update_Param_val(var_x,sample_x[i])
                source.Update_Param_val(var_y, sample_y[j])
            elif recalculate_noise == 'y':
                source.Update_Param_val(var_x,sample_x[j])
                instrument.Update_Param_val(var_y,sample_y[i])
            elif recalculate_noise == 'both':
                instrument.Update_Param_val(var_x,sample_x[i])
                instrument.Update_Param_val(var_y,sample_y[j])
            elif recalculate_noise == 'neither':
                source.Update_Param_val(var_x,sample_x[i])
                source.Update_Param_val(var_y,sample_y[j])

            if recalculate_noise != 'neither':
                #Recalculate noise curves if something is varied
                instrument.Get_ASD()
                instrument.Set_f_opt()
                source.Set_Instrument(instrument)

            source.checkFreqEvol()
            if source.ismono: #Monochromatic Source and not diff EOB SNR
                SNRMatrix[j,i] = calcMonoSNR(source,instrument)
            else: #Chirping Source
                if recalculate_strain == True: #If we need to calculate the waveform everytime
                    #Get PhenomD waveform
                    [phenomD_f,phenomD_h] = source.Get_PhenomD_Strain()
                source.StrainConv(phenomD_f,phenomD_h)
                SNRMatrix[j,i] = calcChirpSNR(source,instrument)

    return [sample_x,sample_y,SNRMatrix]

def Get_Samples(obj,var_x,sampleRate_x,var_y,sampleRate_y):
    ''' Takes in a object (either for the instrument or source), and the variables
        and sample rates desired in the SNR matrix. The function uses that to create a
        sample space for the variable either in linear space (for q,x1,or x2) or logspace
        for everything else.'''
    sample_x = []
    sample_y = []
    recalculate_strain = False

    var_x_dict = obj.Get_Param_Dict(var_x)
    var_y_dict = obj.Get_Param_Dict(var_y)

    if var_x_dict['min'] != None and var_x_dict['max'] != None: #If the variable has non-None 'min',and 'max' dictionary attributes
        if var_x == 'q' or var_x == 'chi1' or var_x == 'chi2':
            #Sample in linear space for mass ratio and spins
            sample_x = np.linspace(var_x_dict['min'],var_x_dict['max'],sampleRate_x)
            recalculate_strain = True #Must recalculate the waveform at each point
        elif var_x == 'T_obs':
            #sample in linear space for instrument variables
            try:
                sample_x = np.linspace(var_x_dict['min'].value,var_x_dict['max'].value,sampleRate_x)
            except:
                sample_x = np.linspace(var_x_dict['min'],var_x_dict['max'],sampleRate_x)
        else:
            #Sample in log space for any other variables
            #Need exception for astropy variables
            try:
                sample_x = np.logspace(np.log10(var_x_dict['min'].value),np.log10(var_x_dict['max'].value),sampleRate_x)
            except:
                sample_x = np.logspace(np.log10(var_x_dict['min']),np.log10(var_x_dict['max']),sampleRate_x)

    if var_y_dict['min'] != None and var_y_dict['max'] != None: #If the variable has non-None 'min',and 'max' dictionary attributes
        if var_y == 'q' or var_y == 'chi1' or var_y == 'chi2':
            #Sample in linear space for mass ratio and spins
            sample_y = np.linspace(var_y_dict['min'],var_y_dict['max'],sampleRate_y)
            recalculate_strain = True #Must recalculate the waveform at each point
        elif var_y == 'T_obs':
            #sample in linear space for instrument variables
            try:
                sample_y = np.linspace(var_y_dict['min'].value,var_y_dict['max'].value,sampleRate_y)
            except:
                sample_y = np.linspace(var_y_dict['min'],var_y_dict['max'],sampleRate_y)
        else:
            #Sample in log space for any other variables
            #Need exception for astropy variables
            try:
                sample_y = np.logspace(np.log10(var_y_dict['min'].value),np.log10(var_y_dict['max'].value),sampleRate_y)
            except:
                sample_y = np.logspace(np.log10(var_y_dict['min']),np.log10(var_y_dict['max']),sampleRate_y)

    return sample_x,sample_y,recalculate_strain

def Handle_Units(sample_x,var_x,sample_y,var_y):
    '''Since I am using astropy units, I need to update units on the selected samples'''
    #Handle x variables
    if var_x == 'L' or var_x == 'A_IMS':
        sample_x = sample_x*u.m
    elif var_x == 'T_obs' or var_x == 'rms':
        sample_x = sample_x*u.s
    elif var_x == 'A_acc':
        sample_x = sample_x*u.m/u.s/u.s
    elif var_x == 'f_acc_break_high' or var_x == 'f_acc_break_low' or var_x == 'cadence':
        sample_x = sample_x/u.s
    #Handle y variables
    if var_y == 'L' or var_y == 'A_IMS':
        sample_y = sample_y*u.m
    elif var_y == 'T_obs' or var_y == 'rms':
        sample_y = sample_y*u.s
    elif var_y == 'A_acc':
        sample_y = sample_y*u.m/u.s/u.s
    elif var_y == 'f_acc_break_high' or var_y == 'f_acc_break_low' or var_y == 'cadence':
        sample_y = sample_y/u.s

    return sample_x,sample_y

def calcMonoSNR(source,instrument):
    #SNR for a monochromatic source in a PTA
    #From Moore,Taylor,and Gair 2015 https://arxiv.org/abs/1406.5199

    if instrument.name == 'NANOGrav' or instrument.name == 'SKA':
        source.Get_MonoStrain(strain_const='Rosado')
    else:
        source.Get_MonoStrain(strain_const='Cornish')

    indxfgw = np.abs(instrument.fT-source.f_init).argmin()

    SNR = source.h_gw/instrument.S_n_f_sqrt[indxfgw]
    return SNR

def calcChirpSNR(source,instrument):
    #Calculates evolving source using non-precessing binary black hole waveform model IMRPhenomD
    #See Husa et al. 2016 (https://arxiv.org/abs/1508.07250) and Khan et al. 2016 (https://arxiv.org/abs/1508.07253)
    #Uses an interpolated method to align waveform and instrument noise, then integrates
    # over the overlapping region. See eqn 18 from Robson,Cornish,and Liu 2018 https://arxiv.org/abs/1803.01944
    # Values outside of the sensitivity curve are arbitrarily set to 1e30 so the SNR is effectively 0
    S_n_f = instrument.S_n_f_sqrt**2 #Amplitude Spectral Density

    #Only want to integrate from observed frequency (f(T_obs_before_merger)) till merger
    indxfgw = np.abs(source.f-source.f_T_obs).argmin()
    if indxfgw >= len(source.f)-1:
        #If the SMBH has already merged set the SNR to ~0
        return 1e-30
    else:
        f_cut = source.f[indxfgw:]
        h_cut = source.h_f[indxfgw:]

    #################################
    #Interpolate the Strain Noise Spectral Density to only the frequencies the
    #strain runs over
    #Set Noise to 1e30 outside of signal frequencies
    S_n_f_interp_old = interp.interp1d(np.log10(instrument.fT.value),np.log10(S_n_f.value),kind='cubic',fill_value=30.0, bounds_error=False)
    S_n_f_interp_new = S_n_f_interp_old(np.log10(f_cut.value))
    S_n_f_interp = 10**S_n_f_interp_new

    #CALCULATE SNR FOR BOTH NOISE CURVES
    denom = S_n_f_interp #Sky Averaged Noise Spectral Density
    numer = f_cut*h_cut**2

    integral_consts = 16/5 # 4 or(4*4/5) from sky/inclination/polarization averaging

    integrand = numer/denom
    SNRsqrd = integral_consts*np.trapz(integrand.value,np.log10(f_cut.value),axis=0) #SNR**2
    SNR = np.sqrt(SNRsqrd)
    return SNR

def calcDiffSNR(source_var_dict,fT,S_n_f_sqrt,diff_f,diff_h,f_init):
    #Calculates the SNR loss from the difference in EOB waveforms and numerical relativity.
    # The strain is from Sean McWilliams in a private communication.
    #Uses an interpolated method to align waveform and instrument noise, then integrates
    # over the overlapping region. See eqn 18 from Robson,Cornish,and Liu 2018 https://arxiv.org/abs/1803.01944
    # Values outside of the sensitivity curve are arbitrarily set to 1e30 so the SNR is effectively 0
    Vars = []
    for name,sub_dict in source_var_dict.items():
        Vars.append(source_var_dict[name]['val'])

    h_n_f = np.sqrt(fT)*S_n_f_sqrt #characteristic strain Amplitude

    diff_f,diff_h = SnN.StrainConv(Vars,diff_f,diff_h)

    #Only want to integrate from observed frequency (f(T_obs_before_merger)) till merger
    indxfgw = np.abs(diff_f-f_init).argmin()
    if indxfgw >= len(diff_f)-1:
        #If the SMBH has already merged set the SNR to ~0
        return 1e-30
    else:
        f_cut = diff_f[indxfgw:]
        h_cut = diff_h[indxfgw:]

    #################################
    #Interpolate the Strain Noise Spectral Density to only the frequencies the
    #strain runs over
    #Set Noise to 1e30 outside of signal frequencies
    h_n_f_interp_old = interp.interp1d(np.log10(fT.value),np.log10(h_n_f.value),kind='cubic',fill_value=30.0, bounds_error=False)
    h_n_f_interp_new = h_n_f_interp_old(np.log10(f_cut.value))
    h_n_f_interp = 10**h_n_f_interp_new

    #CALCULATE SNR FOR BOTH NOISE CURVES
    denom = h_n_f_interp**2 #Sky Averaged Noise Spectral Density
    numer = h_cut**2

    integral_consts = 4 # 4or(4*4/5) from sky/inclination/polarization averaging

    integrand = numer/denom

    SNRsqrd = integral_consts*np.trapz(integrand.value,np.log10(f_cut.value),axis=0) #SNR**2
    SNR = np.sqrt(SNRsqrd)
    return SNR

def plotSNR(source_var_dict,inst_var_dict,var_x,sample_x,var_y,sample_y,SNRMatrix,display=True,dl_axis=False,figloc=None):
    '''Plots the SNR contours from calcSNR'''
    #Selects contour levels to separate sections into
    contLevels = np.array([5,10, 1e2, 1e3, 1e4, 1e5, 1e6])
    #contLevels = np.array([1,10, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7,1e8,1e9])
    #contLevels = np.array([1,10, 1e2, 1e3, 1e4])
    logLevels = np.log10(contLevels)
    axissize = 14
    labelsize = 16
    textsize = 11
    textcolor = 'k'
    linesize = 4
    figsize = (10,8)

    contourcolorPresent = 'plasma'
    transparencyPresent = 1.0
    contourcolorFuture = 'plasma'
    transparencyFuture = 1.0
    colornorm = colors.Normalize(vmin=0.0, vmax=5.0)
    colormap = 'viridis'
    logSNR = np.log10(SNRMatrix)

    inst_name = list(inst_var_dict.keys())[0]
    try:
        xlabel_min = source_var_dict[var_x]['min']
        xlabel_max = source_var_dict[var_x]['max']
    except:
        xlabel_min = inst_var_dict[inst_name][var_x]['min'].value
        xlabel_max = inst_var_dict[inst_name][var_x]['max'].value
    try:
        ylabel_min = source_var_dict[var_y]['min']
        ylabel_max = source_var_dict[var_y]['max']
    except:
        ylabel_min = inst_var_dict[inst_name][var_y]['min'].value
        ylabel_max = inst_var_dict[inst_name][var_y]['max'].value

    #Can't take log of astropy variables
    try:
        tmp_sample_x = np.log10(sample_x)
    except:
        sample_x = sample_x.value
    try:
        tmp_sample_y = np.log10(sample_y)
    except:
        sample_y = sample_y.value

    #Set whether log or linearly spaced axes
    if var_x == 'q' or var_x == 'chi1' or var_x == 'chi2':
        xaxis_type = 'lin'
    else:
        xaxis_type = 'log'

    if var_y == 'q' or var_y == 'chi1' or var_y == 'chi2':
        yaxis_type = 'lin'
    else:
        yaxis_type = 'log'

    #########################
    #Make the Contour Plots
    fig, ax = plt.subplots(figsize=figsize)

    #Set axis scales based on what data sampling we used
    if yaxis_type == 'lin' and xaxis_type == 'log':
        CS1 = ax.contourf(np.log10(sample_x),sample_y,logSNR,logLevels,cmap = colormap)
        ax.set_xlim(np.log10(xlabel_min),np.log10(xlabel_max))
        ax.set_ylim(ylabel_min,ylabel_max)
        x_labels = np.logspace(np.log10(xlabel_min),np.log10(xlabel_max),np.log10(xlabel_max)-np.log10(xlabel_min)+1)
        y_labels = np.linspace(ylabel_min,ylabel_max,ylabel_max-ylabel_min+1)
        ax.set_yticks(y_labels)
        ax.set_xticks(np.log10(x_labels))
    elif yaxis_type == 'log' and xaxis_type == 'lin':
        CS1 = ax.contourf(sample_x,np.log10(sample_y),logSNR,logLevels,cmap = colormap)
        ax.set_xlim(xlabel_min,xlabel_max)
        ax.set_ylim(np.log10(ylabel_min),np.log10(ylabel_max))
        x_labels = np.linspace(xlabel_min,xlabel_max,xlabel_max-xlabel_min+1)
        y_labels = np.logspace(np.log10(ylabel_min),np.log10(ylabel_max),np.log10(ylabel_max)-np.log10(ylabel_min)+1)
        ax.set_xticks(x_labels)
        ax.set_yticks(np.log10(y_labels))
    else:
        CS1 = ax.contourf(np.log10(sample_x),np.log10(sample_y),logSNR,logLevels,cmap = colormap)
        ax.set_xlim(np.log10(xlabel_min),np.log10(xlabel_max))
        ax.set_ylim(np.log10(ylabel_min),np.log10(ylabel_max))
        x_labels = np.logspace(np.log10(xlabel_min),np.log10(xlabel_max),np.log10(xlabel_max)-np.log10(xlabel_min)+1)
        y_labels = np.logspace(np.log10(ylabel_min),np.log10(ylabel_max),np.log10(ylabel_max)-np.log10(ylabel_min)+1)
        ax.set_yticks(np.log10(y_labels))
        ax.set_xticks(np.log10(x_labels))

    #Set axes labels and whether log or linearly spaced
    if var_x == 'M':
        ax.set_xlabel(r'$M_{\rm tot}$ $[M_{\odot}]$',fontsize = labelsize)
        ax.set_xticklabels([r'$10^{%i}$' %x if int(x) > 1 else r'$%i$' %(10**x) for x in np.log10(x_labels)],fontsize = axissize)
    elif var_x == 'q':
        ax.set_xlabel(r'$q$',fontsize = labelsize)
        ax.set_xticklabels(x_labels,fontsize = axissize)
    elif var_x == 'z':
        ax.set_xlabel(r'${\rm Redshift}$',fontsize = labelsize)
        ax.set_xticklabels([x if int(x) < 1 else int(x) for x in x_labels],fontsize = axissize)
    elif var_x == 'chi1' or var_x == 'chi2':
        ax.set_xlabel(r'${\rm Spin}$',fontsize = labelsize)
        ax.set_xticklabels(x_labels,fontsize = axissize)
    elif var_x == 'L':
        ax.set_xlabel(r'${\rm Armlength}$ $[m]$',fontsize = labelsize)
        ax.set_xticklabels([r'$10^{%i}$' %x if int(x) > 1 else r'$%i$' %(10**x) for x in np.log10(x_labels)],fontsize = axissize)

    if var_y == 'M':
        ax.set_ylabel(r'$M_{\rm tot}$ $[M_{\odot}]$',fontsize = labelsize)
        ax.set_yticklabels([r'$10^{%i}$' %y if int(y) > 1 else r'$%i$' %(10**y) for y in np.log10(y_labels)],fontsize = axissize)
    elif var_y == 'q':
        ax.set_ylabel(r'$q$',fontsize = labelsize)
        ax.set_yticklabels(y_labels,fontsize = axissize)
    elif var_y == 'z':
        ax.set_ylabel(r'${\rm Redshift}$',fontsize = labelsize)
        ax.set_yticklabels([y if int(y) < 1 else int(y) for y in y_labels],\
            fontsize = axissize)
    elif var_y == 'chi1' or var_x == 'chi2':
        ax.set_ylabel(r'${\rm Spin}$',fontsize = labelsize)
        ax.set_yticklabels(y_labels,fontsize = axissize)
    elif var_y == 'L':
        ax.set_ylabel(r'${\rm Armlength}$ $[m]$',fontsize = labelsize)
        ax.set_yticklabels([r'$10^{%i}$' %y if int(y) > 1 else r'$%i$' %(10**y) for y in np.log10(y_labels)],fontsize = axissize)

    ax.yaxis.set_label_coords(-.10,.5)

    #If true, display luminosity distance on right side of plot
    if dl_axis:
        dists = np.array([1e-1,1,10,1e2,1e3,1e4])*u.Gpc
        #dists = np.array([1e-1,1,10,1e2])*u.Gpc
        distticks = [z_at_value(cosmo.luminosity_distance,dist) for dist in dists]
        #Set other side y-axis for lookback time scalings
        ax2 = ax.twinx()
        ax2.contour(np.log10(sample_x),np.log10(sample_y),logSNR,logLevels,colors = 'k',alpha=0.0)
        ax2.set_yticks(np.log10(distticks))
        #ax2.set_yticklabels(['%f' %dist for dist in distticks],fontsize = axissize)
        ax2.set_yticklabels([r'$10^{%i}$' %np.log10(dist) if np.abs(int(np.log10(dist))) > 1 else '{:g}'.format(dist) for dist in dists.value],fontsize = axissize)
        ax2.set_ylabel(r'$D_{L}$ [Gpc]',fontsize=labelsize)
        cbar = fig.colorbar(CS1,ax=(ax,ax2),pad=0.01)
    else:
        #########################
        #Make colorbar
        cbar = fig.colorbar(CS1)
    cbar.set_label(r'$SNR$',fontsize = labelsize)
    cbar.ax.tick_params(labelsize = axissize)
    cbar.ax.set_yticklabels([r'$10^{%i}$' %x if int(x) > 1 else r'$%i$' %(10**x) for x in logLevels])

    if display:
        plt.show()

    #########################
    #Save Figure to File
    if figloc != None:
        fig.savefig(figloc,bbox_inches='tight')

def saveSNR(sample_x,sample_y,SNRMatrix,save_location,SNR_filename,sample_filename):
    #Save SNR Matrix
    np.savetxt(save_location+SNR_filename,SNRMatrix)
    #Save Samples
    np.savetxt(save_location+sample_filename,[sample_x,sample_y])

reason: import numpy,import scipy,import astropy,from astropy


text:
from configparser import ConfigParser
from ETL import get_db, get_collection
from keras.models import load_model
import matplotlib.pyplot as plt
import numpy as np

CONFIG = 'config.ini'
DB = 'ETL'
SECTION = 'THRESHOLD'

cfg = ConfigParser()
cfg.read(CONFIG)

addr = cfg.get(DB, 'address')
port = int(cfg.get(DB, 'port'))
database = cfg.get(DB, 'database')
collection = cfg.get(DB, 'collection')

db = get_db(addr, port, database)
col = get_collection(db, collection)

query_result = [l for l in col.find()]
indx = list(set([l['label'] for l in query_result]))
data = [[] for i in indx]
for l in query_result:
    if l['label'] in indx:
        p = indx.index(l['label'])
        data[p].append(l['sentence'])
# tmpl = [l['label'] for l in query_result]

model = load_model('mymodel.h5')
x1 = range(len(data[0]))
y1 = []
x2 = range(len(data[1]))
y2 = []
for l in data[0]:
    t = np.array(l)
    t = t.reshape((1,t.shape[0]))
    rate = float(model.predict(t, verbose=0)[0][0])
    y1.append(rate)
for l in data[1]:
    t = np.array(l)
    t = t.reshape((1,t.shape[0]))
    rate = float(model.predict(t, verbose=0)[0][0])
    y2.append(rate)

max_stars_repo_name: Hochikong/DumbFinance

max_issues_repo_name: Hochikong/DumbFinance

max_forks_repo_name: Hochikong/DumbFinance


text:
import geopandas as gpd
import matplotlib.font_manager as fm
from tqdm import tqdm
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from Analyst import Network
from Agents import Environment, Commuter
from geopy.distance import distance
from matplotlib import rc
from mpl_toolkits.axes_grid1 import make_axes_locatable
from SB0_Variables import *

class TripDemand:
    def __init__(self, origin_gdf, destination_gdf):
        return

    def get_trip_demand(self, destination_gdf, origin_gdf):
        # Load potential destinations from bc assessment
        print("> Estimating travel demand")
        d_gdf = destination_gdf.reset_index()

        # Convert to WGS84
        d_gdf_4326 = d_gdf.to_crs(4326)
        origin_gdf_4326 = origin_gdf.to_crs(4326)

        # Get destinations within the sandbox
        parcel_gdf = gpd.read_file(proxy.gpkg, layer=f'land_parcels_{exp.lower()}')
        sb_dst = parcel_gdf[(parcel_gdf['Landuse'] == 'CM') | (parcel_gdf['Landuse'] == 'MX')].to_crs(4326)
        final_dst = pd.concat([d_gdf_4326, sb_dst])

        # Get distance from blocks to all destinations
        if f'{exp}_td' not in origin_gdf.columns:
            for j, pt0 in tqdm(enumerate(origin_gdf_4326['geometry'].centroid)):
                distances = []
                for pt1 in list(final_dst['geometry'].centroid):
                    distances.append(distance(pt0.coords[0][::-1], pt1.coords[0][::-1]).km)
                origin_gdf.at[j, f'{exp}_td'] = sum(distances) / len(distances)

        origin_gdf[f'{exp}_car_em'] = (origin_gdf[f"pop_drive_{exp}"] * origin_gdf[f'{exp}_td'] * 0.16 * 3 * 365)
        origin_gdf[f'{exp}_car_em_pc'] = origin_gdf[f'{exp}_car_em'] / origin_gdf[f'population_{exp}']
        origin_gdf[f'{exp}_bus_em'] = (origin_gdf[f"pop_bus_{exp}"] * origin_gdf[f'{exp}_td'] * 0.07 * 3 * 365)
        origin_gdf[f'{exp}_bus_em_pc'] = origin_gdf[f'{exp}_bus_em'] / origin_gdf[f'population_{exp}']
        origin_gdf[f'{exp}_total_em'] = origin_gdf[f'{exp}_car_em'] + origin_gdf[f'{exp}_bus_em']
        origin_gdf[f'{exp}_total_em_pc'] = origin_gdf[f'{exp}_total_em'] / origin_gdf[f'population_{exp}']
        origin_gdf = origin_gdf.to_crs(26910)
        return blocks_gdf

pd.set_option('display.width', 700)
fm.fontManager.ttflist += fm.createFontList(['/Volumes/Samsung_T5/Fonts/roboto/Roboto-Light.ttf'])
rc('font', family='Roboto', weight='light')
name = 'Hillside Quadra'
region = 'Capital Regional District, British Columbia'
directory = f'/Volumes/Samsung_T5/Databases/Sandbox/{name}'
experiments = ['e2', 'e3']# ['e0', 'e1', 'e2', 'e3']
em_modes = ['car', 'bus', 'total']
exp_df = pd.DataFrame()
macc = pd.DataFrame()
tf_years = 20
GeoPackage = f'/Volumes/Samsung_T5/Databases/Capital Regional District, British Columbia.gpkg'

for infra, values in {'bus': 'Frequent transit', 'bike': 'Cycling lanes'}.items():

    # Calculate differences from E0
    proxy_files = {
        'E0': gpd.read_file(f'{directory}/{name} Sandbox_mob_{infra}_e0_na.geojson_{infra}_s0.geojson'),
        'E1': gpd.read_file(f'{directory}/{name} Sandbox_mob_{infra}_e1_na.geojson_{infra}_s0.geojson'),
        'E2': gpd.read_file(f'{directory}/{name} Sandbox_mob_{infra}_e2_na.geojson_{infra}_s0.geojson'),
        'E3': gpd.read_file(f'{directory}/{name} Sandbox_mob_{infra}_e3_na.geojson_{infra}_s0.geojson'),
    }

    # Load blocks layer
    blocks_gdf = gpd.read_file(f'{directory}/Mode Shifts - {infra.title()} - Urban Blocks.geojson')
    blocks_gdf.crs = 26910

    # Spatial join from parcels to grid
    proxy = Network('Hillside Quadra Sandbox', crs=26910, directory=directory)

    # Load parcels and assign block ids
    parcels = gpd.read_file(f"{directory}/Hillside Quadra Sandbox.gpkg", layer='land_parcels_e0').reset_index(drop=True)
    parcels['id'] = parcels.index
    parcels = gpd.sjoin(parcels, blocks_gdf.loc[:, ['geometry']], how="left", lsuffix='', rsuffix='block').groupby('id', as_index=False).first()

    # Create environment for analysis
    env = Environment(
        streets=gpd.read_file(GeoPackage, layer='network_links'),
        trip_diary=pd.read_csv('/Volumes/Samsung_T5/Databases/TransLink/TripDiary2017.csv', index_col='time'),
        origins_gdf=parcels,
        destinations_gdf=gpd.read_file(GeoPackage, layer='indeed_employment')
    )

    # Get trip length of commuters based on income level
    income_levels = {
        'Under $10;000 (including loss)': (0, 10000),
        '$10;000 to $19;999': (10000, 19999),
        '$20;000 to $29;999': (20000, 29999),
        '$30;000 to $39;999': (30000, 39999),
        '$40;000 to $49;999': (40000, 49999),
        '$50;000 to $59;999': (50000, 59999),
        '$60;000 to $69;999': (60000, 69999),
        '$70;000 to $79;999': (70000, 79999),
        '$80;000 to $89;999': (80000, 89999),
        '$90;000 to $99;999': (90000, 99999),
        '$100;000 and over': (100000, 1000000),
    }
    income_ratios = [f'r_{level}' for level in income_levels.keys()]

    # Join income level from census dissemination area to sandbox blocks
    census = gpd.read_file(GeoPackage, layer='land_dissemination_area')
    census.columns = [col.strip() for col in census.columns]
    for level in income_levels.keys():
        census[f"r_{level}"] = census[level] / census.loc[:, income_levels].sum(axis=1)
    blocks_gdf_geom = blocks_gdf['geometry']
    blocks_gdf = gpd.sjoin(blocks_gdf, census.loc[:, income_ratios + ['geometry']], how='left').groupby('id', as_index=False).mean()
    blocks_gdf['geometry'] = blocks_gdf_geom
    blocks_gdf = blocks_gdf.set_geometry('geometry')
    # for col in blocks_gdf.columns: blocks_gdf[col] = blocks_gdf[col].fillna(value=None, method='ffill')
    blocks_gdf = blocks_gdf.fillna(method='ffill', axis=0)

    # Get population on each income level at each experiment
    for level in income_levels:
        for exp in experiments:
            blocks_gdf[f"pop_{level}_{exp}"] = (blocks_gdf[f'population_{exp}'] * blocks_gdf[f"r_{level}"]).fillna(0).astype(int)

    # Get trip length of commuters based on mode choice
    commuters = pd.DataFrame()
    paths = gpd.GeoDataFrame()
    for exp in experiments:
        for block in tqdm(blocks_gdf.index):
            for mode in modes:
                for individual in range(int(blocks_gdf.at[block, f"pop_{mode}_{exp}"])):
                    i = len(commuters)

                    # Randomly choose income level based on probabilities of each level in the block
                    comm = Commuter(
                        env, mode=mode, block_id=block,
                        income=income_levels[np.random.choice(list(income_levels.keys()), p=blocks_gdf.loc[block, income_ratios].dropna())],
                    )

                    commuters.at[i, 'experiment'] = exp
                    commuters.at[i, 'block'] = block
                    commuters.at[i, 'mode'] = mode

                    if mode in ['walk', 'bike']:
                        commuters.at[i, 'length'] = comm.trip_length_from_mode() # (data from TransLink survey 2017)

                    elif mode in ['bus', 'drive']:
                        path = comm.path_from_income()
                        if path is not None:
                            commuters.at[i, 'length'] = path.length.sum()  # (data from Census and Indeed.com)
                            paths = pd.concat([paths, path])
        paths.to_file(f'{directory}{exp.title()}_Paths.geojson', driver='GeoJSON')

    # Estimate emissions based on number of drivers and riders
    demand = False
    if demand:
        crd_gdf = gpd.read_file(GeoPackage, layer='land_assessment_fabric')
        d_gdf = crd_gdf[(crd_gdf['n_use'] == 'CM') | (crd_gdf['n_use'] == 'MX')]
        blocks_gdf = trip_demand(d_gdf, blocks_gdf)
        blocks_gdf.to_file(f'{directory}/UrbanBlocks - {infra.title()}.geojson', driver='GeoJSON')

    if not demand: blocks_gdf = gpd.read_file(f'{directory}/UrbanBlocks - {infra.title()}.geojson')

    # Estimate costs per mode on each experiment
    km_per_liter = 9
    price_per_liter = 1.2
    mtn_cost_per_km = 0.07 #https://driving.ca/auto-news/news/part-2-how-much-should-standard-car-maintenance-cost
    annual_insurance = 1400 #https://www.ratelab.ca/how-much-does-it-cost-to-own-and-operate-a-car-in-canada/
    cost_car_year_bc = 700 * 12 #City of Victoria, 2009. Go Victoria: Our urban-mobility Future
    transit_fare = 2.5
    bike_mtn_per_yr = 120 #https://bicycles.stackexchange.com/questions/15874/cost-of-maintenance-of-a-roadbike
    for i, exp in enumerate(experiments):

        # Calculate emissions/person saved from E0
        exp_df.at[i, f'{infra}_em_saved'] = blocks_gdf[f'e0_total_em'].sum() - blocks_gdf[f'{exp}_total_em'].sum()
        exp_df.at[i, f'{infra}_em_saved_per_inh'] = (blocks_gdf[f'e0_total_em'].sum() / blocks_gdf[f'population_e0'].sum()) - (
                    blocks_gdf[f'{exp}_total_em'].sum() / blocks_gdf[f'population_{exp}'].sum())

        for mode in modes:
            if mode == 'drive':
                # Cost based on individual estimates (Old version)
                blocks_gdf[f"cost_drive_{exp}"] = \
                    blocks_gdf[f"pop_drive_{exp}"] * annual_insurance * \
                    ((blocks_gdf[f'{exp}_td']/km_per_liter) * price_per_liter) * \
                        (blocks_gdf[f'{exp}_td'] * mtn_cost_per_km)
                # Cost based on aggregated data (City of Victoria, 2009)
                blocks_gdf[f"cost_drive_{exp}"] = cost_car_year_bc * blocks_gdf[f"pop_drive_{exp}"]
            elif mode == 'bus':
                blocks_gdf[f"cost_bus_{exp}"] = blocks_gdf[f"pop_bus_{exp}"] * transit_fare * 2 * 365
            elif mode == 'bike':
                blocks_gdf[f"cost_bike_{exp}"] = blocks_gdf[f"pop_bike_{exp}"] * bike_mtn_per_yr
            elif mode == 'walk':
                blocks_gdf[f"cost_walk_{exp}"] = 0
            exp_df.at[i, f"cost_{mode}"] = blocks_gdf[f"cost_{mode}_{exp}"].sum()

        # Aggregate costs per experiment
        exp_df["total_cost"] = exp_df[f"cost_bike"] + exp_df[f"cost_drive"] + exp_df[f"cost_bus"]
        exp_df.at[i, "pop"] = blocks_gdf[f'population_{exp}'].sum()
        exp_df["cost_per_cap"] = exp_df["total_cost"]/exp_df["pop"]

    ### Infrastructure costs ###

    # https://www.victoria.ca/EN/main/residents/transportation/cycling/new-cycling-projects/funding.html
    bike_infra_cost_km = 100000

    # https://www.bctransit.com/documents/1529706193497
    vendor_pay = 69262676 + 11577349
    employee_net_pay = 43554814
    operations = 242546000
    maintenance = 55699000
    administration = 30419000
    total_fleet = 1+38+29+128+309+3+87+75+67+12+33+10+6+338+4 # https://bctransit.com/about/fleet
    amortization = 61047000
    debt = 9106000

    # https://bctransit.com/about/facts/victoria
    rev_fares = 37500000
    rev_pass = 4700000
    province = 4760000
    prop_tax = 30100000
    fuel_tax = 20400000
    adv = 700000
    total_trips = 27000000

    # Not sourced
    route_in_boundary = 1.6 # km

    local_gbd = Network('Hillside Quadra Sandbox', crs=26910, directory='/Volumes/Samsung_T5/Databases/Sandbox/Hillside Quadra')

    # Get average length of routes that crosses the neighbourhood
    region = Network('Capital Regional District, British Columbia')
    routes = gpd.read_file(region.gpkg, layer='network_routes')

    # Get routes that intersect with sandbox
    boundary = gpd.read_file(local_gbd.gpkg, layer='land_municipal_boundary')
    routes_overlay = gpd.overlay(routes, boundary)
    routes_overlay['length'] = [geom.length for geom in routes_overlay['geometry']]

    # Calculate the ratio of routes that intersects with the sandbox boundary
    overlay_ratio = routes_overlay['length']/routes[routes['route'].isin(routes_overlay['route'])].reset_index(drop=True)['length']
    ave_ratio_route_in_boundary = overlay_ratio.mean()

    # Calculate cost per trip
    total_cost = operations + maintenance + administration
    total_revenue = rev_fares + rev_pass + province + prop_tax + fuel_tax + adv
    cost = total_cost - total_revenue
    cost_per_trip = (cost/total_trips) * ave_ratio_route_in_boundary

    links = gpd.read_file(local_gbd.gpkg, layer='network_links')
    stops = gpd.read_file(local_gbd.gpkg, layer='network_stops')
    for i, exp in enumerate(experiments):
        if exp == 'e0':
            yr = 2020
            exp_df.at[i, "cycling_cost"] = (links[links[f'cycle_2020'] == 1].length.sum() / 1000) * bike_infra_cost_km
        else:
            yr = 2040
            exp_df.at[i, "cycling_cost"] = (links[(links[f'cycle_2040'] == 1) & (links[f'cycle_2020'] == 0)].length.sum()/1000) * bike_infra_cost_km
        exp_df.at[i, "transit_cost"] = stops[f'frequency_{yr}'].sum() * cost_per_trip * 365 * tf_years

    if infra == 'bus': infra_cost = exp_df["transit_cost"]
    elif infra == 'bike': infra_cost = exp_df["cycling_cost"]
    else: infra_cost = 0

    exp_df[f"{infra}_infra_cost"] = infra_cost
    exp_df[f"{infra}_infra_cost_per_cap"] = exp_df[f"{infra}_infra_cost"] / exp_df["pop"]

    # Plot results
    fig_size = (15, 10)
    fig1 = plt.figure(constrained_layout=False, figsize=fig_size)
    ax1 = {}
    main = 12
    widths = [main, main, main, main, 1]
    heights = [main, main, main]
    gs = fig1.add_gridspec(nrows=len(em_modes), ncols=len(experiments)+1, width_ratios=widths, height_ratios=heights)

    for i, (mode, cmap) in enumerate(zip(em_modes, ['Reds', 'Blues', 'viridis_r'])):
        ax1[i] = {}
        for j, exp in enumerate(experiments):
            print(f"\nPlotting results on {exp}(i: {i}) for {mode}(j: {j})")

            # Create non-NaN GeoDataFrame
            blocks_valid = blocks_gdf.dropna()

            # Calculate mean and median
            mean = blocks_valid[f'{exp}_{mode}_em_pc'].mean()

            # Get minimum and maximum values on each experiment
            cols = [f'{e}_{mode}_em_pc' for e in experiments]
            vmin = min(blocks_valid.loc[:, cols].min())
            vmax = max(blocks_valid.loc[:, cols].max())
            print(f"{cols} min: {vmin}, max: {vmax}, mean: {mean}")

            # Plot block maps with emissions per each mode per capita
            ax1[i][j] = fig1.add_subplot(gs[i, j])
            blocks_valid.plot(f"{exp}_{mode}_em_pc", ax=ax1[i][j], legend=False, vmin=vmin, vmax=vmax, cmap=cmap)
            ax1[i][j].set_title(f"{exp.upper()}, {mode.upper()} | MEAN: {round(mean/1000, 2)} tCO2/p/y")
            ax1[i][j].set_axis_off()

            if j == len(experiments) - 1:
                ax1[i][j + 1] = fig1.add_subplot(gs[i, j + 1])
                divider = make_axes_locatable(ax1[i][j + 1])
                leg_ax = divider.append_axes(position="bottom", size="100%", pad="0%", add_to_figure=False)
                array = np.arange(vmin, vmax)
                show = leg_ax.imshow([array], cmap=cmap, aspect='auto')
                cb = fig1.colorbar(show, cax=ax1[i][j + 1], orientation="vertical")
                cb.set_label('Emissions (tCO2/person/year)')

    # Export plots and maps to files
    plt.tight_layout()
    fig1.savefig(f'{directory}/Mode Shifts - Emissions per Capita - {infra}.png')

    # Export data on MACC DataFrame
    for i, exp in enumerate(experiments):
        if exp != 'e0':
            j = len(macc)
            macc.at[j, 'Abatement Measure'] = f'Densification {exp.title()} + {values}'
            macc.at[j, 'GHGs Abated (tCO2/per capita)'] = (exp_df.loc[i, f"{infra}_em_saved_per_inh"]/1000) * tf_years
            macc.at[j, 'GHGs Abated (tCO2)'] = (macc.at[j, 'GHGs Abated (tCO2/per capita)']) * exp_df.at[i, "pop"]
            macc.at[j, 'Net Cost of Abatement Measure ($)'] = exp_df.loc[i, f"{infra}_infra_cost"]

    macc['Marginal Cost ($/tCO2)'] = macc['Net Cost of Abatement Measure ($)']/macc['GHGs Abated (tCO2)']

print(macc)
print("Exporting MACC to excel")
macc.to_excel(f'{directory}/Mobility MACC.xlsx')
print(f"End")

max_stars_repo_name: nicholas-martino/UrbanMobility

max_issues_repo_name: nicholas-martino/UrbanMobility

max_forks_repo_name: nicholas-martino/UrbanMobility


text:
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import sys
import warnings
import numpy as np
from macros import colors, markers

def operators(dim, operator_location):
    # Given 1, 2, or 4 coordinates like x-y-z-w, depending on dim and operators
    num_operators = len(operator_location.split('-')) / dim
    if num_operators == 1:
        return 'One Operator'
    elif num_operators == 2:
        return 'Two Operators'

def fabric(dim, kary):
    return '{}D {}:1'.format(dim, kary)

def algorithm_font(algorithm):
    return r'\textsf{{{}}}'.format(algorithm)

def plot_mera_treewidth_comparison(data_filename, plot_filename, verbose):
    # Use latex font
    plt.rc('text', usetex=True)
    plt.rc('font', family='serif')

    # Set up Seaborn style
    sns.set(style="darkgrid")

    # Import the dataframe
    dataframe = pd.read_csv(data_filename)
    dataframe['operators'] =\
        np.vectorize(operators)(dataframe['dim'],
                                dataframe['operator_position'])
    dataframe['fabric'] =\
        np.vectorize(fabric)(dataframe['dim'], dataframe['kary'])
    dataframe['algorithm'] =\
        np.vectorize(algorithm_font)(dataframe['algorithm'])

    dataframe = dataframe.rename(
        columns={'time': 'Run Time (sec)',
                 'treewidth': 'Contraction Complexity'})
    if verbose:
        print(dataframe)

    # # Make facet grid of scatter plots
    # grid = sns.FacetGrid(dataframe,
    #                      col="operators",
    #                      row="fabric",
    #                      hue="algorithm",
    #                      palette=[colors[x] for x in ['freetdi', 'meiji', 'netcon']],
    #                      hue_kws=dict(marker=[markers[x] for x in ['freetdi', 'meiji', 'netcon']]),
    #                      margin_titles=True,
    #                      hue_order=['\\textsf{freetdi}',
    #                                 '\\textsf{meiji-e}',
    #                                 '\\textsf{netcon}'])
    # plot = grid.map(plt.scatter,
    #                 "Contraction Complexity",
    #                 "Run Time (sec)",
    #                 alpha=0.7,
    #                 linewidth=0.1,
    #                 edgecolor='black')

    facet_kws = dict()
    warnings.simplefilter(action='ignore', category=FutureWarning)
    plot = sns.factorplot(data=dataframe,
                          col='operators',
                          row='fabric',
                          x="Contraction Complexity",
                          y="Run Time (sec)",
                          hue="algorithm",
                          palette=[colors[x] for x in ['freetdi',
                                                       'meiji',
                                                       'netcon']],
                          facet_kws=facet_kws,
                          kind="strip",
                          dodge=True,
                          jitter=True,
                          alpha=0.7,
                          linewidth=0.1,
                          aspect=1.7,
                          size=2.5,
                          hue_order=['\\textsf{freetdi}',
                                     '\\textsf{meiji-e}',
                                     '\\textsf{netcon}'],
                          legend=False)

    # Manually add dashed lines to facets
    for axis in plot.fig.get_axes():
        axis.axhline(y=1200, color='black', dashes=[3, 3])
        for i in range(len(dataframe["Contraction Complexity"]) - 1):
            axis.axvline(x=i+.5, c="white", dashes=(2, 1))

    # Set axis lengths and format
    plot.set(ylim=(.001, 10000), yscale='log')

    # Set axis labels
    plot.set_titles(col_template="{col_name}", row_template="{row_name} MERA")

    # Add legend
    plot.fig.get_axes()[0].legend(loc="lower right")

    # Save figure
    for extension in ['.pdf', '.png']:
        plot.savefig(plot_filename + extension)

if __name__ == '__main__':
    data_filename = sys.argv[1]
    plot_filename = sys.argv[2]
    plot_mera_treewidth_comparison(data_filename, plot_filename, False)

max_stars_repo_name: TheoryInPractice/ConSequences

max_issues_repo_name: TheoryInPractice/ConSequences

max_forks_repo_name: TheoryInPractice/ConSequences


text:
import datetime
import numpy as np
from scipy import sparse as sp

import torch
import torch.nn as nn
import torch.nn.functional as F
from torch import optim
from torch.utils.data import DataLoader

from tqdm import tqdm

from ..dataset import RecFeatDataset, collate_triplets_with_feature
from ..optimizer import MultipleOptimizer
from ..utils import scisp2tchsp
from ..metric import compute_ndcg

# TODO: make decorator for checking the model is already fitted
class UserItemFM(nn.Module):
    def __init__(self, k, init=0.001, n_iters=10, learn_rate=0.001, l2=0.0001,
                 n_negs=10, use_gpu=False, loss='sgns', loss_agg="sum",
                 alpha=None, no_item_factor=False, warm_start=False):
        """"""
        super().__init__()
        self.k = k
        self.init = init
        self.learn_rate = learn_rate
        self.l2 = l2
        self.n_iters = n_iters
        self.n_negs = n_negs
        # TODO: generalization of the device selection (not only for cuda)
        self.target_device = 'cuda' if use_gpu else 'cpu'
        self.no_item_factor = no_item_factor
        self.warm_start = warm_start

        # setup the loss
        if loss not in {'sgns', 'bpr', 'bce', 'kl'}:
            raise ValueError(
                '[ERROR] Loss should be one of {"sgns", "bpr", "bce", "kl"}'
            )
        self.loss = loss  # {'sgns', 'bpr', 'bce', 'kl', 'mse'}
        self.full_batch = True if loss in {'bce', 'kl'} else False
        if loss == 'bpr' and n_negs > 1:
            print(
                '[Warning] BPR loss only samples 1 negative ' +
                'sample! fall back to n_negs=1...'
            )
            self.n_negs = 1

        self.loss_agg = loss_agg
        self.alpha = alpha  # for weighted loss

    @property
    def device(self):
        return next(self.parameters()).device

    def _init_embeddings(self, user_item, user_feature=None, item_feature=None):
        """"""
        # initialize the factors
        self.is_features_ = (
            False
            if (user_feature is None) and (item_feature is None)
            else True
        )
        item_sparse = True if self.loss in {'sgns', 'bpr'} else False
        self.register_parameter('w0', nn.Parameter(torch.FloatTensor([0])))
        self.embeddings_ = nn.ModuleDict({
            'user': nn.Embedding(user_item.shape[0], self.k + 1, sparse=True),
            'item': nn.Embedding(user_item.shape[1],
                                 self.k + 1, sparse=item_sparse),
        })
        if user_feature is not None:
            self.embeddings_['feat_user'] = nn.Embedding(
                user_feature.shape[1], self.k + 1,
                sparse=sp.issparse(user_feature)
            )
        elif item_feature is not None:
            self.embeddings_['feat_item'] = nn.Embedding(
                item_feature.shape[1], self.k + 1,
                sparse=sp.issparse(item_feature)
            )

        for key, layer in self.embeddings_.items():
            if self.no_item_factor and key == 'item':
                layer.weight.data.zero_()
                layer.weight.requires_grad = False
            else:
                layer.weight.data[:,:-1].normal_(std=self.init)
                layer.weight.data[:,-1].zero_()

    def _init_optimizer(self):
        """
        """
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        params = {'sparse':[], 'dense':[]}
        for lyr in self.embeddings_.children():
            if lyr.sparse:
                params['sparse'].append(lyr.weight)
            else:
                params['dense'].append(lyr.weight)

        # determine if sparse optimizer needed
        if len(params['sparse']) == 0:
            self.opt = optim.Adam(
                params['dense'], lr=self.learn_rate, weight_decay=self.l2
            )
        elif len(params['dense']) == 0:
            self.opt = optim.SparseAdam(params['sparse'], lr=self.learn_rate)
        else:
            # init multi-optimizers
            # register it to the instance
            self.opt = MultipleOptimizer(
                optim.SparseAdam(params['sparse'], lr=self.learn_rate),
                optim.Adam(params['dense'],
                           lr=self.learn_rate,
                           weight_decay=self.l2)
            )

    def _retrieve_factors(self, u, i, feats=None):
        """"""
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        # retrieve embeddins
        user = self.embeddings_['user'](u)
        item = self.embeddings_['item'](i)
        # post-process (devide factor and weight)
        w = torch.cat([user[..., -1], item[..., -1]], dim=1)
        v = torch.cat([user[..., :-1], item[..., :-1]], dim=1)
        for entity, feat in feats.items():
            k = 'feat_{}'.format(entity)
            if k in self.embeddings_:
                if self.embeddings_[k].sparse:
                    feat = self.embeddings_[k](feats)
                else:
                    feat = (
                        self.embeddings_[k].weight[None] * feat[..., None]
                    )
                # post-process (devide factor and weight)
                w = torch.cat([w, feat[..., -1]], dim=1)
                v = torch.cat([v, feat[..., :-1]], dim=1)
        return w, v

    def _update_z(self, features):
        """
        this method pre-compute & cache item-tag factor for prediction
        """
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        self.cpu()  # since below process eats up lots of memory

        # update z
        for entity in ['user', 'item']:
        # for entity, feat in features.items():
            v = self.embeddings_[entity].weight.detach()
            k = 'feat_{}'.format(entity)

            if entity in features:
            # if k in self.embeddings_:
                feat = features[entity].to('cpu')
                vf = self.embeddings_[k].weight.detach()
                zf = feat @ vf
                zf2 = feat**2 @ vf**2
                z = (v + zf)
                z2 = (v**2 + zf2)
            else:
                z, z2 = v, (v**2)

            # register
            self.register_buffer('z{}'.format(entity[0]), z)
            self.register_buffer('z{}2'.format(entity[0]), z2)
        self.to(self.target_device)

    def predict_user(self, users):
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        if np.isscalar(users):
            users = [users]

        if not self.training:
            s, _ = self.forward_batch(users, None)
            return s
        else:
            raise ValueError('You should change to the .eval() mode first!')

    def forward(self, u, i, feats):
        """
        u (torch.LongTensor): user ids
        i (torch.LongTensor): item ids
        feats (torch.FloatTensor): item feature tensor
        """
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        w, v = self._retrieve_factors(u, i, feats)
        w_ = w.sum(1)
        v_ = (v.sum(1)**2 - (v**2).sum(1)).sum(-1) * .5
        s = self.w0 + w_ + v_
        return s, [w, v]

    def forward_batch(self, u, X=None):
        """"""
        if not hasattr(self, 'embeddings_'):
            raise ValueError('You should call .fit first!')

        # retrieve user factors
        u = torch.LongTensor(u).to(self.device)
        vu = self.embeddings_['user'](u)

        if self.training:
            # update item z
            vi = self.embeddings_['item'].weight

            if 'feat_item' in self.embeddings_:
                vfi = self.embeddings_['feat_item'].weight
                zfi = X['item'] @ vfi
                zfi2 = X['item']**2 @ vfi**2

                zi = (vi + zfi)[None]  # (None, n_items, k+1)
                zi2 = (vi**2 + zfi2)[None]
            else:
                zi, zi2 = vi[None], (vi**2)[None]

            if 'feat_user' in self.embeddings_:
                vfu = self.embeddings_['feat_user'].weight
                zfu = X['user'][u] @ vfu
                zfu2 = X['user'][u]**2 @ vfu**2

                zu = (vu + zfu)[:, None]
                zu2 = (vu**2 + zfu2)[:, None]
            else:
                zu, zu2 = vu[:, None], (vu**2)[:, None]

        else:
            zi, zi2 = self.zi[None], self.zi2[None]
            if 'feat_user' in self.embeddings_:
                zu, zu2 = self.zu[u][:, None], self.zu2[u][:, None]
            else:
                zu, zu2 = vu[:, None], (vu**2)[:, None]

        # forward & compute loss
        w = zi[..., -1] + zu[..., -1]
        v = (zi[..., :-1] + zu[..., :-1])**2
        v -= (zi2[..., :-1] + zu2[..., :-1])
        v = v.sum(-1) * .5
        s = self.w0 + w + v
        return s, [vu]

    def _compute_loss(self, s, y, c=None, loss='sgns',
                      aggregate="mean", alpha=1., weights=None):
        """
        s (torch.Tensor): scores
        y (torch.Tensor): target binary indicator (either {0, 1} or {-1, 1})
        c (torch.Tensor): confidence wrt to the y (optional)
        loss (str): type of the loss
        aggregate (str): aggregation method {'mean', 'sum', 'batchmean'}
        alpha (float): 1.
        weights (list of torch.Parameter):
                    weight whose l2 contraint to be explicitly computed
        """
        if loss in {'sgns', 'bpr'}:
            s = s.reshape(-1, self.n_negs + 1)
            y = y.reshape(-1, self.n_negs + 1)
            if loss == 'sgns':
                # requires y \in {-1, 1}
                l = -F.logsigmoid(s * y)

            elif loss == 'bpr':
                # assuming the first columns of scores are the postive
                # and the second is the negative
                # check the shape (should by (batch_sz, 2))
                if s.ndim != 2 or s.shape[-1] != 2:
                    raise ValueError('Should be 2 dimentional and only' +
                        'sample 1 negative sample')

                l = -F.logsigmoid((s * y).sum(1))[:, None]

            else:
                raise ValueError('[ERROR] only support {"sgns", "bpr", "cce", "kl"}')

            if aggregate == 'mean':
                l = l.mean()
            elif aggregate == 'sum':
                l = l.sum()
            elif aggregate == 'batchmean':
                l = l.sum(1).mean(0)

        else:
            # prepare the confidence matrix
            if c is not None:
                c = c.to_dense()
                c = c * alpha + 1

            y = y.to_dense()
            if loss == 'bce':
                if aggregate == 'batchmean':
                    aggregate = 'mean'

                l = F.binary_cross_entropy_with_logits(
                    s, y, reduction=aggregate, weight=c
                )

            elif loss == 'kl':
                l = F.kl_div(
                    torch.log_softmax(s, dim=-1), y / y.sum(1)[:, None],
                    reduction=aggregate
                )

            else:
                raise ValueError('[ERROR] only support {"sgns", "bpr", "cce", "kl"}')

        if (weights is not None) and (len(weights) > 0) and (self.l2 > 0):
            # adding L2 regularization term for sparse entries (user/item)
            l += self.l2 * sum([torch.sum(w**2) for w in weights])

        return l

    def fit(self, user_item, user_feature=None, item_feature=None,
            valid_user_item=None, feature_sparse=False, batch_sz=512,
            n_tests=1000, topk=100, valid_callback=None, verbose=False,
            keep_best=False, n_jobs=4):
        """"""
        # set some variables
        feats = {}
        do_valid = valid_user_item is not None
        n_users, n_items = user_item.shape
        desc_tmp = '[tloss=0.0000]'

        if not self.warm_start:
            self._init_embeddings(user_item, user_feature, item_feature)

        if self.target_device != 'cpu':
            self.to(self.target_device)

        if user_feature is not None:
            feats['user'] = torch.Tensor(user_feature).to(self.device)
        if item_feature is not None:
            feats['item'] = torch.Tensor(item_feature).to(self.device)

        # init optimizer
        self._init_optimizer()

        # prepare dataset
        if self.full_batch:
            dl = np.arange(0, n_users, batch_sz)
        else:
            dl = DataLoader(
                RecFeatDataset(user_item, n_negs=self.n_negs,
                               item_feature=item_feature,
                               user_feature=user_feature),
                batch_sz, collate_fn=collate_triplets_with_feature,
                shuffle=True, num_workers=n_jobs, drop_last=True
            )

        # do the optimization
        c = None
        val_score = 0
        self.best_val = 0
        self.best_model = None
        try:
            for _ in range(self.n_iters):
                with tqdm(total=len(dl), desc=desc_tmp,
                          ncols=80, disable=not verbose) as prog:

                    iter_loss = 0
                    val_timer = datetime.datetime.now()
                    rnd_usr = np.random.permutation(n_users)
                    for sample in dl:
                        # parse the output
                        if self.full_batch:
                            u = rnd_usr[sample:sample + batch_sz]
                            y = user_item[u].copy().tocoo()

                            # for loss weight
                            if (self.alpha is not None) and (self.alpha > 0):
                                c = y.copy()
                                c = scisp2tchsp(c).to(self.device)

                            y.data[:] = 1
                            y = scisp2tchsp(y).to(self.device)

                        else:
                            sample = tuple(d.to(self.device) for d in sample)
                            u, i, y, xu, xi = sample

                        # flush grads
                        self.opt.zero_grad()

                        # forward the model
                        if self.loss in {'bce', 'kl'}:
                            s, w = self.forward_batch(u, feats)
                        elif self.loss in {'sgns', 'bpr'}:
                            s, w = self.forward(u, i, {'user':xu, 'item':xi})

                        # compute the main loss
                        loss = self._compute_loss(
                            s, y, c, self.loss,
                            aggregate=self.loss_agg, weights=w
                        )

                        # backward
                        loss.backward()

                        # update params
                        self.opt.step()

                        # update loss and counter
                        iter_loss += loss.item()
                        now = datetime.datetime.now()
                        elapsed = now - val_timer
                        if do_valid and (elapsed.total_seconds() > 60):
                            val_timer = now
                            val_score = self.validate(
                                user_item, valid_user_item, feats,
                                n_tests, topk, valid_callback
                            )
                        prog.set_description(
                            '[tloss={:.4f} / vacc={:.4f}]'
                            .format(loss.item(), val_score)
                        )
                        prog.update(1)

                    if do_valid:
                        val_score = self.validate(
                            user_item, valid_user_item, feats,
                            n_tests, topk, valid_callback
                        )
                    prog.set_description(
                        '[tloss={:.4f} / vacc={:.4f}]'
                        .format(iter_loss / len(dl), val_score)
                    )

        except KeyboardInterrupt as e:
            print('User stopped the training!...')
        finally:
            # update the cached factors for faster inference
            if do_valid and keep_best:
                self.load_state_dict(self.best_model)
            self._update_z(feats)

    def validate(self, user_item, valid_user_item, feats, n_tests, topk,
                 callback=None):
        """"""
        if callback is None:
            scores = self._validate(
                user_item, valid_user_item, feats,
                n_tests, topk
            )
        else:
            scores = callback(
                self, user_item, valid_user_item, feats,
                n_tests, topk
            )

        # keep the best morel
        val_score = np.mean(scores)
        if val_score > self.best_val:
            self.best_val = val_score
            self.best_model = self.state_dict()

        return val_score

    def _validate(self, user_item, valid_user_item, feats,
                  n_tests, topk):
        self._update_z(feats)
        self.eval()
        scores = []
        n_users, n_items = user_item.shape
        test_users = np.random.choice(n_users, n_tests, False)
        for u in test_users:
            s = self.predict_user(u)[0]
            true = valid_user_item[u]
            train = user_item[u]
            score = compute_ndcg(s, true, train, topk)
            if score is None: continue
            scores.append(score)
        self.train()
        return scores

max_stars_repo_name: eldrin/factormach

max_issues_repo_name: eldrin/factormach

max_forks_repo_name: eldrin/factormach

reason: import numpy,from scipy


text:
#!/usr/bin/env python

# Copyright 2017 Verily Life Sciences Inc.
#
# Use of this source code is governed by a BSD-style
# license that can be found in the LICENSE file.

r"""Convert a particular hdf5 file to sparse, long format.

This code is based loosely on
http://cf.10xgenomics.com/supp/cell-exp/megacell_tutorial.html
and is hardcoded to load
https://support.10xgenomics.com/single-cell/datasets/1M_neurons
"""

import sys
import argparse
import collections
import numpy as np
import os
import scipy.sparse as sp_sparse
import tables

GENOME = 'mm10'

# This dsub-compatible script will read configuration from the environment,
# if available.
DEFAULT_BEGIN_IDX = os.getenv('BEGIN_IDX', 0)  # First cell in the file.
DEFAULT_END_IDX = os.getenv('END_IDX', 1306127)  # Last cell in the file.
DEFAULT_INPUT_FILE = os.getenv('INPUT_FILE',
                               '1M_neurons_filtered_gene_bc_matrices_h5.h5')
DEFAULT_OUTPUT_FILE = os.getenv('OUTPUT_FILE')

np.random.seed(0)

GeneBCMatrix = collections.namedtuple(
    'GeneBCMatrix',
    ['gene_ids', 'gene_names', 'barcodes', 'matrix'])

def get_matrix_from_h5(filename, genome):
  """Load the matrix from the HDF5 file.

  Code is from http://cf.10xgenomics.com/supp/cell-exp/megacell_tutorial.html

  Args:
    filename: HDF5 filename
    genome: Genome of data in the file.

  Returns:
    the sparse matrix of data
  """
  with tables.open_file(filename, 'r') as f:
    try:
      dsets = {}
      for node in f.walk_nodes('/' + genome, 'Array'):
        dsets[node.name] = node.read()
      matrix = sp_sparse.csc_matrix(
          (dsets['data'], dsets['indices'], dsets['indptr']),
          shape=dsets['shape'])
      return GeneBCMatrix(
          dsets['genes'], dsets['gene_names'], dsets['barcodes'], matrix)
    except tables.NoSuchNodeError:
      raise Exception('Genome %s does not exist in %s.' % (genome, filename))
    except KeyError:
      raise Exception('File %s missing one or more required datasets.'
                      % filename)

def run(argv=None):
  """Runs the variant preprocess pipeline.

  Args:
    argv: Pipeline options as a list of arguments.
  """
  parser = argparse.ArgumentParser()
  parser.add_argument(
      '--begin-idx',
      type=int,
      default=DEFAULT_BEGIN_IDX,
      help='Index with which to start reading data (inclusive).')
  parser.add_argument(
      '--end-idx',
      type=int,
      default=DEFAULT_END_IDX,
      help='Index at which to stop reading data (exclusive).')
  parser.add_argument(
      '--input-file',
      default=DEFAULT_INPUT_FILE,
      help='Input file path.')
  parser.add_argument(
      '--output-file',
      default=DEFAULT_OUTPUT_FILE,
      help='Output file path. If None, stdout will be used.')
  args = parser.parse_args(argv)

  sys.stderr.write('Processing cells [%d,%d) from file %s'
                   % (args.begin_idx, args.end_idx, args.input_file))
  mm10_gbm = get_matrix_from_h5(args.input_file, GENOME)

  handle = open(args.output_file, 'w') if args.output_file else sys.stdout

  # Emit the output CSV file header.
  handle.write(','.join(['gene_id', 'gene', 'cell', 'trans_cnt']) + '\n')

  for cell in range(args.begin_idx, args.end_idx):
    for sparse_idx in range(mm10_gbm.matrix.indptr[cell],
                            mm10_gbm.matrix.indptr[cell+1]):
      dense_idx = mm10_gbm.matrix.indices[sparse_idx]
      handle.write(','.join([mm10_gbm.gene_ids[dense_idx],
                             mm10_gbm.gene_names[dense_idx],
                             mm10_gbm.barcodes[cell],
                             str(mm10_gbm.matrix.data[sparse_idx])])
                   + '\n')

  if handle is not sys.stdout:
    handle.close()

if __name__ == '__main__':
  run()