text:
''' Some tests for filters '''

import numpy as np

from numpy.testing import assert_equal, assert_raises

import scipy.ndimage as sndi

def test_ticket_701():
    # Test generic filter sizes
    arr = np.arange(4).reshape((2,2))
    func = lambda x: np.min(x)
    res = sndi.generic_filter(arr, func, size=(1,1))
    # The following raises an error unless ticket 701 is fixed
    res2 = sndi.generic_filter(arr, func, size=1)
    assert_equal(res, res2)

def test_orders_gauss():
    # Check order inputs to Gaussians
    arr = np.zeros((1,))
    yield assert_equal, 0, sndi.gaussian_filter(arr, 1, order=0)
    yield assert_equal, 0, sndi.gaussian_filter(arr, 1, order=3)
    yield assert_raises, ValueError, sndi.gaussian_filter, arr, 1, -1
    yield assert_raises, ValueError, sndi.gaussian_filter, arr, 1, 4
    yield assert_equal, 0, sndi.gaussian_filter1d(arr, 1, axis=-1, order=0)
    yield assert_equal, 0, sndi.gaussian_filter1d(arr, 1, axis=-1, order=3)
    yield assert_raises, ValueError, sndi.gaussian_filter1d, arr, 1, -1, -1
    yield assert_raises, ValueError, sndi.gaussian_filter1d, arr, 1, -1, 4

def test_valid_origins():
    """Regression test for #1311."""
    func = lambda x: np.mean(x)
    data = np.array([1,2,3,4,5], dtype=np.float64)
    assert_raises(ValueError, sndi.generic_filter, data, func, size=3,
                  origin=2)
    func2 = lambda x, y: np.mean(x + y)
    assert_raises(ValueError, sndi.generic_filter1d, data, func,
                  filter_size=3, origin=2)
    assert_raises(ValueError, sndi.percentile_filter, data, 0.2, size=3,
                  origin=2)

    for filter in [sndi.uniform_filter, sndi.minimum_filter,
                   sndi.maximum_filter, sndi.maximum_filter1d,
                   sndi.median_filter, sndi.minimum_filter1d]:
        # This should work, since for size == 3, the valid range for origin is
        # -1 to 1.
        filter(data, 3, origin=-1)
        filter(data, 3, origin=1)
        # Just check this raises an error instead of silently accepting or
        # segfaulting.
        assert_raises(ValueError, filter, data, 3, origin=2)

max_stars_repo_name: seberg/scipy

max_issues_repo_name: seberg/scipy

max_forks_repo_name: seberg/scipy

reason: import numpy,from numpy,import scipy


text:
from PJLink.HelperClasses import *

class ImageData:
    """Holder class for Image data

    """

    def __init__(self, dims, cs, itype, data):
        self.dimensions = dims[:2]
        self.data = data
        self.color_space = cs
        self.image_type = itype
        self.channels = dims[2] if len(dims) > 2 else 1

    @staticmethod
    def _detmode(cs, itype, channels):

        mode = cs
        dtype = None
        if itype == "Bit":
            mode = "1"
            dtype = "uint8"
        elif itype == "Byte":
            dtype = "uint8"
            if cs == "Grayscale":
                mode = "L"
            elif cs == "RGB":
                mode = "RGB"
                if channels == 4:
                    mode = "RGBA"
            if cs == "HSB":
                mode = "HSV"
        elif itype == "Bit16":
            dtype = "int32"
            mode = "I"
        elif itype == "Real" or itype == "Real32":
            mode = "F"
            dtype = "float32"

        return mode, dtype

    _mode_cs_map = {
        "1" : "Grayscale",
        "L" : "Grayscale",
        "RGB" : "RGB",
        "RGBA" : "RGB",
        "HSV"  : "HSB",
        "CMYK" : "CMYK",
        "LAB"  : "LAB",
        "F"    : "RGB",
        "I"    : "RGB"
    }

    _mode_it_map = {
        "1" : "Bit",
        "L" : "Byte",
        "RGB" : "Byte",
        "RGBA" : "Byte",
        "HSV"  : "Byte",
        "CMYK" : "Byte",
        "LAB"  : "Byte",
        "F"    : "Real32",
        "I"    : "Bit16"
    }

    @classmethod
    def _invmode(cls, mode):
        cs = cls._mode_cs_map[mode]
        itype = cls._mode_it_map[mode]

        return cs, itype

    @classmethod
    def _topil(cls, dims, cs, itype, data, channels):
        from PIL import Image
        mode, dtype = cls._detmode(cs, itype, channels)

        try:
            cast = data.astype(dtype)
        except AttributeError:
            cast = data

        try:
            dbuf = cast.data.tobytes()
        except AttributeError:
            try:
                dbuf = cast.tobytes()
            except AttributeError:
                dbuf = cast.data

        return Image.frombuffer(mode, dims, dbuf, "raw", mode, 0, 1)

    @property
    def pil_mode(self):
        return self._detmode(self.color_space, self.image_type, self.channels)

    @classmethod
    def frompil(cls, img):
        """frompil constructs a ImageData object from a PIL image

        :param img:
        :return:
        """

        cs, itype = cls._invmode(img.mode)

        if Env.HAS_NUMPY:
            import numpy as np
            dats = np.asarray(img)
        else:
            import array

            array_data = img.__array_interface__
            # to_type_code = array_data["typestr"]
            # bom = to_type_code[0]
            # form = to_type_code[1]
            # num_bits = to_type_code[2]
            if itype == "Bit" or itype == "Byte":
                tc = "B"
            elif itype == "Bit16":
                tc = "h"
            elif itype == "Real32":
                tc = "f"
            dbuf = array.array(tc, array_data["data"])
            dats = BufferedNDArray(dbuf, array_data["shape"])

        return cls(dats.shape, cs, itype, dats)

    def topil(self):
        return self._topil(tuple(self.dimensions), self.color_space, self.image_type, self.data, self.channels)

    @property
    def expr(self):
        return MPackage.Image(
            self.data,
            self.image_type,
            ColorSpace_ = self.color_space
        )

def _get_type(head, link, stack):
    return stack["dtype"]
def _get_dims(head, link, stack):
    return stack["dims"]
def _get_image(name, *items, ImageData=ImageData):
    odict = dict(items)
    return ImageData(odict["dims"], odict["color_space"], odict["bit_size"], odict["data"])

ImageArrayInfoDecoder = (
    "image",             #object name
    "ImageArrayInfo",  #object head
    # decoder tuples
    # they come as (key_name, head, typename, dim_list )
    ("dims",        (None, "Integer", [ 1 ])), #_getArray works by depth, not true dimension, if the list isn't [ 0 ]
    ("color_space", (None, "String", None)),
    ("bit_size",    (None, "String", None)),
    ("dtype",       (None, "String", None)),
    ("data",        (None, _get_type, _get_dims)),
    _get_image
)

_decoder = ImageArrayInfoDecoder

ext: py

reason: import numpy


text:
# Import common libraries
import sys
import gym
import numpy as np
import matplotlib.pyplot as plt

# Set plotting options
%matplotlib inline
plt.style.use('ggplot')
np.set_printoptions(precision=3, linewidth=120)
# Create an environment
env = gym.make('Acrobot-v1')
env.seed(505);

# Explore state (observation) space
print("State space:", env.observation_space)
print("- low:", env.observation_space.low)
print("- high:", env.observation_space.high)

# Explore action space
print("Action space:", env.action_space)
def create_tiling_grid(low, high, bins=(10, 10), offsets=(0.0, 0.0)):
    """Define a uniformly-spaced grid that can be used for tile-coding a space.

    Parameters
    ----------
    low : array_like
        Lower bounds for each dimension of the continuous space.
    high : array_like
        Upper bounds for each dimension of the continuous space.
    bins : tuple
        Number of bins or tiles along each corresponding dimension.
    offsets : tuple
        Split points for each dimension should be offset by these values.

    Returns
    -------
    grid : list of array_like
        A list of arrays containing split points for each dimension.
    """
    # TODO: Implement this
    pass

low = [-1.0, -5.0]
high = [1.0, 5.0]
create_tiling_grid(low, high, bins=(10, 10), offsets=(-0.1, 0.5))  # [test]
def create_tilings(low, high, tiling_specs):
    """Define multiple tilings using the provided specifications.

    Parameters
    ----------
    low : array_like
        Lower bounds for each dimension of the continuous space.
    high : array_like
        Upper bounds for each dimension of the continuous space.
    tiling_specs : list of tuples
        A sequence of (bins, offsets) to be passed to create_tiling_grid().

    Returns
    -------
    tilings : list
        A list of tilings (grids), each produced by create_tiling_grid().
    """
    # TODO: Implement this
    pass

# Tiling specs: [(<bins>, <offsets>), ...]
tiling_specs = [((10, 10), (-0.066, -0.33)),
                ((10, 10), (0.0, 0.0)),
                ((10, 10), (0.066, 0.33))]
tilings = create_tilings(low, high, tiling_specs)
from matplotlib.lines import Line2D

def visualize_tilings(tilings):
    """Plot each tiling as a grid."""
    prop_cycle = plt.rcParams['axes.prop_cycle']
    colors = prop_cycle.by_key()['color']
    linestyles = ['-', '--', ':']
    legend_lines = []

    fig, ax = plt.subplots(figsize=(10, 10))
    for i, grid in enumerate(tilings):
        for x in grid[0]:
            l = ax.axvline(x=x, color=colors[i % len(colors)], linestyle=linestyles[i % len(linestyles)], label=i)
        for y in grid[1]:
            l = ax.axhline(y=y, color=colors[i % len(colors)], linestyle=linestyles[i % len(linestyles)])
        legend_lines.append(l)
    ax.grid('off')
    ax.legend(legend_lines, ["Tiling #{}".format(t) for t in range(len(legend_lines))], facecolor='white', framealpha=0.9)
    ax.set_title("Tilings")
    return ax  # return Axis object to draw on later, if needed

visualize_tilings(tilings);
def discretize(sample, grid):
    """Discretize a sample as per given grid.

    Parameters
    ----------
    sample : array_like
        A single sample from the (original) continuous space.
    grid : list of array_like
        A list of arrays containing split points for each dimension.

    Returns
    -------
    discretized_sample : array_like
        A sequence of integers with the same number of dimensions as sample.
    """
    # TODO: Implement this
    pass

def tile_encode(sample, tilings, flatten=False):
    """Encode given sample using tile-coding.

    Parameters
    ----------
    sample : array_like
        A single sample from the (original) continuous space.
    tilings : list
        A list of tilings (grids), each produced by create_tiling_grid().
    flatten : bool
        If true, flatten the resulting binary arrays into a single long vector.

    Returns
    -------
    encoded_sample : list or array_like
        A list of binary vectors, one for each tiling, or flattened into one.
    """
    # TODO: Implement this
    pass

# Test with some sample values
samples = [(-1.2 , -5.1 ),
           (-0.75,  3.25),
           (-0.5 ,  0.0 ),
           ( 0.25, -1.9 ),
           ( 0.15, -1.75),
           ( 0.75,  2.5 ),
           ( 0.7 , -3.7 ),
           ( 1.0 ,  5.0 )]
encoded_samples = [tile_encode(sample, tilings) for sample in samples]
print("\nSamples:", repr(samples), sep="\n")
print("\nEncoded samples:", repr(encoded_samples), sep="\n")
from matplotlib.patches import Rectangle

def visualize_encoded_samples(samples, encoded_samples, tilings, low=None, high=None):
    """Visualize samples by activating the respective tiles."""
    samples = np.array(samples)  # for ease of indexing

    # Show tiling grids
    ax = visualize_tilings(tilings)

    # If bounds (low, high) are specified, use them to set axis limits
    if low is not None and high is not None:
        ax.set_xlim(low[0], high[0])
        ax.set_ylim(low[1], high[1])
    else:
        # Pre-render (invisible) samples to automatically set reasonable axis limits, and use them as (low, high)
        ax.plot(samples[:, 0], samples[:, 1], 'o', alpha=0.0)
        low = [ax.get_xlim()[0], ax.get_ylim()[0]]
        high = [ax.get_xlim()[1], ax.get_ylim()[1]]

    # Map each encoded sample (which is really a list of indices) to the corresponding tiles it belongs to
    tilings_extended = [np.hstack((np.array([low]).T, grid, np.array([high]).T)) for grid in tilings]  # add low and high ends
    tile_centers = [(grid_extended[:, 1:] + grid_extended[:, :-1]) / 2 for grid_extended in tilings_extended]  # compute center of each tile
    tile_toplefts = [grid_extended[:, :-1] for grid_extended in tilings_extended]  # compute topleft of each tile
    tile_bottomrights = [grid_extended[:, 1:] for grid_extended in tilings_extended]  # compute bottomright of each tile

    prop_cycle = plt.rcParams['axes.prop_cycle']
    colors = prop_cycle.by_key()['color']
    for sample, encoded_sample in zip(samples, encoded_samples):
        for i, tile in enumerate(encoded_sample):
            # Shade the entire tile with a rectangle
            topleft = tile_toplefts[i][0][tile[0]], tile_toplefts[i][1][tile[1]]
            bottomright = tile_bottomrights[i][0][tile[0]], tile_bottomrights[i][1][tile[1]]
            ax.add_patch(Rectangle(topleft, bottomright[0] - topleft[0], bottomright[1] - topleft[1],
                                   color=colors[i], alpha=0.33))

            # In case sample is outside tile bounds, it may not have been highlighted properly
            if any(sample < topleft) or any(sample > bottomright):
                # So plot a point in the center of the tile and draw a connecting line
                cx, cy = tile_centers[i][0][tile[0]], tile_centers[i][1][tile[1]]
                ax.add_line(Line2D([sample[0], cx], [sample[1], cy], color=colors[i]))
                ax.plot(cx, cy, 's', color=colors[i])

    # Finally, plot original samples
    ax.plot(samples[:, 0], samples[:, 1], 'o', color='r')

    ax.margins(x=0, y=0)  # remove unnecessary margins
    ax.set_title("Tile-encoded samples")
    return ax

visualize_encoded_samples(samples, encoded_samples, tilings);
class QTable:
    """Simple Q-table."""

    def __init__(self, state_size, action_size):
        """Initialize Q-table.

        Parameters
        ----------
        state_size : tuple
            Number of discrete values along each dimension of state space.
        action_size : int
            Number of discrete actions in action space.
        """
        self.state_size = state_size
        self.action_size = action_size

        # TODO: Create Q-table, initialize all Q-values to zero
        # Note: If state_size = (9, 9), action_size = 2, q_table.shape should be (9, 9, 2)

        print("QTable(): size =", self.q_table.shape)

class TiledQTable:
    """Composite Q-table with an internal tile coding scheme."""

    def __init__(self, low, high, tiling_specs, action_size):
        """Create tilings and initialize internal Q-table(s).

        Parameters
        ----------
        low : array_like
            Lower bounds for each dimension of state space.
        high : array_like
            Upper bounds for each dimension of state space.
        tiling_specs : list of tuples
            A sequence of (bins, offsets) to be passed to create_tilings() along with low, high.
        action_size : int
            Number of discrete actions in action space.
        """
        self.tilings = create_tilings(low, high, tiling_specs)
        self.state_sizes = [tuple(len(splits)+1 for splits in tiling_grid) for tiling_grid in self.tilings]
        self.action_size = action_size
        self.q_tables = [QTable(state_size, self.action_size) for state_size in self.state_sizes]
        print("TiledQTable(): no. of internal tables = ", len(self.q_tables))

    def get(self, state, action):
        """Get Q-value for given <state, action> pair.

        Parameters
        ----------
        state : array_like
            Vector representing the state in the original continuous space.
        action : int
            Index of desired action.

        Returns
        -------
        value : float
            Q-value of given <state, action> pair, averaged from all internal Q-tables.
        """
        # TODO: Encode state to get tile indices

        # TODO: Retrieve q-value for each tiling, and return their average
        pass

    def update(self, state, action, value, alpha=0.1):
        """Soft-update Q-value for given <state, action> pair to value.

        Instead of overwriting Q(state, action) with value, perform soft-update:
            Q(state, action) = alpha * value + (1.0 - alpha) * Q(state, action)

        Parameters
        ----------
        state : array_like
            Vector representing the state in the original continuous space.
        action : int
            Index of desired action.
        value : float
            Desired Q-value for <state, action> pair.
        alpha : float
            Update factor to perform soft-update, in [0.0, 1.0] range.
        """
        # TODO: Encode state to get tile indices

        # TODO: Update q-value for each tiling by update factor alpha
        pass

# Test with a sample Q-table
tq = TiledQTable(low, high, tiling_specs, 2)
s1 = 3; s2 = 4; a = 0; q = 1.0
print("[GET]    Q({}, {}) = {}".format(samples[s1], a, tq.get(samples[s1], a)))  # check value at sample = s1, action = a
print("[UPDATE] Q({}, {}) = {}".format(samples[s2], a, q)); tq.update(samples[s2], a, q)  # update value for sample with some common tile(s)
print("[GET]    Q({}, {}) = {}".format(samples[s1], a, tq.get(samples[s1], a)))  # check value again, should be slightly updated

max_stars_repo_name: odellus/deep-reinforcement-learning

max_issues_repo_name: odellus/deep-reinforcement-learning

max_forks_repo_name: odellus/deep-reinforcement-learning


text:
"""
This is a collection of helper functions to yt.load_sample
"""
import json
import re
from functools import lru_cache
from itertools import chain
from pathlib import Path
from typing import Optional, Union

import pkg_resources

from yt.funcs import mylog
from yt.utilities.on_demand_imports import (
    _pandas as pd,
    _pooch as pooch,
    _requests as requests,
)

num_exp = re.compile(r"\d*(\.\d*)?")
byte_unit_exp = re.compile(r"[KMGT]?B")

def _parse_byte_size(s: str):
    """
    Convert a string size specification to integer byte size.

    This function should be insensitive to case and whitespace.

    It doesn't always return an int, as a temporary measure to deal with missing
    or corrupted data in the registry. This should be fixed in the future.

    Examples
    --------
    # most of the following examples are adapted from
    # https://stackoverflow.com/a/31631711/5489622
    >>> _parse_byte_size(None)
    <NA>
    >>> from numpy import nan
    >>> _parse_byte_size(nan)
    <NA>
    >>> _parse_byte_size("1B")
    1
    >>> _parse_byte_size("1.00 KB")
    1024
    >>> _parse_byte_size("488.28 KB")
    500000
    >>> _parse_byte_size("1.00 MB")
    1049000
    >>> _parse_byte_size("47.68 MB")
    50000000
    >>> _parse_byte_size("1.00 GB")
    1074000000
    >>> _parse_byte_size("4.66 GB")
    5004000000
    >>> _parse_byte_size("1.00 TB")
    1100000000000
    >>> _parse_byte_size("4.55 TB")
    5003000000000
    """
    try:
        s = s.upper()
    except AttributeError:
        # input is not a string (likely a np.nan)
        return pd.NA

    val = float(re.search(num_exp, s).group())
    unit = re.search(byte_unit_exp, s).group()
    prefixes = ["B", "K", "M", "G", "T"]
    raw_res = val * 1024 ** prefixes.index(unit[0])
    return int(float(f"{raw_res:.3e}"))

@lru_cache(maxsize=128)
def get_data_registry_table():
    """
    Load the sample data registry as a pandas.Dataframe instance.

    This function is considered experimental and is exposed for exploratory purposed.
    The output format is subject to change.

    The output of this function is cached so it will only generate one request per session.
    """

    # it would be nicer to have an actual api on the yt website server,
    # but this will do for now
    api_url = "https://raw.githubusercontent.com/yt-project/website/master/data/datafiles.json"  # noqa

    response = requests.get(api_url)

    if not response.ok:
        raise RuntimeError(
            "Could not retrieve registry data. Please check your network setup."
        )

    website_json = response.json()
    # this dict follows this schema: {frontend_name: {flat dataframe-like}}

    columns = ["code", "filename", "size", "url", "description"]
    website_table = pd.concat(pd.DataFrame(d) for d in website_json.values())[columns]

    # add a int-type byte size column
    # note that we cast to pandas specific type "Int64" because we expect missing values
    # see https://pandas.pydata.org/pandas-docs/stable/user_guide/missing_data.html#integer-dtypes-and-missing-data
    website_table["byte_size"] = (
        website_table["size"].apply(_parse_byte_size).astype("Int64")
    )

    # normalize urls to match the local json
    website_table["url"] = website_table["url"].apply(
        lambda u: u.replace("http:", "https:")
    )

    # load local data
    pooch_json = json.load(
        pkg_resources.resource_stream("yt", "sample_data_registry.json")
    )
    pooch_table = pd.DataFrame(pooch_json.values())

    # merge tables
    unified_table = website_table.merge(pooch_table, on="url", how="outer")

    # PR 3089
    # ideally we should be able to do this, but it's not possible
    # at the time of writing because fhe "filename" is incomplete
    # see the companion comment in load_sample
    # unified_table.set_index("filename", inplace=True)
    # unified_table.index.rename("id", inplace=True)
    return unified_table

def _get_test_data_dir_path():
    from yt.config import ytcfg

    return Path(ytcfg.get("yt", "test_data_dir"))

def lookup_on_disk_data(fn):
    """
    Look for data file/dir on disk.

    Returns
    -------

    pathlib.Path to a file/dir matching fn if any

    Raises
    ------
    FileNotFoundError
    """

    path = Path(fn).expanduser()

    if path.exists():
        return path

    alt_path = _get_test_data_dir_path() / fn
    err_msg = f"No such file or directory: '{fn}'."
    if not alt_path.parent.is_dir():
        raise FileNotFoundError(err_msg)

    if not alt_path.exists():
        err_msg += f"\n(Also tried '{alt_path}')."
        raise FileNotFoundError(err_msg)

    return alt_path

@lru_cache(maxsize=128)
def _get_pooch_instance():
    data_registry = get_data_registry_table()
    cache_storage = _get_test_data_dir_path() / "yt_download_cache"

    sample_data_registry = json.load(
        pkg_resources.resource_stream("yt", "sample_data_registry.json")
    )
    registry = {k: v["hash"] for k, v in sample_data_registry.items()}
    return pooch.create(
        path=cache_storage, base_url="https://yt-project.org/data/", registry=registry
    )

def _download_sample_data_file(
    filename, progressbar=True, timeout: Optional[int] = None
):
    """
    Download a file by url.

    Returns
    -------
    storage_filename : location of the downloaded file

    """
    downloader = pooch.HTTPDownloader(progressbar=progressbar, timeout=timeout)

    poochie = _get_pooch_instance()
    poochie.fetch(filename, downloader=downloader)
    return Path.joinpath(poochie.path, filename)

max_stars_repo_name: munkm/yt

max_issues_repo_name: munkm/yt

max_forks_repo_name: munkm/yt

reason: from numpy


text:
"""
Compare commercial models with mean, std, CI over the original test set split.
"""

from ser_evaluation.evaluate import get_test_set_filenames, preprocess_crema, scores
from pathlib import Path
import numpy as np
import pandas as pd
from tqdm import tqdm

def get_metrics(pred_path):
    #bootstrap_path = 'predictions_bootstrap'
    #pred_path = "predictions/commercial_results/results_csv"
    target_intended_observed_ls = []
    model_ls = []
    subset_ls = []

    f1_score = []
    f1_score_anger = []
    f1_score_happy = []
    f1_score_neutral = []
    f1_score_sad = []

    recall = []
    recall_anger = []
    recall_happy = []
    recall_neutral = []
    recall_sad = []

    data_path = Path(pred_path)
    paths = [path for path in data_path.glob('**/*') if path.is_file()]
    paths.sort()

    for path in tqdm(paths):
        model = path.stem
        emotions = [
            'Filename',
            'Anger',
            'Happy',
            'Neutral',
            'Sad']

        pred_df = pd.read_csv(
            path,
            usecols=(emotions)
            )

        col_names = {k: "{}_pred".format(k[:3].upper()) if k != 'Filename'
                     else k.lower() for k in pred_df.columns}
        pred_df.rename(columns = col_names, inplace = True)

        #Convert probabilities to 0s and 1s
        for col in pred_df.drop(['filename'], axis = 1).columns:
            pred_df[col] = np.where(pred_df[col] >= 0.5, 1, 0)

        test_file_name = get_test_set_filenames(path)
        crema = preprocess_crema(demographics_path = 'Data/VideoDemographics.csv',
            ratings_path = 'Data/processedResults/summaryTable.csv',
            test_file_names= test_file_name)

        for key, cat in crema.items():

            subset = key.split("_")[-1]

            if 'observed' in key:
                target_intended_observed = 'observed'
            else:
                target_intended_observed = 'intended'

            cat_pred_df = cat.merge(pred_df, on = 'filename')
            cat_pred_df.columns

            y_true = cat_pred_df[cat_pred_df.filter(regex="{}_".format(target_intended_observed))
                                 .columns]. \
                                drop(columns = ['{}_DIS'.format(target_intended_observed),
                                                '{}_FEA'.format(target_intended_observed)]). \
                                to_numpy()
            y_pred = cat_pred_df[cat_pred_df.filter(regex='pred').columns].to_numpy()

            target_intended_observed_ls.append(target_intended_observed)
            model_ls.append(model)
            subset_ls.append(subset)

            scores_dict = scores(y_true, y_pred)

            f1_score.append(scores_dict['macro avg']['f1-score'])
            f1_score_anger.append(scores_dict['0']['f1-score'])
            f1_score_happy.append(scores_dict['1']['f1-score'])
            f1_score_neutral.append(scores_dict['2']['f1-score'])
            f1_score_sad.append(scores_dict['3']['f1-score'])

            recall.append(scores_dict['macro avg']['recall'])
            recall_anger.append(scores_dict['0']['recall'])
            recall_happy.append(scores_dict['1']['recall'])
            recall_neutral.append(scores_dict['2']['recall'])
            recall_sad.append(scores_dict['3']['recall'])

    f1_df = pd.DataFrame(
        data=zip(
            target_intended_observed_ls,
            model_ls,
            subset_ls,
            f1_score,
            f1_score_anger,
            f1_score_happy,
            f1_score_neutral,
            f1_score_sad,
            recall,
            recall_anger,
            recall_happy,
            recall_neutral,
            recall_sad),
        columns=[
            'target_intended_observed',
            'model',
            'subset',
            'f1_score',
            'f1_score_anger',
            'f1_score_happy',
            'f1_score_neutral',
            'f1_score_sad',
            'recall',
            'recall_anger',
            'recall_happy',
            'recall_neutral',
            'recall_sad'])
    return f1_df

if __name__ == '__main__':
    df = get_metrics("predictions/commercial_results/results_csv")

    output_path = Path('bias_results/commercial/commercial_models.csv')
    output_path.parent.mkdir(parents=True, exist_ok=True)

    df.to_csv(output_path, index = False, float_format='%.5f')

max_stars_repo_name: wwongwk/speech-recognition-bias

max_issues_repo_name: kingsleykuan/speech-recognition-bias

max_forks_repo_name: kingsleykuan/speech-recognition-bias


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

""" Plot probabilities of entity occurrence.

Usage: python plot_fig6_neg_binomial.py
Input data files: ../data/[app_name]_out/complete_user_[app_name].txt, ../data/[app_name]_out/user_[app_name]_all.txt
Time: ~5M
"""

import sys, os, platform
from collections import defaultdict, Counter
import numpy as np
from scipy.special import comb
from scipy import stats

import matplotlib as mpl
if platform.system() == 'Linux':
    mpl.use('Agg')  # no UI backend

import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter

sys.path.append(os.path.join(os.path.dirname(__file__), '../'))
from utils.helper import Timer
from utils.plot_conf import ColorPalette, hide_spines

def negative_binomial(n, k, rho):
    # n: number of trials, k: number of success (being sampled)
    return comb(n-1, k-1, exact=True) * (rho ** k) * ((1 - rho) ** (n - k))

def main():
    timer = Timer()
    timer.start()

    cc4 = ColorPalette.CC4
    blue = cc4[0]

    app_name = 'cyberbullying'
    rho = 0.5272
    entity = 'user'

    fig, axes = plt.subplots(1, 2, figsize=(10, 3.3))

    print('for entity: {0}'.format(entity))
    sample_entity_freq_dict = defaultdict(int)
    with open('../data/{1}_out/{0}_{1}_all.txt'.format(entity, app_name), 'r') as sample_datefile:
        for line in sample_datefile:
            sample_entity_freq_dict[line.rstrip().split(',')[1]] += 1

    complete_entity_freq_dict = defaultdict(int)
    with open('../data/{1}_out/complete_{0}_{1}.txt'.format(entity, app_name), 'r') as complete_datefile:
        for line in complete_datefile:
            complete_entity_freq_dict[line.rstrip().split(',')[1]] += 1

    complete_to_sample_freq_dict = defaultdict(list)
    sample_to_complete_freq_dict = defaultdict(list)

    for item, complete_vol in complete_entity_freq_dict.items():
        if item in sample_entity_freq_dict:
            complete_to_sample_freq_dict[complete_vol].append(sample_entity_freq_dict[item])
        else:
            complete_to_sample_freq_dict[complete_vol].append(0)

    for item, sample_vol in sample_entity_freq_dict.items():
        sample_to_complete_freq_dict[sample_vol].append(complete_entity_freq_dict[item])

    for item in set(complete_entity_freq_dict.keys()) - set(sample_entity_freq_dict.keys()):
        sample_to_complete_freq_dict[0].append(complete_entity_freq_dict[item])

    ax1_x_axis = range(1, 101)

    ax1_y_axis = []
    empirical_mean_list = []
    expected_mean_list = []
    for num_sample in ax1_x_axis:
        # compute sample to complete
        empirical_cnt_dist = sample_to_complete_freq_dict[num_sample]
        neg_binomial_cnt_dist = []
        for x in range(num_sample, max(30, 3 * num_sample + 1)):
            neg_binomial_cnt_dist.extend([x] * int(negative_binomial(x, num_sample, rho) * len(empirical_cnt_dist)))
        ks_test = stats.ks_2samp(empirical_cnt_dist, neg_binomial_cnt_dist)
        empirical_mean = sum(empirical_cnt_dist) / len(empirical_cnt_dist)
        empirical_mean_list.append(empirical_mean)
        expected_mean = sum(neg_binomial_cnt_dist) / len(neg_binomial_cnt_dist)
        expected_mean_list.append(expected_mean)
        print('num_sample: {0}, number of Bernoulli trials: {1}, d_statistic: {2:.4f}, p: {3:.4f}, expected mean: {4:.2f}, empirical mean: {5:.2f}'
              .format(num_sample, len(empirical_cnt_dist), ks_test[0], ks_test[1], expected_mean, empirical_mean))
        ax1_y_axis.append(ks_test[0])

    axes[0].plot(ax1_x_axis, ax1_y_axis, c='k', lw=1.5, ls='-')

    axes[0].set_xlabel(r'sample frequency $n_s$', fontsize=16)
    axes[0].set_ylabel('D-statistic', fontsize=16)
    axes[0].set_xlim([-2, 102])
    axes[0].set_xticks([0, 25, 50, 75, 100])
    axes[0].set_ylim([0, 0.17])
    axes[0].yaxis.set_major_formatter(FuncFormatter(lambda x, _: '{0:.2f}'.format(x)))
    axes[0].tick_params(axis='both', which='major', labelsize=16)
    axes[0].set_title('(a)', fontsize=18, pad=-3*72, y=1.0001)

    # show an example
    num_sample = np.argmin(ax1_y_axis) + 1

    axes[0].scatter(num_sample, ax1_y_axis[num_sample - 1], s=40, c=blue, zorder=30)
    axes[0].set_yticks([0, ax1_y_axis[num_sample - 1], 0.05, 0.1, 0.15])
    axes[0].plot([axes[0].get_xlim()[0], num_sample], [ax1_y_axis[num_sample - 1], ax1_y_axis[num_sample - 1]], color=blue, ls='--', lw=1)
    axes[0].plot([num_sample, num_sample], [axes[0].get_ylim()[0], ax1_y_axis[num_sample - 1]], color=blue, ls='--', lw=1)

    # plot sample to complete
    ax2_x_axis = range(num_sample, max(30, 3 * num_sample + 1))
    num_items = len(sample_to_complete_freq_dict[num_sample])
    sample_to_complete_cnt = Counter(sample_to_complete_freq_dict[num_sample])
    ax2_y_axis = [sample_to_complete_cnt[x] / num_items for x in ax2_x_axis]
    ax2_neg_binomial_axis = [negative_binomial(x, num_sample, rho) for x in ax2_x_axis]

    axes[1].plot(ax2_x_axis, ax2_y_axis, c=blue, lw=1.5, ls='-', marker='o', zorder=20, label='empirical')
    axes[1].plot(ax2_x_axis, ax2_neg_binomial_axis, c='k', lw=1.5, ls='-', marker='x', zorder=10, label='negative binomial')

    axes[1].set_xlabel(r'complete frequency $n_c$', fontsize=16)
    axes[1].set_ylabel(r'Pr($n_c$|$n_s$={0})'.format(num_sample), fontsize=16)
    axes[1].set_xticks([num_sample, 2 * num_sample, 3 * num_sample])
    axes[1].set_ylim([-0.005, 0.15])
    axes[1].set_yticks([0, 0.05, 0.1])
    axes[1].tick_params(axis='both', which='major', labelsize=16)
    axes[1].legend(frameon=False, fontsize=16, ncol=1, fancybox=False, shadow=True, loc='upper left')
    axes[1].set_title('(b)', fontsize=18, pad=-3*72, y=1.0001)

    axes[1].plot([empirical_mean_list[num_sample - 1], empirical_mean_list[num_sample - 1]], [axes[1].get_ylim()[0], 0.1], color=blue, ls='--', lw=1)
    axes[1].plot([expected_mean_list[num_sample - 1], expected_mean_list[num_sample - 1]], [axes[1].get_ylim()[0], 0.1], color='k', ls='--', lw=1)

    hide_spines(axes)

    timer.stop()

    plt.tight_layout(rect=[0, 0.05, 1, 1])
    plt.savefig('../images/entity_negative_binomial.pdf', bbox_inches='tight')
    if not platform.system() == 'Linux':
        plt.show()

if __name__ == '__main__':
    main()

max_stars_repo_name: avalanchesiqi/twitter-sampling

max_issues_repo_name: avalanchesiqi/twitter-sampling

max_forks_repo_name: avalanchesiqi/twitter-sampling

reason: import numpy,from scipy


text:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Constrain the evolution of the atmospheres of the TRAPPIST-1 Planets. Mostly
built of David Fleming's modeling of the stellar evolution.

Rory Barnes
rory@astro.washington.edu
22 Jan 2020

"""

import os
from approxposterior import approx, utility, gpUtils
from scipy.stats import norm
import emcee
import numpy as np
import george
import random
import re
import subprocess
import vplot as vpl

"""

Model parameters (daStateVector) [prior-type, (prior-properties), (bounds), ref]:
1. Stellar mass (dStarMass) [uniform, (0.07,0.11), (0.07,0.11), Fleming20]
2. Log saturation fraction (dFSat) [normal, (-2.92, 0.26), (-5,-1), Wright11,Chadney15]
3. Saturation time (dTSat) [uniform, (0.1,12), (0.1,12), Fleming20]
4. Age (dAge) [normal, (7.6,2.2), (0.1,12), Burgasser17]
5. XUV Exponential decay term (dBeta) [normal, (-1.18,0.31), (-2,0), Jackson12]
6. Fraction of solid body in silicates (dFRock) [uniform, (0.05,0.9), (0.05,0.9)]
7. Planet's initial water mass (dWaterMass) [uniform, (0,0.5), (0,0.5)]
8. Planet's initial envelope mass (dEnvMass) [uniform, (0,10), (0,10)]
9. Efficiency of hydrogen escape (dEscCoeffH) [uniform, (0.1,0.5), (0.1,0.5)]
10. Efficiecny of water photolysis and H escape (dEscCoeffH20) [uniform, (0.06,0.13), (0.06,0.13), Bolmont17]
11. Pressure of XUV absorption (dPressXUV) [uniform, (0.1,10), (0.1,10), LC17]
12. Planetary Albedo (dAlbedo) [uniform, (0.05,0.8), (0.05,0.8)]
13. Initial planetary mass (dPlanetMass) [uniform, (0.5,10), (0.5,10)]
14. Orbital period (dOrbPer) [normal, (6.099043,1.5e-5), (6.098,6.1), Delrez17]

Observational constraints [type, (limits), ref]:
1. Stellar luminosity (dLum) [normal, (5.22e-4,1.9e-5), vanGrootel18]
2. L_bol/L_XUV (dLRatio) [normal, (7.5e-4,1.5e-4), Wheatley17]
3. Stellar Effective Temperature (dTEff)[normal, (2511,37), Delrez17]
4. Current mass of planet e (dPlanetMass) [normal, (0.772,0.077), Grimm18]
5. Radius of planet e (dPlanetRad) [normal, (0.91,0.027), Grimm18]

Additional Parameters (daParams)
1. Luminosity (dLumTrial)
2. L_XUV (dLumXUVTrial)
3. Eff. Temp (dTEffTrial)
4. Period of e (dPerETrial)
5. Mass of e (dMassETrial)
6. Radius of e (dRadETrial)
7. Stellar radius (dStarRad)
8. Planetary equilibrium temp. (dEqTemp)
9. Surface pressure (dSurfPress)
10. Final water mass (dWaterMass)
11. Final H mass fraction (dEnvMass)
12. Atmospheric oxygen pressure (dOxygenPress)
13. Semi-major axis (dSemi)

"""

# Constants
LSUN = 3.846e26         # Solar luminosity in ergs/s
YEARSEC = 3.15576e7     # seconds per year
BIGG = 6.67428e-11      # Universal Gravitational Constant in cgs
DAYSEC = 86400.0        # seconds per day
MSUN = 1.988416e30      # mass of sun in g
AUCM = 1.49598e11       # cm per AU
MTO = 1.4e24            # Mass of all of Earth's water in g
MEarth = 5.972e27       # Mass of Earth in g

def SampleStateVector(iSize=1, **kwargs):

    daStateVector = []
    for iSample in range(iSize):
        while True:
            dStarMassGuess = np.random.uniform(low=dStarMassMin, high=dStarMassMax)
            dFSatGuess = norm.rvs(loc=dFSatMean, scale=dFSatSig, size=1)[0]
            dTSatGuess = np.random.uniform(low=dTSatMin, high=dTSatMax)
            dAgeGuess = norm.rvs(loc=dAgeMean, scale=dAgeSig, size=1)[0]
            dBetaGuess = norm.rvs(loc=dBetaMean, scale=dBetaSig, size=1)[0]

            dFRockGuess = np.random.uniform(low=dFRockMin, high=dFRockMax)
            dWaterMassGuess = np.random.uniform(low=dWaterMassMin, high=dWaterMassMax)
            dEnvMassGuess = np.random.uniform(low=dEnvMassMin, high=dEnvMassMax)
            dEscCoeffHGuess = np.random.uniform(low=dEscCoeffHMin, high=dEscCoeffHMax)
            dEscCoeffH2OGuess = np.random.uniform(low=dEscCoeffH2OMin, high=dEscCoeffH2OMax)
            dPressXUVGuess =  np.random.uniform(low=dPressXUVMin, high=dPressXUVMax)
            dAlbedoGuess = np.random.uniform(low=dAlbedoMin, high=dAlbedoMax)
            dPlanetMassGuess = np.random.uniform(low=dPlanetMassMin, high=dPlanetMassMax)
            dOrbPerGuess = norm.rvs(loc=dOrbPerMean, scale=dOrbPerSig, size=1)[0]

            daGuess = [dStarMassGuess,
                       dFSatGuess,
                       dTSatGuess,
                       dAgeGuess,
                       dBetaGuess,

                       dFRockGuess,
                       dWaterMassGuess,
                       dEnvMassGuess,
                       dEscCoeffHGuess,
                       dEscCoeffH2OGuess,
                       dPressXUVGuess,
                       dAlbedoGuess,
                       dPlanetMassGuess,
                       dOrbPerGuess
                       ]
            if not np.isinf(LnPrior(daGuess, **kwargs)):
                daStateVector.append(daGuess)
                break

    if iSize > 1:
        return daStateVector
    else:
        return daStateVector[0]
# end function

def LnPrior(daStateVector, **kwargs):
    """
    log prior
    """

    dLnPrior = 0

    # Get the current vector
    dMass,dFSat,dTSat,dAge,dBeta,dFRock,dWaterMass,dEnvMass,dEscCoeffH,dEscCoeffH2O,dPressXUV,dAlbedo,dPlanetMass,dOrbPer = daStateVector

    # Uniform priors need no scaling, but must be in the limits
    if (dMass < dStarMassMin) or (dMass > dStarMassMax):
        return -np.inf
    if (dFSat < dFSatMin) or (dFSat > dFSatMax):
        return -np.inf
    if (dTSat < dTSatMin) or (dTSat > dTSatMax):
        return -np.inf
    if (dAge < dAgeMin) or (dAge > dAgeMax):
        return -np.inf
    if (dBeta < dBetaMin) or (dBeta > dBetaMax):
        return -np.inf
    if (dFRock < dFRockMin) or (dFRock > dFRockMax):
        return -np.inf
    if (dWaterMass < dWaterMassMin) or (dWaterMass > dWaterMassMax):
        return -np.inf
    if (dEnvMass < dEnvMassMin) or (dEnvMass > dEnvMassMax):
        return -np.inf
    if (dEscCoeffH < dEscCoeffHMin) or (dEscCoeffH > dEscCoeffHMax):
        return -np.inf
    if (dEscCoeffH2O < dEscCoeffH2OMin) or (dEscCoeffH2O > dEscCoeffH2OMax):
        return -np.inf
    if (dPressXUV < dPressXUVMin) or (dPressXUV > dPressXUVMax):
        return -np.inf
    if (dAlbedo < dAlbedoMin) or (dAlbedo > dAlbedoMax):
        return -np.inf
    if (dPlanetMass < dPlanetMassMin) or (dPlanetMass > dPlanetMassMax):
        return -np.inf
    if (dOrbPer < dOrbPerMin) or (dOrbPer > dOrbPerMax):
        return -np.inf

    # Scale the normally distributed priors
    dLnPrior += norm.logpdf(dAge, dAgeMean, dAgeSig)
    dLnPrior += norm.logpdf(dBeta, dBetaMean, dBetaSig)
    dLnPrior += norm.logpdf(dFSat, dFSatMean, dFSatSig)

    return dLnPrior
# end function

### Loglikelihood and MCMC functions ###

def LnLike(daStateVector, **kwargs):
    """
    loglikelihood function: runs VPLanet simulation
    """

    # Get the current state vector
    dMass,dFSat,dTSat,dAge,dBeta,dFRock,dWaterMass,dEnvMass,dEscCoeffH,dEscCoeffH2O,dPressXUV,dAlbedo,dPlanetMass,dOrbPer = daStateVector

    # Convert to VPLanet input
    dFSat = 10 ** dFSat # Unlog
    dStopTime = dAge*1.e9 # Convert from Gyr -> yr
    dOutputTime = dStopTime # Output only at the end of the simulation

    # Initialize the output parameter vector
    daParams = np.zeros(iNumOutputPrms)

    # Get the prior probability to ignore unphysical state vectors
    # Do this to prevent errors stemming from VPLanet not finishing
    dLnPrior = kwargs["LnPrior"](daStateVector, **kwargs)
    if np.isinf(dLnPrior):
        for iPrm in range(iNumOutputPrms):
            daParams[iPrm] = np.nan
        return -np.inf, daParams

    # Get strings containing VPLanet input files (they must be provided!)
    try:
        sStarFileIn = kwargs.get("STARIN")
        sPrimaryFileIn = kwargs.get("VPLIN")
        sPlanetFileIn = kwargs.get("PLANETIN")
    except KeyError as err:
        print("ERROR: Must supply VPLIN, STARIN, and EIN in LnLike.")
        raise

    # Get PATH
    try:
        PATH = kwargs.get("PATH")
    except KeyError as err:
        print("ERROR: Must supply PATH.")
        raise

    # Randomize file names to prevent overwrites
    sVPLName = 'vpl%012x' % random.randrange(16**12)
    sStarName = 'st%012x' % random.randrange(16**12)
    sPlanetName = 'e%012x' % random.randrange(16**12)
    sPrimaryFile = sVPLName + '.in'
    sStarFile = sStarName + '.in'
    sPlanetFile = sPlanetName + '.in'
    sLogFile = sVPLName + '.log'
    sStarFwFile = '%s.star.forward' % sVPLName
    sPlanetFwFile = '%s.e.forward' % sPlanetName

# mkdir output!

    # Populate the input files
    sStarFileIn = re.sub("%s(.*?)#" % "dMass", "%s %.6e #" % ("dMass", dMass), sStarFileIn)
    sStarFileIn = re.sub("%s(.*?)#" % "dSatXUVFrac", "%s %.6e #" % ("dSatXUVFrac", dFSat), sStarFileIn)
    sStarFileIn = re.sub("%s(.*?)#" % "dSatXUVTime", "%s %.6e #" % ("dSatXUVTime", -dTSat), sStarFileIn)
    sStarFileIn = re.sub("%s(.*?)#" % "dXUVBeta", "%s %.6e #" % ("dXUVBeta", -dBeta), sStarFileIn)

    sPlanetFileIn = re.sub("%s(.*?)#" % "dFracRock", "%s %.6e #" % ("dFracRock", dFRock), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dSurfaceWaterMass", "%s %.6e #" % ("dSurfaceWaterMass", dWaterMass), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dEnvelopeMass", "%s %.6e #" % ("dEnvelopeMass", -dEnvMass), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dAtmXAbsEffH2O", "%s %.6e #" % ("dAtmXAbsEffH2O", dEscCoeffH2O), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dAtmXAbsEffH ", "%s %.6e #" % ("dAtmXAbsEffH", dEscCoeffH), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dPresXUV", "%s %.6e #" % ("dPresXUV", dPressXUV), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dAlbedoGlobal", "%s %.6e #" % ("dAlbedoGlobal", dAlbedo), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dMass", "%s %.6e #" % ("dMass", dPlanetMass), sPlanetFileIn)
    sPlanetFileIn = re.sub("%s(.*?)#" % "dOrbPeriod", "%s %.6e #" % ("dOrbPeriod", dOrbPer), sPlanetFileIn)

    with open(os.path.join(PATH, "output", sStarFile), 'w') as f:
        print(sStarFileIn, file = f)

    with open(os.path.join(PATH, "output", sPlanetFile), 'w') as f:
        print(sPlanetFileIn, file = f)

    # Populate the primary input file
    # Populate list of planets
    saBodyFiles = str(sStarFile) + " " + str(sPlanetFile) +" #"

    saBodyFiles = saBodyFiles.strip()

    sPrimaryFileIn = re.sub('%s(.*?)#' % "dStopTime", '%s %.6e #' % ("dStopTime", dStopTime), sPrimaryFileIn)
    sPrimaryFileIn = re.sub('%s(.*?)#' % "dOutputTime", '%s %.6e #' % ("dOutputTime", dOutputTime), sPrimaryFileIn)
    sPrimaryFileIn = re.sub('sSystemName(.*?)#', 'sSystemName %s #' % sVPLName, sPrimaryFileIn)
    sPrimaryFileIn = re.sub('saBodyFiles(.*?)#', 'saBodyFiles %s #' % saBodyFiles, sPrimaryFileIn)
    with open(os.path.join(PATH, "output", sPrimaryFile), 'w') as f:
        print(sPrimaryFileIn, file = f)

    #exit()

    # Run VPLANET and get the output, then delete the output files
    subprocess.call(["vplanet", '-q', sPrimaryFile], cwd = os.path.join(PATH, "output"))
    output = vpl.GetOutput(os.path.join(PATH, "output"), logfile = sLogFile)

    try:
        os.remove(os.path.join(PATH, "output", sStarFile))
        os.remove(os.path.join(PATH, "output", sPlanetFile))
        os.remove(os.path.join(PATH, "output", sPrimaryFile))
        os.remove(os.path.join(PATH, "output", sStarFwFile))
        os.remove(os.path.join(PATH, "output", sLogFile))
    except FileNotFoundError:
        # Run failed!
        for iPrm in range(iNumOutputPrms):
            daParams[iPrm] = np.nan
        return -np.inf, daParams

    # Ensure we ran for as long as we set out to
    # XXX Why not == pr < dEpsilon?
    if not output.log.final.system.Age / YEARSEC >= dStopTime:
        for iPrm in range(iNumOutputPrms):
            daParams[iPrm] = np.nan
        return -np.inf, daParams

    # Get final values of observed parameters
    dLumTrial = float(output.log.final.star.Luminosity)
    dLumXUVTrial = float(output.log.final.star.LXUVStellar)
    dTEffTrial = float(output.log.final.star.Temperature)
    dMassETrial = float(output.log.final.e.Mass)
    dRadETrial = float(output.log.final.e.Radius)

    # Compute ratio of XUV to bolometric luminosity
    dLumXUVRatioTrial = dLumXUVTrial / dLumTrial

    # Get extra parameters
    dStarRad = float(output.log.final.star.Radius)
    dEqTemp = float(output.log.final.e.ThermTemp)
    dSurfPress = float(output.log.final.e.PresSurf)
    dWaterMass = float(output.log.final.e.SurfWaterMass)
    dEnvMass = float(output.log.final.e.EnvelopeMass)
    dOxygenPress = float(output.log.final.e.OxygenMass)
    dSemi = float(output.log.final.e.SemiMajorAxis)

    # Extract constraints
    # Must at least have luminosity, err for star
    dLum = kwargs.get("LUM")
    dLumSig = kwargs.get("LUMSIG")

    dLumXUVRatio = kwargs.get("LUMXUVRATIO")
    dLumXUVRatioSig = kwargs.get("LUMXUVRATIOSIG")

    dTEff = kwargs.get("TEFF")
    dTEffSig = kwargs.get("TEFFSIG")

    dMassE = kwargs.get("MASSE")
    dMassESig = kwargs.get("MASSESIG")

    dRadE = kwargs.get("RADE")
    dRadESig = kwargs.get("RADESIG")

    # Compute the likelihood using provided constraints, assuming we have
    # luminosity constraints for host star
    dLnLike = 0

    dLnLike += ((dLum - dLumTrial) / dLumSig) ** 2
    dLnLike += ((dLumXUVRatio - dLumXUVRatioTrial) / dLumXUVRatioSig) ** 2
    dLnLike += ((dTEff - dTEffTrial) / dTEffSig) ** 2
    dLnLike += ((dMassE - dMassETrial) / dMassESig) ** 2
    dLnLike += ((dRadE - dRadETrial) / dRadESig) ** 2

    dLnLike = -0.5 * dLnLike

    print("dLnLike: ",dLnLike)
    # Return likelihood and diognostic parameters

    daParams = np.array([dLumTrial,
                         dLumXUVTrial,
                         dTEffTrial,
                         dMassETrial,
                         dRadETrial,
                         dStarRad,
                         dEqTemp,
                         dSurfPress,
                         dWaterMass,
                         dEnvMass,
                         dOxygenPress,
                         dSemi
                         ])
    return dLnLike, daParams
#end function

######## Main code begins here ###########

# Observational constraints
iNumObs = 6

dLumMean = 0.000522               # Van Grootel et al. (2018) [Lsun]
dLumSig = 0.000019            # Van Grootel et al. (2018) [Lsun]

dTEffMean = 2511
dTEffSig = 37

dLRatioMean = 7.5e-4              # Wheatley et al. (2017)
dLRatioSig = 1.5e-4           # Wheatley et al. (2017)

dMassEMean = 0.772
dMassESig = 0.077

dRadEMean = 0.91
dRadESig = 0.027

# Model parameters with normally distributed priors
iNumModelPrms = 14
iNumOutputPrms = 12

dBetaMean = -1.18                 # Jackson et al. (2012)
dBetaSig = 0.31               # Jackson et al. (2012)

dAgeMean = 7.6                    # Burgasser et al. (2017) [Gyr]
dAgeSig = 2.2                 # Burgasser et al. (2017) [Gyr]

dFSatMean = -2.92                 # Wright et al. (2011) and Chadney et al. (2015)
dFSatSig = 0.26               # Wright et al. (2011) and Chadney et al. (2015)

dOrbPerMean = 6.099043
dOrbPerSig = 1.5e-5

# Define bounds of the sampled parameter space for the model parameters
dStarMassMin = 0.07
dStarMassMax = 0.11

dTSatMin = 0.1
dTSatMax = 12

dFSatMin = -5
dFSatMax = -1

dAgeMin = 0.1
dAgeMax = 12

dBetaMin = -2
dBetaMax = 0

dTSatMin = -5
dTSatMax = -1

dFRockMin = 0.05
dFRockMax = 0.9

# Earth masses
dWaterMassMin = 0
dWaterMassMax = 0.5

dEnvMassMin = 0.0
dEnvMassMax = 10

dEscCoeffHMin = 0.1
dEscCoeffHMax = 0.5

dEscCoeffH2OMin = 0.06
dEscCoeffH2OMax = 0.13

dPressXUVMin = 0.1
dPressXUVMax = 10

dAlbedoMin = 0.05
dAlbedoMax = 0.9

dPlanetMassMin = 0.5
dPlanetMassMax = 10

dOrbPerMin = 6.098
dOrbPerMax = 6.1

# Dictionary to hold all constraints
kwargs = {"PATH" : ".",                          # Path to all files
          "LnPrior" : LnPrior,          # Function for priors
          "PriorSample" : SampleStateVector,  # Function to sample priors
          "LUM" : dLumMean,                  # Best fit luminosity constraint
          "LUMSIG" : dLumSig,            # Luminosity uncertainty (Gaussian)
          "LUMXUVRATIO" : dLRatioMean,       # L_bol/L_XUV best fit
          "LUMXUVRATIOSIG" : dLRatioSig, # L_bol/L_XUV uncertainty (Gaussian)
          "TEFF" : dTEffMean,
          "TEFFSIG" : dTEffSig,
          "MASSE" : dMassEMean,
          "MASSESIG" : dMassESig,
          "RADE" : dRadEMean,
          "RADESIG" : dRadESig
          }
# Define approxposterior parameters
iTrainInit = 10                         # Initial size of training set
iNewPoints = 10                          # Number of new points to find each iteration
iMaxIter = 10                        # Maximum number of iterations
iSeed = 77                        # RNG seed
iGPRestarts = 10                 # Number of times to restart GP hyperparameter optimizations
iMinObjRestarts = 10             # Number of times to restart objective fn minimization
iGPOptInterval = 25                 # Optimize GP hyperparameters even this many iterations
dConvergenceLimit = 0.1
iNumConverged = 3

# Prior bounds
daBounds = (
          (dStarMassMin, dStarMassMax),
          (dFSatMin, dFSatMax),
          (dTSatMin, dTSatMax),
          (dAgeMin, dAgeMax),
          (dBetaMin, dBetaMax),
          (dFRockMin, dFRockMax),
          (dWaterMassMin, dWaterMassMax),
          (dEnvMassMin, dEnvMassMax),
          (dEscCoeffHMin, dEscCoeffHMax),
          (dEscCoeffH2OMin, dEscCoeffH2OMax),
          (dPressXUVMin, dPressXUVMax),
          (dAlbedoMin, dAlbedoMax),
          (dPlanetMassMin, dPlanetMassMax),
          (dOrbPerMin,dOrbPerMax)
          )

sAlgorithm = "BAPE"              # Kandasamy et al. (2015) formalism

# Set RNG seed
np.random.seed(iSeed)

# emcee.EnsembleSampler, emcee.EnsembleSampler.run_mcmc and GMM parameters
samplerKwargs = {"nwalkers" : 100}
mcmcKwargs = {"iterations" : int(1.0e4)}

# Loglikelihood function setup required to run VPLanet simulations
#kwargs = trappist1.kwargsTRAPPIST1
#PATH = os.path.dirname(os.path.abspath(__file__))
#kwargs["PATH"] = PATH

# Extract path
PATH = kwargs["PATH"]

print (PATH)
# Get the input files, save them as strings
with open(os.path.join(PATH, "star.in"), 'r') as f:
    sStarFile = f.read()
    kwargs["STARIN"] = sStarFile
with open(os.path.join(PATH, "e.in"), 'r') as f:
    sPlanetFile = f.read()
    kwargs["PLANETIN"] = sPlanetFile
with open(os.path.join(PATH, "vpl.in"), 'r') as f:
    sPrimaryFile = f.read()
    kwargs["VPLIN"] = sPrimaryFile
#with open(os.path.join(PATH, "e.in"), 'r') as f:
#    sEFile = f.read()
#    kwargs["EIN"] = sEFile

# Make sure directory output/ exists, if not make it!
if not os.path.exists("output"):
    subprocess.call(["mkdir","output"])

# Generate initial training set using latin hypercube sampling over parameter bounds
# Evaluate forward model log likelihood + lnprior for each theta
if not os.path.exists("apRunAPFModelCache.npz"):
    y = np.zeros(iTrainInit)
    theta = np.zeros((iTrainInit,iNumModelPrms))
    for iTrial in range(iTrainInit):
        print("Training simulation: %d" % iTrial)
        theta[iTrial,:] = SampleStateVector()
        y[iTrial] = LnLike(theta[iTrial], **kwargs)[0] + LnPrior(theta[iTrial], **kwargs)
    #np.savez("apRunAPFModelCache.npz", theta=theta, y=y)

else:
    print("Loading in cached simulations...")
    sims = np.load("apRunAPFModelCache.npz")
    theta = sims["theta"]
    y = sims["y"]

### Initialize GP ###

# Use ExpSquared kernel, the approxposterior default option
gp = gpUtils.defaultGP(theta, y, order=None, white_noise=-12)

# Initialize approxposterior
ap = approx.ApproxPosterior(theta=theta,
                            y=y,
                            gp=gp,
                            lnprior=LnPrior,
                            lnlike=LnLike,
                            priorSample=SampleStateVector,
                            bounds=daBounds,
                            algorithm=sAlgorithm)

# Run!
ap.run(m=iNewPoints, nmax=iMaxIter, estBurnin=True, mcmcKwargs=mcmcKwargs,
       thinChains=True,samplerKwargs=samplerKwargs, verbose=True,
       nGPRestarts=iGPRestarts,nMinObjRestarts=iMinObjRestarts,
       optGPEveryN=iGPOptInterval,seed=iSeed, cache=True, convergenceCheck=True,
       eps=dConvergenceLimit,kmax=iNumConverged,**kwargs)
# Done!


text:
#!/usr/bin/env python

# Import modules
import numpy as np
import sklearn
from sklearn.preprocessing import LabelEncoder
import pickle
from sensor_stick.srv import GetNormals
from sensor_stick.features import compute_color_histograms
from sensor_stick.features import compute_normal_histograms
from visualization_msgs.msg import Marker
from sensor_stick.marker_tools import *
from sensor_stick.msg import DetectedObjectsArray
from sensor_stick.msg import DetectedObject
from sensor_stick.pcl_helper import *

import rospy
import tf
from geometry_msgs.msg import Pose
from std_msgs.msg import Float64
from std_msgs.msg import Int32
from std_msgs.msg import String
from pr2_robot.srv import *
from rospy_message_converter import message_converter
import yaml

# Helper function to get surface normals
def get_normals(cloud):
    get_normals_prox = rospy.ServiceProxy('/feature_extractor/get_normals', GetNormals)
    return get_normals_prox(cloud).cluster

# Helper function to create a yaml friendly dictionary from ROS messages
def make_yaml_dict(test_scene_num, arm_name, object_name, pick_pose, place_pose):
    yaml_dict = {}
    yaml_dict["test_scene_num"] = test_scene_num.data
    yaml_dict["arm_name"]  = arm_name.data
    yaml_dict["object_name"] = object_name.data
    yaml_dict["pick_pose"] = message_converter.convert_ros_message_to_dictionary(pick_pose)
    yaml_dict["place_pose"] = message_converter.convert_ros_message_to_dictionary(place_pose)
    return yaml_dict

# Helper function to output to yaml file
def send_to_yaml(yaml_filename, dict_list):
    data_dict = {"object_list": dict_list}
    with open(yaml_filename, 'w') as outfile:
        yaml.dump(data_dict, outfile, default_flow_style=False)

# Callback function for your Point Cloud Subscriber
def pcl_callback(pcl_msg):

    cloud = ros_to_pcl(pcl_msg)

    vox = cloud.make_voxel_grid_filter()

    # Choose a voxel (also known as leaf) size
    LEAF_SIZE = 0.003

    # Set the voxel (or leaf) size
    vox.set_leaf_size(LEAF_SIZE, LEAF_SIZE, LEAF_SIZE)

    # Call the filter function to obtain the resultant downsampled point cloud
    cloud_filtered = vox.filter()

### Statistical Outlier Removal Filter
    # Much like the previous filters, we start by creating a filter object:
    outlier_filter = cloud_filtered.make_statistical_outlier_filter()

    # Set the number of neighboring points to analyze for any given point
    outlier_filter.set_mean_k(50)

    # Set threshold scale factor
    x = 1.0

    # Any point with a mean distance larger than global (mean distance+x*std_dev) will be considered outlier
    outlier_filter.set_std_dev_mul_thresh(x)
    cloud_filtered = outlier_filter.filter()

### PassThrough Filter
    # Create a PassThrough filter object first across the Z axis
    passThroughZ = cloud_filtered.make_passthrough_filter()
    filter_axis = 'z'
    passThroughZ.set_filter_field_name(filter_axis)
    axis_min = 0.6
    axis_max = 1.1
    passThroughZ.set_filter_limits(axis_min, axis_max)
    cloud_filtered = passThroughZ.filter()

## Now, Create a PassThrough filter object across the Y axis
    passThroughY = cloud_filtered.make_passthrough_filter()
    filter_axis = 'y'
    passThroughY.set_filter_field_name(filter_axis)
    axis_min = -0.5
    axis_max = 0.5
    passThroughY.set_filter_limits(axis_min, axis_max)
    cloud_filtered = passThroughY.filter()

### RANSAC Plane Segmentation
    # Create the segmentation object
    seg = cloud_filtered.make_segmenter()
    seg.set_model_type(pcl.SACMODEL_PLANE)
    seg.set_method_type(pcl.SAC_RANSAC)

    max_distance = 0.01
    seg.set_distance_threshold(max_distance)
    inliers, coefficients = seg.segment()

    ### Extract inliers and outliers
    pcl_table_cloud = cloud_filtered.extract(inliers, negative=False)
    pcl_object_cloud = cloud_filtered.extract(inliers, negative=True)

### Euclidean Clustering
    # Go from XYZRGB to RGB since to build the k-d tree we only needs spatial data
    white_cloud = XYZRGB_to_XYZ(pcl_object_cloud)
    # Apply function to convert XYZRGB to XYZ
    tree = white_cloud.make_kdtree()

    ec = white_cloud.make_EuclideanClusterExtraction()
    ec.set_ClusterTolerance(0.005)
    ec.set_MinClusterSize(150)
    ec.set_MaxClusterSize(50000)
    ec.set_SearchMethod(tree)
    cluster_indices = ec.Extract()

    #Create Cluster-Mask Point Cloud to visualize each cluster separately
    cluster_color = get_color_list(len(cluster_indices))
    color_cluster_point_list = []

    for j, indices in enumerate(cluster_indices):
        for i, indice in enumerate(indices):
             color_cluster_point_list.append([white_cloud[indice][0],
                                              white_cloud[indice][1],
                                              white_cloud[indice][2],
                                     rgb_to_float(cluster_color[j])])

    #Create new cloud containing all clusters, each with unique color
    cluster_cloud = pcl.PointCloud_PointXYZRGB()
    cluster_cloud.from_list(color_cluster_point_list)

    #Convert PCL data to ROS messages
    ros_cluster_cloud = pcl_to_ros(cluster_cloud)
    ros_object_cloud = pcl_to_ros(pcl_object_cloud)
    ros_table_cloud = pcl_to_ros(pcl_table_cloud)

    #Publish ROS messages
    pcl_objects_pub.publish(ros_object_cloud)
    pcl_table_pub.publish(ros_table_cloud)
    pcl_cluster_pub.publish(ros_cluster_cloud)

    det_obj_label = []
    det_obj = []
        # Grab the points for the cluster
    for index, pts_list in enumerate(cluster_indices):
	pcl_cluster = pcl_object_cloud.extract(pts_list)
	ros_cluster = pcl_to_ros(pcl_cluster)

        # Compute associated feature vector
	color_hist = compute_color_histograms(ros_cluster, using_hsv=True)

	norms = get_normals(ros_cluster)
	norm_hist = compute_normal_histograms(norms)

        # Make prediction
	features = np.concatenate((color_hist, norm_hist))
	predict = clf.predict(scaler.transform(features.reshape(1,-1)))
	label = encoder.inverse_transform(predict)[0]
	det_obj_label.append(label)

        # Publish a label into RViz
	label_pos = list(white_cloud[pts_list[0]])
	label_pos[2] += .2
	object_markers_pub.publish(make_label(label,label_pos,index))

	do = DetectedObject()
	do.label = label
	do.cloud = ros_cluster
	det_obj.append(do)

    rospy.loginfo('Detected {} objects: {}'.format(len(det_obj_label), det_obj_label))

    # Publish the list of detected objects
    detected_objects_pub.publish(det_obj)

    # Suggested location for where to invoke your pr2_mover() function within pcl_callback()
    true_object_list = rospy.get_param('/object_list')
    dropbox_param = rospy.get_param('/dropbox')
    true_list = []

    for i in range(len(true_object_list)):
	true_list.append(true_object_list[i]['name'])

    print "\n"
    print "Need to find: "
    print true_list
    print "\n"
    trueset = set(true_list)
    detset = set(det_obj_label)

    # Could add some logic to determine whether or not your object detections are robust
    # before calling pr2_mover()
    if detset == trueset:
	try:
            pr2_mover(det_obj)
	except rospy.ROSInterruptException:
            pass
    else:
	rospy.loginfo("Wrong object(s) detected")

# function to load parameters and request PickPlace service
def pr2_mover(object_list):

    #Initialize variables
    test_scene_num = Int32()
    object_name = String()
    pick_pose = Pose()
    place_pose = Pose()
    arm_name = String()
    yaml_dict_list = []
    test_scene_num.data = 2

    dropbox_pose = []
    dropbox_name = []
    dropbox_group = []

    true_object_list = rospy.get_param('/object_list')
    dropbox_param = rospy.get_param('/dropbox')

    for i in range(0, len(dropbox_param)):
	dropbox_pose.append(dropbox_param[i]['position'])
	dropbox_group.append(dropbox_param[i]['group'])

    centroids = [] # to be list of tuples (x, y, z)
    labels = []

    for i in range(0, len(true_object_list)):
        object_name.data = true_object_list[i]['name' ]
        object_group     = true_object_list[i]['group']
	for j in object_list:
	    labels.append(j.label)
	    point_arr = ros_to_pcl(j.cloud).to_array()
	    avg = np.mean(point_arr, axis=0)[:3]
	    centroids.append(avg)

	if object_group == 'red':
	    arm_name.data = 'left'
	    place_pose.position.x = dropbox_pose[0][0] - random.randint(.05,.1)
	    place_pose.position.y = dropbox_pose[0][1]
	    place_pose.position.z = dropbox_pose[0][2]
	else:
	    arm_name.data = 'right'
	    place_pose.position.x = dropbox_pose[1][0] - random.randint(.05,.1)
	    place_pose.position.y = dropbox_pose[1][1]
	    place_pose.position.z = dropbox_pose[1][2]
        try:
            index = labels.index(object_name.data)
            pick_pose.position.x = np.asscalar(centroids[index][0])
            pick_pose.position.y = np.asscalar(centroids[index][1])
            pick_pose.position.z = np.asscalar(centroids[index][2])
        except ValueError:
            continue

        # Create a list of dictionaries for later output to yaml format
        yaml_dict = make_yaml_dict(test_scene_num, arm_name, object_name, pick_pose, place_pose)
        yaml_dict_list.append(yaml_dict)
        rospy.wait_for_service('pick_place_routine')

        try:
            pick_place_routine = rospy.ServiceProxy('pick_place_routine', PickPlace)

            #Insert your message variables to be sent as a service requestc
            resp = pick_place_routine(test_scene_num, object_name, arm_name, pick_pose, place_pose)

            print ("Response: ",resp.success)

        except rospy.ServiceException, e:
            print "Service call failed: %s"%e

    #Output your request parameters into output yaml file
    yaml_filename = 'output_'+str(test_scene_num.data)+'.yaml'
    send_to_yaml(yaml_filename, yaml_dict_list)

    return

if __name__ == '__main__':

    # ROS node initialization
    rospy.init_node('clustering', anonymous=True)

    # Create Subscribers
    pcl_sub = rospy.Subscriber("/pr2/world/points", pc2.PointCloud2, pcl_callback,queue_size=1)

    # Create Publishers
    pcl_objects_pub = rospy.Publisher("/pcl_objects", PointCloud2,queue_size=1)
    pcl_table_pub = rospy.Publisher("/pcl_table", PointCloud2, queue_size=1)
    pcl_cluster_pub = rospy.Publisher("/pcl_cluster", PointCloud2, queue_size=1)
    object_markers_pub = rospy.Publisher("/object_markers", Marker, queue_size=1)
    detected_objects_pub = rospy.Publisher("/detected_objects", DetectedObjectsArray, queue_size=1)
    pr2_base_mover_pub = rospy.Publisher("/pr2/world_joint_controller/command", Float64, queue_size=10)

    #Load Model From disk
    model = pickle.load(open('model.sav','rb'))
    clf = model['classifier']
    encoder = LabelEncoder()
    encoder.classes_ = model['classes']
    scaler = model['scaler']

    # Initialize color_list
    get_color_list.color_list = []

    # Spin while node is not shutdown
    while not rospy.is_shutdown():
	rospy.spin()