text:
import warnings
from functools import lru_cache
from typing import Optional, List, Dict, Union

import numpy as np
import pandas as pd
from nba_api.stats.endpoints import boxscoretraditionalv2, commonallplayers
from requests import ConnectionError

warnings.filterwarnings('ignore')

def convert_timestring_to_second(data: pd.DataFrame, column: str) -> pd.Series:
    """Converting a string of game time into a numeric of seconds from the start of the game.

    Args:
        data (pd.DataFrame): play-by-play data from stats.nba.com
        column (str): name of column with game time

    Returns:
        pd.Series: seconds from start of game
    """
    data.columns = [x.upper() for x in data.columns]
    s = data.loc[:, column].str.split(pat=":", expand=True).astype('int32')
    s.columns = ["MIN", "SEC"]
    timegame = pd.Series([np.abs((minute * 60 + sec) - 720 * period) if period < 5 \
                              else np.abs((minute * 60 + sec) - (2880 + 300 * (period - 4))) \
                          for (minute, sec, period) in zip(s["MIN"], s["SEC"], data["PERIOD"])])
    return timegame

def player_on_court(data: pd.DataFrame, all_id: Optional[np.ndarray] = None) -> np.ndarray:
    """Getting array of PLAYER_ID players who were on court at start of quarter.

    Args:
        data (pd.Dataframe): play-by-play data from stats.nba.com
        all_id (np.array): Array of PLAYER_ID of all players who appeared on court in quarter from nba boxscore.
        Used only when 10 players cannot be retrieved from play-by-play data. By default, None.

    Returns:
        np.array: PLAYER_ID of players on court at beginning of quarter.
    """
    data.columns = [x.upper() for x in data.columns]
    if all_id is None:
        pl1_id = data.loc[(~data["EVENTMSGTYPE"].isin([9, 18])) & (~data["PERSON1TYPE"].isin([6, 7])) & \
                          (~pd.isna(data["PLAYER1_NAME"])), "PLAYER1_ID"].unique()
        pl2_id = data.loc[(~data["EVENTMSGTYPE"].isin([9, 18])) & (~data["PERSON2TYPE"].isin([6, 7])) & \
                          (~pd.isna(data["PLAYER2_NAME"])), "PLAYER1_ID"].unique()
        pl3_id = data.loc[(~data["EVENTMSGTYPE"].isin([9, 18])) & (~data["PERSON3TYPE"].isin([6, 7])) & \
                          (~pd.isna(data["PLAYER3_NAME"])), "PLAYER3_ID"].unique()
        all_id = np.unique(np.concatenate((pl1_id, pl2_id, pl3_id)))
        all_id = all_id[(all_id != 0) & (all_id < 1610612737)]

    sub_off = data.loc[data["EVENTMSGTYPE"] == 8, "PLAYER1_ID"].unique()
    sub_on = data.loc[data["EVENTMSGTYPE"] == 8, "PLAYER2_ID"].unique()
    all_id = all_id[~np.in1d(all_id, sub_on[~np.in1d(sub_on, sub_off)])]
    sub_on_off = sub_on[np.in1d(sub_on, sub_off)]

    for i in sub_on_off:
        on = np.min(data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER2_ID"] == i), "PCTIMESTRING"])
        off = np.min(data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER1_ID"] == i), "PCTIMESTRING"])
        if off > on:
            all_id = all_id[~np.in1d(all_id, i)]
        elif off == on:
            on_event = np.min(data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER2_ID"] == i), "EVENTNUM"])
            off_event = np.min(data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER1_ID"] == i), "EVENTNUM"])
            if off_event > on_event:
                all_id = all_id[~np.in1d(all_id, i)]
    return all_id

def sort_players(data: pd.DataFrame, all_id: np.ndarray) -> List:
    """Sorting players on court by teams (first away team, then home team)

    Args:
        data (pd.Dataframe): play-by-play data from stats.nba.com
        all_id: unsorted by teams array of PLAYER_ID

    Returns:
`       list: sorted by teams (first away team, then home team) array of PLAYER_ID of players
        on court at beginning of quarter.
    """
    data.columns = [x.upper() for x in data.columns]
    pl1_id = data.loc[(data["PLAYER1_ID"] != 0) & (data["PERSON1TYPE"].isin([4, 5])), ["PLAYER1_ID", "PERSON1TYPE"]] \
        .rename(columns={"PLAYER1_ID": "PLAYER_ID", "PERSON1TYPE": "PERSONTYPE"})
    pl2_id = data.loc[(data["PLAYER2_ID"] != 0) & (data["PERSON2TYPE"].isin([4, 5])), ["PLAYER2_ID", "PERSON2TYPE"]] \
        .rename(columns={"PLAYER2_ID": "PLAYER_ID", "PERSON2TYPE": "PERSONTYPE"})
    pl3_id = data.loc[(data["PLAYER3_ID"] != 0) & (data["PERSON3TYPE"].isin([4, 5])), ["PLAYER3_ID", "PERSON3TYPE"]] \
        .rename(columns={"PLAYER3_ID": "PLAYER_ID", "PERSON3TYPE": "PERSONTYPE"})
    pl_id = np.array(pd.concat([pl1_id, pl2_id, pl3_id], axis=0, ignore_index=True).drop_duplicates() \
                     .sort_values(by="PERSONTYPE", ascending=False).loc[:, "PLAYER_ID"])
    pl_id = [x for x in pl_id if x in all_id]
    return pl_id

def fill_columns_data(data: pd.DataFrame, all_id: Union[np.ndarray, List]) -> pd.DataFrame:
    """Adding columns with PLAYER_ID on court

    Args:
        data (pd.DataFrame): play-by-play data from stats.nba.com
        all_id (np.array): sorted by teams (before away team, after home team) array of PLAYER_ID of players
        on court at beginning of quarter.

    Returns:
        pd.Dataframe: play-by-play data with PLAYER_ID on court
    """
    data.columns = [x.upper() for x in data.columns]
    for i in range(1, 11):
        data.loc[:, "PLAYER" + str(i)] = all_id[i - 1]
        subs = data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER1_ID"] == data["PLAYER" + str(i)])]
        nsub = subs.shape[0]
        while nsub > 0:
            n = subs.index[0]
            player_on = subs.iloc[0]["PLAYER2_ID"]
            data.loc[n:, "PLAYER" + str(i)] = player_on
            subs = data.loc[(data["EVENTMSGTYPE"] == 8) & (data["PLAYER1_ID"] == data["PLAYER" + str(i)])]
            nsub = subs.shape[0]
    return data

def adding_player_on_court(data: pd.DataFrame, **kwargs: Dict[str, float]) -> pd.DataFrame:
    """Adding players on court to play-by-play data from stats.nba.com

    Args:
        data (pd.DataFrame): play-by-play data from stats.nba.com
        **kwargs (Dict): Arbitrary keyword arguments.

    Returns:
        pd.DataFrame: play-by-play data with players on court
    """
    args = kwargs
    data.columns = [x.upper() for x in data.columns]
    if isinstance(data["PCTIMESTRING"][0], str):
        data["PCTIMESTRING"] = convert_timestring_to_second(data, "PCTIMESTRING")
    d = dict()
    for period in data["PERIOD"].unique():
        df = data.loc[data["PERIOD"] == period]
        all_id = player_on_court(df)
        if len(all_id) == 10:
            all_id = sort_players(df, all_id)
            d[period] = fill_columns_data(df, all_id)
        else:
            retry = 0
            bx = ""
            while retry < args.get("retry", 5):
                try:
                    bx = boxscoretraditionalv2.BoxScoreTraditionalV2(
                        game_id=args.get("game_id", "00" + str(df.GAME_ID.unique()[0])),
                        start_period=args.get("start_period", str(period)),
                        end_period=args.get("end_period", str(period)),
                        range_type=args.get("range_type", "1"),
                        timeout=args.get("timeout", 10)
                    )
                    break
                except ConnectionError as e:
                    bx = e
                    retry += 1
            if not isinstance(bx, boxscoretraditionalv2.BoxScoreTraditionalV2):
                raise ConnectionError(bx)
            player_stats = bx.player_stats.get_data_frame()
            all_id = np.array([player_stats.PLAYER_ID]).ravel()
            all_id = player_on_court(df, all_id)
            d[period] = fill_columns_data(df, all_id)

    return pd.concat(d, axis=0, ignore_index=True)

@lru_cache()
def _cache_player_info(**kwargs: Dict[str, float]) -> commonallplayers.CommonAllPlayers:
    """Getting data on all players in history of NBA

    Args:
        **kwargs (Dict): Arbitrary keyword arguments.

    Returns:
        commonallplayers.CommonAllPlayers: class nba_api with information about all players in NBA history
    """
    args = kwargs
    info = commonallplayers.CommonAllPlayers(
        is_only_current_season=args.get("is_only_current_season", 0),
        league_id=args.get("league_id", "00"),
        season=args.get("season", "2021-22"),
        timeout=args.get("timeout", 10)
    )
    return info

def replace_id_on_name(id_players: Union[pd.Series, np.ndarray, List],
                       player_data: Optional[pd.DataFrame] = None, **kwargs: Dict[str, float]) -> List:
    """Replacing player's PLAYER_ID with his full name

    Args:
        id_players: List of player IDs on court
        player_data: Data on all players in history of NBA. If None, data is requested from stats.nba.com.
            By default, None.
        **kwargs: Arbitrary keyword arguments.

    Returns:
        list: List of names of players on court
    """
    args = kwargs
    if player_data is None:
        retry = 0
        info = ""
        while retry < args.get("retry", 5):
            try:
                info = _cache_player_info()
                break
            except ConnectionError as e:
                info = e
                retry += 1
        if not isinstance(info, commonallplayers.CommonAllPlayers):
            raise ConnectionError(info)
        player_data = info.common_all_players.get_data_frame()

    player_data.columns = [x.upper() for x in player_data.columns]
    name_players = []
    for player in id_players:
        pl_id = player_data.loc[player_data["PERSON_ID"] == player, "DISPLAY_FIRST_LAST"].reset_index(drop=True)[0]
        name_players.append(pl_id)

    return name_players


text:
# Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
# SPDX-License-Identifier: Apache-2.0

import gc
import matplotlib.pyplot as plt
import multiprocessing as mp
from multiprocessing import Pool
import os
import pickle
import psutil
import time
import argparse

import gpflow
import gpflowopt
import numpy as np
import scipy as sp
import tensorflow as tf

from dpareto.utils.object_io import save_object
from dpareto.utils.random_seed_setter import set_random_seed

os.environ['TF_CPP_MIN_LOG_LEVEL'] = '2'  # disables TensorFlow debugging info dump that appears on Windows machines
#os.environ['CUDA_VISIBLE_DEVICES'] = '-1'  # forces TensorFlow to use CPU instead of defaulting to GPU

class HypervolumeImprovementHarness:
    def __init__(self, model_class, problem_name, initial_data_options, anti_ideal_point, optimization_domain, fixed_hyperparams, instance_options={}, num_instances=0, num_replications=1, num_workers=1, plot_options={}, output_base_dir=''):
        mp.set_start_method('spawn', force=True)  # workaround for weird matplotlib+multiprocessing compatibility issue

        self._model_class = model_class
        self._problem_name = problem_name
        self._num_instances = num_instances
        self._num_replications = num_replications

        # This is used as an ordering for consistently converting the hyperparams dict of
        # hyperparam-name-string -> hyperparam-value that we use into a numpy array of hyperparam
        # values that GPFlowOpt uses.
        self._optimizable_hyperparams_names = list(optimization_domain.keys())

        parser = argparse.ArgumentParser()
        parser.add_argument('--init-data-file')
        args = parser.parse_args()

        # Loads initial data into numpy arrays from file or generates it.
        if args.init_data_file is not None:
            initial_data_options["filename"] = args.init_data_file
        initial_data = self._get_initial_points(initial_data_options)
        self._initial_hyperparams_vals = initial_data[0]
        self._initial_privacy_vals = initial_data[1]
        self._initial_utility_vals = initial_data[2]
        self._max_initial_utility_variance = initial_data[3]

        # Transforms the initial data into a format that GPFlowOpt can use.
        transformed_data = self._create_initial_transformations(initial_data)
        self._transformed_initial_hyperparams_vals = transformed_data[0]
        self._transformed_initial_privacy_vals = transformed_data[1]
        self._transformed_initial_utility_vals = transformed_data[2]

        # Creates separate GPs for modeling privacy and utility.
        self._privacy_gp = self._create_privacy_gp()
        self._utility_gp = self._create_utility_gp()

        # Runs a hard-coded sanity check against data from the loaded file.
        self._print_sanity_check()

        # Sets the information for the reference (aka anti-ideal) point and its transformation.
        self._anti_ideal_point = anti_ideal_point
        transformed_anti_ideal_eps = self._transform_privacy(self._anti_ideal_point[0])
        transformed_anti_ideal_err = self._transform_utility(self._anti_ideal_point[1])
        self._transformed_anti_ideal_point = [transformed_anti_ideal_eps, transformed_anti_ideal_err]

        # Optimization domain dictates the range of each to-be-optimized-over hyperparamter. Fixed hyperparams
        # are the hyperparams that aren't optimized over.
        self._optimization_domain = optimization_domain
        self._gpflowopt_domain = self._construct_gpflowopt_domain(optimization_domain)
        self._fixed_hyperparams = fixed_hyperparams
        self._instance_options = instance_options

        # Create acquisition function, set it's anti-ideal point, and create it's optimizer.
        self._acquisition = gpflowopt.acquisition.HVProbabilityOfImprovement([self._privacy_gp, self._utility_gp])
        self._acquisition.reference = np.array([self._transformed_anti_ideal_point])
        self._acquisition_optimizer = gpflowopt.optim.StagedOptimizer([gpflowopt.optim.MCOptimizer(self._gpflowopt_domain, 5000),
                                                                       gpflowopt.optim.SciPyOptimizer(self._gpflowopt_domain)])

        # Create the Bayesian Optimizer.
        self._bayesian_optimizer = gpflowopt.BayesianOptimizer(self._gpflowopt_domain,
                                                               self._acquisition,
                                                               optimizer=self._acquisition_optimizer,
                                                               scaling=False)

        # Multiprocessing worker info
        self._num_workers = num_workers

        # Misc fields
        self._saved_hypervolumes = []
        self._all_results = []
        self._output_base_dir = output_base_dir
        self._run_name = int(time.time() * 1000)
        self._save_results = True
        self._make_results_directory()
        self._show_plot = plot_options.get('show_plot', False)
        self._plot_options = plot_options

    ##################################################################################
    # Initial data-related methods
    ##################################################################################
    @staticmethod
    def _process_replications(data):
        aggregated_points = {}
        for point in data:
            hyperparams_dict = point[0]
            frozen_hyperparams = frozenset(hyperparams_dict.items())
            if frozen_hyperparams not in aggregated_points:
                aggregated_points[frozen_hyperparams] = ([], [])
            aggregated_points[frozen_hyperparams][0].append(point[1])
            aggregated_points[frozen_hyperparams][1].append(point[2])

        processed_data = []
        seen_points = set()
        for point in data:
            hyperparams_dict = point[0]
            frozen_hyperparams = frozenset(hyperparams_dict.items())
            if frozen_hyperparams not in seen_points:
                seen_points.add(frozen_hyperparams)
                # compute privacy as max value
                privacy = np.max(aggregated_points[frozen_hyperparams][0])
                # compute utility as mean value
                utility = np.mean(aggregated_points[frozen_hyperparams][1])
                # add variance of utility values
                utility_variance = np.var(aggregated_points[frozen_hyperparams][1])
                # add everything back into the list of processed data
                processed_data.append((hyperparams_dict, privacy, utility, utility_variance))
        return processed_data

    def _get_initial_points(self, initial_data_options):
        num_initial_points = initial_data_options['num_points']

        # get initial data from file or generate it with random sampling
        if 'filename' in initial_data_options:
            filename = initial_data_options['filename']
            raw_data = self._load_points_from_file(filename)
        else:
            raise ValueError("Please provide path to the file full_results.pkl that contains initial data points from random sampling run.")

        # Process out the replications and extract the proper number of initial points
        processed_data = self._process_replications(raw_data)[:num_initial_points]

        # Set input hyperparams and output privacy & utility
        initial_hyperparams_dicts = np.array([instance[0] for instance in processed_data])

        initial_hyperparams_vals = np.array([np.array([instance[name] for name in self._optimizable_hyperparams_names]) for instance in initial_hyperparams_dicts])
        initial_privacy_vals = np.expand_dims(np.array([instance[1] for instance in processed_data]), axis=1)
        initial_utility_vals = 1-np.expand_dims(np.array([instance[2] for instance in processed_data]), axis=1)
        max_initial_utility_variance = np.max([instance[3] for instance in processed_data])
        return initial_hyperparams_vals, initial_privacy_vals, initial_utility_vals, max_initial_utility_variance

    def _load_points_from_file(self, filename):
        with open(filename, 'rb') as f:
            raw_data = pickle.load(f)
        return raw_data

    def _create_initial_transformations(self, initial_data):
        transformed_initial_hyperparams_vals = self._transform_initial_hyperparameters(initial_data[0])
        transformed_initial_privacy_vals = self._transform_initial_privacy(initial_data[1])
        transformed_initial_utility_vals = self._transform_utility(initial_data[2])
        return transformed_initial_hyperparams_vals, transformed_initial_privacy_vals, transformed_initial_utility_vals

    ##################################################################################
    # Hyperparameter-related methods
    ##################################################################################
    def _transform_initial_hyperparameters(self, initial_hyperparams_vals):
        self._hyperparam_normalization_factor = initial_hyperparams_vals.max(axis=0)
        return self._transform_hyperparameters(initial_hyperparams_vals)

    def _transform_hyperparameters(self, hyperparams):
        return hyperparams / self._hyperparam_normalization_factor

    def _untransform_hyperparameters(self, transformed_hyperparams):
        return transformed_hyperparams * self._hyperparam_normalization_factor

    ##################################################################################
    # Privacy-related methods
    ##################################################################################
    def _transform_initial_privacy(self, initial_privacy_vals):
        log_initial_privacy_vals = np.log(initial_privacy_vals)
        self._privacy_normalization_factor = np.asscalar(log_initial_privacy_vals.max(axis=0))
        return self._transform_privacy(initial_privacy_vals)

    def _transform_privacy(self, val):
        return np.log(val) / self._privacy_normalization_factor

    def _untransform_privacy(self, transformed_val):
        return np.exp(transformed_val * self._privacy_normalization_factor)

    def _create_privacy_gp(self):
        privacy_gp = gpflow.gpr.GPR(self._transformed_initial_hyperparams_vals.copy(),
                                    self._transformed_initial_privacy_vals.copy(),
                                    self.SafeMatern52(input_dim=5, ARD=True), name='privacy')
        privacy_gp.likelihood.variance = 0
        privacy_gp.likelihood.variance.fixed = True
        privacy_gp.optimize()
        return privacy_gp

    def predict_privacy(self, hyperparams, transform_hyperparams=False, untransform_prediction=True):
        gp_input = self._transform_hyperparameters(hyperparams) if transform_hyperparams else hyperparams
        pred = self._privacy_gp.predict_f(gp_input)[0]
        if untransform_prediction:
            pred = self._untransform_privacy(pred)
        return pred

    ##################################################################################
    # Utility-related methods
    ##################################################################################
    def _transform_initial_utility(self, initial_utility_vals):
        return self._transform_utility(initial_utility_vals)

    def _transform_utility(self, val):
        return sp.special.logit(val)

    def _untransform_utility(self, transformed_val):
        return sp.special.expit(transformed_val)

    def _create_utility_gp(self):
        utility_gp = gpflow.gpr.GPR(self._transformed_initial_hyperparams_vals.copy(),
                                    self._transformed_initial_utility_vals.copy(),
                                    self.SafeMatern52(input_dim=5, ARD=True), name='utility')
        utility_gp.likelihood.variance = 0.01
        utility_gp.optimize()
        return utility_gp

    def predict_utility(self, hyperparams, transform_hyperparams=False, untransform_prediction=True):
        gp_input = self._transform_hyperparameters(hyperparams) if transform_hyperparams else hyperparams
        pred = self._utility_gp.predict_f(gp_input)[0]
        if untransform_prediction:
            pred = self._untransform_utility(pred)
        return pred

    ##################################################################################
    # Multi-objective optimization-related methods
    ##################################################################################
    def predict(self, hyperparams, transform_hyperparams=False, untransform_prediction=True):
        privacy_pred = self.predict_privacy(hyperparams, transform_hyperparams, untransform_prediction)
        utility_pred = self.predict_utility(hyperparams, transform_hyperparams, untransform_prediction)
        return privacy_pred, utility_pred

    def _construct_gpflowopt_domain(self, optimization_domain):
        gpflowopt_domain = None
        for idx, name in enumerate(self._optimizable_hyperparams_names):
            transformed_domain = optimization_domain[name] / self._hyperparam_normalization_factor[idx]
            new_dim = gpflowopt.domain.ContinuousParameter(name, *transformed_domain)
            if not gpflowopt_domain:
                gpflowopt_domain = new_dim
            else:
                gpflowopt_domain += new_dim
        return gpflowopt_domain

    def get_transformed_hypervolume(self):
        hv = self._acquisition.pareto.hypervolume(self._transformed_anti_ideal_point)
        return hv

    def get_untransformed_hypervolume(self):
        # Get untransformed pf points
        pf, dom = gpflowopt.pareto.non_dominated_sort(self._acquisition.data[1])
        pf_eps = self._untransform_privacy(pf[:, 0])
        pf_err = self._untransform_utility(pf[:, 1])
        # Sort both lists by increasing x-axis (eps)
        pf_eps, pf_err = (list(t) for t in zip(*sorted(zip(pf_eps, pf_err))))

        # TODO: Do this before untransforming (in case the untransform doesn't preserve ordering).
        # Beginning with the right-most point on the pf, discard out of both lists if it is not in the region-of-interest
        for i in range(len(pf_eps)-1, -1):
            if pf_eps > self._anti_ideal_point[0] or pf_err > self._anti_ideal_point[1]:
                del pf_eps[i]
                del pf_err[i]

        hv = 0
        if len(pf_eps) > 0:
            for i in range(len(pf_eps)-1, 0):
                # Accumulate the rectangle formed by: 1) current point to the anti-ideal privacy point, 2) current point up to the next point.
                hv += ((self._anti_ideal_point[0] - pf_eps[i]) * (pf_err[i-1] - pf_err[i]))
            # Final point on the pf has rectangle formed by: 1) current point to the anti-ideal privacy point, 2) current point to the anti-ideal error point.
            hv += ((self._anti_ideal_point[0] - pf_eps[0]) * (self._anti_ideal_point[1] - pf_err[0]))
        return hv

    def run(self):
        self._previous_hyperparams = None
        self._optimization_point_number = 0

        optimization_result = self._bayesian_optimizer.optimize([self._optimization_handler], n_iter=self._num_instances)

        # Post final run tasks
        hypervolume = self.get_untransformed_hypervolume()
        self._saved_hypervolumes.append(hypervolume)
        print("Hypervolume is: {}".format(hypervolume))
        self._estimate_last_point()
        self._plot_and_save()

    ##################################################################################
    # Methods relating to the expensive function being optimized over.
    ##################################################################################
    def _convert_model_input_to_problem_instance_input(self, transformed_hyperparams_arr):
        hyperparams_arr = self._untransform_hyperparameters(transformed_hyperparams_arr)
        hyperparams_dict = self._fixed_hyperparams.copy()
        # TODO: make better
        for idx, name in enumerate(self._optimizable_hyperparams_names):
            val = hyperparams_arr[0][idx]
            if name in ['epochs', 'lot_size']:
                val = int(round(val))
            hyperparams_dict[name] = val
        return hyperparams_dict

    def _run_pre_iteration_tasks(self, transformed_hyperparams):
        # Print and store hypervolume from previous run
        hypervolume = self.get_untransformed_hypervolume()
        self._saved_hypervolumes.append(hypervolume)
        print("Hypervolume is: {}".format(hypervolume))

        # Run check to ensure previously acquired point was learned properly, then save current hyperparams for next run.
        self._estimate_last_point()
        self._previous_hyperparams = transformed_hyperparams

        # Untransform hyperparams and put into our expected format (a dictionary of hyperparam names to vals).
        hyperparams_dict = self._convert_model_input_to_problem_instance_input(transformed_hyperparams)

        # Predict privacy and utility before training model.
        estimated_privacy, estimated_utility = self.predict(transformed_hyperparams, transform_hyperparams=False)
        print("curr point: {}".format(hyperparams_dict))
        print("curr pts estimated (eps, err): {}".format((estimated_privacy, estimated_utility)))

        self._plot_and_save()

        self._optimization_point_number += 1

        save_object('{}/optimization-hyperparameters'.format(self._results_dir),
                    hyperparams_dict,
                    'round-{}'.format(self._optimization_point_number))

        return hyperparams_dict

    def _run_post_iteration_tasks(self, eps, acc):
        self._print_sanity_check()

        transformed_eps = self._transform_privacy(eps)
        transformed_err = self._transform_utility(1-acc)

        return transformed_eps, transformed_err

    def _run_replications(self, hyperparams_dict, options):
        print("Beginning run of {} replications using {} workers.".format(self._num_replications, self._num_workers))

        # Create a list of inputs for the to-be-spawned process
        training_instances_inputs = []
        for i in range(self._num_replications):
            process_hyperparams = hyperparams_dict.copy()
            process_options = options.copy()
            # Compute privacy only for the first instance of its kind, not for any subsequent replications. Determine
            # this simply using the ordering of the instances in the hyperparam_instances list
            if i > 0:
                process_options['accumulate_privacy'] = False
                process_options['verbose'] = False
            training_instance_input = {'hyperparams': process_hyperparams,
                                       'options': process_options,
                                       'model_class': self._model_class}
            training_instances_inputs.append(training_instance_input)

        # Begin multiprocessing across the instances
        self._id_queue = mp.Manager().Queue()
        for pid in range(self._num_workers):
            self._id_queue.put(pid)
        with Pool(self._num_workers, self._process_id_init, (self._id_queue,)) as pool:
            results = pool.map_async(self._run_single_instance, training_instances_inputs, chunksize=1)
            # while not results.ready():
            #     completed_instances = self._num_replications - results._number_left
            #     print("Status: {} of {} complete.".format(completed_instances, self._num_replications), end='\r')
            #     time.sleep(1)
            self._all_results += results.get()

        # process results
        aggregated_privacy = np.max([result[0] for result in results.get()])
        aggregated_utility = np.mean([result[1] for result in results.get()])
        return aggregated_privacy, aggregated_utility

    # Given some hyperparameters and options, runs a single instance on its own compute device
    @staticmethod
    def _run_single_instance(inputs):

        hyperparams = inputs['hyperparams']
        options = inputs['options']
        model_class = inputs['model_class']

        global process_id
        options.update({'device_id': process_id})
        # Ensures the process is running on the CPU identified by process_id.
        proc = psutil.Process()
        try:
            proc.cpu_affinity([process_id])
        except:
            pass  # doesn't work on Macs...but Macs seem to handle CPU assignment fine without this, so just continue.

        # Change this process' random seed
        custom_seed = int(time.time() + hash(frozenset(hyperparams.items())) + 547*float(process_id)) # don't ask
        set_random_seed(custom_seed)

        print("Worker {} evaluating hyperparameters {}.".format(process_id, hyperparams))

        # Run instance with random hyperparams
        priv, acc = -1, -1
        try:
            instance = model_class(hyperparams, options)
            priv, acc = instance.run()
            gc.collect()  # necessary to avoid out-of-memory errors
            if 'accumulate_privacy' in options and options['accumulate_privacy'] == False:
                priv = -2
        # If an instance crashes, just return (-1,-1) as a failure code
        except Exception as e:
            print("Worker {} running with hyperparameters \n{}\n threw the following exception:\n{}".format(process_id, hyperparams, e))
        finally:
            print("Worker {} returning result: ({}, {})".format(process_id, priv, acc))
            return priv, acc

    # Gets the unique process id (to be used to assign a unique compute device)
    @staticmethod
    def _process_id_init(queue):
        global process_id
        process_id = queue.get()

    def _optimization_handler(self, transformed_hyperparams):
        # GPFlowOpt discards output on Windows; the following lines correct this.
        import sys
        sys.stdout = sys.__stdout__

        # Do any computation before testing a point on the objective function.
        hyperparams_dict = self._run_pre_iteration_tasks(transformed_hyperparams)

        options = self._instance_options
        options['output_base_dir'] = self._results_dir

        eps, acc = self._run_replications(hyperparams_dict, options)
        print("actual (eps, err): {}".format((eps, 1-acc)))

        # Do any computation after testing a point on the objective function before updating the model.
        transformed_eps, transformed_err = self._run_post_iteration_tasks(eps, acc)

        return np.array([(transformed_eps, transformed_err)])

    def _estimate_last_point(self):
        if self._previous_hyperparams is not None:
            estimated_privacy, estimated_utility = self.predict(self._previous_hyperparams, transform_hyperparams=False)
            print("#######################################################################################")
            print("Checking last point acquired in optimization loop.")
            print("estimated (eps, err) of point: {}".format((estimated_privacy, estimated_utility)))
            print("#######################################################################################")

    ##################################################################################
    # Miscellaneous code.
    ##################################################################################
    def _plot_and_save(self):
        self._plot_empirical_pareto_frontier(suffix=str(self._optimization_point_number))

        # save_object('{}/saved_gp_models'.format(self._results_dir),
        #             [self._privacy_gp, self._utility_gp],
        #             'round-{}'.format(self._optimization_point_number),
        #             save_text=False)

        save_object('{}/saved_optimization_state'.format(self._results_dir),
                    self,
                    'round-{}'.format(self._optimization_point_number),
                    save_text=False)

    def _plot_empirical_pareto_frontier(self, suffix=''):
        pf, dom = gpflowopt.pareto.non_dominated_sort(self._acquisition.data[1])
        pf_eps = self._untransform_privacy(pf[:, 0])
        pf_err = self._untransform_utility(pf[:, 1])
        # Sort both lists by increasing x-axis (eps)
        pf_eps, pf_err = (list(t) for t in zip(*sorted(zip(pf_eps, pf_err))))

        all_eps = self._untransform_privacy(self._acquisition.data[1][:, 0])
        all_err = self._untransform_utility(self._acquisition.data[1][:, 1])

        plt.figure(figsize=(9, 7))

        # Plot all points, color by dominance
        plt.scatter(all_eps, all_err, c=dom, marker='.')

        # Mark all points on the pf with red X and build a stepwise dominance line
        plt.scatter(pf_eps, pf_err, color='r', marker='x')
        boundary_pt = ((pf_eps[0], 2*self._anti_ideal_point[1]), (2*self._anti_ideal_point[0], pf_eps[-1]))
        pf_eps.insert(0, boundary_pt[0][0])
        pf_err.insert(0, boundary_pt[0][1])
        pf_eps.insert(len(pf_eps), boundary_pt[1][0])
        pf_err.insert(len(pf_err), boundary_pt[1][1])
        plt.step(pf_eps, pf_err, 'k-', where='post')

        # Plot options
        self._set_plot_properties(plt, self._plot_options)
        if self._save_results:
            plots_dir = '{}/pareto-front-plots'.format(self._results_dir)
            os.makedirs(plots_dir, exist_ok=True)
            plt.savefig('{}/round-{}.png'.format(plots_dir, suffix))
        if self._show_plot:
            # TODO: resolve blocking issue when showing plot on Windows
            plt.show()

    def _set_plot_properties(self, plt, plot_options):
        title = plot_options.get('title', 'Empirical Pareto Frontier')
        xlabel = plot_options.get('xlabel', 'epsilon')
        ylabel = plot_options.get('ylabel', 'classification error')
        plt.title(title, size=20)
        plt.xlabel(xlabel, size=18)
        plt.ylabel(ylabel, size=18)

        left = plot_options.get('left', 0)
        right = plot_options.get('right', self._anti_ideal_point[0])
        bottom = plot_options.get('bottom', 0)
        top = plot_options.get('top', self._anti_ideal_point[1])
        plt.xlim(left=left, right=right)
        plt.ylim(bottom=bottom, top=top)

    def _print_sanity_check(self):
        print("#######################################################################################")
        print("Running GP sanity check on a training point.")

        # first point in loaded data
        first_point = np.expand_dims(self._initial_hyperparams_vals[0], axis=0)

        actual_privacy = self._initial_privacy_vals[0]
        actual_utility = self._initial_utility_vals[0]
        print("actual (eps, err): {}".format((actual_privacy, actual_utility)))

        # gp estimate of first point
        estimated_privacy, estimated_utility = self.predict(first_point, transform_hyperparams=True)
        print("estimated (eps, err): {}".format((estimated_privacy, estimated_utility)))

        # perturbed first point in loaded data (lot size increase by 1)
        almost_first_point = np.copy(first_point)
        almost_first_point[0][self._optimizable_hyperparams_names.index('lot_size')] += 1
        estimated_privacy, estimated_utility = self.predict(almost_first_point, transform_hyperparams=True)
        print("estimated (eps, err) of perturbed: {}".format((estimated_privacy, estimated_utility)))

        print("#######################################################################################")

    def _make_results_directory(self):
        if not self._output_base_dir:
            self._results_dir = 'results/hvpoi_results/{}/{}'.format(self._problem_name, self._run_name)
        else:
            self._results_dir = '{}/hvpoi_results/{}'.format(self._output_base_dir, self._run_name)
        os.makedirs(self._results_dir, exist_ok=True)

    # Matern52 kernel isn't numerically safe, so this work-around is necessary.
    class SafeMatern52(gpflow.kernels.Matern52):
        def euclid_dist(self, X, X2):
            r2 = self.square_dist(X, X2)
            return tf.sqrt(r2 + 1e-7)

ext: py

reason: import numpy,import scipy


text:
import cv2
import os
import pickle
import numpy as np
from openVO import StereoCamera

config_path = "../Data/Calibration/competition1/"
assert(os.path.exists(config_path))
intr_left = pickle.load(open(config_path + "intrinsics_left.p", "rb"))
K_left, dist_left = intr_left["K"], intr_left["dist"]
intr_right = pickle.load(open(config_path + "intrinsics_right.p", "rb"))
K_right, dist_right = intr_right["K"], intr_right["dist"]
rect_params = pickle.load(open(config_path + "rectParams.p", "rb"))
print("rect params", rect_params)

# Generated with script
# {'minDisparity': 0, 'numDisparities': 80, 'blockSize': 3, 'P1': 0, 'P2': 675, 'disp12MaxDiff': 2, 'preFilterCap': 0, 'uniquenessRatio': 6, 'speckleWindowSize': 50, 'speckleRange': 2}
# Existing params
# "minDisparity": 0,
# "numDisparities": 48,
# "blockSize": 3,
# "P1": 27,
# "P2": 675,
# "disp12MaxDiff": 50,
# "preFilterCap": 0,
# "uniquenessRatio": 10,
# "speckleWindowSize": 0,
# "speckleRange": 4}

sgbm_params = {"minDisparity": 0,
               "numDisparities": 80,
               "blockSize": 3,
               "P1": 27,
               "P2": 675,
               "disp12MaxDiff": 2,
               "preFilterCap": 0,
               "uniquenessRatio": 6,
               "speckleWindowSize": 50,
               "speckleRange": 2}

cam_left = cv2.VideoCapture(2)
cam_right = cv2.VideoCapture(4)

# Adjust to overlay images (depends on distance to object)
shift = 100 # pixels

while True:
    got_left, left = cam_left.read()
    got_right, right = cam_right.read()
    if not (got_left and got_right):
        break
    combined = np.concatenate((left, right), axis=1)
    cv2.imshow("original", combined)
    undist_left = cv2.undistort(left, K_left, dist_left)
    undist_right = cv2.undistort(right, K_right, dist_right)
    undist = np.concatenate((undist_left, undist_right), axis=1)
    cv2.imshow("Undistorted", undist)

    stereo = StereoCamera(K_left, dist_left, K_right, dist_right, rect_params, sgbm_params, (left.shape[1], left.shape[0]))
    rect_left = stereo.undistort_rectify_left(left)
    rect_right = stereo.undistort_rectify_right(right)
    rect = np.concatenate((rect_left, rect_right), axis=1)
    cv2.imshow("Rectified", rect)
    gray_left = cv2.cvtColor(left, cv2.COLOR_BGR2GRAY)
    gray_right = cv2.cvtColor(right, cv2.COLOR_BGR2GRAY)
    gray_undist_left = cv2.cvtColor(undist_left, cv2.COLOR_BGR2GRAY)
    gray_undist_right = cv2.cvtColor(undist_right, cv2.COLOR_BGR2GRAY)
    gray_rect_left = cv2.cvtColor(rect_left, cv2.COLOR_BGR2GRAY)
    gray_rect_right = cv2.cvtColor(rect_right, cv2.COLOR_BGR2GRAY)
    h, w = gray_left.shape
    shift_block = np.zeros((h, shift))
    red_channel = np.zeros((h, w+shift))
    gray_left = np.concatenate((gray_left, shift_block), axis = 1)
    gray_undist_left = np.concatenate((gray_undist_left, shift_block), axis = 1)
    gray_rect_left = np.concatenate((gray_rect_left, shift_block), axis = 1)
    gray_right = np.concatenate((shift_block, gray_right), axis = 1)
    gray_undist_right = np.concatenate((shift_block, gray_undist_right), axis = 1)
    gray_rect_right = np.concatenate((shift_block, gray_rect_right), axis = 1)
    stacked = np.stack((gray_left, gray_right, red_channel), axis=2).astype(np.uint8)
    stacked_undist = np.stack((gray_undist_left, gray_undist_right, red_channel), axis=2).astype(np.uint8)
    stacked_rect = np.stack((gray_rect_left, gray_rect_right, red_channel), axis=2).astype(np.uint8)

    for x in range(h//10 - 1):
        cv2.line(stacked, (0, 10 * x), (w+shift - 1, 10 * x), (0,0,255))
        cv2.line(stacked_undist, (0, 10 * x), (w+shift - 1, 10 * x), (0,0,255))
        cv2.line(stacked_rect, (0, 10 * x), (w+shift - 1, 10 * x), (0,0,255))
    cv2.imshow("Stacked Original", stacked)
    cv2.imshow("Stacked Undistoreted", stacked_undist)
    cv2.imshow("Stacked Rectified", stacked_rect)

    _, disparity, _ = stereo.compute_3d(left, right)
    max, min = disparity.max(), disparity.min()
    normed_disp = (disparity - min) / (max - min)
    cv2.imshow('disparity', normed_disp)
    if cv2.waitKey(40) == 27:
        break

reason: import numpy


text:
#!/usr/bin/env python

# Reader for CINE files produced by Vision Research Phantom Software

################################################################################

# This code was originally developed (w/o a license) for the Irvine lab at the
#  university of Chicago.  Moved to an Apache License (with permision) when
#  ported to Python 3 in 2019 by Dustin Kleckner (dkleckner@ucmerced.edu).

################################################################################

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

################################################################################

import sys
import os
import time
import struct
import numpy as np
from threading import Lock
import datetime
import hashlib

FRACTION_MASK = (2**32-1)
MAX_INT = 2**32

#Define types for fields
BYTE = 'B'
WORD = 'H'
INT16 = 'h'
SHORT = 'h'
BOOL = 'i'
DWORD = 'I'
UINT = 'I'
LONG = 'l'
INT = 'l'
FLOAT = 'f'
DOUBLE = 'd'
TIME64 = 'Q'
RECT = '4i'
WBGAIN = '2f'
IMFILTER = '28i'

#Define Header fields
TAGGED_FIELDS = {
    1000: ('ang_dig_sigs', ''),
    1001: ('image_time_total', TIME64),
    1002: ('image_time_only', TIME64),
    1003: ('exposure_only', DWORD),
    1004: ('range_data', ''),
    1005: ('binsig', ''),
    1006: ('anasig', ''),
    1007: ('undocumented', '')}  # 1007 exists in my files, but is not in documentation I can find

HEADER_FIELDS = [
    ('type', '2s'),
    ('header_size', WORD),
    ('compression', WORD),
    ('version', WORD),
    ('first_movie_image', LONG),
    ('total_image_count', DWORD),
    ('first_image_no', LONG),
    ('image_count', DWORD),
    ('off_image_header', DWORD),
    ('off_setup', DWORD),
    ('off_image_offsets', DWORD),
    ('trigger_time', TIME64),
]

BITMAP_INFO_FIELDS = [
    ('bi_size', DWORD),
    ('bi_width', LONG),
    ('bi_height', LONG),
    ('bi_planes', WORD),
    ('bi_bit_count', WORD),
    ('bi_compression', DWORD),
    ('bi_image_size', DWORD),
    ('bi_x_pels_per_meter', LONG),
    ('bi_y_pels_per_meter', LONG),
    ('bi_clr_used', DWORD),
    ('bi_clr_important', DWORD),
]

SETUP_FIELDS = [
    ('frame_rate_16', WORD),
    ('shutter_16', WORD),
    ('post_trigger_16', WORD),
    ('frame_delay_16', WORD),
    ('aspect_ratio', WORD),
    ('contrast_16', WORD),
    ('bright_16', WORD),
    ('rotate_16', BYTE),
    ('time_annotation', BYTE),
    ('trig_cine', BYTE),
    ('trig_frame', BYTE),
    ('shutter_on', BYTE),
    ('description_old', '121s'),  # Guessed at length... because it isn't documented!  This seems to work.
    ('mark', '2s'),
    ('length', WORD),
    ('binning', WORD),
    ('sig_option', WORD),
    ('bin_channels', SHORT),
    ('samples_per_image', BYTE)] + \
    [('bin_name%d' % i, '11s') for i in range(8)] + [
        ('ana_option', WORD),
        ('ana_channels', SHORT),
        ('res_6', BYTE),
        ('ana_board', BYTE)] + \
    [('ch_option%d' % i, SHORT) for i in range(8)] + \
    [('ana_gain%d' % i, FLOAT) for i in range(8)] + \
    [('ana_unit%d' % i, '6s') for i in range(8)] + \
    [('ana_name%d' % i, '11s') for i in range(8)] + [
    ('i_first_image', LONG),
    ('dw_image_count', DWORD),
    ('n_q_factor', SHORT),
    ('w_cine_file_type', WORD)] + \
    [('sz_cine_path%d' % i, '65s') for i in range(4)] + [
    ('b_mains_freq', WORD),
    ('b_time_code', BYTE),
    ('b_priority', BYTE),
    ('w_leap_sec_dy', DOUBLE),
    ('d_delay_tc', DOUBLE),
    ('d_delay_pps', DOUBLE),
    ('gen_bits', WORD),
    ('res_1', INT16),  # Manual says INT, but this is clearly wrong!
    ('res_2', INT16),
    ('res_3', INT16),
    ('im_width', WORD),
    ('im_height', WORD),
    ('edr_shutter_16', WORD),
    ('serial', UINT),
    ('saturation', INT),
    ('res_5', BYTE),
    ('auto_exposure', UINT),
    ('b_flip_h', BOOL),
    ('b_flip_v', BOOL),
    ('grid', UINT),
    ('frame_rate', UINT),
    ('shutter', UINT),
    ('edr_shutter', UINT),
    ('post_trigger', UINT),
    ('frame_delay', UINT),
    ('b_enable_color', BOOL),
    ('camera_version', UINT),
    ('firmware_version', UINT),
    ('software_version', UINT),
    ('recording_time_zone', INT),
    ('cfa', UINT),
    ('bright', INT),
    ('contrast', INT),
    ('gamma', INT),
    ('reserved1', UINT),
    ('auto_exp_level', UINT),
    ('auto_exp_speed', UINT),
    ('auto_exp_rect', RECT),
    ('wb_gain', '8f'),
    ('rotate', INT),
    ('wb_view', WBGAIN),
    ('real_bpp', UINT),
    ('conv_8_min', UINT),
    ('conv_8_max', UINT),
    ('filter_code', INT),
    ('filter_param', INT),
    ('uf', IMFILTER),
    ('black_cal_sver', UINT),
    ('white_cal_sver', UINT),
    ('gray_cal_sver', UINT),
    ('b_stamp_time', BOOL),
    ('sound_dest', UINT),
    ('frp_steps', UINT),
    ] + [('frp_img_nr%d' % i, INT) for i in range(16)] + \
        [('frp_rate%d' % i, UINT) for i in range(16)] + \
        [('frp_exp%d' % i, UINT) for i in range(16)] + [
    ('mc_cnt', INT),
    ] + [('mc_percent%d' % i, FLOAT) for i in range(64)] + [
    ('ci_calib', UINT),
    ('calib_width', UINT),
    ('calib_height', UINT),
    ('calib_rate', UINT),
    ('calib_exp', UINT),
    ('calib_edr', UINT),
    ('calib_temp', UINT),
    ] + [('header_serial%d' % i, UINT) for i in range(4)] + [
    ('range_code', UINT),
    ('range_size', UINT),
    ('decimation', UINT),
    ('master_serial', UINT),
    ('sensor', UINT),
    ('shutter_ns', UINT),
    ('edr_shutter_ns', UINT),
    ('frame_delay_ns', UINT),
    ('im_pos_xacq', UINT),
    ('im_pos_yacq', UINT),
    ('im_width_acq', UINT),
    ('im_height_acq', UINT),
    ('description', '4096s')
]

T64_F = lambda x: int(x) / 2.**32
T64_F_ms = lambda x: '%.3f' % (float(x.rstrip('L')) / 2.**32)
T64_S = lambda s: lambda t: time.strftime(s, time.localtime(float(t.rstrip('L'))/2.**32))

#Processing the data in chunks keeps it in the L2 catch of the processor, increasing speed for large arrays by ~50%
CHUNK_SIZE = 6 * 10**5 #Should be divisible by 3, 4 and 5!  This seems to be near-optimal.

def ten2sixteen(a):
    b = np.zeros(a.size//5*4, dtype='u2')

    for j in range(0, len(a), CHUNK_SIZE):
        (a0, a1, a2, a3, a4) = [a[j+i:j+CHUNK_SIZE:5].astype('u2') for i in range(5)]

        k = j//5 * 4
        k2 = k + CHUNK_SIZE//5 * 4

        b[k+0:k2:4] = ((a0 & 0b11111111) << 2) + ((a1 & 0b11000000) >> 6)
        b[k+1:k2:4] = ((a1 & 0b00111111) << 4) + ((a2 & 0b11110000) >> 4)
        b[k+2:k2:4] = ((a2 & 0b00001111) << 6) + ((a3 & 0b11111100) >> 2)
        b[k+3:k2:4] = ((a3 & 0b00000011) << 8) + ((a4 & 0b11111111) >> 0)

    return b

def sixteen2ten(b):
    a = np.zeros(b.size//4*5, dtype='u1')

    for j in range(0, len(a), CHUNK_SIZE):
        (b0, b1, b2, b3) = [b[j+i:j+CHUNK_SIZE:4] for i in range(4)]

        k = j//4 * 5
        k2 = k + CHUNK_SIZE//4 * 5

        a[k+0:k2:5] =                              ((b0 & 0b1111111100) >> 2)
        a[k+1:k2:5] = ((b0 & 0b0000000011) << 6) + ((b1 & 0b1111110000) >> 4)
        a[k+2:k2:5] = ((b1 & 0b0000001111) << 4) + ((b2 & 0b1111000000) >> 6)
        a[k+3:k2:5] = ((b2 & 0b0000111111) << 2) + ((b3 & 0b1100000000) >> 8)
        a[k+4:k2:5] = ((b3 & 0b0011111111) << 0)

    return a

def twelve2sixteen(a):
    b = np.zeros(a.size//3*2, dtype='u2')

    for j in range(0, len(a), CHUNK_SIZE):
        (a0, a1, a2) = [a[j+i:j+CHUNK_SIZE:3].astype('u2') for i in range(3)]

        k = j//3 * 2
        k2 = k + CHUNK_SIZE//3 * 2

        b[k+0:k2:2] = ((a0 & 0xFF) << 4) + ((a1 & 0xF0) >> 4)
        b[k+1:k2:2] = ((a1 & 0x0F) << 8) + ((a2 & 0xFF) >> 0)

    return b

def sixteen2twelve(b):
    a = np.zeros(b.size//2*3, dtype='u1')

    for j in range(0, len(a), CHUNK_SIZE):
        (b0, b1) = [b[j+i:j+CHUNK_SIZE:2] for i in range(2)]

        k = j//2 * 3
        k2 = k + CHUNK_SIZE//2 * 3

        a[k+0:k2:3] =                       ((b0 & 0xFF0) >> 4)
        a[k+1:k2:3] = ((b0 & 0x00F) << 4) + ((b1 & 0xF00) >> 8)
        a[k+2:k2:3] = ((b1 & 0x0FF) << 0)

    return a

class Cine(object):
    '''Class for reading Vision Research CINE files, e.g. from Phantom cameras.

    Supports indexing, so frame can be accessed like list items, and ``len``
    returns the number of frames.  Iteration is also supported.

    Cine objects also use locks for file reading, allowing cine objects to
    be shared safely by several threads.

    Parameters
    ---------
    filename : string
        Source filename
    '''

    def __init__(self, fn):
        self.f = open(fn, 'rb')
        self.fn = fn

        self.read_header(HEADER_FIELDS)
        self.read_header(BITMAP_INFO_FIELDS, self.off_image_header)
        self.read_header(SETUP_FIELDS, self.off_setup)
        self.image_locations = self.unpack('%dQ' % self.image_count, self.off_image_offsets)
        if type(self.image_locations) not in (list, tuple):
            self.image_locations = [self.image_locations]

        self.width = self.bi_width
        self.height = self.bi_height

        self.file_lock = Lock()  # Allows Cine object to be accessed from multiple threads!

        self._hash = None

    def gamma_corrected_frame(self, frame_number, bottom_clip=0, top_clip=None, gamma=2.2):
        '''Return a frame as a gamma corrected 'u1' array, suitable for saving
        to a standard image.

        Output is equal to: ``255 * ((original - bottom_clip) / (top_clip - bottom_clip))**(1/gamma)``

        Parameters
        ----------
        frame_number : integer
        gamma : float (default: 2.2)
            The gamma correction to apply
        top_clip : integer (default: 0)
        bottom_clip : integer (default: 2**real_bpp)

        Returns
        -------
        frame : numpy array (dtype='u1')

        '''
        if top_clip is None: top_clip = 2**self.real_bpp

        return (255 * np.clip(((self.get_frame(frame_number) - bottom_clip) / float(top_clip - bottom_clip)), 0, 1)**(1.0/gamma)).astype('u1')

    def get_frame(self, frame_number):
        '''Get a frame from the cine file.

        Parameters
        ----------
        frame_number : integer

        Returns
        -------
        frame : numpy array (dtype='u1' or 'u2', depending on bit depth)
        '''

        self.file_lock.acquire()

        image_start = self.image_locations[frame_number]
        annotation_size = self.unpack(DWORD, image_start)
        annotation = self.unpack('%db' % (annotation_size - 8))
        image_size = self.unpack(DWORD)

        #self.f.seek(image_start + annotation_size-8)
        data_type = 'u1' if self.bi_bit_count in (8, 24) else 'u2'

        actual_bits = image_size * 8 // (self.width * self.height)

        if actual_bits in (10, 12):
            data_type = 'u1'

        self.f.seek(image_start + annotation_size)

        frame = np.frombuffer(self.f.read(image_size), data_type)

        if (actual_bits == 10):
           frame = ten2sixteen(frame)
        elif (actual_bits == 12):
           frame = twelve2sixteen(frame)

        # if (actual_bits % 8):
        #     raise ValueError('Data should be byte aligned, packed data not supported' % actual_bits)

        frame = frame.reshape(self.height, self.width)[::-1]

        if actual_bits in (10, 12):
           frame = frame[::-1, :]  # Don't know why it works this way, but it does...

        self.file_lock.release()

        return frame

    def __len__(self):
        return self.image_count

    len = __len__

    def __getitem__(self, key):
        if type(key) == slice:
            return map(self.get_frame, range(self.image_count)[key])

        return self.get_frame(key)

    def get_time(self, frame_number):
        '''Get the time of a specific frame.

        Parameters
        ----------
        frame_number : integer

        Returns
        -------
        time : float
            Time from start in seconds.'''
        return float(frame_number) / self.frame_rate

    def get_fps(self):
        '''Get the frames per second of the movie.

        Returns
        -------
        fps : int
        '''
        return self.frame_rate

    def __iter__(self):
        self._iter_current_frame = -1
        return self

    def next(self):
        self._iter_current_frame += 1
        if self._iter_current_frame >= self.image_count:
            raise StopIteration
        else:
            return self.get_frame(self._iter_current_frame)

    def close(self):
        '''Closes the cine file.'''
        self.f.close()

    #---------------------------------------------------------------------------
    #These functions are not meant to be used externally, and so are
    #   undocumented.
    #---------------------------------------------------------------------------
    def unpack(self, fs, offset=None):
        if offset is not None:
            self.f.seek(offset)
        s = struct.Struct('<' + fs)
        vals = s.unpack(self.f.read(s.size))
        if len(vals) == 1:
            return vals[0]
        else:
            return vals

    def read_tagged_blocks(self):
        if not self.off_setup + self.length < self.off_image_offsets:
            return
        next_tag_exists = True
        next_tag_offset = 0
        while next_tag_exists:
            block_size, next_tag_exists = self._read_tag_block(next_tag_offset)
            next_tag_offset += block_size

    def _read_tag_block(self, off_set):
        self.file_lock.acquire()
        self.f.seek(self.off_setup + self.length + off_set)
        block_size = self.unpack(DWORD)
        b_type = self.unpack(WORD)
        more_tags = self.unpack(WORD)

        if b_type == 1004:
            # docs say to ignore range data
            # it seems to be a poison flag, if see this, give up tag parsing
            return block_size, 0

        try:
            d_name, d_type = TAGGED_FIELDS[b_type]

        except KeyError:
            #            print 'unknown type, find an updated version of file spec', b_type
            return block_size, more_tags

        if d_type == '':
            #            print "can't deal with  <" + d_name + "> tagged data"
            return block_size, more_tags

        s_tmp = struct.Struct('<' + d_type)
        if (block_size-8) % s_tmp.size != 0:
            #            print 'something is wrong with your data types'
            return block_size, more_tags

        d_count = (block_size-8)//(s_tmp.size)

        data = self.unpack('%d' % d_count + d_type)
        if not isinstance(data, tuple):
            # fix up data due to design choice in self.unpack
            data = (data, )

        # parse time
        if b_type == 1002 or b_type == 1001:
            data = [(datetime.datetime.fromtimestamp(d >> 32), float((FRACTION_MASK & d))/MAX_INT) for d in data]
        # convert exposure to seconds
        if b_type == 1003:
            data = [float(d)/(MAX_INT) for d in data]

        setattr(self, d_name, data)

        self.file_lock.release()
        return block_size, more_tags

    def read_header(self, fields, offset=0):
        self.f.seek(offset)
        for name, format in fields:
            setattr(self, name, self.unpack(format))

    def __unicode__(self):
        return self.fn

    def __str__(self):
        return unicode(self).encode('utf-8')

    __repr__ = __unicode__

    @property
    def trigger_time_p(self):
        '''Returns the time of the trigger, tuple of (datatime_object, fraction_in_ns)'''
        return datetime.datetime.fromtimestamp(self.trigger_time >> 32), float(FRACTION_MASK & self.trigger_time)/(MAX_INT)

    @property
    def hash(self):
        if self._hash is None:
            self._hash_fun()
        return self._hash

    def __hash__(self):
        return int(self.hash, base=16)

    def _hash_fun(self):
        """
        Generates the md5 hash of the header of the file.  Here the
        header is defined as everything before the first image starts.

        This includes all of the meta-data (including the plethora of
        time stamps) so this will be unique.
        """
        # get the file lock (so we don't screw up any other reads)
        self.file_lock.acquire()

        self.f.seek(0)
        max_loc = self.image_locations[0]
        md5 = hashlib.md5()

        chunk_size = 128*md5.block_size
        chunk_count = (max_loc//chunk_size) + 1

        for j in range(chunk_count):
            md5.update(self.f.read(128*md5.block_size))

        self._hash = md5.hexdigest()

        self.file_lock.release()

    def __eq__(self, other):
        return self.hash == other.hash

    def __ne__(self, other):
        return not self == other

max_stars_repo_name: klecknerlab/muvi

max_issues_repo_name: klecknerlab/muvi

max_forks_repo_name: klecknerlab/muvi


text:
#!/usr/bin/env python2
# -*- coding: utf-8 -*-

import numpy as np
import statsmodels.stats.power as pw

def cohen_d(orig_accuracy_list, accuracy_list):
    nx = len(orig_accuracy_list)
    ny = len(accuracy_list)
    dof = nx + ny - 2
    pooled_std = np.sqrt(((nx-1)*np.std(orig_accuracy_list, ddof=1) ** 2 + (ny-1)*np.std(accuracy_list, ddof=1) ** 2) / dof)
    result = (np.mean(orig_accuracy_list) - np.mean(accuracy_list)) / pooled_std
    return result

def power(orig_accuracy_list, mutation_accuracy_list):
    eff_size = cohen_d(orig_accuracy_list, mutation_accuracy_list)
    pow = pw.FTestAnovaPower().solve_power(effect_size=eff_size, nobs=len(orig_accuracy_list) + len(mutation_accuracy_list), alpha=0.05)
    return pow

max_stars_repo_name: dlfaults/DCReplication

max_issues_repo_name: dlfaults/DCReplication

max_forks_repo_name: dlfaults/DCReplication

reason: import numpy,import statsmodels


text:
# Some helper functions
import numpy as np
import scipy.special
from scipy.signal import butter, lfilter

# Return real part of exponential integral, same as matlab expint()
def expint(v):
    return np.real(-scipy.special.expi(-v)-np.pi*1j)

# Circular shift of an array
def circular_shift(x,t):
    return [x[t:len(x)], x[0:t]]

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

def bandpass(raw_input,preprocess,high_cut,fs):
    if (preprocess == "butter"):
        input = butter_bandpass_filter(raw_input,20,high_cut,fs)
        input = input / max(abs(input))
    elif (preprocess == "ellip"):
        sos = scipy.signal.ellip(4,5,40,[20/(fs/2),high_cut/(fs/2)],btype='bandpass', output='sos')
        input = scipy.signal.sosfilt(sos, raw_input)
        input = input / max(abs(input))
    else:
        input = raw_input
    return input

reason: import numpy,import scipy,from scipy


text:
import tensorflow as tf
#from libs import batch_norm
#from libs import data_pros

import pandas as pd
import numpy as np

# Loading the data from our  dataset

#dirname = ...
#dirname_1 = "/home/stanlee321/Desktop/SIN_TODO/newton/1.0/data.csv"
#dirname_2 = "/home/stanlee321/Desktop/SIN_TODO/newton/1.0/newton.csv"

dirname_1 = "/notebooks/Power/newton/1.0/data.csv"
dirname_2 = "/notebooks/Power/newton/1.0/newton.csv"

busd   = pd.read_csv(dirname_1)
solbus = pd.read_csv(dirname_2)
busd.fillna(0.0)
solbus.fillna(0.0)

busd["Bus2"] = busd['Bus']*busd['Bus']
busd["Type2"] = busd['Type']*busd['Type']
busd["Vsp2"] = busd['Vsp']*busd['Vsp']
busd["Pgi2"] = busd['Pgi']*busd['Pgi']
busd["Qgi2"] = busd['Qgi']*busd['Qgi']

busd["Qgi/Pgi"] = busd["Qgi"]/(busd["Pgi"]+1.0)

busd["Vsp/Pgi"] = busd["Vsp"]/(busd["Pgi"]+1.0)
busd["Vsp/Qgi"] = busd["Vsp"]/(busd["Qgi"]+1.0)
busd["Vsp/Pli"] = busd["Vsp"]/busd["Pli"]
busd["Vsp/Qli"] = busd["Vsp"]/busd["Qli"]

columns_busd = list(busd.columns.values)
columns_solbus = list(solbus.columns.values)

df1 = busd
df2 = solbus

# Pre-processing the data

from sklearn import preprocessing
from sklearn.cross_validation import train_test_split

import matplotlib.pyplot as plt

# Function for Normalize the Datas and Split into Training and CV set...
def process_data(df1,df2,size,r_state):

    # Para normalizar los DataFrame, se extrae las columnas que no se quieren normalizar,
    # en nuestro caso para df1 es Index=[0,1] y para df2 es Index = [0]
    # Nombres de las columnas originales
    columns_features = list(df1.columns.values)
    columns_labels   = list(df2.columns.values)

    ####Normalizando Busdata
    x = df1.values
    x_scaled = preprocessing.scale(x)
    df1_norm_all = pd.DataFrame(x_scaled)

    df1_norm_all = df1_norm_all.drop(df1_norm_all.columns[[3]], axis=1)

    # For Solution data set
    y = df2.values
    y_scaled = preprocessing.scale(y)
    df2_norm_all = pd.DataFrame(y_scaled)

    # Spliting
    X_train, X_test, y_train, y_test = train_test_split(df1_norm_all, df2_norm_all, test_size= size, random_state = r_state)

    X_cv, X_test, y_cv, y_test = train_test_split(X_test,y_test, test_size = 0.50, random_state = r_state)

    X_train = np.array(X_train)
    y_train = np.array(y_train)

    X_test = np.array(X_test)
    y_test = np.array(y_test)

    X_cv = np.array(X_cv)
    y_cv = np.array(y_cv)

    return (X_train, y_train, X_test, y_test, X_cv, y_cv, df1_norm_all,df2_norm_all)

X_train, y_train, X_test, y_test, X_cv, y_cv, df1_norm_all, df2_norm_all\
 = process_data(df1,df2, 0.30, r_state = None )

learning_rate = 0.001
batch_size = 1000
#n_iterations,\
total_batches = int(X_train.shape[0]/batch_size)
n_neurons = 500
n_layers = 4
activation_fn = tf.nn.relu
final_activation_fn = tf.nn.softmax
cost_type = 'l2_norm'
hm_epochs = 1000

logs_path = '/notebooks/DeepLearningTasks/proyect/tensorboard_logs'

# Linear model

def linear(x, n_output, name=None, activation=None, reuse=None):
    """Fully connected layer.

    Parameters
    ----------
    x : tf.Tensor
        Input tensor to connect
    n_output : int
        Number of output neurons
    name : None, optional
        Scope to apply

    Returns
    -------
    h, W : tf.Tensor, tf.Tensor
        Output of fully connected layer and the weight matrix
    """
    if len(x.get_shape()) != 2:
        x = flatten(x, reuse=reuse)

    n_input = x.get_shape().as_list()[1]

    with tf.variable_scope(name or "fc", reuse=reuse):
        W = tf.get_variable(
            name='W',
            shape=[n_input, n_output],
            dtype=tf.float32,
            initializer=tf.contrib.layers.xavier_initializer())

        b = tf.get_variable(
            name='b',
            shape=[n_output],
            dtype=tf.float32,
            initializer=tf.constant_initializer(0.0))

        h = tf.nn.bias_add(
            name='h',
            value=tf.matmul(x, W),
            bias=b)

        if activation:
            h = activation(h)

        return h, W

def build_model(Xs,ys,n_neurons, n_layers, activation_fn, final_activation_fn, cost_type):

	xs = Xs
	ys = ys

	n_xs = xs.shape[1]
	n_ys = ys.shape[1]

	X = tf.placeholder(name='X', shape=[None, n_xs], dtype=tf.float32)
	Y = tf.placeholder(name='Y', shape=[None, n_ys], dtype=tf.float32)

	current_input = X

	for layer_i in range(n_layers):
		current_input = linear(current_input, n_neurons,\
		 activation=activation_fn, name='layer{}'.format(layer_i))[0]

	Y_pred = linear(current_input, n_ys, activation=final_activation_fn, name='prediction')[0]

	if cost_type == 'l1_norm':
		cost = tf.reduce_mean(tf.reduce_sum(tf.abs(Y - Y_pred), 1))
	elif cost_type == 'l2_norm':
		cost = tf.reduce_mean(tf.reduce_sum(tf.squared_difference(Y,Y_pred), 1))
	else:
		raise ValuError("Unkonw value, chosse between l1 or l2")

	saver = tf.train.Saver()

	tf_log = 'tf_log'

	return {'X': X,'Y': Y, 'Y_pred': Y_pred, 'cost': cost, 'saver': saver, 'tf_log': tf_log}
	#return {'X': X,'Y': Y, 'Y_pred': Y_pred, 'cost': cost}

def train_nn(X,y,X_cv,y_cv,X_test,y_test,\
		learning_rate,\
		n_neurons,\
		n_layers,\
		activation_fn,\
		final_activation_fn,\
		cost_type,\
		batch_size,\
		total_batches,\
		hm_epochs):

	Xs = np.array(X).reshape(-1,X.shape[1])
	ys = np.array(y).reshape(-1,y.shape[1])

	Xs_cv = np.array(X_cv).reshape(-1,X_cv.shape[1])
	ys_cv = np.array(y_cv).reshape(-1,y_cv.shape[1])

	Xs_test = np.array(X_test).reshape(-1,X_test.shape[1])
	ys_test = np.array(y_test).reshape(-1,y_test.shape[1])

	model = build_model(Xs, ys, n_neurons, n_layers,\
		activation_fn, final_activation_fn, cost_type)
	with tf.name_scope('Model'):
		prediction = model['Y_pred']
	with tf.name_scope('Loss'):
		cost = model['cost']
	with tf.name_scope('SGD'):
		optimizer = tf.train.AdamOptimizer(learning_rate=learning_rate).minimize(cost)
	with tf.name_scope('Accuracy'):
		correct = tf.equal(tf.argmax(prediction, 1), tf.argmax(model['Y'],1))
		acc = tf.reduce_mean(tf.cast(correct, 'float'))
	# Tensorboard
	tf.scalar_summary('loss',cost)
	tf.scalar_summary('accuracy', acc)

	merged_summary_op = tf.merge_all_summaries()

	saver = model['saver']
	tf_log = model['tf_log']

	with tf.Session() as sess:

		sess.run(tf.global_variables_initializer())

		summary_writer = tf.train.SummaryWriter(logs_path, graph = tf.get_default_graph())

		try:
			epoch = int(open(tf_log,'r').read().split('\n')[-2]) +1
			print('STARTING:', epoch)
		except:
			epoch = 1

		for epoch in range(hm_epochs):
			if epoch != 1:
				saver.restore(sess,"/notebooks/DeepLearningTasks/proyect/model.ckpt")

			epoch_loss = 0

			i = 0

			while i < (Xs.shape[0]):
				start = i
				end = i + batch_size
				batch_x = np.array(Xs[start:end])
				batch_y = np.array(ys[start:end])

				_, c, summary = sess.run([optimizer, cost, merged_summary_op], feed_dict={model['X']: batch_x, model['Y'] : batch_y})
				summary_writer.add_summary(summary, epoch * total_batches + 1)
				#epoch_loss += c
				epoch_loss = c / total_batches

				i +=batch_size
			saver.save(sess, "model.ckpt")

			with open(tf_log, 'a') as f:
				f.write(str(epoch)+'\n')

			print('Epoch', epoch + 1, 'completed out of ', hm_epochs,'loss:', epoch_loss)

			correct = tf.equal(tf.argmax(prediction, 1), tf.argmax(model['Y'],1))
			accuracy = tf.reduce_mean(tf.cast(correct, 'float'))
			print('Accuracy: ', accuracy.eval({model['X']:Xs_test, model['Y']:ys_test}))
		print("Run teh comand: \n" \
			"--> tensorboard -- logdir=/tmp/tensorflow_logs"\
			"\nThen open http://0.0.0.0:6666/ into your web browser")

train_nn(X_train,y_train,X_cv,y_cv,X_test,y_test,\
		learning_rate,\
		n_neurons,\
		n_layers,\
		activation_fn,\
		final_activation_fn,\
		cost_type,\
		batch_size,\
		total_batches,\
		hm_epochs)


text:
import numpy as np

def step_function(x):
    y = x > 0
    return y.astype(np.int)

max_stars_repo_name: BlueWay-KU/Study

max_issues_repo_name: BlueWay-KU/Study

max_forks_repo_name: BlueWay-KU/Study


text:
#
# Copyright (c) 2017, UT-BATTELLE, LLC
# All rights reserved.
#
# This software is released under the BSD license detailed
# in the LICENSE file in the top level a-prime directory
#
import numpy

def get_season_months_index(begin_month, end_month):

    if end_month < begin_month:
            n_months_season_yr1 = 12 - begin_month
            n_months_season_yr2 = end_month + 1
            n_months_season = n_months_season_yr1 + n_months_season_yr2

            index_months = numpy.zeros(n_months_season, dtype = numpy.int)
            index_months[0:n_months_season_yr1] = numpy.arange(begin_month,12)
            index_months[n_months_season_yr1:n_months_season] = numpy.arange(0, end_month+1)

    else:
            n_months_season = end_month - begin_month + 1
            index_months = numpy.arange(begin_month, end_month+1)

    return (index_months, n_months_season)

max_issues_repo_name: ACME-Climate/a-prime

max_forks_repo_name: ACME-Climate/a-prime


text:
from tool.runners.python import SubmissionPy
import numpy as np

class ThomasSubmission(SubmissionPy):
    def run(self, s):
        puzzle_input = [int(c) for c in s]
        recipes = [3, 7]
        elves = [0, 1]
        found_from = None
        while found_from is None:
            if np.array_equal(recipes[-len(puzzle_input) - 1 : -1], puzzle_input):
                found_from = -len(puzzle_input) - 1
                break
            if np.array_equal(recipes[-len(puzzle_input) :], puzzle_input):
                found_from = -len(puzzle_input)
                break
            new = recipes[elves[0]] + recipes[elves[1]]
            recipes += [int(c) for c in str(new)]
            elves = [
                (elves[0] + 1 + recipes[elves[0]]) % len(recipes),
                (elves[1] + 1 + recipes[elves[1]]) % len(recipes),
            ]
        return str(len(recipes) + found_from)


text:
# Blackhurst_IO.py (scripts)
# !/usr/bin/env python3
# coding=utf-8
"""
Scrapes data from Blackhurst paper 'Direct and Indirect Water Withdrawals
for US Industrial Sectors' (Supplemental info)
Includes supporting functions for Blackhurst paper data.
"""

import io
import tabula
import pandas as pd
import numpy as np
import flowsa
from flowsa.common import US_FIPS, flow_by_activity_fields, fba_fill_na_dict
from flowsa.flowbyfunctions import assign_fips_location_system, \
    proportional_allocation_by_location_and_activity, filter_by_geoscale
from flowsa.dataclean import harmonize_units, clean_df
from flowsa.mapping import add_sectors_to_flowbyactivity
from flowsa.data_source_scripts.BLS_QCEW import clean_bls_qcew_fba

# Read pdf into list of DataFrame
def bh_call(url, bh_response, args):
    """
    Load the data in pdf format and produce a response
    :param url: url where data found
    :param bh_response: response to reading url
    :param args: arguments as specified in flowbyactivity.py ('year' and 'source')
    :return: dataframe to format into a flowbyactivty data set
    """

    pages = range(5, 13)
    bh_df_list = []
    for x in pages:
        bh_df = tabula.read_pdf(io.BytesIO(bh_response.content), pages=x, stream=True)[0]
        bh_df_list.append(bh_df)

    bh_df = pd.concat(bh_df_list, sort=False)

    return bh_df

def bh_parse(dataframe_list, args):
    """
    Functions to being parsing and formatting data into flowbyactivity format
    :param dataframe_list: list of dataframes to concat and format
    :param args: arguments as specified in flowbyactivity.py ('year' and 'source')
    :return: dataframe parsed and partially formatted to flowbyactivity specifications
    """
    # concat list of dataframes (info on each page)
    df = pd.concat(dataframe_list, sort=False)
    df = df.rename(columns={"I-O code": "ActivityConsumedBy",
                            "I-O description": "Description",
                            "gal/$M": "FlowAmount",
                            })
    # hardcode
    df.loc[:, 'FlowAmount'] = df['FlowAmount'] / 1000000  # original data in gal/million usd
    df['Unit'] = 'gal/USD'
    df['SourceName'] = 'Blackhurst_IO'
    df['Class'] = 'Water'
    df['FlowName'] = 'Water Withdrawals IO Vector'
    df['Location'] = US_FIPS
    df = assign_fips_location_system(df, '2002')
    df['Year'] = '2002'

    return df

def convert_blackhurst_data_to_gal_per_year(df, attr):
    """
    Load BEA Make After Redefinition data to convert dataframe units
    to gallon per year
    :param df: fba ataframe
    :param attr: attribute data from method yaml
    :return: transformed fba df
    """
    # import flowsa
    # from flowsa.common import fba_fill_na_dict
    # from flowsa.dataclean import harmonize_units
    # from flowsa.dataclean import clean_df

    # load the bea make table
    bmt = flowsa.getFlowByActivity(datasource='BEA_Make_AR',
                                   year=2002, flowclass='Money')
    # clean df
    bmt = clean_df(bmt, flow_by_activity_fields, fba_fill_na_dict)
    bmt = harmonize_units(bmt)
    # drop rows with flowamount = 0
    bmt = bmt[bmt['FlowAmount'] != 0]

    bh_df_revised = pd.merge(df, bmt[['FlowAmount', 'ActivityProducedBy', 'Location']],
                             left_on=['ActivityConsumedBy', 'Location'],
                             right_on=['ActivityProducedBy', 'Location']
                             )

    bh_df_revised.loc[:, 'FlowAmount'] = ((bh_df_revised['FlowAmount_x']) *
                                          (bh_df_revised['FlowAmount_y']))
    bh_df_revised.loc[:, 'Unit'] = 'gal'
    # drop columns
    bh_df_revised = bh_df_revised.drop(columns=["FlowAmount_x", "FlowAmount_y",
                                                'ActivityProducedBy_y'])
    bh_df_revised = bh_df_revised.rename(columns={"ActivityProducedBy_x":
                                                      "ActivityProducedBy"})

    return bh_df_revised

def convert_blackhurst_data_to_gal_per_employee(df_wsec, attr, method):
    """
    Load BLS employment data and use to transporm original units to gallons per employee
    :param df_wsec: dataframe that includes sector columns
    :param attr: attribute data from fba method yaml
    :param method: The name of the fba method yaml
    :return: transformed fba dataframe with sector columns
    """

    # import flowsa
    # from flowsa.mapping import add_sectors_to_flowbyactivity
    # from flowsa.flowbyfunctions import filter_by_geoscale
    # from flowsa.common import fba_fill_na_dict
    # from flowsa.dataclean import harmonize_units
    # from flowsa.dataclean import clean_df
    # from flowsa.data_source_scripts.BLS_QCEW import clean_bls_qcew_fba

    bls = flowsa.getFlowByActivity(datasource='BLS_QCEW', year=2002, flowclass='Employment')

    bls = filter_by_geoscale(bls, 'national')

    # clean df
    bls = clean_df(bls, flow_by_activity_fields, fba_fill_na_dict)
    bls = harmonize_units(bls)
    bls = clean_bls_qcew_fba(bls, attr=attr)

    # assign naics to allocation dataset
    bls_wsec = add_sectors_to_flowbyactivity(bls, sectorsourcename=method['target_sector_source'])
    # drop rows where sector = None ( does not occur with mining)
    bls_wsec = bls_wsec[~bls_wsec['SectorProducedBy'].isnull()]
    bls_wsec = bls_wsec.rename(columns={'SectorProducedBy': 'Sector',
                                        'FlowAmount': 'HelperFlow'})

    # merge the two dfs
    df = pd.merge(df_wsec,
                  bls_wsec[['Location', 'Sector', 'HelperFlow']],
                  how='left',
                  left_on=['Location', 'SectorConsumedBy'],
                  right_on=['Location', 'Sector'])
    # drop any rows where sector is None
    df = df[~df['Sector'].isnull()]
    # fill helperflow values with 0
    df['HelperFlow'] = df['HelperFlow'].fillna(0)

    # calculate proportional ratios
    df_wratio = proportional_allocation_by_location_and_activity(df, 'Sector')

    df_wratio = df_wratio.rename(columns={'FlowAmountRatio': 'EmployeeRatio',
                                          'HelperFlow': 'Employees'})

    # drop rows where helperflow = 0
    df_wratio = df_wratio[df_wratio['Employees'] != 0]

    # calculate gal/employee in 2002
    df_wratio.loc[:, 'FlowAmount'] = (df_wratio['FlowAmount'] *
                                      df_wratio['EmployeeRatio']) / df_wratio['Employees']
    df_wratio.loc[:, 'Unit'] = 'gal/employee'

    # drop cols
    df_wratio = df_wratio.drop(columns=['Sector', 'Employees', 'EmployeeRatio'])

    return df_wratio

def scale_blackhurst_results_to_usgs_values(df_to_scale, attr):
    """
    Scale the initial estimates for Blackhurst-based mining estimates to USGS values.
    Oil-based sectors are allocated a larger percentage of the difference between initial
    water withdrawal estimates and published USGS values.

    This method is based off the Water Satellite Table created by Yang and Ingwersen, 2017
    :param df_to_scale: fba dataframe to be modified
    :param attr: attributes from the fba method yaml
    :return: scaled fba results
    """

    # import flowsa
    # from flowsa.dataclean import harmonize_units

    # determine national level published withdrawal data for usgs mining in FBS method year
    pv_load = flowsa.getFlowByActivity(datasource="USGS_NWIS_WU",
                                       year=str(attr['helper_source_year']),
                                       flowclass='Water'
                                       )
    pv_load = harmonize_units(pv_load)
    pv_sub = pv_load[(pv_load['Location'] == str(US_FIPS)) &
                     (pv_load['ActivityConsumedBy'] == 'Mining')].reset_index(drop=True)
    pv = pv_sub['FlowAmount'].loc[0] * 1000000  # usgs unit is Mgal, blackhurst unit is gal

    # sum quantity of water withdrawals already allocated to sectors
    av = df_to_scale['FlowAmount'].sum()

    # calculate the difference between published value and allocated value
    vd = pv - av

    # subset df to scale into oil and non-oil sectors
    df_to_scale['sector_label'] = np.where(
        df_to_scale['SectorConsumedBy'].apply(lambda x: x[0:5] == '21111'), 'oil','nonoil')
    df_to_scale['ratio'] = np.where(df_to_scale['sector_label'] == 'oil', 2 / 3, 1 / 3)
    df_to_scale['label_sum'] = df_to_scale.groupby(['Location',
                                                    'sector_label'])['FlowAmount'].transform('sum')
    df_to_scale.loc[:, 'value_difference'] = vd.astype(float)

    # calculate revised water withdrawal allocation
    df_scaled = df_to_scale.copy()
    df_scaled.loc[:, 'FlowAmount'] = df_scaled['FlowAmount'] + \
                                     (df_scaled['FlowAmount'] / df_scaled['label_sum']) * \
                                     (df_scaled['ratio'] * df_scaled['value_difference'])
    df_scaled = df_scaled.drop(columns=['sector_label', 'ratio', 'label_sum', 'value_difference'])

    return df_scaled

max_stars_repo_name: JohnAndrewTaylor/flowsa

max_issues_repo_name: JohnAndrewTaylor/flowsa

max_forks_repo_name: JohnAndrewTaylor/flowsa


text:
# The MIT License (MIT)
# Copyright (c) 2019 by the xcube development team and contributors
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of
# this software and associated documentation files (the "Software"), to deal in
# the Software without restriction, including without limitation the rights to
# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
# of the Software, and to permit persons to whom the Software is furnished to do
# so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
import datetime
from abc import ABCMeta
from typing import Tuple, Dict, Any

import numpy as np
import pandas as pd
import xarray as xr

from xcube.constants import CRS_WKT_EPSG_4326
from xcube.core.gen.iproc import ReprojectionInfo, XYInputProcessor, _normalize_lon_360
from xcube.core.timecoord import to_time_in_days_since_1970
from .transexpr import translate_snap_expr_attributes
from .vectorize import new_band_coord_var, vectorize_wavebands

class SnapNetcdfInputProcessor(XYInputProcessor, metaclass=ABCMeta):
    """
    Input processor for SNAP L2 NetCDF inputs.
    """

    @property
    def default_parameters(self) -> Dict[str, Any]:
        default_parameters = super().default_parameters
        default_parameters.update(input_reader='netcdf4',
                                  input_reader_params=dict(decode_cf=True,
                                                           decode_coords=True,
                                                           decode_times=False),
                                  xy_names=('lon', 'lat'),
                                  xy_tp_names=('TP_longitude', 'TP_latitude'),
                                  xy_crs=CRS_WKT_EPSG_4326)
        return default_parameters

    def get_time_range(self, dataset: xr.Dataset) -> Tuple[float, float]:

        if "time_coverage_start" in dataset.attrs:
            t1 = str(dataset.attrs["time_coverage_start"])
            t2 = str(dataset.attrs.get("time_coverage_end", t1))
        else:
            t1 = dataset.attrs.get('start_date')
            t2 = dataset.attrs.get('stop_date', t1)
        if t1 is None or t2 is None:
            raise ValueError('illegal L2 input: missing start/stop time')
        t1 = to_time_in_days_since_1970(t1)
        t2 = to_time_in_days_since_1970(t2)
        return t1, t2

    def pre_process(self, dataset: xr.Dataset) -> xr.Dataset:
        """ Do any pre-processing before reprojection. """
        return translate_snap_expr_attributes(dataset)

    def post_process(self, dataset: xr.Dataset) -> xr.Dataset:
        def new_band_coord_var_ex(band_dim_name: str, band_values: np.ndarray) -> xr.DataArray:
            # Bug in HIGHROC OLCI L2 data: both bands 20 and 21 have wavelengths at 940 nm
            if band_values[-2] == band_values[-1] and band_values[-1] == 940.:
                band_values[-1] = 1020.
            return new_band_coord_var(band_dim_name, band_values)

        return vectorize_wavebands(dataset, new_band_coord_var_ex)

# noinspection PyAbstractClass
class SnapOlciHighrocL2InputProcessor(SnapNetcdfInputProcessor):
    """
    Input processor for SNAP Sentinel-3 OLCI HIGHROC Level-2 NetCDF inputs.
    """

    def __init__(self, **parameters):
        super().__init__('snap-olci-highroc-l2', **parameters)

# noinspection PyAbstractClass
class SnapOlciCyanoAlertL2InputProcessor(SnapNetcdfInputProcessor):
    """
    Input processor for SNAP Sentinel-3 OLCI CyanoAlert Level-2 NetCDF inputs.
    """

    def __init__(self, **parameters):
        super().__init__('snap-olci-cyanoalert-l2', **parameters)

class CMEMSInputProcessor(XYInputProcessor):
    """
    CMEMS input processor that expects input datasets that do not have time bounds:

    * Have dimensions ``lat``, ``lon`` and ``time`` of length 1;
    * have coordinate variables ``lat[lat]``, ``lon[lat]``, ``time[time]``;
    * have coordinate variables ``lat[lat]``, ``lon[lat]`` as decimal degrees on WGS84 ellipsoid,
      both linearly increasing with same constant delta;
    * have coordinate variable ``time[time]`` representing a date+time value object in cftime.DatetimeGregorian;
    * have any data variables of form ``<var>[time, lat, lon]``;

    The CMEMS input processor can be configured by the following parameters:

    * ``input_reader`` the input reader identifier, default is "netcdf4".

    """

    def __init__(self, **parameters):
        super().__init__('cmems', **parameters)

    @property
    def default_parameters(self) -> Dict[str, Any]:
        default_parameters = super().default_parameters
        default_parameters.update(input_reader='netcdf4',
                                  xy_crs=CRS_WKT_EPSG_4326)
        return default_parameters

    def get_reprojection_info(self, dataset: xr.Dataset) -> ReprojectionInfo:
        return super().get_reprojection_info(dataset).derive(
            xy_gcp_step=(max(1, len(dataset.lon) // 4),
                         max(1, len(dataset.lat) // 4))
        )

    def pre_process(self, dataset: xr.Dataset) -> xr.Dataset:
        if 'longitude' in dataset.dims:
            dataset = dataset.rename(({'longitude': 'lon'}))

        if 'latitude' in dataset.dims:
            dataset = dataset.rename(({'latitude': 'lat'}))

        self._validate(dataset)

        if "time" in dataset:
            # Remove time dimension of length 1.
            dataset = dataset.squeeze("time")

        return _normalize_lon_360(dataset)

    def get_time_range(self, dataset: xr.Dataset) -> Tuple[float, float]:
        time_coverage_start, time_coverage_end = None, None
        date_format = '%Y-%m-%d %H:%M:%S'
        if "time" in dataset:
            time_coverage_start = str(dataset.time[0].values)
            date = pd.to_datetime(str(dataset.time[0].values), utc=True)
            if 'daily' in str(dataset.attrs.values()):
                time_coverage_end = datetime.datetime.strftime((date + datetime.timedelta(days=1)), date_format)
            if 'hourly' in str(dataset.attrs.values()):
                time_coverage_end = datetime.datetime.strftime((date + datetime.timedelta(hours=1)), date_format)

        return to_time_in_days_since_1970(time_coverage_start), to_time_in_days_since_1970(time_coverage_end)

    def _validate(self, dataset):
        self._check_coordinate_var(dataset, "lon", min_length=2)
        self._check_coordinate_var(dataset, "lat", min_length=2)
        if "time" in dataset.dims:
            self._check_coordinate_var(dataset, "time", max_length=1)
            required_dims = ("time", "lat", "lon")
        else:
            required_dims = ("lat", "lon")
        count = 0
        for var_name in dataset.data_vars:
            var = dataset.data_vars[var_name]
            if var.dims == required_dims:
                count += 1
        if count == 0:
            raise ValueError(f"dataset has no variables with required dimensions {required_dims!r}")

    # noinspection PyMethodMayBeStatic
    def _check_coordinate_var(self, dataset: xr.Dataset, coord_var_name: str,
                              min_length: int = None, max_length: int = None):
        if coord_var_name not in dataset.coords:
            raise ValueError(f'missing coordinate variable "{coord_var_name}"')
        coord_var = dataset.coords[coord_var_name]
        if len(coord_var.shape) != 1:
            raise ValueError('coordinate variable "lon" must be 1D')
        coord_var_bnds_name = coord_var.attrs.get("bounds", coord_var_name + "_bnds")
        if coord_var_bnds_name in dataset:
            coord_bnds_var = dataset[coord_var_bnds_name]
            expected_shape = (len(coord_var), 2)
            if coord_bnds_var.shape != expected_shape:
                raise ValueError(f'coordinate bounds variable "{coord_bnds_var}" must have shape {expected_shape!r}')
        else:
            if min_length is not None and len(coord_var) < min_length:
                raise ValueError(f'coordinate variable "{coord_var_name}" must have at least {min_length} value(s)')
            if max_length is not None and len(coord_var) > max_length:
                raise ValueError(f'coordinate variable "{coord_var_name}" must have no more than {max_length} value(s)')


text:
"""
 Author: Matt Hanson
 Created: 18/12/2020 9:30 AM
 """

from Pasture_Growth_Modelling.initialisation_support.inital_long_term_runs import run_past_basgra_irrigated, \
    run_past_basgra_dryland, plot_multiple_results
import pandas as pd
import numpy as np

if __name__ == '__main__':
    reseed=True
    irr_ox = run_past_basgra_irrigated(site='oxford',reseed=reseed)
    eyrewell = run_past_basgra_irrigated(reseed=reseed)
    dry_ox = run_past_basgra_dryland(site='oxford', reseed=reseed)
    plot_multiple_results({'irr_ox':irr_ox, 'irr_eyre':eyrewell, 'dry_ox':dry_ox},
                          out_vars=['BASAL'], rolling=30, main_kwargs={'alpha': 0.2}, label_main=False,
                          label_rolling=True, show=True)
    t = dry_ox
    t.loc[t.doy == 152].BASAL.describe(percentiles=np.arange(5, 55, 5) / 100)


text:
"""
Copyright 2020 The OneFlow Authors. All rights reserved.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""
from __future__ import absolute_import

import os
import unittest

import numpy as np
import oneflow as flow
import oneflow.typing as oft

func_config = flow.FunctionConfig()
func_config.default_data_type(flow.float)

def multi_square_sum(
    x, name=None,
):

    return (
        flow.user_op_builder(name if name is not None else "MultiSquareSum")
        .Op("multi_square_sum")
        .Input("x", x)
        .Output("y")
        .Build()
        .InferAndTryRun()
        .RemoteBlobList()[0]
    )

def _check(test_case, xs, y):
    ref_y = np.sum(np.array([np.sum(x ** 2) for x in xs]))
    test_case.assertTrue(np.allclose(y, ref_y))

def _run_test(test_case, x, n, dtype, device):
    flow.clear_default_session()

    @flow.global_function(function_config=func_config)
    def multi_square_sum_job(x: oft.Numpy.Placeholder(x.shape, dtype=dtype)):
        with flow.scope.placement(device, "0:0"):
            xs = [x + 0.1 * i for i in range(n)]
            return multi_square_sum(xs)

    y = multi_square_sum_job(x).get()
    _check(test_case, [(x + 0.1 * i).astype(np.float32) for i in range(n)], y.numpy())

@flow.unittest.skip_unless_1n1d()
class TestMultiSquareSum(flow.unittest.TestCase):
    @unittest.skipIf(os.getenv("ONEFLOW_TEST_CPU_ONLY"), "only test cpu cases")
    def test_multi_square_sum_random_gpu(test_case):
        x = np.random.rand(3, 4, 5).astype(np.float32)
        _run_test(test_case, x, 5, flow.float32, "gpu")
        _run_test(test_case, x, 5, flow.float32, "gpu")
        _run_test(test_case, x, 88, flow.float32, "gpu")
        _run_test(test_case, x, 64, flow.float32, "gpu")

    def test_multi_square_sum_random_cpu(test_case):
        x = np.random.rand(3, 4, 5).astype(np.float32)
        _run_test(test_case, x, 5, flow.float32, "cpu")

if __name__ == "__main__":
    unittest.main()

max_issues_repo_name: duijiudanggecl/oneflow

max_forks_repo_name: duijiudanggecl/oneflow


text:
import numpy as np
import scipy
import matplotlib.pyplot as plt

def plot_gaussians(ax, s, mu):
    """
        :param X            : 2D array of our dataset
        :param y            : 1D array of the groundtruth labels of the dataset
        :param ax          : existing subplot, to draw on it
    """

    k = s.shape[0]
    color_list = ['r', 'g', 'b', 'c', 'm', 'y']
    N_sides = 11

    # Iterate over all gaussians
    for k_cnt in range(k):

        # pick the covariance matrix and mean values of each gaussian
        s_k = s[k_cnt]
        mu_k = mu[k_cnt]

        if s_k.shape[0]==2:
            # dataset with 2 features
            side_range = np.arange(0, N_sides, 1)
            theta = side_range/(N_sides-1)*2*np.pi
            matrix = np.array([np.cos(theta), np.sin(theta)])
            mu_k = np.expand_dims(mu_k, axis=1)
            mu_repeated = np.repeat(mu_k, N_sides, axis=1)
            epoints = np.matmul(scipy.linalg.sqrtm(s_k), matrix) + mu_repeated
            plt.plot(epoints[0,:], epoints[1,:], color=color_list[k_cnt], linewidth=2)
            plt.scatter(mu_k[0], mu_k[1], color=color_list[k_cnt], marker='o', linewidth=2)
        elif s_k.shape[0]==3:
            # dataset with 3 features
            side_range = np.arange(0, N_sides, 1)
            theta = side_range/(N_sides-1)*np.pi
            phi = side_range/(N_sides-1)*2*np.pi
            sin_theta = np.expand_dims(np.sin(theta), axis=1)
            cos_theta = np.expand_dims(np.cos(theta), axis=1)
            sin_phi = np.expand_dims(np.sin(phi), axis=1)
            cos_phi = np.expand_dims(np.cos(phi), axis=1)
            sx = np.matmul(sin_theta, cos_phi.transpose())
            sy = np.matmul(sin_theta, sin_phi.transpose())
            sz = np.matmul(cos_theta, np.ones((1, N_sides)))
            svect = np.array([sx.reshape((N_sides*N_sides)), sy.reshape((N_sides*N_sides)), sz.reshape((N_sides*N_sides))])

            mu_k = np.expand_dims(mu_k, axis=1)
            mu_repeated = np.repeat(mu_k, N_sides*N_sides, axis=1)
            epoints = np.matmul(scipy.linalg.sqrtm(s_k), svect) + mu_repeated

            ex = epoints[0]
            ey = epoints[1]
            ez = epoints[2]

            ax.plot3D(epoints[0,:], epoints[1,:], epoints[2,:], color=color_list[k_cnt])
            ax.plot3D(ex, ey, ez, color=color_list[k_cnt])
        else:
            print('Each dataset sample should have either 2 or 3 features. Not plotting this one.')

        # Pause for a short time, to allow the plotted points to be shown in the figure
        plt.pause(0.001)

max_stars_repo_name: ShellySrivastava/Machine-Learning

max_issues_repo_name: ShellySrivastava/Machine-Learning

max_forks_repo_name: ShellySrivastava/Machine-Learning