text:
#!/usr/bin/python

"""
Brief python script to calculate multiyear seasonal means of TROPOMI UT NO2.

This is done by weighting single year seasonal means by the number of
observations.

"""
# Import packages:
import glob
import sys
import os
from os import path
from netCDF4 import Dataset

import numpy as np

# Steps to follow:
# (1) Define input:
# Season of interest (eventually change to argparse)
SEAS = "jja"
# Data product of interest:
PRODUCT = "fresco"
# Data directory:
IODIR = "/home/ucfaea1/python/uptrop/Data/"
# Years of interest:
YEARS = ['2019', '2020']

# (2) Get files of interest:
tfile1 = path.join(IODIR, 'tropomi-ut-no2-fresco-07-1x1-' + SEAS + '-' + YEARS[0] + '-v6.nc')
#tfile1 = glob.glob(tfile_glob_string)

tfile2 = path.join(IODIR, 'tropomi-ut-no2-fresco-07-1x1-' + SEAS + '-' + YEARS[1] + '-v6.nc')
#tfile2 = glob.glob(tfile_glob_string)

utno2_files = [tfile1, tfile2]

# Get number of files:
nfiles = len(utno2_files)

# Initialize:
first = 0

# Loop over files:
for f in range(nfiles):

    # (3) Read data from files:
    fh = Dataset(utno2_files[f], mode='r')

    # latitude:
    tlat = fh.variables['lat'][:]

    # longitude:
    tlon = fh.variables['lon'][:]

    # TROPOMI UT NO2:
    tutno2 = fh.variables['utno2'][:]

    # TROPOMI UT NO2 error:
    tutno2err = fh.variables['utno2err'][:]

    # Cloud top pressre range:
    tcld_top_p_range = fh.variables['cld_top_p_range'][:]

    # Cloud top pressure ceiling:
    tcld_top_p_ceil = fh.variables['cld_top_p_ceil'][:]

    # Number of observations:
    tnobs = fh.variables['nobs'][:]

    # Close file:
    fh.close()

    if first==0:
        # Get dimensions:
        nx = len(tlon)
        ny = len(tlat)

        # Define arrays:
        my_utno2 = np.zeros((nx, ny))
        my_utno2err = np.zeros((nx, ny))
        my_cld_top_p_range = np.zeros((nx, ny))
        my_cld_top_p_ceil = np.zeros((nx, ny))
        my_nobs = np.zeros((nx, ny))

        # Redefine:
        first = 1

    # Add individual years:
    indices = np.where(~np.isnan(tnobs))
    my_utno2[indices] += tutno2[indices] * tnobs[indices]
    my_utno2err[indices] += tutno2err[indices] * tnobs[indices]
    my_cld_top_p_range[indices] += tcld_top_p_range[indices] * tnobs[indices]
    my_cld_top_p_ceil[indices] += tcld_top_p_ceil[indices] * tnobs[indices]
    my_nobs[indices] += tnobs[indices]

# Get multiyear mean:
my_utno2 = np.where(~np.isnan(my_nobs), np.divide(my_utno2, my_nobs), np.nan)
my_utno2err = np.where(~np.isnan(my_nobs), np.divide(my_utno2err, my_nobs), np.nan)
my_cld_top_p_range = np.where(~np.isnan(my_nobs), np.divide(my_cld_top_p_range, my_nobs), np.nan)
my_cld_top_p_ceil = np.where(~np.isnan(my_nobs), np.divide(my_cld_top_p_ceil, my_nobs), np.nan)

# (5) Save data:
outfile = path.join(IODIR, 'tropomi-ut-no2-fresco-07-1x1-' + SEAS + '-multiyr-mean-2019-2020-v6.nc')

ncout = Dataset(outfile, mode='w',format="NETCDF4")

ncout.createDimension('lat', ny)
ncout.createDimension('lon', nx)

# create longitude axis:
lon = ncout.createVariable('lon', np.float32, ('lon',))
lon.units = 'degrees_east'
lon.long_name = 'longiitude'
lon[:] = tlon

# Create latitude axis:
lat = ncout.createVariable('lat', np.float32, ('lat',))
lat.units = 'degrees_north'
lat.long_name = 'latitude'
lat[:] = tlat

# Save data: UT NO2 VMR (gno2vmr), UT NO2 error (gerr), No. of data points (gcnt)
utno2 = ncout.createVariable('utno2', np.float32, ('lon', 'lat'))
utno2.units = 'pptv'
utno2.long_name = 'NO2 mixing ratios in the UT (180-450 hPa) obtained using cloud-slicing'
utno2[:] = my_utno2

utno2err = ncout.createVariable('utno2err', np.float32, ('lon', 'lat'))
utno2err.units = 'pptv'
utno2err.long_name = 'Standard error of the NO2 mixing ratios in the UT (180-450 hPa) obtained using cloud-slicing'
utno2err[:] = my_utno2err

utcldrange = ncout.createVariable('cld_top_p_range', np.float32, ('lon', 'lat'))
utcldrange.units = 'hPa'
utcldrange.long_name = 'Gridded mean range in cloud top pressures used to cloud-slice TROPOMI NO2'
utcldrange[:] = my_cld_top_p_range

utcldceil = ncout.createVariable('cld_top_p_ceil', np.float32, ('lon', 'lat'))
utcldceil.units = 'hPa'
utcldceil.long_name = 'Gridded mean ceiling of cloud top pressures used to cloud-slice TROPOMI NO2'
utcldceil[:] = my_cld_top_p_ceil

nobs = ncout.createVariable('nobs', np.float32, ('lon', 'lat'))
nobs.units = 'unitless'
nobs.long_name = 'Number of observations in each gridsquare used to obtain cloud-sliced UT NO2 mixing ratios'
nobs[:] = my_nobs

ncout.close()

max_stars_repo_name: eamarais/erc-uptrop

max_issues_repo_name: eamarais/erc-uptrop

max_forks_repo_name: eamarais/erc-uptrop


text:
"""
Implements deep maximum entropy inverse reinforcement learning based on
Ziebart et al., 2008 and Wulfmeier et al., 2015, using symbolic methods with
Theano.

Matthew Alger, 2015
matthew.alger@anu.edu.au
"""

from itertools import product

import numpy as np
import numpy.random as rn
import theano as th
import theano.tensor as T
import time
import maxent
import pickle

FLOAT = th.config.floatX

def find_svf(n_states, trajectories):
    """
    Find the state vistiation frequency from trajectories.

    n_states: Number of states. int.
    trajectories: 3D array of state/action pairs. States are ints, actions
        are ints. NumPy array with shape (T, L, 2) where T is the number of
        trajectories and L is the trajectory length.
    -> State visitation frequencies vector with shape (N,).
    """

    svf = np.zeros(n_states)

    for trajectory in trajectories:
        for state, _, _ in trajectory:
            svf[state] += 1

    svf /= trajectories.shape[0]

    return th.shared(svf, "svf", borrow=True)

def optimal_value(n_states, n_actions, transition_probabilities, reward,
                  discount, threshold=1e-2):
    """
    Find the optimal value function.

    n_states: Number of states. int.
    n_actions: Number of actions. int.
    transition_probabilities: Function taking (state, action, state) to
        transition probabilities.
    reward: Vector of rewards for each state.
    discount: MDP discount factor. float.
    threshold: Convergence threshold, default 1e-2. float.
    -> Array of values for each state
    """

    v = T.zeros(n_states, dtype=FLOAT)

    def update(s, prev_diff, v, reward, tps):
        max_v = float("-inf")
        v_template = T.zeros_like(v)
        for a in range(n_actions):
            tp = tps[s, a, :]
            max_v = T.largest(max_v, T.dot(tp, reward + discount*v))
        new_diff = abs(v[s] - max_v)
        if T.lt(prev_diff, new_diff):
            diff = new_diff
        else:
            diff = prev_diff
        return (diff, T.set_subtensor(v_template[s], max_v)), {}

    def until_converged(diff, v):
        (diff, vs), _ = th.scan(
                fn=update,
                outputs_info=[{"initial": diff, "taps": [-1]},
                              None],
                sequences=[T.arange(n_states)],
                non_sequences=[v, reward, transition_probabilities])
        return ((diff[-1], vs.sum(axis=0)), {},
                th.scan_module.until(diff[-1] < threshold))

    (_, vs), _ = th.scan(fn = until_converged,
                         outputs_info=[
                            # Need to force an inf into the right Theano
                            # data type and this seems to be the only way that
                            # works.
                            {"initial": getattr(np, FLOAT)(float("inf")),
                             "taps": [-1]},
                            {"initial": v,
                             "taps": [-1]}],
                         n_steps=1000)

    return vs[-1]

def find_policy(n_states, n_actions, transition_probabilities, reward, discount,
                threshold=1e-2, v=None):
    """
    Find the optimal policy.

    n_states: Number of states. int.
    n_actions: Number of actions. int.
    transition_probabilities: Function taking (state, action, state) to
        transition probabilities.
    reward: Vector of rewards for each state.
    discount: MDP discount factor. float.
    threshold: Convergence threshold, default 1e-2. float.
    v: Optimal value array (if known). Default None.
    -> Action probabilities for each state.
    """

    if v is None:
        v = optimal_value(n_states, n_actions, transition_probabilities, reward,
                          discount, threshold)

    # Get Q using equation 9.2 from Ziebart's thesis.
    Q = T.zeros((n_states, n_actions))
    def make_Q(i, j, tps, Q, reward, v):
        Q_template = T.zeros_like(Q)
        tp = transition_probabilities[i, j, :]
        return T.set_subtensor(Q_template[i, j], tp.dot(reward + discount*v)),{}

    prod = np.array(list(product(range(n_states), range(n_actions))))
    state_range = th.shared(prod[:, 0])
    action_range = th.shared(prod[:, 1])
    Qs, _ = th.scan(fn=make_Q,
                    outputs_info=None,
                    sequences=[state_range, action_range],
                    non_sequences=[transition_probabilities, Q, reward, v])
    Q = Qs.sum(axis=0)
    Q -= Q.max(axis=1).reshape((n_states, 1))  # For numerical stability.
    Q = T.exp(Q)/T.exp(Q).sum(axis=1).reshape((n_states, 1))
    return Q

def find_expected_svf(n_states, r, n_actions, discount,
                      transition_probability, trajectories):
    """
    Find the expected state visitation frequencies using algorithm 1 from
    Ziebart et al. 2008.

    n_states: Number of states N. int.
    alpha: Reward. NumPy array with shape (N,).
    n_actions: Number of actions A. int.
    discount: Discount factor of the MDP. float.
    transition_probability: NumPy array mapping (state_i, action, state_k) to
        the probability of transitioning from state_i to state_k under action.
        Shape (N, A, N).
    trajectories: 3D array of state/action pairs. States are ints, actions
        are ints. NumPy array with shape (T, L, 2) where T is the number of
        trajectories and L is the trajectory length.
    -> Expected state visitation frequencies vector with shape (N,).
    """

    n_trajectories = trajectories.shape[0]
    trajectory_length = trajectories.shape[1]

    policy = find_policy(n_states, n_actions,
                         transition_probability, r, discount)

    start_state_count = T.extra_ops.bincount(trajectories[:, 0, 0],
                                             minlength=n_states)
    p_start_state = start_state_count.astype(FLOAT)/n_trajectories

    def state_visitation_step(i, j, prev_svf, policy, tps):
        """
        The sum of the outputs of a scan over this will be a row of the svf.
        """

        svf = prev_svf[i] * policy[i, j] * tps[i, j, :]
        return svf, {}

    prod = np.array(list(product(range(n_states), range(n_actions))))
    state_range = th.shared(prod[:, 0])
    action_range = th.shared(prod[:, 1])
    def state_visitation_row(prev_svf, policy, tps, state_range, action_range):
        svf_t, _ = th.scan(fn=state_visitation_step,
                           sequences=[state_range, action_range],
                           non_sequences=[prev_svf, policy, tps])
        svf_t = svf_t.sum(axis=0)
        return svf_t, {}

    svf, _ = th.scan(fn=state_visitation_row,
                     outputs_info=[{"initial": p_start_state, "taps": [-1]}],
                     n_steps=trajectories.shape[1]-1,
                     non_sequences=[policy, transition_probability, state_range,
                                 action_range])

    return svf.sum(axis=0) + p_start_state

def saveNN(structure, file_name, weights, biases, α, e):
    N_layers = len(structure)
    for l in range(N_layers-1):
        # save weights
        fo_NNweights = file_name+'_NNpara_l'+str(l+1)+'W.csv'
        weights_temp = weights[l].get_value()
        np.savetxt(fo_NNweights, weights_temp, delimiter = ',')
        # save biases
        fo_NNbiases = file_name+'_NNpara_l'+str(l+1)+'B.csv'
        biases_temp = biases[l].get_value()
        np.savetxt(fo_NNbiases, biases_temp, delimiter = ',')
    # save alpha
    fo_alpha = file_name+'_NNpara_alpha.csv'
    alpha_temp = α.get_value()
    np.savetxt(fo_alpha, alpha_temp, delimiter = ',')

    print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' NN parameters of round '+str(e)+' have been saved.')

def irl(structure, feature_matrix, n_actions, discount, transition_probability,
        trajectories, epochs,  learning_rate, file_name, initialisation="normal", l1=0.1,
        l2=0.1):
    #
    """
    Find the reward function for the given trajectories.

    structure: Neural network structure tuple, e.g. (10, 3, 3) would be a
        3-layer neural network with 10 inputs.
    feature_matrix: Matrix with the nth row representing the nth state. NumPy
        array with shape (N, D) where N is the number of states and D is the
        dimensionality of the state.
    n_actions: Number of actions A. int.
    discount: Discount factor of the MDP. float.
    transition_probability: NumPy array mapping (state_i, action, state_k) to
        the probability of transitioning from state_i to state_k under action.
        Shape (N, A, N).
    trajectories: 3D array of state/action pairs. States are ints, actions
        are ints. NumPy array with shape (T, L, 2) where T is the number of
        trajectories and L is the trajectory length.
    epochs: Number of gradient descent steps. int.
    learning_rate: Gradient descent learning rate. float.
    initialisation: What distribution to use. str in {normal, uniform}. Default
        normal.
    l1: L1 regularisation. Default 0.1. float.
    l2: L2 regularisation. Default 0.1. float.
    -> Reward vector with shape (N,).
    """

    n_states, d_states = feature_matrix.shape
    transition_probability = th.shared(transition_probability, borrow=True)
    trajectories = th.shared(trajectories, borrow=True)

    # Initialise W matrices; b biases.
    n_layers = len(structure)-1
    weights = []
    hist_w_grads = []  # For AdaGrad.
    biases = []
    hist_b_grads = []  # For AdaGrad.
    for i in range(n_layers):
        # W
        shape = (structure[i+1], structure[i])
        if initialisation == "normal":
            matrix = th.shared(rn.normal(size=shape), name="W", borrow=True)
        else:
            matrix = th.shared(rn.uniform(size=shape), name="W", borrow=True)
        weights.append(matrix)
        hist_w_grads.append(th.shared(np.zeros(shape), name="hdW", borrow=True))

        # b
        shape = (structure[i+1], 1)
        if initialisation == "normal":
            matrix = th.shared(rn.normal(size=shape), name="b", borrow=True)
        else:
            matrix = th.shared(rn.uniform(size=shape), name="b", borrow=True)
        biases.append(matrix)
        hist_b_grads.append(th.shared(np.zeros(shape), name="hdb", borrow=True))

    # Initialise α weight, β bias.
    if initialisation == "normal":
        α = th.shared(rn.normal(size=(1, structure[-1])), name="alpha",
                      borrow=True)
    else:
        α = th.shared(rn.uniform(size=(1, structure[-1])), name="alpha",
                      borrow=True)
    hist_α_grad = T.zeros(α.shape)  # For AdaGrad.

    adagrad_epsilon = 1e-6  # AdaGrad numerical stability.

    #### Theano symbolic setup. ####

    # Symbolic input.
    s_feature_matrix = T.matrix("x")
    # Feature matrices.
    # All dimensions of the form (d_layer, n_states).
    φs = [s_feature_matrix.T]
    # Forward propagation.
    for W, b in zip(weights, biases):
        φ = T.nnet.sigmoid(th.compile.ops.Rebroadcast((0, False), (1, True))(b)+ W.dot(φs[-1]))
        φs.append(φ)
        # φs[1] = φ1 etc.
    # Reward.
    r = α.dot(φs[-1]).reshape((n_states,))
    # Engineering hack: z-score the reward.
    r = (r - r.mean())/r.std()
    # Associated feature expectations.
    expected_svf = find_expected_svf(n_states, r,
                                     n_actions, discount,
                                     transition_probability,
                                     trajectories)
    svf = maxent.find_svf(n_states, trajectories.get_value())
    # Derivatives (backward propagation).
    updates = []
    α_grad = φs[-1].dot(svf - expected_svf).T
    hist_α_grad += α_grad**2
    adj_α_grad = α_grad/(adagrad_epsilon + T.sqrt(hist_α_grad))
    updates.append((α, α + adj_α_grad*learning_rate))

    def grad_for_state(s, theta, svf_diff, r):
        """
        Calculate the gradient with respect to theta for one state.
        """

        regularisation = abs(theta).sum()*l1 + (theta**2).sum()*l2
        return svf_diff[s] * T.grad(r[s], theta) - regularisation, {}

    #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' enumerate(weights) starts.')
    for i, W in enumerate(weights):
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' enumerate(weights) '+str(i)+'.')
        w_grads, _ = th.scan(fn=grad_for_state,
                             sequences=[T.arange(n_states)],
                             non_sequences=[W, svf - expected_svf, r])
        w_grad = w_grads.sum(axis=0)
        hist_w_grads[i] += w_grad**2
        adj_w_grad = w_grad/(adagrad_epsilon + T.sqrt(hist_w_grads[i]))
        updates.append((W, W + adj_w_grad*learning_rate))

    #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' enumerate(biases) starts.')
    for i, b in enumerate(biases):
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' enumerate(biases) '+str(i)+'.')
        b_grads, _ = th.scan(fn=grad_for_state,
                             sequences=[T.arange(n_states)],
                             non_sequences=[b, svf - expected_svf, r])
        b_grad = b_grads.sum(axis=0)
        hist_b_grads[i] += b_grad**2
        adj_b_grad = b_grad/(adagrad_epsilon + T.sqrt(hist_b_grads[i]))
        updates.append((b, b + adj_b_grad*learning_rate))

    train = th.function([s_feature_matrix], updates=updates, outputs=r)
    run = th.function([s_feature_matrix], outputs=r)

    # build result file
    fo_reward = open(file_name+'_reward_record.csv',"w")
    for n in range(n_states):
        fo_reward.write('S'+str(n)+', ')
    fo_reward.write('\n')
    fo_reward.close()

    for e in range(epochs):
        if (e+1)%5 == 0:
            print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' reward NN training round '+str(e+1)+ ' starts.')
        reward = train(feature_matrix)
        # reward write
        fo_reward = open(file_name+'_reward_record.csv',"a")
        for n in range(n_states):
            fo_reward.write(str(reward[n])+', ')
        fo_reward.write('\n')
        fo_reward.close()
        if (e+1)%5== 0:
            print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' reward: '+str(reward[[0,5,10,15,20,25,30]])+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' reward: '+str(r.mean()+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' weights1: '+str(weights[0].get_value() )+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' weights2: '+str(weights[1].get_value() )+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' biases1: '+str(biases[0].get_value() )+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' biases2: '+str(biases[1].get_value() )+ '.')
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' alpha: '+str(α.get_value() )+ '.')
        #rewardtest = run(feature_matrix)
        #print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' rewardtest: '+str(rewardtest[[0,5,10,15,20,25,30]])+ '.')

        if (e+1)%20 == 0:
            saveNN(structure, file_name, weights, biases, α, e)

    print(str(time.strftime("%Y%m%d %X", time.localtime()) )+' NN training is done.')

    return reward.reshape((n_states,)), weights, biases

ext: py

reason: import numpy,import theano


text:
from __future__ import print_function, division, unicode_literals
import numpy as np
from psy import McmcHoDina
from psy.utils import r4beta

attrs = np.random.binomial(1, 0.5, (5, 60))

g = r4beta(1, 2, 0, 0.6, (1, 60))
no_s = r4beta(2, 1, 0.4, 1, (1, 60))

theta = np.random.normal(0, 1, (1000, 1))
lam00 = np.random.normal(0, 1, 5)
lam11 = np.random.uniform(0.5, 3, 5)

ho_dina = McmcHoDina(attrs=attrs)
skills_p = McmcHoDina.get_skills_p(lam0=lam00, lam1=lam11, theta=theta)
skills = np.random.binomial(1, skills_p)

yita = ho_dina.get_yita(skills)
p_val = ho_dina.get_p(yita, guess=g, no_slip=no_s)
score = np.random.binomial(1, p_val)

ho_dina_est = McmcHoDina(attrs=attrs, score=score, max_iter=10000, burn=5000)
est_lam0, est_lam1, est_theta, est_skills, est_no_s, est_g = ho_dina_est.mcmc()

max_stars_repo_name: cegfdb/IRT

max_issues_repo_name: glan-wxl/IRT

max_forks_repo_name: glan-wxl/IRT

reason: import numpy


text:
import os
import json

import numpy as np
import torch
from tqdm.auto import tqdm

from models import CNN
from configs import CnnConfig
from utils import load_json
import cv2
import glob

def replace_nan(x, y):
    x = 0.0 if np.isnan(x) else x
    y = 0.0 if np.isnan(y) else y
    return x, y

def draw_hands(
    image,
    hand_x,
    hand_y,
    connections,
    connection_color,
    thickness,
    point_color,
    frame_length,
    frame_width,
):
    for connection in connections:
        x0 = hand_x[connection[0]] * frame_width
        y0 = hand_y[connection[0]] * frame_length
        x1 = hand_x[connection[1]] * frame_width
        y1 = hand_y[connection[1]] * frame_length

        x0, y0 = replace_nan(x0, y0)
        x1, y1 = replace_nan(x1, y1)

        cv2.line(
            image,
            (int(x0), int(y0)),
            (int(x1), int(y1)),
            connection_color,
            thickness,
        )

    for x, y in zip(hand_x, hand_y):
        x, y = replace_nan(x, y)
        cv2.circle(image, (int(x), int(y)), thickness, point_color, thickness)

    return image

def draw_pose(
    image,
    pose_x,
    pose_y,
    links,
    connection_color,
    thickness,
    point_color,
    frame_length,
    frame_width,
):
    for link in links:
        x0 = pose_x[link[0]] * frame_width
        y0 = pose_y[link[0]] * frame_length
        x1 = pose_x[link[1]] * frame_width
        y1 = pose_y[link[1]] * frame_length

        x0, y0 = replace_nan(x0, y0)
        x1, y1 = replace_nan(x1, y1)

        cv2.line(
            image,
            (int(x0), int(y0)),
            (int(x1), int(y1)),
            connection_color,
            thickness,
        )
        cv2.circle(image, (int(x0), int(y0)), thickness, point_color, thickness)
        cv2.circle(image, (int(x1), int(y1)), thickness, point_color, thickness)

    return image

def cnn_feat(file_path, save_dir):
    video_record = load_json(file_path)

    FRAME_LENGTH = 1080
    FRAME_WIDTH = 1920
    POINT_COLOR = (255, 0, 0)
    CONNECTION_COLOR = (0, 255, 0)
    THICKNESS = 2

    connections = [
        (0, 1),
        (1, 2),
        (2, 3),
        (3, 4),
        (5, 6),
        (6, 7),
        (7, 8),
        (9, 10),
        (10, 11),
        (11, 12),
        (13, 14),
        (14, 15),
        (15, 16),
        (17, 18),
        (18, 19),
        (19, 20),
        (0, 5),
        (5, 9),
        (9, 13),
        (13, 17),
        (0, 17),
    ]

    links = [
        (11, 12),
        (11, 23),
        (12, 24),
        (23, 24),
        (11, 13),
        (13, 15),
        (12, 14),
        (14, 16),
        (15, 21),
        (15, 17),
        (17, 19),
        (19, 15),
        (22, 16),
        (16, 18),
        (18, 20),
        (16, 20),
    ]

    model = CNN(CnnConfig)
    features = np.empty((0, CnnConfig.output_dim))

    assert video_record["n_frames"] > 0, "Number of frames should be greater than zero"
    for i in range(video_record["n_frames"]):
        image = np.zeros((FRAME_LENGTH, FRAME_WIDTH, 3), np.uint8)
        pose_x = video_record["pose_x"][i]
        pose_y = video_record["pose_y"][i]
        hand1_x = video_record["hand1_x"][i]
        hand1_y = video_record["hand1_y"][i]
        hand2_x = video_record["hand2_x"][i]
        hand2_y = video_record["hand2_y"][i]

        if hand1_x[0] != 0:
            image = draw_hands(
                image,
                hand1_x,
                hand1_y,
                connections,
                CONNECTION_COLOR,
                THICKNESS,
                POINT_COLOR,
                FRAME_LENGTH,
                FRAME_WIDTH,
            )

        if hand2_x[0] != 0:
            image = draw_hands(
                image,
                hand2_x,
                hand2_y,
                connections,
                CONNECTION_COLOR,
                THICKNESS,
                POINT_COLOR,
                FRAME_LENGTH,
                FRAME_WIDTH,
            )

        image = draw_pose(
            image,
            pose_x,
            pose_y,
            links,
            CONNECTION_COLOR,
            THICKNESS,
            POINT_COLOR,
            FRAME_LENGTH,
            FRAME_WIDTH,
        )
        image = image.astype(np.float32) / 255
        image = cv2.resize(image, (224, 224))
        feat = model(torch.FloatTensor(image).permute(2, 0, 1).unsqueeze(0))
        features = np.vstack([features, feat.numpy()])

    save_path = os.path.join(save_dir, video_record["uid"] + ".npy")
    np.save(save_path, features)

def runner(args, mode):
    json_files = sorted(
        glob.glob(
            os.path.join(args.data_dir, f"{args.dataset}_{mode}_keypoints", "*.json")
        )
    )

    save_dir = os.path.join(args.data_dir, f"{args.dataset}_{mode}_features")
    if not os.path.exists(save_dir):
        os.mkdir(save_dir)

    saved_files = glob.glob(os.path.join(save_dir, "*.npy"))
    if len(saved_files) != len(json_files):
        for file_path in tqdm(json_files, desc=f"Saving CNN features - {mode} files"):
            cnn_feat(file_path, save_dir)
    else:
        print(mode, "CNN features already exist!")

def save_cnn_features(args):
    runner(args, mode="train")
    runner(args, mode="val")
    runner(args, mode="test")

max_stars_repo_name: munishm/INCLUDE-teams-expt

max_issues_repo_name: munishm/INCLUDE-teams-expt

max_forks_repo_name: munishm/INCLUDE-teams-expt


text:
import cv2
import numpy as np
#from math import sqrt

windowName = 'tic tac toe'
img = np.zeros((512, 512, 3), np.uint8)
cv2.namedWindow(windowName)
count =1

a = [[2,2,2],[2,2,2],[2,2,2]]

for i in range (106,406):

    cv2.line(img,(306,i),(306,i+1),(255,0,0),2)
    cv2.line(img,(206,i),(206,i+1),(255,0,0),2)
    cv2.waitKey(5)
    cv2.imshow(windowName, img)

for j in range (106,406):

    cv2.line(img,(j,206),(j+1,206),(255,0,0),2)
    cv2.line(img,(j,306),(j+1,306),(255,0,0),2)
    cv2.waitKey(5)
    cv2.imshow(windowName, img)

def win_check(b, m):
    l = [[[0,0],[0,1],[0,2]],[[1,0],[1,1],[1,2]],[[2,0],[2,1],[2,2]],[[0,0],[1,0],[2,0]],[[0,1],[1,1],[2,1]],[[0,2],[1,2],[2,2]],[[0,0],[1,1],[2,2]],[[2,0],[1,1],[0,2]]]
    k = 0
    for i in l:

        p = 0
        for j in i:
            c = j[0]
            d = j[1]

            if b[c][d] == m:
                p+=1
        if p==3:
            start = i[0]
            end = i[2]
            k+=1
    if k>0:

        x = [156,256,356]
        y = [156,256,356]
        x_1 = start[0]
        x_2 = start[1]
        y_1 = end[0]
        y_2 = end[1]
        x1 = x[x_1]
        y1 = y[y_1]
        x2 = x[x_2]
        y2 = y[y_2]

        if m == 1:
            cv2.line(img, (x1,x2),(y1,y2), (0, 255, 0), 3)
        else:
            cv2.line(img, (x1,x2),(y1,y2), (0, 0, 255), 3)

        return True

    else:
        return False

    pass

def full_board_check(a):
    k = 0
    for i in range(0,3):
      for j in range (0,3):
        if a[i][j] == 2:
            k+=1
    if k == 0:
        return True
    else:
        return False

    pass

def win_text(win):

    if win == 1:

        image = cv2.rectangle(img, (106,406), (406,512), (0,0,0), -1)
        cv2.imshow(windowName, image)
        font = cv2.FONT_HERSHEY_SIMPLEX
        org = (106, 459)
        image = cv2.putText(img, 'Player 1 has won!!', org, font, 1, (0, 0, 255), 2, cv2.LINE_AA)
        cv2.imshow(windowName, image)

    if win == 2:

        image = cv2.rectangle(img, (106,406), (406,512), (0,0,0), -1)
        cv2.imshow(windowName, image)
        font = cv2.FONT_HERSHEY_SIMPLEX
        org = (106, 459)
        image = cv2.putText(img, 'Player 2 has won!!', org, font, 1, (0, 255, 0), 2, cv2.LINE_AA)
        cv2.imshow(windowName, image)

    if win == 3:

        image = cv2.rectangle(img, (106,406), (406,512), (0,0,0), -1)
        cv2.imshow(windowName, image)
        font = cv2.FONT_HERSHEY_SIMPLEX
        org = (106, 459)
        image = cv2.putText(img, 'Match is Drawn', org, font, 1, (0, 255, 255), 2, cv2.LINE_AA)
        cv2.imshow(windowName, image)

def draw_shape(event, x, y, flags, param):

    if event == cv2.EVENT_LBUTTONDOWN :

        def isInside(circle_x, circle_y, rad, x, y):

            if ((x - circle_x) * (x - circle_x) + (y - circle_y) * (y - circle_y) <= rad * rad):
                return True;
            else:
                return False;

        circlex = [156,256,356]
        circley = [156,256,356]
        rad = 50;
        global a
        for i in range (0,3):
            for j in range(0,3):
                circle1_x = circlex[i]
                circle1_y = circley[j]
                if(isInside(circle1_x, circle1_y, rad, x, y)):
                    global count
                    if(count%2 == 1):
                        if a[i][j] == 2:
                            cv2.line(img, (circle1_x-20,circle1_y-20),(circle1_x+20,circle1_y+20), (0, 0, 255), 3)
                            cv2.line(img, (circle1_x+20,circle1_y-20),(circle1_x-20,circle1_y+20), (0, 0, 255), 3)
                            a[i][j] = 1

                            win1  = win_check(a, 1)
                            if win1 == True:
                                victory = 1
                                win_text(victory)
                            count+=1

                    else:
                        if a[i][j] == 2:
                            cv2.circle(img, (circle1_x, circle1_y), 25, (0, 255, 0), 3)
                            a[i][j] = 0

                            win2 = win_check(a, 0)
                            if win2 == True:
                                victory = 2
                                win_text(victory)
                            count+=1

        if count == 10:
            victory =3
            win_text(victory)

cv2.setMouseCallback(windowName, draw_shape)

def main():

    while(True):
        cv2.imshow(windowName, img)

        k = cv2.waitKey(1)

        if k == ord('Q') :

            break
        if k == ord('q') :

            break
        if k == 27:
            break

    cv2.destroyAllWindows()

if __name__ == "__main__":

    main()

max_stars_repo_path: tic tac toe.py

max_stars_repo_name: dukedaffy/tic-tac-toe-opencv

max_issues_repo_path: tic tac toe.py

max_issues_repo_name: dukedaffy/tic-tac-toe-opencv

max_forks_repo_path: tic tac toe.py

max_forks_repo_name: dukedaffy/tic-tac-toe-opencv


text:
import logging
log = logging.getLogger(__name__)

import numpy
import operator

import westpa, west
from westpa.yamlcfg import check_bool
from westpa.kinetics import RateAverager
from westext.weed.ProbAdjustEquil import probAdjustEquil

EPS = numpy.finfo(numpy.float64).eps

class WEEDDriver:
    def __init__(self, sim_manager, plugin_config):
        if not sim_manager.work_manager.is_master:
            return

        self.sim_manager = sim_manager
        self.data_manager = sim_manager.data_manager
        self.system = sim_manager.system
        self.work_manager = sim_manager.work_manager

        self.do_reweight = check_bool(plugin_config.get('do_equilibrium_reweighting', False)
                           or plugin_config.get('do_reweighting', False))
        self.windowsize = 0.5
        self.windowtype = 'fraction'

        windowsize = plugin_config.get('window_size')
        if windowsize is not None:
            if isinstance(windowsize, float):
                self.windowsize = windowsize
                self.windowtype = 'fraction'
                if self.windowsize <= 0 or self.windowsize > 1:
                    raise ValueError('WEED parameter error -- fractional window size must be in (0,1]')
            elif isinstance(windowsize,int):
                self.windowsize = int(windowsize)
                self.windowtype = 'fixed'
            else:
                raise ValueError('WEED parameter error -- invalid window size {!r}'.format(windowsize))
        log.info('using window size of {!r} ({})'.format(self.windowsize, self.windowtype))

        self.max_windowsize = plugin_config.get('max_window_size')
        if self.max_windowsize is not None:
            log.info('Using max windowsize of {:d}'.format(self.max_windowsize))

        self.reweight_period = plugin_config.get('reweight_period', 0)
        self.priority = plugin_config.get('priority', 0)

        self.rate_calc_queue_size = plugin_config.get('rate_calc_queue_size', 1)
        self.rate_calc_n_blocks = plugin_config.get('rate_calc_n_blocks', 1)

        if self.do_reweight:
            sim_manager.register_callback(sim_manager.prepare_new_iteration,self.prepare_new_iteration, self.priority)

    def get_rates(self, n_iter, mapper):
        '''Get rates and associated uncertainties as of n_iter, according to the window size the user
        has selected (self.windowsize)'''

        if self.windowtype == 'fraction':
            if self.max_windowsize is not None:
                eff_windowsize = min(self.max_windowsize,int(n_iter * self.windowsize))
            else:
                eff_windowsize = int(n_iter * self.windowsize)
        else: # self.windowtype == 'fixed':
            eff_windowsize = min(n_iter, self.windowsize or 0)

        averager = RateAverager(mapper, self.system, self.data_manager, self.work_manager)
        averager.calculate(max(1, n_iter-eff_windowsize), n_iter+1, self.rate_calc_n_blocks, self.rate_calc_queue_size)
        self.eff_windowsize = eff_windowsize

        return averager

    def prepare_new_iteration(self):
        n_iter = self.sim_manager.n_iter
        we_driver = self.sim_manager.we_driver

        if we_driver.target_states and self.do_reweight:
            log.warning('equilibrium reweighting requested but target states (sinks) present; reweighting disabled')
            return

        if not self.do_reweight:
            # Reweighting not requested
            log.debug('equilibrium reweighting not enabled')
            return

        with self.data_manager.lock:
            weed_global_group = self.data_manager.we_h5file.require_group('weed')
            last_reweighting = int(weed_global_group.attrs.get('last_reweighting', 0))

        if n_iter - last_reweighting < self.reweight_period:
            # Not time to reweight yet
            log.debug('not reweighting')
            return
        else:
            log.debug('reweighting')

        mapper = we_driver.bin_mapper
        bins = we_driver.next_iter_binning
        n_bins = len(bins)

        # Create storage for ourselves
        with self.data_manager.lock:
            iter_group = self.data_manager.get_iter_group(n_iter)
            try:
                del iter_group['weed']
            except KeyError:
                pass

            weed_iter_group = iter_group.create_group('weed')
            avg_populations_ds = weed_iter_group.create_dataset('avg_populations', shape=(n_bins,), dtype=numpy.float64)
            unc_populations_ds = weed_iter_group.create_dataset('unc_populations', shape=(n_bins,), dtype=numpy.float64)
            avg_flux_ds = weed_iter_group.create_dataset('avg_fluxes', shape=(n_bins,n_bins), dtype=numpy.float64)
            unc_flux_ds = weed_iter_group.create_dataset('unc_fluxes', shape=(n_bins,n_bins), dtype=numpy.float64)
            avg_rates_ds = weed_iter_group.create_dataset('avg_rates', shape=(n_bins,n_bins), dtype=numpy.float64)
            unc_rates_ds = weed_iter_group.create_dataset('unc_rates', shape=(n_bins,n_bins), dtype=numpy.float64)

        averager = self.get_rates(n_iter, mapper)

        with self.data_manager.flushing_lock():
            avg_populations_ds[...] = averager.average_populations
            unc_populations_ds[...] = averager.stderr_populations
            avg_flux_ds[...] = averager.average_flux
            unc_flux_ds[...] = averager.stderr_flux
            avg_rates_ds[...] = averager.average_rate
            unc_rates_ds[...] = averager.stderr_rate

            binprobs = numpy.fromiter(map(operator.attrgetter('weight'),bins), dtype=numpy.float64, count=n_bins)
            orig_binprobs = binprobs.copy()

        westpa.rc.pstatus('Calculating equilibrium reweighting using window size of {:d}'.format(self.eff_windowsize))
        westpa.rc.pstatus('\nBin probabilities prior to reweighting:\n{!s}'.format(binprobs))
        westpa.rc.pflush()

        probAdjustEquil(binprobs, averager.average_rate, averager.stderr_rate)

        # Check to see if reweighting has set non-zero bins to zero probability (should never happen)
        assert (~((orig_binprobs > 0) & (binprobs == 0))).all(), 'populated bin reweighted to zero probability'

        # Check to see if reweighting has set zero bins to nonzero probability (may happen)
        z2nz_mask = (orig_binprobs == 0) & (binprobs > 0)
        if (z2nz_mask).any():
            westpa.rc.pstatus('Reweighting would assign nonzero probability to an empty bin; not reweighting this iteration.')
            westpa.rc.pstatus('Empty bins assigned nonzero probability: {!s}.'
                                .format(numpy.array_str(numpy.arange(n_bins)[z2nz_mask])))
        else:
            westpa.rc.pstatus('\nBin populations after reweighting:\n{!s}'.format(binprobs))
            for (bin, newprob) in zip(bins, binprobs):
                bin.reweight(newprob)

            weed_global_group.attrs['last_reweighting'] = n_iter

        assert (abs(1 - numpy.fromiter(map(operator.attrgetter('weight'),bins), dtype=numpy.float64, count=n_bins).sum())
                < EPS * numpy.fromiter(map(len,bins), dtype=numpy.int, count=n_bins).sum())

        westpa.rc.pflush()

max_issues_repo_name: burntyellow/westpa

max_forks_repo_name: burntyellow/westpa


text:
from __future__ import absolute_import, division, print_function, unicode_literals

import logging
import numpy as np
from itertools import product

from ..utils.various import mdot, less_logging, math_commands
from ..ml import ScoreEstimator, Ensemble, load_estimator
from ..utils.histo import Histo
from ..sampling import SampleAugmenter
from .. import sampling
from .base import BaseLikelihood

logger = logging.getLogger(__name__)

class HistoLikelihood(BaseLikelihood):
    def create_negative_log_likelihood(
        self,
        x_observed,
        observables=None,
        score_components=None,
        n_observed=None,
        x_observed_weights=None,
        include_xsec=True,
        luminosity=300000.0,
        mode="sampled",
        n_histo_toys=100000,
        model_file=None,
        hist_bins=None,
        thetas_binning=None,
        test_split=None,
    ):
        """
        Returns a function which calculates the negative log likelihood for a given
        parameter point, evaulated with a dataset (x_observed,n_observed,x_observed_weights).

        Parameters
        ----------

        x_observed : list of ndarray
            Set of event observables with shape `(n_events, n_observables)`.

        observables : list of str or None , optional
            Kinematic variables used in the histograms. The names are the same as
            used for instance in `DelphesReader`.

        score_components : None or list of int, optional
            Defines the score components used. Default value: None.

        n_observed : int or None , optional
            If int, number of observed events. If None, n_observed is defined by
            the length of x_observed. Default: None.

        x_observed_weights : list of float or None , optional
            List of event weights with shape `(n_events)`. If None, all events have equal
            weights. Default: None.

        include_xsec : bool, optional
            Whether the Poisson likelihood representing the total number of events is
            included in the analysis. Default value: True.

        luminosity : float, optional
            Integrated luminosity in pb^{-1} assumed in the analysis. Default value: 300000.

        mode : {"weighted" , "sampled", "histo"} , optional
            If "sampled", for each evaulation of the likelihood function, a separate
            set of events are sampled and histogram is created to construct the
            likelihood function. If "weighted", first a set of weighted events is
            sampled which is then used to create histograms. Default value: "sampled"

        n_histo_toys : int or None, optional
            Number of events drawn to construct the histograms used. If None and weighted_histo
            is True, all events in the training fraction of the MadMiner file are used. If None
            and weighted_histo is False, 100000 events are used. Default value: 100000.

        model_file : str or None, optional
            Filename of a saved neural network estimating the score. Required if
            score_components is not None. Default value: None.

        hist_bins : int or list of (int or ndarray) or None, optional
            Defines the histogram binning. If int, gives the number of bins automatically
            chosen for each summary statistic. If list, each entry corresponds to one
            summary statistic (e.g. kinematic variable specified by hist_vars); an int
            entry corresponds to the number of automatically chosen bins, an ndarray specifies
            the bin edges along this dimension explicitly. If None, the bins are chosen according
            to the defaults: for one summary statistic the default is 25 bins, for 2 it's 8 bins
            along each direction, for more it's 5 per dimension. Default value: None.

        thetas_binning : list of ndarray or None
            Specifies the parameter points used to determine the optimal binning.
            This is requires if hist_bins doesn't already fully specify the
            binning of the histogram. Default value : None

        test_split :

        Returns
        -------
        negative_log_likelihood : likelihood
            Function that evaluates the negative log likelihood for a given parameter point

        """

        # Check input and join observables and score components - nothing interesting
        if observables is None:
            observables = list([])
        if score_components is None:
            score_components = list([])
        if observables == [] and score_components == []:
            logger.info("No observables and scores provided. Calculate LLR due to rate and set include_xsec=True.")
            include_xsec = True
        observables = list(observables) + list(score_components)

        if n_observed is None:
            n_observed = len(x_observed)

        supported_modes = ["sampled", "weighted", "histo"]
        if mode not in supported_modes:
            raise ValueError("Mode %s unknown. Choose one of the following methods: %s", mode, supported_modes)

        # Load model - nothing interesting
        if score_components != []:
            assert all([isinstance(score_component, int) for score_component in score_components])
            if model_file is None:
                raise ValueError("You need to provide a model_file!")
            model = load_estimator(model_file)
        else:
            model = None

        # Create summary function
        logger.info("Setting up standard summary statistics")
        summary_function = None
        if observables != []:
            summary_function = self._make_summary_statistic_function(observables=observables, model=model)

        # Weighted sampled
        data, summary_stats, weights_benchmarks = None, None, None
        if mode == "weighted" or mode == "histo":
            logger.info("Getting weighted data")
            data, weights_benchmarks = self._make_histo_data_weighted(
                summary_function=summary_function, n_toys=n_histo_toys, test_split=test_split
            )
        if (mode == "weighted" or mode == "histo") and observables != []:
            summary_stats = summary_function(x_observed)

        # find binning
        logger.info("Setting up binning")
        if observables != [] and (
            hist_bins is None or not all([hasattr(hist_bin, "__len__") for hist_bin in hist_bins])
        ):
            if thetas_binning is None:
                raise ValueError("Your input requires adaptive binning: thetas_binning can not be None.")
            hist_bins = self._find_bins(hist_bins=hist_bins, n_summary_stats=len(observables))
            hist_bins = self._fixed_adaptive_binning(
                thetas_binning=thetas_binning,
                x_bins=hist_bins,
                data=data,
                weights_benchmarks=weights_benchmarks,
                n_toys=n_histo_toys,
                summary_function=summary_function,
            )
        logger.info("Use binning: %s", hist_bins)

        if mode == "histo":
            if hist_bins is not None:
                benchmark_histograms, _ = self._get_benchmark_histograms(data, weights_benchmarks, hist_bins)
                total_weights = np.array(
                    [sum(benchmark_histogram.flatten()) for benchmark_histogram in benchmark_histograms]
                )
            else:
                benchmark_histograms = None
                total_weights = np.array([sum(weights_benchmark) for weights_benchmark in weights_benchmarks.T])
        else:
            total_weights, benchmark_histograms = None, None

        # define negative likelihood function
        def nll(params):
            # Just return the expected Length
            if params is None:
                return self.n_nuisance_parameters + self.n_parameters

            # Process input
            if len(params) != self.n_nuisance_parameters + self.n_parameters:
                logger.warning(
                    "Number of parameters is %s, expected %s physical parameters and %s nuisance paramaters",
                    len(params),
                    self.n_parameters,
                    self.n_nuisance_parameters,
                )
            theta = params[: self.n_parameters]
            nu = params[self.n_parameters :]
            if len(nu) == 0:
                nu = None

            # Compute Log Likelihood
            log_likelihood = self._log_likelihood(
                n_events=n_observed,
                xs=x_observed,
                theta=theta,
                nu=nu,
                include_xsec=include_xsec,
                luminosity=luminosity,
                x_weights=x_observed_weights,
                mode=mode,
                n_histo_toys=n_histo_toys,
                hist_bins=hist_bins,
                summary_function=summary_function,
                data=data,
                summary_stats=summary_stats,
                weights_benchmarks=weights_benchmarks,
                benchmark_histograms=benchmark_histograms,
                total_weights=total_weights,
            )
            return -log_likelihood

        return nll

    def create_expected_negative_log_likelihood(
        self,
        theta_true,
        nu_true,
        observables=None,
        score_components=None,
        include_xsec=True,
        luminosity=300000.0,
        n_asimov=None,
        mode="sampled",
        n_histo_toys=100000,
        model_file=None,
        hist_bins=None,
        thetas_binning=None,
        test_split=None,
    ):
        """
        Returns a function which calculates the expected negative log likelihood for a given
        parameter point, evaulated with test data sampled according to theta_true.

        Parameters
        ----------
        theta_true : ndarray
            Specifies the physical paramaters according to which the test data is sampled.

        nu_true : ndarray
            Specifies the nuisance paramaters according to which the test data is sampled.

        observables : list of str or None , optional
            Kinematic variables used in the histograms. The names are the same as
            used for instance in `DelphesReader`.

        score_components : None or list of int, optional
            Defines the score components used. Default value: None.

        include_xsec : bool, optional
            Whether the Poisson likelihood representing the total number of events is
            included in the analysis. Default value: True.

        luminosity : float, optional
            Integrated luminosity in pb^{-1} assumed in the analysis. Default value: 300000.

        n_asimov : int or None, optional
            Size of the Asimov sample. If None, all weighted events in the MadMiner
            file are used. Default value: None.

        mode : {"weighted" , "sampled"} , optional
            If "sampled", for each evaulation of the likelihood function, a separate
            set of events are sampled and histogram is created to construct the
            likelihood function. If "weighted", first a set of weighted events is
            sampled which is then used to create histograms. Default value: "sampled"

        n_histo_toys : int or None, optional
            Number of events drawn to construct the histograms used. If None and weighted_histo
            is True, all events in the training fraction of the MadMiner file are used. If None
            and weighted_histo is False, 100000 events are used. Default value: 100000.

        model_file : str or None, optional
            Filename of a saved neural network estimating the score. Required if
            score_components is not None. Default value: None.

        hist_bins : int or list of (int or ndarray) or None, optional
            Defines the histogram binning. If int, gives the number of bins automatically
            chosen for each summary statistic. If list, each entry corresponds to one
            summary statistic (e.g. kinematic variable specified by hist_vars); an int
            entry corresponds to the number of automatically chosen bins, an ndarray specifies
            the bin edges along this dimension explicitly. If None, the bins are chosen according
            to the defaults: for one summary statistic the default is 25 bins, for 2 it's 8 bins
            along each direction, for more it's 5 per dimension. Default value: None.

        thetas_binning : list of ndarray or None
            Specifies the parameter points used to determine the optimal binning.
            If none, theta_true will be used. Default value : None

        test_split :

        Returns
        -------
        negative_log_likelihood : likelihood
            Function that evaluates the negative log likelihood for a given parameter point

        """

        if thetas_binning is None:
            thetas_binning = [theta_true]

        x_asimov, x_weights = self._asimov_data(theta_true, n_asimov=n_asimov)
        n_observed = luminosity * self.xsecs([theta_true], [nu_true])[0]

        return self.create_negative_log_likelihood(
            observables=observables,
            score_components=score_components,
            x_observed=x_asimov,
            n_observed=n_observed,
            x_observed_weights=x_weights,
            include_xsec=include_xsec,
            luminosity=luminosity,
            mode=mode,
            n_histo_toys=n_histo_toys,
            model_file=model_file,
            hist_bins=hist_bins,
            thetas_binning=thetas_binning,
        )

    def _log_likelihood(
        self,
        n_events,
        xs,
        theta,
        nu,
        include_xsec=True,
        luminosity=300000.0,
        x_weights=None,
        mode="sampled",
        n_histo_toys=100000,
        hist_bins=None,
        summary_function=None,
        data=None,
        summary_stats=None,
        weights_benchmarks=None,
        benchmark_histograms=None,
        total_weights=None,
    ):
        """
        Low-level function which calculates the value of the log-likelihood ratio.
        See create_negative_log_likelihood for options.
        """

        log_likelihood = 0.0
        if include_xsec:
            log_likelihood = log_likelihood + self._log_likelihood_poisson(
                n_events, theta, nu, luminosity, weights_benchmarks, total_weights
            )

        if summary_function is not None:
            if x_weights is None:
                x_weights = n_events / float(len(xs)) * np.ones(len(xs))
            else:
                x_weights = x_weights * n_events / np.sum(x_weights)
            log_likelihood_events = self._log_likelihood_kinematic(
                xs,
                theta,
                nu,
                mode,
                n_histo_toys,
                hist_bins,
                summary_function,
                data,
                summary_stats,
                weights_benchmarks,
                benchmark_histograms,
            )
            log_likelihood = log_likelihood + np.dot(x_weights, log_likelihood_events)

        if nu is not None:
            log_likelihood = log_likelihood + self._log_likelihood_constraint(nu)

        logger.debug("Total log likelihood: %s", log_likelihood)
        return log_likelihood

    def _log_likelihood_kinematic(
        self,
        xs,
        theta,
        nu,
        mode="sampled",
        n_histo_toys=100000,
        hist_bins=None,
        summary_function=None,
        data=None,
        summary_stats=None,
        weights_benchmarks=None,
        benchmark_histograms=None,
    ):
        """
        Low-level function which calculates the value of the kinematic part of the
        log-likelihood. See create_negative_log_likelihood for options.
        """
        # shape of theta
        if nu is not None:
            theta = np.concatenate((theta, nu), axis=0)

        # Calculate summary statistics
        if summary_stats is None:
            summary_stats = summary_function(xs)

        # Make histograms
        if mode == "sampled":
            data = self._make_histo_data_sampled(
                summary_function=summary_function, theta=theta, n_histo_toys=n_histo_toys
            )
            histo = Histo(data, weights=None, bins=hist_bins, epsilon=1.0e-12)
        elif mode == "weighted":
            weights = self._weights([theta], [nu], weights_benchmarks)[0]
            histo = Histo(data, weights=weights, bins=hist_bins, epsilon=1.0e-12)
        elif mode == "histo":
            bin_centers = [np.array([(bins[i] + bins[i + 1]) / 2 for i in range(len(bins) - 1)]) for bins in hist_bins]
            bin_centers = np.array(list(product(*bin_centers)))
            histo = self._histogram_morphing(theta, benchmark_histograms, hist_bins, bin_centers)

        # calculate log-likelihood from histogram
        log_p = histo.log_likelihood(summary_stats)

        return log_p

    def _make_summary_statistic_function(self, observables=None, model=None):
        """
        Low-level function that returns a function "summary_function" which
        evaluates the summary statistic for an event.
        """
        variables = math_commands()
        x_indices = self._find_x_indices(observables)

        def summary_function(xs):
            # only prefined observables - very fast
            if not "score" in x_indices and not "function" in x_indices:
                return xs[:, x_indices]

            # evaulate some observables using eval() - more slow
            data_events = []
            for x in xs:
                data_event = []
                if "function" in x_indices:
                    for observable_name, observable_value in zip(self.observables, x):
                        variables[observable_name] = observable_value
                if "score" in x_indices:
                    if isinstance(model, ScoreEstimator):
                        score = model.evaluate_score(x=np.array([x]))[0]
                    elif isinstance(model, Ensemble) and model.estimator_type == "score":
                        score, _ = model.evaluate_score(x=np.array([x]), calculate_covariance=False)[0]
                    else:
                        raise RuntimeError("Model has to be 'ScoreEstimator' or Ensemble thereof.")

                for observable, x_indice in zip(observables, x_indices):
                    if x_indice == "function":
                        data_event.append(float(eval(observable, variables)))
                    elif x_indice == "score":
                        data_event.append(score[observable])
                    else:
                        data_event.append(x[x_indice])
                data_events.append(data_event)
            return np.array(data_events)

        return summary_function

    def _find_x_indices(self, observables):
        """
        Low-level function that finds the indices corresponding to the observables
        and returns them as a list.
        """
        x_names = list(self.observables.keys())
        x_indices = []
        for obs in observables:
            if isinstance(obs, int):
                x_indices.append("score")
            else:
                try:
                    x_indices.append(x_names.index(obs))
                except:
                    x_indices.append("function")
        logger.debug("Using x indices %s", x_indices)
        return x_indices

    def _make_histo_data_sampled(self, summary_function, theta, n_histo_toys=1000):
        """
        Low-level function that creates histogram data sampled from one benchmark
        """
        # Get unweighted events
        with less_logging():
            sampler = SampleAugmenter(self.madminer_filename, include_nuisance_parameters=True)
            x, theta, _ = sampler.sample_train_plain(
                theta=sampling.morphing_point(theta),
                n_samples=n_histo_toys,
                test_split=False,
                filename=None,
                folder=None,
            )

        # Calculate summary stats
        data = summary_function(x)

        return data

    def _make_histo_data_weighted(self, summary_function, n_toys, test_split=None):
        """
        Low-level function that creates weighted histogram data
        """
        # Get weighted events
        start_event, end_event, correction_factor = self._train_test_split(True, test_split)
        x, weights_benchmarks = self.weighted_events(start_event=start_event, end_event=end_event, n_draws=n_toys)
        weights_benchmarks *= self.n_samples / n_toys

        # Calculate summary stats
        if summary_function is not None:
            data = summary_function(x)
        else:
            data = None

        return data, weights_benchmarks

    def _find_bins(self, hist_bins, n_summary_stats):
        """
        Low-level function that sets up the binning of the histograms (I)
        """
        if hist_bins is None:
            if n_summary_stats == 1:
                hist_bins = [25]
            elif n_summary_stats == 2:
                hist_bins = [8, 8]
            else:
                hist_bins = [5 for _ in range(n_summary_stats)]
        elif isinstance(hist_bins, int):
            # hist_bins = tuple([hist_bins] * n_summary_stats)
            hist_bins = [hist_bins for _ in range(n_summary_stats)]
        return hist_bins

    def _fixed_adaptive_binning(
        self,
        thetas_binning,
        x_bins,
        data=None,
        weights_benchmarks=None,
        n_toys=None,
        test_split=None,
        summary_function=None,
    ):
        """
        Low-level function that sets up the binning of the histograms (II)
        """

        # Get weighted data
        if data is None:
            data, weights_benchmarks = self._make_histo_data_weighted(
                summary_function=summary_function, n_toys=n_toys, test_split=test_split
            )

        # Calculate weights for thetas
        weights = self._weights(thetas_binning, None, weights_benchmarks)
        weights = np.asarray(weights)
        weights = np.mean(weights, axis=0)

        # Histogram
        histo = Histo(data, weights, x_bins, epsilon=1.0e-12)
        x_bins = histo.edges
        return x_bins

    def _get_benchmark_histograms(self, data, weights_benchmarks, hist_bins, epsilon=1.0e-12):
        """
        Low-level function that returns histograms for morphing benchmarks
        """
        # get histogram bins
        histo_benchmarks = []
        for weights in weights_benchmarks.T:
            # histo = Histo(data, weights=weights, bins=hist_bins, epsilon=1.0e-12)
            ranges = [(edges[0], edges[-1]) for edges in hist_bins]
            histo, _ = np.histogramdd(data, bins=hist_bins, range=ranges, normed=False, weights=weights)
            histo[:] += epsilon
            histo_benchmarks.append(histo)

        # get bin centers
        bin_centers = [np.array([(bins[i] + bins[i + 1]) / 2 for i in range(len(bins) - 1)]) for bins in hist_bins]
        bin_centers = np.array(list(product(*bin_centers)))
        return histo_benchmarks, bin_centers

    def _histogram_morphing(self, theta, histogram_benchmarks, hist_bins, bin_centers):
        """
        Low-level function that morphes histograms
        """
        # get binning
        hist_nbins = [len(bins) - 1 for bins in hist_bins]

        # get array of flattened histograms
        flattened_histo_weights = []
        for histo in histogram_benchmarks:
            flattened_histo_weights.append(histo.flatten())
        flattened_histo_weights = np.array(flattened_histo_weights).T

        # calculate dot product
        theta_matrix = self._get_theta_benchmark_matrix(theta)
        histo_weights_theta = mdot(theta_matrix, flattened_histo_weights)

        # create histogram
        histo = Histo(bin_centers, weights=histo_weights_theta, bins=hist_bins, epsilon=1.0e-12)
        return histo

    def _clean_nans(self, array):
        not_finite = np.any(~np.isfinite(array), axis=0)
        if np.sum(not_finite) > 0:
            logger.warning("Removing %s inf / nan results from calculation")
            array[:, not_finite] = 0.0
        return array

max_stars_repo_name: irenecortinovis/madminer

max_issues_repo_name: irenecortinovis/madminer

max_forks_repo_name: irenecortinovis/madminer


text:
import numpy as np
from qcodes import VisaInstrument, validators as vals
from qcodes.utils.validators import Numbers
from qcodes.utils.helpers import create_on_off_val_mapping

from qcodes.utils.deprecate import deprecate_moved_to_qcd

def parse_on_off(stat):
    if stat.startswith('0'):
        stat = 'Off'
    elif stat.startswith('1'):
        stat = 'On'
    return stat

frequency_mode_docstring = """
This command sets the frequency mode of the signal generator.

*FIXed* and *CW* These choices are synonymous. Any currently running frequency sweeps are
turned off, and the current CW frequency settings are used to control the output
frequency.

*SWEep* The effects of this choice are determined by the sweep generation type selected. In analog sweep generation, the ramp
sweep frequency settings (start, stop, center, and span) control the output
frequency. In step sweep generation, the current step sweep frequency settings
control the output frequency. In both cases, this selection also activates the sweep.
This choice is available with Option 007 only.

*LIST* This choice selects the swept frequency mode. If sweep triggering is set to
immediate along with continuous sweep mode, executing the command starts the LIST or STEP frequency sweep.
"""

IQsource_docstring = """
This command selects the I/Q modulator source for one of the two possible paths.

*EXTernal* This choice selects an external 50 ohm source as the I/Q input to I/Q modulator.
*INTernal* This choice is for backward compatibility with ESG E44xxB models and performs
the same function as the BBG1 selection.
*BBG1* This choice selects the baseband generator as the source for the I/Q modulator.
*EXT600* This choice selects a 600 ohm impedance for the I and Q input connectors and
routes the applied signals to the I/Q modulator.
*OFF* This choice disables the I/Q input.
"""

deprecate_moved_to_qcd(alternative="qcodes_contrib_drivers.drivers.Keysight.Keysight_E8267D.Keysight_E8267D")
class Keysight_E8267D(VisaInstrument):
    """
    This is the qcodes driver for the Keysight_E8267D signal generator

    Status: beta-version.
        TODO:
        - Add all parameters that are in the manual

    This driver will most likely work for multiple Agilent sources.

    This driver does not contain all commands available for the E8527D but
    only the ones most commonly used.
    """

    def __init__(self, name, address, step_attenuator=False, **kwargs):
        super().__init__(name, address, **kwargs)

        # Only listed most common spellings idealy want a
        # .upper val for Enum or string
        on_off_validator = vals.Enum('on', 'On', 'ON',
                                     'off', 'Off', 'OFF')
        on_off_mapping = create_on_off_val_mapping(1, 0)

        self.add_parameter(name='frequency',
                           label='Frequency',
                           unit='Hz',
                           get_cmd='FREQ:CW?',
                           set_cmd='FREQ:CW' + ' {:.4f}',
                           get_parser=float,
                           set_parser=float,
                           vals=vals.Numbers(1e5, 20e9),
                           docstring='Adjust the RF output frequency')
        self.add_parameter(name='frequency_offset',
                           label='Frequency offset',
                           unit='Hz',
                           get_cmd='FREQ:OFFS?',
                           set_cmd='FREQ:OFFS {}',
                           get_parser=float,
                           vals=Numbers(min_value=-200e9,
                                        max_value=200e9))
        self.add_parameter('frequency_mode',
                           label='Frequency mode',
                           set_cmd='FREQ:MODE {}',
                           get_cmd='FREQ:MODE?',
                           get_parser=lambda s: s.strip(),
                           vals=vals.Enum('FIX', 'CW', 'SWE', 'LIST'))
        self.add_parameter(name='phase',
                           label='Phase',
                           unit='deg',
                           get_cmd='PHASE?',
                           set_cmd='PHASE' + ' {:.8f}',
                           get_parser=self.rad_to_deg,
                           set_parser=self.deg_to_rad,
                           vals=vals.Numbers(-180, 180))
        self.add_parameter(name='power',
                           label='Power',
                           unit='dBm',
                           get_cmd='POW:AMPL?',
                           set_cmd='POW:AMPL' + ' {:.4f}',
                           get_parser=float,
                           set_parser=float,
                           vals=vals.Numbers(-130, 25))
        # self.add_parameter('status',
        self.add_parameter('rf_output', ### MLu edit 2020-04-07 Instrument driver looks for key 'rf_output'
                           get_cmd=':OUTP?',
                           set_cmd='OUTP {}',
                           get_parser=parse_on_off,
                           vals=on_off_validator)
        self.add_parameter(name='modulation_rf_enabled',
                           get_cmd='OUTP:MOD?',
                           set_cmd='OUTP:MOD {}',
                           val_mapping={'on': 1, 'off': 0})
                           # val_mapping=on_off_mapping)
        self.add_parameter('IQmodulator_enabled',
                           get_cmd='DM:STATe?',
                           set_cmd='DM:STATe {}',
                           val_mapping=on_off_mapping,
                           docstring='Enables or disables the internal I/Q modulator. Source can be external or internal.')

        for source in [1, 2]:
            self.add_parameter(f'IQsource{source}',
                               get_cmd=f'DM:SOUR{source}?',
                               set_cmd=f'DM:SOUR{source} {{}}',
                               get_parser=lambda s: s.strip(),
                               vals=vals.Enum('OFF', 'EXT', 'EXT600', 'INT'),
                               docstring=IQsource_docstring)

        self.add_parameter(f'IQadjustments_enabled', get_cmd=f'DM:IQAD?', set_cmd=f'DM:IQAD {{}}', val_mapping=on_off_mapping, docstring='Enable or disable IQ adjustments')

        IQoffset_parameters = dict(get_parser=float, set_parser=float, vals=vals.Numbers(-100,100))
        self.add_parameter(f'I_offset', get_cmd=f'DM:IQAD:IOFF?', set_cmd=f'DM:IQAD:IOFF {{}}', **IQoffset_parameters, docstring='I channel offset in percentage')
        self.add_parameter(f'Q_offset', get_cmd=f'DM:IQAD:QOFF?', set_cmd=f'DM:IQAD:QOFF {{}}',  **IQoffset_parameters, docstring='Q channel offset in percentage')
        self.add_parameter(f'IQ_quadrature', get_cmd=f'DM:IQAD:QSK?', set_cmd=f'DM:IQAD:QSK {{}}', get_parser=float, set_parser=float, docstring='IQ quadrature offset', unit='deg')

        self.add_parameter(f'pulse_modulation_enabled', get_cmd=f'PULM:STATe?', set_cmd=f'PULM:STATe {{}}', val_mapping=on_off_mapping, docstring='Enable or disable pulse modulation path')
        self.add_parameter(f'pulse_modulation_source', get_cmd=f'PULM:SOURce?', set_cmd=f'PULM:SOURce {{}}', get_parser=lambda s: s.strip(), vals=vals.Enum('EXT', 'INT', 'SCAL'))

        self.add_parameter(f'wideband_amplitude_modulation_enabled', get_cmd=f'AM:WID:STATe?', set_cmd=f'AM:WID:STATe {{}}', val_mapping=on_off_mapping, docstring='This command enables or disables wideband amplitude modulation')

        self.connect_message()

    def on(self):
        self.set('status', 'on')

    def off(self):
        self.set('status', 'off')

    @staticmethod
    def deg_to_rad(angle_deg):
        return np.deg2rad(float(angle_deg))

    @staticmethod
    def rad_to_deg(angle_rad):
        return np.rad2deg(float(angle_rad))

max_stars_repo_name: qnl/Qcodes

max_issues_repo_name: qnl/Qcodes

max_forks_repo_name: qnl/Qcodes


text:
import shutil

import progressbar
import scipy.stats as stats

from alexnet_data import *

def generate_two_correlated_time_series(size, rho):
    num_samples = size
    num_variables = 2
    cov = [[1.0, rho], [rho, 1.0]]

    L = np.linalg.cholesky(cov)

    uncorrelated = np.random.standard_normal((num_variables, num_samples))
    correlated = np.dot(L, uncorrelated)
    x, y = correlated
    rho, p_val = stats.pearsonr(x, y)
    return x, y, rho

def generate():
    shutil.rmtree(DATA_FOLDER, ignore_errors=True)
    size = 1024
    total_num_images = 100000
    bar = progressbar.ProgressBar()
    for i in bar(range(total_num_images)):
        if i % 2 == 0:
            class_name = 'CORRELATED'
            rho = 0.8
        else:
            class_name = 'UNCORRELATED'
            rho = 0.0

        # more fun to consider non-stationary time series.
        # still the correlation holds almost surely.
        x_tr, y_tr, _ = generate_two_correlated_time_series(size, rho)
        x_tr = np.cumsum(x_tr)
        y_tr = np.cumsum(y_tr)

        x_te, y_te, _ = generate_two_correlated_time_series(size, rho)
        x_te = np.cumsum(x_te)
        y_te = np.cumsum(y_te)

        train_output_dir = os.path.join(DATA_FOLDER, 'train', class_name)
        mkdir_p(train_output_dir)
        generate_multi_time_series([x_tr, y_tr], os.path.join(train_output_dir, '{}.png'.format(i)))

        test_output_dir = os.path.join(DATA_FOLDER, 'test', class_name)
        mkdir_p(test_output_dir)
        generate_multi_time_series([x_te, y_te], os.path.join(test_output_dir, '{}.png'.format(i)))

if __name__ == '__main__':
    generate()

max_stars_repo_name: philipperemy/tensorflow-cnn-time-series

max_issues_repo_name: afcarl/tensorflow-cnn-time-series

max_forks_repo_name: afcarl/tensorflow-cnn-time-series

reason: import scipy


text:
import decode

import cv2
import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt

from tensorflow import keras

def main():
    config = tf.compat.v1.ConfigProto(gpu_options=tf.compat.v1.GPUOptions(per_process_gpu_memory_fraction=0.8))
    config.gpu_options.allow_growth = True
    session = tf.compat.v1.Session(config=config)
    tf.compat.v1.keras.backend.set_session(session)

    classes_names = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', 'p', 'q', 'r',
                     's', 't', 'u', 'v', 'w', 'x', 'y', 'z']
    path = "./test/normal/q.png"  # real test image path

    model = keras.models.load_model("LeNet5_model")
    # print(model.summary())

    # test_data_validate(model, classes_names)
    real_write_test(model, classes_names, path)

def test_data_validate(model, classes_names):
    """checkout test images recognition result

    choose ten of them, show images in plt and result in console

    :param model: trained LeNet5 model
    :param classes_names: readable classes names
    :return: none
    """
    test_images = decode.decode_idx3_ubyte("raw_data_files/emnist-letters-train-images-idx3-ubyte", 50)
    test_images = test_images / 255.0
    for i in range(10):
        test_image = test_images[i]
        plt.figure()
        plt.imshow(test_image)
        plt.colorbar()
        test_image = test_image.reshape(28, 28, 1)
        test_image = np.expand_dims(test_image, 0)
        predictions = model.predict(test_image)
        plt.xlabel("This is " + classes_names[np.argmax(predictions[0])])
        plt.xticks([])
        plt.yticks([])
        plt.show(block=True)

def real_write_test(model, classes_names, path):
    """test the real write characters by hand

    :param model: trained LeNet5 model
    :param classes_names: readable classes names
    :param path: test image path
    :return: none
    """
    test_img = cv2.imread(path)

    # preprocess
    test_img = cv2.resize(test_img, (28, 28))
    plt.figure()
    plt.imshow(test_img)
    plt.colorbar()
    plt.yticks([])
    plt.xticks([])
    plt.show()
    test_img = cv2.cvtColor(test_img, cv2.COLOR_RGB2GRAY)
    plt.figure()
    plt.imshow(test_img)
    plt.xticks([])
    plt.yticks([])
    plt.colorbar()
    plt.show()
    test_img = np.array(test_img)

    test_img = 255 - test_img
    test_img = (test_img - np.min(test_img)) / (np.max(test_img) - np.min(test_img))

    # simple filter, enhance signal
    for i in range(28):
        for j in range(28):
            if test_img[i][j] < 0.2:
                test_img[i][j] = 0
            elif test_img[i][j] >= 0.8:
                test_img[i][j] = 1

    plt.figure()
    plt.imshow(test_img)
    test_img = test_img.reshape(28, 28, 1)
    test_img = np.expand_dims(test_img, 0)
    predictions = model.predict(test_img)
    plt.yticks([])
    plt.xticks([])
    plt.colorbar()
    plt.xlabel("This is " + classes_names[np.argmax(predictions[0])])
    plt.show()

if __name__ == '__main__':
    main()


text:
# Set Matplotlib to not require X server (for cron use)
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import argparse
import glob
import re
from collections import OrderedDict

# Import global variables
SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
sys.path.append(os.path.dirname(SCRIPT_DIR))
from global_variables import *

def get_sorted_progress_files(files):
    '''
    Sort files by progress log number. Assume files include basenames progress.*, progress_1.*, progress_2.*, ...,
    with no gaps, but may be out of order
    :param files: A list of progress files.
    :return:
    '''
    files_arr = np.array(files)
    order = []
    for file in files:
        # Check file name starts with "progress"
        basename = os.path.basename(file)
        m = re.match(r'progress(_\d+)*', basename)
        if m and not m.group(1):
            # match is 'progress'
            order.append(0)
        elif m and m.group(1):
            # match is 'progress_*', first group is _*
            order.append(int(m.group(1)[1:]))
    files_arr = files_arr[np.argsort(order)]
    return files_arr.tolist()

def lines_to_acc_loss(lines):
    '''
    Extract iteration numbers as well as accuracy and loss for three rotation angles
    :param lines: The lines containing the above information
    :return:
    '''
    iters = []
    az_acc = []
    el_acc = []
    ti_acc = []
    az_loss = []
    el_loss = []
    ti_loss = []
    for line in lines:
        # Get iterations
        m = re.search('Iteration ([0-9]+)', line)
        if m:
            iters.append(int(m.group(1)))
        # Get azimuths
        m = re.search('accuracy_azimuth = (([0-9]|\.)+)', line)
        if m:
            az_acc.append(float(m.group(1)))
        m = re.search('loss_azimuth = (([0-9]|\.)+)', line)
        if m:
            az_loss.append(float(m.group(1)))
        # Get elevations
        m = re.search('accuracy_elevation = (([0-9]|\.)+)', line)
        if m:
            el_acc.append(float(m.group(1)))
        m = re.search('loss_elevation = (([0-9]|\.)+)', line)
        if m:
            el_loss.append(float(m.group(1)))
        # Get tilts
        m = re.search('accuracy_tilt = (([0-9]|\.)+)', line)
        if m:
            ti_acc.append(float(m.group(1)))
        m = re.search('loss_tilt = (([0-9]|\.)+)', line)
        if m:
            ti_loss.append(float(m.group(1)))

    # Compute sum of accuracies and losses
    acc_sum = [np.sum(t) for t in zip(az_acc, el_acc, ti_acc)]
    loss_sum = [np.sum(t) for t in zip(az_loss, el_loss, ti_loss)]

    return zip(iters, az_acc, el_acc, ti_acc, az_loss, el_loss, ti_loss, acc_sum, loss_sum)

def train_info_to_tsv(save_path, tuples):
    '''
    Save the given list of training info tuples to a file in TSV format
    :param tuples: A list of tuples, where each tuple includes iteration number, accuracies, and losses
    :return:
    '''
    with open(save_path, 'w') as f:
        f.write('#Iters\tAcc_Azimuth\tAcc_Elevation\tAcc_Tilt\tLoss_Azimuth\tLoss_Elevation\tLoss_Tilt\tAcc_Sum\tLoss_Sum\n')
        for t in tuples:
            f.write('%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n' % t)

def plot_loss(ax, train_info):
    '''
    Generates a loss plot for the given data
    :param ax: The pyplot axes to draw on
    :param train_info: A list of tuples, where each tuple includes iteration number, accuracies, and losses
    :return:
    '''
    iters = [t[0] for t in train_info]
    az_loss = [t[4] for t in train_info]
    el_loss = [t[5] for t in train_info]
    ti_loss = [t[6] for t in train_info]
    loss_sum = [t[8] for t in train_info]

    ax.plot(iters, az_loss, color='g', label='Azimuth')
    ax.plot(iters, el_loss, color='r', label='Elevation')
    ax.plot(iters, ti_loss, color='c', label='Tilt')
    ax.plot(iters, loss_sum, color='b', label='Total')

    # Adjust graph so legend appears to the right of the plot
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax.set_title('Training loss vs. iterations')
    ax.set_xlabel('Number of iterations')
    ax.set_ylabel('Loss')
    ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
    ax.set_ylim([0, 120])

def plot_acc(ax, train_info):
    '''
    Generates an accuracy plot for the given data
    :param ax: The pyplot axes to draw on
    :param train_info: A list of tuples, where each tuple includes iteration number, accuracies, and losses
    :return:
    '''
    iters = [t[0] for t in train_info]
    az_acc = [t[1] for t in train_info]
    el_acc = [t[2] for t in train_info]
    ti_acc = [t[3] for t in train_info]
    acc_sum = [t[7] for t in train_info]

    ax.plot(iters, az_acc, color='g', label='Azimuth')
    ax.plot(iters, el_acc, color='r', label='Elevation')
    ax.plot(iters, ti_acc, color='c', label='Tilt')
    ax.plot(iters, acc_sum, color='b', label='Total')

    # Adjust graph so legend appears to the right of the plot
    box = ax.get_position()
    ax.set_position([box.x0, box.y0, box.width * 0.8, box.height])
    ax.set_title('Testing accuracy vs. iterations')
    ax.set_xlabel('Number of iterations')
    ax.set_ylabel('Sum of accuracies')
    ax.set_ylim([0, 3])

def main(exp_num):
    print('Creating training plots for experiment %d...' % exp_num)
    # Get experiment path
    exp_dir = glob.glob(os.path.join(g_experiments_root_folder, '%06d*' % exp_num))
    if len(exp_dir) == 0:
        print('Experiment %d does not exist, quitting' % exp_num)
        return
    exp_dir = exp_dir[0]

    # Get progress log paths
    progress_dir = os.path.join(exp_dir, 'progress')
    progress_log_names = get_sorted_progress_files(os.listdir(progress_dir))

    # Keep track of training info across progress logs, keeping latest versions of values
    all_test_info = OrderedDict()
    all_train_info = OrderedDict()
    for progress_log_name in progress_log_names:
        # Read Caffe log file
        progress_log_contents = open(os.path.join(progress_dir, progress_log_name), 'r').read()
        # Extract info from test lines
        test_lines = re.findall('.*Test.*', progress_log_contents, re.MULTILINE)
        test_info = lines_to_acc_loss(test_lines)
        # Store/overwrite info
        for tuple in test_info:
            all_test_info[tuple[0]] = tuple

        # Extract info from training lines
        train_lines = re.findall('.*Train.*', progress_log_contents, re.MULTILINE)
        train_lines += re.findall('Iteration [0-9]+, loss = ', progress_log_contents)
        train_info = lines_to_acc_loss(train_lines)
        for tuple in train_info:
            all_train_info[tuple[0]] = tuple

    # Convert info back to lists
    all_test_info = all_test_info.values()
    all_train_info = all_train_info.values()

    # Save info to TSV files
    train_info_to_tsv(os.path.join(progress_dir, 'progress_test.log'), all_test_info)
    train_info_to_tsv(os.path.join(progress_dir, 'progress_train.log'), all_train_info)

    # Draw progress plots
    plt.figure(figsize=(21, 6))
    ax = plt.subplot(131)
    plot_loss(ax, all_train_info)
    ax = plt.subplot(132)
    plot_loss(ax, all_test_info)
    ax = plt.subplot(133)
    plot_acc(ax, all_test_info)
    # Make global plot adjustments
    plt.suptitle('Experiment %d' % exp_num)
    plt.tight_layout(pad=3, w_pad=18)
    # Save plots to file
    plt.savefig(os.path.join(progress_dir, 'plots.png'))

    print('Done.')

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('exp_num', type=int, help='The experiment number')
    args = parser.parse_args()

    main(args.exp_num)

max_stars_repo_name: rszeto/click-here-cnn

max_issues_repo_name: rszeto/click-here-cnn

max_forks_repo_name: rszeto/click-here-cnn


text:
########################################
### IMPORT MODULES                   ###
########################################
import pandas as pd
from math import sqrt
# import matplotlib.pyplot as plt
import csv
from random import sample
import numpy as np
import googlemaps  # imports google maps library
from datetime import datetime
from geopy import distance

########################################
### SET CONFIGURATIONS               ###
########################################
gmaps = googlemaps.Client(key='####################################')
sdate = f"{datetime.now():%m_%d_%H_%M}"
cols = [
    'Orig_Zip', 'Dest_Zip', 'Orig_Lon', 'Orig_Lat', 'Dest_Lon', 'Dest_Lat', 'Transit'
]
with open('dataset/routes_{}.csv'.format(sdate), 'w', encoding='utf-8') as csv_file:
    writer = csv.writer(csv_file)
    writer.writerow(cols)

########################################
### LOAD DATA                        ###
########################################
data = pd.read_csv('dataset/London postcodes.csv')
data = data[data['In Use?'] == 'Yes']
keep = ['Postcode', 'Latitude', 'Longitude', 'Population']
df = data[keep]
df['Distance'] = np.zeros(len(df))
df = df[df.Population > 100]
df = df.reset_index()

########################################
### SELECT SAMPLE                    ###
########################################
# Compute distance from centre
for i in range(len(df)):
    df.loc[i, 'Distance'] = distance.distance((-0.115425, 51.512376),(df.loc[i, 'Longitude'],df.loc[i, 'Latitude'])).meters
    print(i / len(df))
cut = df[df['Distance'] < 10000]
cut = cut.reset_index()

########################################
### OBTAIN TRAVEL TIME FROM GOOGLE   ###
########################################
def draw_connection(df):
    orig_draw = sample(list(df.index), 1)[0]
    dest_draw = sample(list(df.index), 1)[0]
    return orig_draw, dest_draw

full_df = pd.DataFrame(index=range(29000), columns=cols)
for idx in range(25000):
    time_of_travel = datetime(2018, 8, 1, 9, 0)
    orig_draw, dest_draw = draw_connection(cut)
    full_df.loc[idx, 'Orig_Zip'] = cut.loc[orig_draw, 'Postcode']
    full_df.loc[idx, 'Dest_Zip'] = cut.loc[dest_draw, 'Postcode']
    full_df.loc[idx, 'Orig_Lon'] = cut.loc[orig_draw, 'Longitude']
    full_df.loc[idx, 'Orig_Lat'] = cut.loc[dest_draw, 'Latitude']
    full_df.loc[idx, 'Dest_Lon'] = cut.loc[orig_draw, 'Longitude']
    full_df.loc[idx, 'Dest_Lat'] = cut.loc[dest_draw, 'Latitude']
    try:
        transit_directions = gmaps.directions(
            (cut.loc[orig_draw, 'Latitude'], cut.loc[orig_draw, 'Longitude']),
            (cut.loc[dest_draw, 'Latitude'], cut.loc[dest_draw, 'Longitude']),
            mode='transit',
            departure_time=time_of_travel,
        )
        travel_time_atob = transit_directions[0].get('legs')[0].get('duration').get(
            'value'
        )
        full_df.loc[idx, 'Transit'] = travel_time_atob
        with open(
            'dataset/routes_{}.csv'.format(sdate), 'a', encoding='utf-8'
        ) as csv_file:
            writer = csv.writer(csv_file)
            writer.writerow(full_df.iloc[idx])
    except BaseException as e:
        print(e)
    print(idx / 32000)

########################################
### SAVE DUPLICATE AS EXCEL          ###
########################################
writer = pd.ExcelWriter('dataset/routelist_{}.xlsx'.format(sdate))
full_df.to_excel(writer, 'Sheet1', index=False)
writer.save()

max_stars_repo_name: wbglaeser/travel-times-london

max_issues_repo_name: wbglaeser/travel-times-london

max_forks_repo_name: wbglaeser/travel-times-london