text:
import gym
import gym_gazebo
from qlearn import QLearn
import time
import numpy
import pandas
from itertools import product
from functools import reduce
import tensorflow as tf
from tensorboard_logging import Logger

# Inspired by Basic Reinforcement Learning Tutorial 4: Q-learning in OpenAI gym
#  https://github.com/vmayoral/basic_reinforcement_learning/blob/master/tutorial4/q-learning-gym-1.py

def build_state(features):
    return int("".join(map(lambda feature: str(int(feature)), features)))

def to_bin(value, bins):
    return numpy.digitize(x=[value], bins=bins)[0]

def step(action, state):
    """
    Implementation of "step" which uses the given action
    as an offset of the existing state. This function overloads
    the environment step method and should be used when actions-as-offsets
    are required for this particular environment.
    """
    # if action == (0.0, 0.0, 0.0):
        # print("\t\trobot decided to stay still!")
    offset_action = [a + s for a,s in zip(list(action), state)]
    # print("step: action: ", action)
    # print("step: state: ", state)
    # print("step: offset_action: ", offset_action)
    observation, reward, done, info = env.step(offset_action)
    return observation, reward, done, info

# Create the environment
env = gym.make('GazeboModularScara3DOF-v2')

# print(env.action_space) # Box(3,)
# print(env.observation_space) # Box(9,)
# print(env.observation_space.shape) # (9,)
# print(env.observation_space.shape[0]) # 9
# print(env.action_space) # 9

goal_average_steps = 2
max_number_of_steps = 20
last_time_steps = numpy.ndarray(0)
n_bins = 10
epsilon_decay_rate = 0.99 ########
save_model_with_prefix = 'qlearn_test' ######
it = 1 ######

number_of_features = env.observation_space.shape[0] # 9 however we're only using 3 of them (joint angles)
last_time_steps = numpy.ndarray(0)

# Number of states is huge so in order to simplify the situation
# typically, we discretize the space to: n_bins ** number_of_features
joint1_bins = pandas.cut([-numpy.pi/2, numpy.pi/2], bins=n_bins, retbins=True)[1][1:-1]
joint2_bins = pandas.cut([-numpy.pi/2, numpy.pi/2], bins=n_bins, retbins=True)[1][1:-1]
joint3_bins = pandas.cut([-numpy.pi/2, numpy.pi/2], bins=n_bins, retbins=True)[1][1:-1]
action_bins = pandas.cut([-numpy.pi/2, numpy.pi/2], bins=n_bins, retbins=True)[1][1:-1]
# print("joint1_bins: ", joint1_bins)

difference_bins = abs(joint1_bins[0] - joint1_bins[1])
# print("difference_bins: ", difference_bins)

# Generate actions
############
# test actions
# actions = [(0.0, 0.0, 0.0), (numpy.pi/2, numpy.pi/2, numpy.pi/2), (0 , 0, numpy.pi/2)]
# actions = [(0.0, 0.0, 0.0), (numpy.pi/2, numpy.pi/2, numpy.pi/2), (0, 0, numpy.pi/5)]

# generate actions so that it can re #####ach any "measurable" point in the space from anywhere
# with 10 bins and combinations of 3 elements: 729 possible actions.
# actions = []
# for x in product(action_bins, repeat=3):
#     actions.append(x)

# the following actions represent an offset over the existing state,
#  modifications over the traditional "step" method of the environment
#  are implemented in the "step" function
actions = [(difference_bins, 0.0, 0.0), (-difference_bins, 0.0, 0.0),
            (0.0, difference_bins, 0.0), (0.0, -difference_bins, 0.0),
            (0.0, 0.0, difference_bins), (0.0, 0.0, -difference_bins),
            (0.0, 0.0, 0.0)]
############

# The Q-learn algorithm
#qlearn = QLearn(actions=actions,
#    alpha=0.2, gamma=0.90, epsilon=0.5, epsilon_decay_rate=0.99)

qlearn = QLearn(actions=actions,
    alpha=0.5, gamma=0.30, epsilon=0.8, epsilon_decay_rate=0.999)

# logs_path = "/tmp/modular_scara_3dof_qlearn/"
logger = Logger("/tmp/rosrl/GazeboModularScara3DOF-v2/q-learn/")

with tf.Session() as sess:
    # op to write logs to Tensorboard
    # summary_writer = tf.summary.FileWriter(logs_path, graph=tf.get_default_graph())

    for i_episode in range(1000): # episodes
        print("I_EPISODE", i_episode)#####
        observation = env.reset()

        joint1_position, joint2_position, joint3_position  = observation[:3]
        state = build_state([to_bin(joint1_position, joint1_bins),
                         to_bin(joint2_position, joint2_bins),
                         to_bin(joint3_position, joint3_bins)])
        episode_reward = 0

        for t in range(max_number_of_steps):
            env.render()
            # print("join1_bins", joint1_bins)
            # print("Number of steps", t)
            # print("\t\tq: ",qlearn.q)
            # Pick an action based on the current f
            # print("state: ",state)
            # Pick an action based on the current state
            action = qlearn.chooseAction(state)
            # print("state:", state)##
            # print("action: ",action)
            # print("state observation: ",observation[:3])###

            # Execute the action and get feedback
            # observation, reward, done, info = env.step(action) # environment step method
            observation, reward, done, info = step(action, observation[:3])   # local step method, applies the action as an offset
                                                                                # to the state
            # print("reward: ",reward)
            episode_reward += reward
            # print("next state observation: ",observation[:3])###

            # Digitize the observation to get a state
            joint1_position, joint2_position, joint3_position  = observation[:3]
            nextState = build_state([to_bin(joint1_position, joint1_bins),
                            to_bin(joint2_position, joint2_bins),
                            to_bin(joint3_position, joint3_bins)])

            # print("nextState", nextState)
            # summary = tf.Summary(value=[tf.Summary.Value(tag="episode_reward",simple_value=episode_reward)])
            # summary_writer.add_summary(summary, i_episode * max_number_of_steps + t)

            if done:
                last_time_steps = numpy.append(last_time_steps, [int(t + 1)])
                break
            else:
                # Q-learn stuff
                #qlearn.learn(state, action, reward, nextState)
                qlearn.learn(state, action, reward, nextState, save_model_with_prefix, it)
                state = nextState

                it += 1 #####

        logger.log_scalar("episode_reward", episode_reward, i_episode)
        print("episode reward: ",episode_reward)

# l = last_time_steps.tolist()
# l.sort()
# print("Overall score: {:0.2f}".format(last_time_steps.mean()))
# print("Best 100 score: {:0.2f}".format(reduce(lambda x, y: x + y, l[-100:]) / len(l[-100:])))


text:
import numpy as np
import os
from utils.data_funcs import compute_ensemble_accuracy, get_conf_matrix, boolean_string, compute_q_statistics
from utils.data_funcs import compute_f1
from argparse import ArgumentParser

parser = ArgumentParser()
parser.add_argument("--directory", type=str, default="cifar10_cnn_part_nary_ecoc_4", help="")
parser.add_argument("--save_plot", type=boolean_string, default=False, help="")
config = parser.parse_args()

if "_full_" in config.directory:
    folder = os.path.join("assets", "results", "full_share", config.directory)
elif "_part_" in config.directory:
    folder = os.path.join("assets", "results", "part_share", config.directory)
else:
    folder = os.path.join("assets", "results", "private", config.directory)

label_folder = os.path.join("assets", "results", "a_golden_labels")

if "cifar10_" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "cifar10_golden_labels.npz"))["embeddings"]
elif "cifar100_" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "cifar100_golden_labels.npz"))["embeddings"]
elif "flower102_" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "flower_golden_labels.npz"))["embeddings"]
elif "mnist" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "mnist_golden_labels.npz"))["embeddings"]
elif "stsa" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "stsa_golden_labels.npz"))["embeddings"]
elif "trec" in config.directory:
    golden_labels = np.load(os.path.join(label_folder, "trec_golden_labels.npz"))["embeddings"]
else:
    raise ValueError("Unknown task!!! Only support [cifar10 | cifar100 | flower102 | mnist | stsa | trec]")

pred_labels = np.load(folder + "/pred_labels.npz")["embeddings"]
num_classifier = pred_labels.shape[1]

if num_classifier == 1:
    golden_labels = np.asarray(golden_labels)
    golden_labels = np.reshape(golden_labels, newshape=golden_labels.shape[0])
    pred_labels = np.reshape(pred_labels, newshape=pred_labels.shape[0])
    accuracy = np.mean(np.equal(pred_labels, golden_labels).astype(np.float32))
    print("Accuracy: {}\t{:4.2f}%".format(accuracy, accuracy * 100))
    if config.plot:
        get_conf_matrix(pred_labels, None, golden_labels, filename="assets/{}_heatmap".format(config.directory))
else:
    nary_ecoc = np.load(folder + "/nary_ecoc.npz")["embeddings"]
    nl = list(range(5, num_classifier + 1, 5))
    accuracies = list()
    for n in nl:
        accuracy = compute_ensemble_accuracy(pred_labels[:, 0:n], nary_ecoc[:, 0:n], golden_labels)
        accuracies.append(float("{:4.2f}".format(accuracy * 100)))
        print("Accuracy: n={}\t{}\t{:4.2f}%".format(n, accuracy, accuracy * 100))
    print(",".join([str(x) for x in nl]))
    print(",".join([str(x) for x in accuracies]))
    q_val = compute_q_statistics(pred_labels, nary_ecoc, golden_labels)
    print("Average Q-statistics: {}".format(q_val))
    micro_f1, macro_f1, weighted_f1 = compute_f1(pred_labels, nary_ecoc, golden_labels)
    print("Micro-F1 Score: {}".format(micro_f1))
    print("Macro-F1 Score: {}".format(macro_f1))
    print("Weighted Macro-F1 Score: {}".format(weighted_f1))
    get_conf_matrix(pred_labels, nary_ecoc, golden_labels, filename="assets/{}_heatmap".format(config.directory),
                    save_file=config.save_plot)

ext: py

max_stars_repo_name: IsaacChanghau/DeepNaryECOC

max_issues_repo_name: IsaacChanghau/DeepNaryECOC

max_forks_repo_name: IsaacChanghau/DeepNaryECOC

reason: import numpy


text:
import numpy as np
from scipy.stats import norm
from al_utilities import safe_div

class AlAgent:
    """
    This class definition specifies the properties of the object that implements the reduced Bayesian model.
    The model infers the mean of the outcome generating distribution according to changepoint probability and
    relative uncertainty.
    """

    def __init__(self, agent_vars):
        # This function creates an agent object of class Agent based on the agent initialization input

        # Set variable task properties based on input
        self.h = agent_vars.h
        self.u = agent_vars.u
        self.sigma = agent_vars.sigma
        self.sigma_t_sq = agent_vars.sigma_0
        self.tau_t = agent_vars.tau_0
        self.omega_t = agent_vars.omega_0
        self.mu_t = agent_vars.mu_0
        self.q = agent_vars.q
        self.s = agent_vars.s
        self.sigma_H = agent_vars.sigma_H
        self.C = agent_vars.C

        # Initialize variables
        self.a_t = np.nan
        self.alpha_t = np.nan
        self.tot_var = np.nan

    def learn(self, delta_t, b_t, v_t, mu_H, high_val):
        """ This function implements the inference of the reduced Bayesian model

        :param delta_t: Current prediction error
        :param b_t: Last prediction of participant
        :param v_t: Helicopter visibility
        :param mu_H: True helicopter location
        :param high_val: High-value index
        """

        # Update variance of predictive distribution
        self.tot_var = self.sigma ** 2 + self.sigma_t_sq

        # Compute changepoint probability
        # -------------------------------

        # Likelihood of prediction error given that no changepoint occurred: (U(delta_t;[0,300]))^s * h
        term_1 = ((1/300) ** self.s) * self.h

        # Likelihood of prediction error given that changepoint occurred: (N(delta_t;0,sigma^2_t+\sigma^2))^s * (1-h)
        term_2 = (norm.pdf(delta_t, 0, np.sqrt(self.tot_var)) ** self.s) * (1 - self.h)

        # Compute changepoint probability
        self.omega_t = safe_div(term_1, (term_2 + term_1))

        # Compute learning rate and update belief
        # ---------------------------------------
        self.alpha_t = self.omega_t + self.tau_t - self.tau_t * self.omega_t

        # Add reward bias to learning rate and correct for learning rates > 1 and < 0
        self.alpha_t = self.alpha_t + self.q * high_val
        if self.alpha_t > 1.0:
            self.alpha_t = 1.0
        elif self.alpha_t < 0.0:
            self.alpha_t = 0.0

        # Set model belief equal to last prediction of participant to estimate model using subjective prediction errors
        self.mu_t = b_t

        # Compute belief update
        # ---------------------

        # hat{a_t} := alpha_t * delta_t
        self.a_t = self.alpha_t * delta_t

        # mu_{t+1} := mu_t + a_t
        self.mu_t = self.mu_t + self.a_t

        # During catch trials, take helicopter location into consideration
        # ----------------------------------------------------------------
        if v_t:

            # Compute  helicopter weight
            # w_t := sigma_t^2 / (sigma_t^2 + sigma_H^2)
            w_t = self.sigma_t_sq / (self.sigma_t_sq + self.sigma_H ** 2)

            # Compute mean of inferred distribution with additional helicopter information
            # mu_t = (1 - w_t ) * mu_{t+1} + w_t * mu_H
            self.mu_t = (1 - w_t) * self.mu_t + w_t * mu_H

            # Recompute the model's update under consideration of the catch-trial information
            self.a_t = self.mu_t - b_t

            # Compute mixture variance of the two distributions...
            # 1 ((1/sigma_t^2) + (1/sigma_H^2))
            term_1 = safe_div(1, self.sigma_t_sq)
            term_2 = safe_div(1, self.sigma_H ** 2)
            self.C = safe_div(1, term_1 + term_2)

            # ...and update relative uncertainty accordingly
            # tau_t = C / (C + sigma^2)
            self.tau_t = safe_div(self.C, self.C + self.sigma ** 2)

        # Update relative uncertainty for the next trial
        # ----------------------------------------------

        # Update estimation uncertainty:
        # sigma_{t+1}^2 := (omega_t * sigma^2
        #                   + (1-omega_t) * tau_t *  sigma^2
        #                   + omega_t * (1 - omega_t) * (delta_t * (1 - tau_t))^2) / exp(u)

        # Note that u is already in exponential form
        term_1 = self.omega_t * (self.sigma ** 2)
        term_2 = (1 - self.omega_t) * self.tau_t * (self.sigma ** 2)
        term_3 = self.omega_t * (1 - self.omega_t) * ((delta_t * (1 - self.tau_t)) ** 2)
        self.sigma_t_sq = safe_div((term_1 + term_2 + term_3), self.u)

        # Update relative uncertainty:
        # tau_{t+1} := sigma_{t+1}^2 / (sigma_{t+1}^2 + sigma^2)
        self.tau_t = safe_div(self.sigma_t_sq, (self.sigma_t_sq + self.sigma ** 2))

reason: import numpy,from scipy


text:
from scipy import spatial
import numpy as np

def gd(reference: list, computed: list, norm='euclidean'):
    """
    :param reference: list of the reference points [(coord_1, coord_2, coord_3,..., coord_n)]
    :return: the sum of the euclidean disctance from the given reference points.

    Possible norms:  ‘braycurtis’, ‘canberra’, ‘chebyshev’, ‘cityblock’, ‘correlation’, ‘cosine’, ‘dice’,
                     ‘euclidean’, ‘hamming’, ‘jaccard’, ‘jensenshannon’, ‘kulsinski’, ‘mahalanobis’, ‘matching’,
                     ‘minkowski’, ‘rogerstanimoto’, ‘russellrao’, ‘seuclidean’, ‘sokalmichener’, ‘sokalsneath’,
                     ‘sqeuclidean’, ‘wminkowski’, ‘yule’.

    Type: Convergence Indicator

    The generational distance performence indicator [1,2,3] measueres the distancee from the
    given reference list. T

    .. math::

        GD(S,P) = \frac{1}{n}(\Sigma_{s \in S} min ||F(s) - F(r)||^2)^{\frac{1}{p}}

    where |S| is the number of the points in a Pareto-set approximation and P a discrete representation
    of the Pareto-front. Generally, p=2. In this case it is equivalent with M1 measure [4].
    When P = 1 it is equivalent with the Gamma-Metric.

    References
    ----------

    .: [1] David A. Van Veldhuizen and David A. Van Veldhuizen.
           Multi-objective evolutionary algorithms: classifications, analyses, and new innovations.
           Technical Report, Evolutionary Computation, 1999.
    .: [2] Audet, Charles, J. Bigeon, D. Cartier, Sébastien Le Digabel, and Ludovic Salomon.
           "Performance indicators in multiobjective optimization." Optimization Online (2018).
    .: [3] https://pymoo.org/misc/performance_indicator.html
    .: [4] E. Zitzler, K. Deb, L. Thiele, Comparison of multiobjective evolutionary algorithms: Empirical
           results, Evolutionary computation 8 (2) (2000) 173–195.

    """

    distances = spatial.distance.cdist(reference, computed, metric=norm)
    minimums = np.nanmin(distances, axis=0)

    return np.sum(minimums) / len(computed)

if __name__ == '__main__':
    # 2d test - euclidean
    ref = [(1., 1.), (2., 1. / 2.), (3., 1. / 3.), (4., 1. / 4.), (5., 1. / 5)]  ### 1/x
    calc = [(1., 1.01), (2., 0.51), (4., 0.26)]

    print(gd(ref, calc))  # result is 0.01

    # 3d test - chebyshev
    ref = [(1., 1., 1.), (2., 1. / 2., 1. / 2.), (3., 1. / 3., 1. / 3.), (4., 1. / 4., 1. / 4.),
           (5., 1. / 5, 1. / 5,)]  ### 1/x
    calc = [(1., 1.01, 1.01), (2., 0.51, 0.51), (4., 0.26, 0.26)]

    print(gd(ref, calc, norm='chebyshev'))  # 0.01

max_stars_repo_name: tamasorosz/pymoqi

max_issues_repo_name: tamasorosz/pymoqi

max_forks_repo_name: tamasorosz/pymoqi


text:
from itertools import product
from typing import overload, Optional, Union, List

from matplotlib.figure import Figure
from mpl_format.figures import FigureFormatter
from numpy.ma import arange
from pandas import DataFrame, Series

from probability.distributions.continuous.normal import Normal
from probability.distributions.mixins.conjugate import ConjugateMixin, \
    PredictiveMixin

from probability.utils import none_are_none, all_are_none, num_format, \
    one_is_none

class NormalNormalConjugate(
    ConjugateMixin,
    PredictiveMixin,
    object
):
    """
    Class for estimating the mean of a normal distribution with
    known variance σ² or known precision τ (where τ = 1 / σ²).

    Prior Hyper-parameters
    ----------------------
    * `μ₀` and one of `σ₀²` or `τ₀` are the hyper-parameters for the value of
    the normally-distributed mean of the distribution.
    * `μ₀` is the mean value of the mean.
    * `σ₀²` is the variance of the mean.
    * `τ₀` is the precision (inverse of the variance) of the mean.
    * Interpretation is that the mean was estimated from prior observations with
        - total precision (sum of all individual precisions) τ₀ = 1 / σ₀²
        - sample mean μ₀.

    Posterior Hyper-parameters
    --------------------------
    * `n` is the number of observations.
    * `x_sum` is the sum of the observations Σx.
    * `σ²` is the variance of the observations.
    * `τ` is the precision (inverse of the variance) of the observations

    Model parameters
    ----------------
    * `P(x)` is the probability that μ = x.
    * -∞ ≤ x ≤ ∞

    Links
    -----
    * https://en.wikipedia.org/wiki/Normal_distribution
    * https://en.wikipedia.org/wiki/Conjugate_prior
    """
    @overload
    def __init__(
            self, n: int,
            x_sum: float,
            sigma_sq: float,
            mu_0: float,
            sigma_0_sq: float
    ):
        """
        Create a new normal-normal conjugate distribution.

        :param n: Number of observations.
        :param x_sum: Sum of observations.
        :param sigma_sq: Variance of observations.
        :param mu_0: Value for the μ₀ (mean) hyper-parameter of the prior Normal
                     distribution describing the mean.
        :param sigma_0_sq: Value for the σ₀² (variance) hyper-parameter of the
                           prior Normal distribution describing the mean.
        """
        pass

    @overload
    def __init__(
            self, n: int,
            x_sum: float,
            tau: float,
            mu_0: float,
            tau_0: float,
    ):
        """
        Create a new normal-normal conjugate distribution.

        :param n: Number of observations.
        :param x_sum: Sum of observations.
        :param tau: Variance of observations.
        :param mu_0: Value for the μ₀ (mean) hyper-parameter of the prior Normal
                     distribution describing the mean.
        :param tau_0: Value for the τ₀ (precision) hyper-parameter of the
                      prior Normal distribution describing the mean.
        """
        pass

    def __init__(
            self, n: int,
            x_sum: float,
            mu_0: float,
            sigma_sq: Optional[float] = None,
            tau: Optional[float] = None,
            sigma_0_sq: Optional[float] = None,
            tau_0: Optional[float] = None
    ):
        """
        Create a new normal-normal conjugate distribution.

        :param n: Number of observations.
        :param x_sum: Sum of observations.
        :param sigma_sq: Variance of observations.
                         Distribution describing the mean.
        :param tau: Precision of observations.
        :param mu_0: Value for the μ₀ (mean) hyper-parameter of the prior Normal
        :param sigma_0_sq: Value for the σ₀² (variance) hyper-parameter of the
                           prior Normal distribution describing the mean.
        :param tau_0: Value for the τ₀ (precision) hyper-parameter of the
                      prior Normal distribution describing the mean.
        """
        if not (
            none_are_none(sigma_0_sq, sigma_sq) and
            all_are_none(tau_0, tau)
        ) and not (
            all_are_none(sigma_0_sq, sigma_sq) and
            none_are_none(tau_0, tau)
        ):
            raise ValueError('Give either σ² and σ₀² or τ and τ₀')

        self._n: int = n
        self._x_sum: float = x_sum
        self._mu_0: float = mu_0
        if sigma_0_sq is not None:
            self._sigma_sq: float = sigma_sq
            self._sigma_0_sq: float = sigma_0_sq
            self._parameterization: str = 'σ²'
        else:
            self._sigma_sq: float = 1 / tau
            self._sigma_0_sq: float = 1 / tau_0
            self._parameterization: str = 'τ'

    @property
    def mu_0(self) -> float:
        return self._mu_0

    @mu_0.setter
    def mu_0(self, value: float):
        self._mu_0 = value

    @property
    def sigma_sq(self) -> float:
        return self._sigma_sq

    @sigma_sq.setter
    def sigma_sq(self, value: float):
        self._sigma_sq = value
        self._parameterization = 'σ²'

    @property
    def sigma_0_sq(self) -> float:
        return self._sigma_0_sq

    @sigma_0_sq.setter
    def sigma_0_sq(self, value: float):
        self._sigma_0_sq = value
        self._parameterization = 'σ²'

    @property
    def tau(self) -> float:
        return 1 / self._sigma_sq

    @tau.setter
    def tau(self, value: float):
        self._sigma_sq = 1 / value
        self._parameterization = 'τ'

    @property
    def tau_0(self) -> float:
        return 1 / self._sigma_0_sq

    @tau_0.setter
    def tau_0(self, value: float):
        self._sigma_0_sq = 1 / value
        self._parameterization = 'τ'

    @property
    def n(self) -> int:
        return self._n

    # region posterior hyper-parameters

    @property
    def mu_0_prime(self) -> float:
        return (
            (self.tau_0 * self.mu_0 + self.tau * self._x_sum) /
            (self.tau_0 + self.n * self.tau)
        )

    @property
    def tau_0_prime(self) -> float:
        return self.tau_0 + self.n * self.tau

    @property
    def sigma_0_sq_prime(self) -> float:
        return 1 / self.tau_0_prime

    # endregion

    def prior(self) -> Normal:
        """
        Return a Normal distribution reflecting the prior belief about the
        distribution of the parameter μ, before seeing any data.
        """
        if self._parameterization == 'σ²':
            y_label = (
                '$P(μ_{Norm}=x|'
                'μ_{0, Norm},'
                'σ²_{0, Norm})$'
            )
        else:
            y_label = (
                '$P(μ_{Norm}=x|'
                'μ_{0, Norm},'
                'τ_{0, Norm})$'
            )
        return Normal(
            mu=self.mu_0,
            sigma_sq=self.sigma_0_sq
        ).with_x_label('μ').with_y_label(y_label).prepend_to_label('Prior: ')

    def likelihood(self, **kwargs) -> Normal:
        """
        Return a distribution reflecting the likelihood of observing
        the data, under a Normal model, independent of the prior belief
        about the distribution of parameter μ.
        """
        return Normal(
            mu=self._x_sum / self._n,
            sigma_sq=self._sigma_sq
        )

    def posterior(self) -> Normal:
        """
        Return a Normal distribution reflecting the posterior belief about the
        distribution of the parameter μ, after observing the data.
        """
        if self._parameterization == 'σ²':
            y_label = (
                '$P(μ_{Norm}=x|'
                'μ_{0, Norm}′,'
                'σ²_{0, Norm}′)$'
            )
        else:
            y_label = (
                '$P(μ_{Norm}=x|'
                'μ_{0, Norm}′,'
                'τ_{0, Norm}′)$'
            )

        return Normal(
            mu=self.mu_0_prime,
            sigma_sq=self.sigma_0_sq_prime
        ).with_x_label('μ').with_y_label(y_label).prepend_to_label(
            'Posterior: '
        )

    # region predictive

    def prior_predictive(self) -> Normal:
        """
        Return a Normal describing the expected distribution of future
        values based on the prior parameter estimates.
        """
        return Normal(
            mu=self.mu_0,
            sigma_sq=self.sigma_0_sq + self.sigma_sq
        ).with_y_label(
            r'$P(\tilde{X}=x|'
            r'μ_{Norm}=μ_{0, Norm},'
            r'σ²=σ²_{0, Norm}+σ²)$'
        )

    def posterior_predictive(self) -> Normal:
        """
        Return a Normal describing the expected distribution of future
        values based on the posterior parameter estimates.
        """
        return Normal(
            mu=self.mu_0_prime,
            sigma_sq=self.sigma_0_sq_prime + self.sigma_sq
        ).with_y_label(
            r'$P(\tilde{X}=x|'
            r'μ_{Norm}=μ_{0, Norm}′,'
            r'σ²=σ²_{0, Norm}′+σ²)$'
        )

    # endregion

    def plot(self, **kwargs) -> Figure:
        """
        Plot a grid of the different components of the Compound Distribution.

        :param kwargs: kwargs for plot methods
        """
        ppf_prior_01 = self.prior().ppf().at(0.01)
        ppf_prior_99 = self.prior().ppf().at(0.99)
        ppf_posterior_01 = self.posterior().ppf().at(0.99)
        ppf_posterior_99 = self.posterior().ppf().at(0.99)
        x_norm_min = int(min(ppf_prior_01, ppf_posterior_01)) - 1
        x_norm_max = int(max(ppf_prior_99, ppf_posterior_99)) + 1
        x_param = arange(x_norm_min, x_norm_max + 0.001, 0.001)
        ff = FigureFormatter(n_rows=2, n_cols=3)
        (
            ax_prior, ax_data, ax_posterior,
            ax_prior_predictive, ax_likelihood, ax_posterior_predictive
        ) = ff.axes.flat
        # plot prior and posterior parameters
        self.prior().plot(x=x_param, ax=ax_prior.axes, **kwargs)
        self.posterior().plot(x=x_param, ax=ax_posterior.axes, **kwargs)
        y_max_params = max(ax_prior.get_y_max(), ax_posterior.get_y_max())
        ax_prior.set_y_lim(0, y_max_params)
        ax_posterior.set_y_lim(0, y_max_params)
        ax_prior.set_title_text('prior').add_legend()
        ax_posterior.set_title_text('posterior').add_legend()
        # plot prior and posterior predictives
        ppf_prior_pred_01 = self.prior_predictive().ppf().at(0.01)
        ppf_prior_pred_99 = self.prior_predictive().ppf().at(0.99)
        ppf_posterior_pred_01 = self.posterior_predictive().ppf().at(0.01)
        ppf_posterior_pred_99 = self.posterior_predictive().ppf().at(0.99)
        x_pred_min = int(min(ppf_prior_pred_01, ppf_posterior_pred_01)) - 1
        x_pred_max = int(max(ppf_prior_pred_99, ppf_posterior_pred_99)) + 1
        x_pred = arange(x_pred_min, x_pred_max + 0.001, 0.001)
        self.prior_predictive().plot(
            x=x_pred, kind='line',
            ax=ax_prior_predictive.axes,
            **kwargs
        )
        self.posterior_predictive().plot(
            x=x_pred, kind='line',
            ax=ax_posterior_predictive.axes,
            **kwargs
        )
        y_max_pred = max(ax_prior_predictive.get_y_max(),
                         ax_posterior_predictive.get_y_max())
        ax_prior_predictive.set_y_lim(0, y_max_pred)
        ax_posterior_predictive.set_y_lim(0, y_max_pred)
        ax_prior_predictive.set_title_text('prior predictive').add_legend()
        ax_posterior_predictive.set_title_text(
            'posterior predictive'
        ).add_legend()
        # plot data
        observations = self.likelihood().rvs(num_samples=self._n)
        observations.plot.bar(ax=ax_data.axes, **kwargs)
        ax_data.set_text(title='data', x_label='i', y_label='$X_i$')
        y_abs_max = max([abs(y_lim) for y_lim in ax_data.get_y_lim()])
        ax_data.set_y_lim(-y_abs_max, y_abs_max)
        # plot likelihood
        x_like_min = int(self.likelihood().ppf().at(0.01)) - 1
        x_like_max = int(self.likelihood().ppf().at(0.99)) + 1
        x_exponential = arange(x_like_min, x_like_max + 0.001, 0.001)
        self.likelihood().plot(x=x_exponential, ax=ax_likelihood.axes)
        ax_likelihood.set_title_text('likelihood')
        ax_likelihood.add_legend()
        return ff.figure

    @staticmethod
    @overload
    def infer_posterior(
            data: Series,
            mu_0: float,
            sigma_0_sq: Optional[float] = None,
    ) -> Normal:
        pass

    @staticmethod
    @overload
    def infer_posterior(
            data: Series,
            mu_0: float,
            tau_0: Optional[float] = None
    ) -> Normal:
        pass

    @staticmethod
    def infer_posterior(
            data: Series,
            mu_0: float,
            sigma_0_sq: Optional[float] = None,
            tau_0: Optional[float] = None
    ) -> Normal:
        """
        Return a new Normal distribution of the posterior most likely to
        generate the given data.

        :param data: Series of float observations.
        :param mu_0: Value for the μ₀ (mean) hyper-parameter of the prior Normal
                     distribution describing the mean.
        :param sigma_0_sq: Value for the σ₀² (variance) hyper-parameter of the
                           prior Normal distribution describing the mean.
        :param tau_0: Value for the τ₀ (precision) hyper-parameter of the
                      prior Normal distribution describing the mean.
        """
        if not one_is_none(sigma_0_sq, tau_0):
            raise ValueError('Give either σ₀² or τ₀')

        n = len(data)
        x_sum = data.sum()
        if sigma_0_sq is None:
            tau = 1 / data.var()
            return NormalNormalConjugate(
                n=n, x_sum=x_sum, mu_0=mu_0,
                tau=tau, tau_0=tau_0
            ).posterior()
        else:
            sigma_sq = data.var()
            return NormalNormalConjugate(
                n=n, x_sum=x_sum, mu_0=mu_0,
                sigma_sq=sigma_sq, sigma_0_sq=sigma_0_sq
            ).posterior()

    @staticmethod
    @overload
    def infer_posteriors(
            data: DataFrame, prob_vars: Union[str, List[str]],
            cond_vars: Union[str, List[str]],
            mu_0: float,
            sigma_0_sq: Optional[float] = None,
            stats: Optional[Union[str, List[str]]] = None
    ) -> DataFrame:

        pass

    @staticmethod
    @overload
    def infer_posteriors(
            data: DataFrame, prob_vars: Union[str, List[str]],
            cond_vars: Union[str, List[str]],
            mu_0: float,
            tau_0: Optional[float] = None,
            stats: Optional[Union[str, List[str]]] = None
    ) -> DataFrame:

        pass

    @staticmethod
    def infer_posteriors(
            data: DataFrame, prob_vars: Union[str, List[str]],
            cond_vars: Union[str, List[str]],
            mu_0: float,
            sigma_0_sq: Optional[float] = None,
            tau_0: Optional[float] = None,
            stats: Optional[Union[str, List[str]]] = None,
    ) -> DataFrame:
        """
        Return a DataFrame mapping probability and conditional variables to Beta
        distributions of posteriors most likely to generate the given data.

        :param data: DataFrame containing discrete data.
        :param prob_vars: Name(s) of poisson variables whose posteriors to
                          find probability of.
        :param cond_vars: Names of discrete variables to condition on.
                          Calculations will be done for the cartesian product
                          of variable values
                          e.g if cA={1,2} and cB={3,4} then
                          cAB = {(1,3), (1, 4), (2, 3), (2, 4)}.
        :param mu_0: Value for the μ₀ (mean) hyper-parameter of the prior Normal
                     distribution describing the mean.
        :param sigma_0_sq: Value for the σ₀² (variance) hyper-parameter of the
                           prior Normal distribution describing the mean.
        :param tau_0: Value for the τ₀ (precision) hyper-parameter of the
                      prior Normal distribution describing the mean.
        :param stats: Optional stats to append to the output e.g. 'alpha',
                      'median'. To pass arguments use a dict mapping stat
                      name to iterable of args.
        :return: DataFrame with columns for each conditioning variable, a
                 'prob_var' column indicating the probability variable, a
                 `prob_val` column indicating the value of the probability
                 variable, and a `Beta` column containing the distribution.
        """
        if isinstance(prob_vars, str):
            prob_vars = [prob_vars]
        if isinstance(cond_vars, str):
            cond_vars = [cond_vars]
        cond_products = product(
            *[data[cond_var].unique() for cond_var in cond_vars]
        )
        if stats is not None:
            if isinstance(stats, str) or isinstance(stats, dict):
                stats = [stats]
        else:
            stats = []
        mus = []
        # iterate over conditions
        for cond_values in cond_products:
            cond_data = data
            cond_dict = {}
            for cond_var, cond_value in zip(cond_vars, cond_values):
                cond_data = cond_data.loc[cond_data[cond_var] == cond_value]
                cond_dict[cond_var] = cond_value
            for prob_var in prob_vars:
                prob_dict = cond_dict.copy()
                prob_dict['prob_var'] = prob_var
                if sigma_0_sq is not None:
                    posterior = NormalNormalConjugate.infer_posterior(
                        data=cond_data[prob_var],
                        mu_0=mu_0, sigma_0_sq=sigma_0_sq
                    )
                else:
                    posterior = NormalNormalConjugate.infer_posterior(
                        data=cond_data[prob_var],
                        mu_0=mu_0, tau_0=tau_0
                    )
                prob_dict['mu'] = posterior
                for stat in stats:
                    prob_dict = {**prob_dict,
                                 **posterior.stat(stat, True)}
                mus.append(prob_dict)

        mus_data = DataFrame(mus)

        return mus_data

    def __str__(self):

        if self._parameterization == 'σ²':
            return (
                f'NormalNormalConjugate('
                f'n={self._n}, '
                f'Σx={self._x_sum}, '
                f'σ²={num_format(self._sigma_sq, 3)})'
                f'μ₀={num_format(self._mu_0, 3)}, '
                f'σ₀²={num_format(self._sigma_0_sq, 3)}, '
            )
        elif self._parameterization == 'τ':
            return (
                f'NormalNormalConjugate('
                f'n={self._n}, '
                f'Σx={self._x_sum}, '
                f'τ={num_format(self.tau, 3)})'
                f'μ₀={num_format(self._mu_0, 3)}, '
                f'τ₀={num_format(self.tau_0, 3)}, '
            )
        else:
            return 'Invalid parameterization!'

    def __repr__(self):

        if self._parameterization == 'σ²':
            return (
                f'NormalNormalConjugate('
                f'n={self._n}, '
                f'x_sum={self._x_sum}, '
                f'sigma_sq={self._sigma_sq})'
                f'mu_0={self._mu_0}, '
                f'sigma_0_sq={self._sigma_0_sq}, '
            )
        elif self._parameterization == 'τ':
            return (
                f'NormalNormalConjugate('
                f'n={self._n}, '
                f'x_sum={self._x_sum}, '
                f'tau={self.tau})'
                f'mu_0={self._mu_0}, '
                f'tau_0={self.tau_0}, '
            )
        else:
            return 'Invalid parameterization!'

max_stars_repo_name: vahndi/probability

max_issues_repo_name: vahndi/probability

max_forks_repo_name: vahndi/probability

reason: from numpy


text:
import numpy as np

def handler_particles_collisions(arr, grid, currents, dt, average, pmax, cross_sections, volume_cell, particles_weight, remains, species_mass = None, monitoring = False, group_fn = None):
    """ Handles the collisions between species in DSMC. It works INPLACE.

    Args:
        arr (list): list of arrays, each is associated with a species and has shape : number of particles in species x 5
        grid (np.ndarray): a 2D-array of (shape number of cells x max size x 2) containing for each cell the particles that are inside.
                           Each particle is describes by its 'container index' (or index in arr, or species index) and its index in the associated species container.
        currents (np.ndarray): 1D-array containing the number of particles in each cell.
        dt (float): time step
        average (np.ndarray): 1D-array containing the average number of particles in each cell
        pmax (np.ndarray): 1D-array of floats of size (number of cells) containing for each cell the normalizing maximum probability.
        cross_sections (np.ndarray): 2D-arrays containing on line i, column j the cross section between species i and j.
        volume_cell (float): volume of a cell (structured mesh)
        particles_weight (int): the particles weight (number of real particles, one simulated one represents)
        remains (np.ndarray): when selecting possible candidates to compute collisions from (using Bird formula), the results is not a integer.
                              This is the float part that should be passed to the next time step.
        species_mass (np.ndarray, optional): 1D-array containing the mass of the species, useful when particles are not the same mass. Defaults to None.
        monitoring (bool, optional): if it should track the process. In which case is computed :
                                        - tracking = ['cell_idx','max_proba','mean_proba','mean_number_of_particles','mean_distance'] - shape : nb_cells x 5
                                        - collisions = ['cell_idx','° couples','quantity'] - shape : nb_cells*nb_species**2 x 3
                                    Defaults to False.
        group_fn (function, optional): Function giving the groups for a colliding couples (in order to be saved as such).
                                       This allows to know what collisions occured (between which species). If monitoring is True, it should not be None.
                                       Defaults to None.

    Returns:
        [np.ndarray (optional), np.ndarray (optional)]: arrays describing the collisions process as described in *monitoring*.

    TODO :
        - it's a little counter intuitive to have average being updated outside this function but pmax inside. Even if there is logical grounds.
    """
    remains[:], cands = candidates(currents, dt, average, pmax, volume_cell, particles_weight, remains)
    nb_cells = grid.shape[0]
    nb_species = len(arr)
    masses = None
    if(monitoring):
        # setting collisions
        nb_groups = (nb_species*(nb_species+1))//2
        n = nb_cells*nb_groups
        collisions = np.zeros((n,3)) # ['cell_idx','° couples','quantity'] # per cell and colliding couples, thus : nb_cells x nb_species**2 lines
        # initializing collisions - we could do it outside ...
        groups_list = np.arange(nb_groups)
        for k in range(nb_cells):
            collisions[k*nb_groups:(k+1)*nb_groups] = k
            collisions[k*nb_groups:(k+1)*nb_groups,1] = groups_list

        # setting tracking
        tracking = np.zeros((nb_cells, 5)) # ['cell_idx','max_proba','mean_proba','mean_number_of_particles','mean_distance'], thus nb_cells lines

    count_collisions = 0
    for idx, k in enumerate(currents): # TODO : parallelize # looping over cells right now
        if(cands[idx]>0):
            choice = index_choosen_couples(currents[idx], cands[idx]) # returns couples of colliding particles indexes in GRID - this is why it works !

            g = grid[idx]
            parts = np.array([[g[c[0]], g[c[1]]] for c in choice], dtype = int)
            array = np.array([[ arr[c[0,0]][c[0,1]] , arr[c[1,0]][c[1,1]] ] for c in parts]) # at this point arrays is an array of couples of [idx_container, idx_in_container]

            # selecting the right cross_sections
            if(np.isscalar(cross_sections)):
                cross_sections_couples = cross_sections
            else:
                cross_sections_couples = np.array([cross_sections[c[0,0], c[1,0]] for c in parts])

            vr_norm = np.linalg.norm((array[:,1,2:]-array[:,0,2:]), axis = 1)
            d = np.linalg.norm((array[:,1,:2]-array[:,0,:2]), axis = 1)

            proba = probability(vr_norm = vr_norm, pmax = pmax[idx], cross_sections = cross_sections_couples)

            max_proba = np.max(proba)

            if(max_proba>1):
                pmax[idx] = max_proba*pmax[idx]

            colliding_couples = is_colliding(proba) # indexes in array of the ACTUALLY colliding couples

            colliding_array = array[colliding_couples]
            count_collisions+= colliding_couples.shape[0]

            if(monitoring):
                if(colliding_couples.shape[0]>0):
                    # a 2D-array with [[[c1,i1],[c2,i2]], [[c3,i3],[c4,i4]], ...] ; c for container index, i for index in container
                    # ideally, we would like to do a groupby on container indexes and then count the number
                    # this is the use of group_fn, defined as a lambda function of 'arr' using first *set_groups* then *get_groups*
                    # and outside this function as it can be setup prior to simulation
                    groups = group_fn(colliding_array.astype(int)) # once we have the groups we can count the quantity for each groups, and save it
                    unique, counts = np.unique(groups, return_counts=True) # so there we may not have all groups in unique
                    collisions[(idx*nb_groups+unique).astype(int),2] = counts

                    # tracking - we could add more than that easily
                    tracking[idx,:] = np.array([idx, pmax[idx], np.mean(proba[colliding_couples]),\
                        average[idx], np.mean(d[colliding_couples])])
                else:
                    tracking[idx,:] = np.array([idx, pmax[idx], np.nan,\
                                            average[idx], np.nan])

            # a very long process that should be avoided whenever possible
            if(species_mass is not None):
                colliding_parts = parts[colliding_couples]
                masses = np.zeros((colliding_parts.shape[0],2))
                for k, part in enumerate(colliding_parts):
                    masses[k,0] = species_mass[part[0,0]]
                    masses[k,1] = species_mass[part[1,0]]

            array[colliding_couples] = reflect(colliding_array, vr_norm[colliding_couples], masses)

            for k in range(len(array)):
                c1, c2 = array[k,0], array[k,1]
                c = parts[k]
                arr[c[0,0]][c[0,1]][:] = c1 # copy
                arr[c[1,0]][c[1,1]][:] = c2

    if(monitoring):
        return count_collisions, collisions, tracking # in theory it is useless to return pmax
    return count_collisions

def candidates(currents, dt, average, pmax, volume_cell, particles_weight, remains):
    """ Returns the number of candidates couples to perform dsmc collisions between particles. (Note that this formula is for one type of particle only => in fact, it is not)

    Args:
        currents (ndarray - 2D - float): number of particles per cell
        dt (float): time step
        average (ndarray - 2D - float): average number of particle in the cell
        pmax (ndarray - 2D - float): max probability per cell ()
        volume_cell (float or ndarray - 2D - float): volume of a cell
        particles_weight (float): "macro-ratio" - ratio of real particles over macro-particles

    Returns:
        (ndarray - 2D - float, ndarray - 2D - int) : the number of candidates to select per cell to perform collisions - fractional part first and then int part.
    """
    # for one type of particle for now
    remains, cands = np.modf(0.5*currents*average*pmax*particles_weight/volume_cell*dt+remains) # (Nc Nc_avg particles_weight (sigma vr)max dt)*/(2V_c)
    return remains, cands.astype(int)

def index_choosen_couples(current, candidates, verbose = False):
    """ Return as many couple indexes as there is candidates couples of particles that should be processed. The index are in [0, number of particles in current cell].

    Args:
        current (int): number of particles in the current cells
        candidates (int): number of couples to select, which mean 2*candidates particles needs to be selected.
        verbose (bool, optional): If verbose is True, then it will print if there is not enough particles compared to the number that have to be selected.
                                  Note : initially,  *np.random.default_rng().choice* with *replace = False* was used and it returned an error. It was then removed
                                  as the command is not available in all versions of NumPy.
                                  Defaults to False.
        TODO : should the selection of couples be as much as possible without repetitions or is it the right way of doing to randomly select all couples with repetitions ?
    Returns:
        [np.ndarray]: 2D-arrays containing the candidates couples (part1, part2) as indexes in the cell.
    """
    return np.random.randint(low = 0, high = current, size = (candidates,2))

        # previous version - all NumPy version do not have it
    # try :
    #     return np.random.default_rng().choice(current, size = (candidates, 2), replace=False, axis = 1, shuffle = False)
    # except ValueError as e:
    #     if(verbose):
    #         print(e)
    #         print(f'{current} < 2 x {candidates} => Some macro-particles will collide several times during the time step.')
    #     return np.random.default_rng().choice(current, size = (candidates, 2), replace=True, axis = 1, shuffle = False) # we dont need shuffling

def probability(vr_norm, pmax, cross_sections): # still per cell
    """ Returns the collisions probability associated with the possibly colliding couples.

    Args:
        vr_norm (np.ndarray): array of size (number of candidates) with the relative speed norm for each couple.
                              Can be computed using something like : np.linalg.norm((arr[choices][:,1,2:]-arr[choices][:,0,2:]), axis = 1)
        pmax (float): maximum probability for the current cell
        cross_sections (np.ndarray): array of size (number of candidates) containing the cross sections for each couple.
                                     Note that a *float* can also be used and it will return the right result.
                                     It is useful when dealing with constant cross sections accross couples.

    Returns:
        [np.ndarray]: 1D-array of size (number of candidates) containing the probability of collisions for each couple.
    """
    return cross_sections/pmax*vr_norm

def is_colliding(proba):
    """ Returns the indexes where there is in fact collisions (only the indexes).

    Args:
        proba (np.ndarray): the colliding probability for each possibly colliding couples.

    Returns:
        np.ndarray: returns a 1D-array with the indexes of the actually colliding couples.
    """
    r = np.random.random(size = proba.shape)
    return np.where(proba>r)[0]

def reflect(arr, vr_norm, masses = None):
    """ Reflect particles following their collisions. Reflections are for now purely specular (no randomness).

    Args:
        arr (np.ndarray): 3D-array of size (number of actually colliding couples x 2 x 5)
        vr_norm (np.ndarray): relative velocity norm for each colliding couple
        masses (np.ndarray, optional): Mass for each particle, for each couple. If mass is the same for every particle, then use None as it simplified the computations.
                                       Defaults to None.

    Returns:
        np.ndarray  : Returns the new array with velocities having been updated accordingly.
                      Note that *arr* is in fact changed in place, however, we still return it in case it is needed  (in the case of
                      the main algorithm *handler_particles_collisions*, *arr* is a copy and not a reference to the initial array).
                      E.g. arr = array[is_colliding(proba)]
    """

    if(masses is None):
        coeff1, coeff2 = 0.5, 0.5 # same mass
    else:
        mass_sum = masses[:,0]+masses[:,1]
        coeff1, coeff2 = masses[:,0]/mass_sum, masses[:,1]/mass_sum # this time they are 1D-array of size the number of particles
        coeff1 = np.expand_dims(coeff1, axis = 1)
        coeff2 = np.expand_dims(coeff2, axis = 1)

    r = np.random.random(size = (2,arr.shape[0]))
    ctheta = 2*r[0,:]-1
    stheta = np.sqrt(1-ctheta*ctheta)
    phi = 2*np.pi*r[1,:]

    v_cm = coeff1*arr[:,0,2:]+coeff2*arr[:,1,2:]
    v_r_ = np.expand_dims(vr_norm, axis = 1)*np.stack((stheta*np.cos(phi), stheta*np.sin(phi), ctheta),  axis = 1) #

    arr[:,0,2:] = v_cm + coeff1*v_r_
    arr[:,1,2:] = v_cm - coeff2*v_r_

    return arr

# -------------------------------- Groups functions ------------------------ #

# default groups functions
def set_groups(n):
    """ Associate for each possibile binart collisions an unique index. This allows the tracking.
        The goal is : when monitoring dsmc, for each actually colliding couples we want to save which species collided.
        Another way, instead of saving which species collided, is to save an identifier of the collision (e.g. for 'I' and 'I-', the id is 3).
        The pros are :
            - less savings
            - no difference between (A,B) and (B,A)
            - it is vectorized : if we were to recognize that (A,B) is the same than (B,A) then, it would not be anymore vectorized.
        Of course, it means that when post-processing the data, we have to know which species are associated to the id.
        Note : it's a little bit more complicated for not much improvement because we could have saved simply the species and post-process everything afterwards.
    Args:
        n (int): number of species

    Returns:
        np.ndarray: symmetric matrix with with an identifier for each collision.
    """
    groups = np.zeros((n, n))
    count = 0
    for i in range(n):
        for j in range(i, n):
            groups[i,j] = count
            groups[j,i] = count
            count +=1
    return groups

def get_groups(arr, groups):
    """ Function that returns, the groups associated with the couples given in arr.

    Args:
        arr (np.ndarray): 3D-array of size (candidates x 2 x 2) where a couple is described as : [part1, part2] where part1 = [species_id_1, index_in_array_1], part2 = [species_id_2, index_in_array_2]
        groups (np.ndarray): symmetric matrix with with an identifier for each collision.

    Returns:
        np.ndarray : 1D-array of size (candidates) containing the identifiers for each collision type.
    """
    return groups[arr[:,0,0],arr[:,1,0]]

max_stars_repo_name: Quettle/lppydsmc

max_issues_repo_name: Quettle/lppydsmc

max_forks_repo_name: Quettle/lppydsmc


text:
#  Copyright (c) 2021, Greg Michael
#  Licensed under BSD 3-Clause License. See LICENSE.txt for details.

import unittest
from unittest.mock import patch
import textwrap

import io
import numpy as np

import craterstats as cst
import craterstats.gm as gm

class TestPlottingClasses(unittest.TestCase):

    root = cst.gm.filename(cst.__file__, 'p')
    file_fns = root + 'config/functions.txt'
    cf = cst.Chronologyfn(file_fns, 'Mars, Hartmann & Neukum (2001)')
    pf = cst.Productionfn(file_fns, 'Mars, Ivanov (2001)')
    cc = cst.Cratercount(root + 'sample/Pickering.scc')

    def test_Craterplotset_CreatePlotSpace(self):
        cps=cst.Craterplotset()
        cps.CreatePlotSpace()
        self.assertTrue(cps.fig)
        self.assertTrue(cps.ax)

    def test_Craterplot_get_data_range(self):
        cp = cst.Craterplot(cratercount=self.cc)
        cps = cst.Craterplotset(craterplot=[cp])

        res = cp.get_data_range(cps)
        self.assertTrue(np.allclose(res,(.1,2.7,6.5e-4,8.2),rtol=.05))

    def test_Craterplotset_autoscale(self):
        cp = cst.Craterplot(cratercount=self.cc)
        cps = cst.Craterplotset(craterplot=[cp])

        cps.autoscale()
        self.assertEqual(list(cps.xrange)+list(cps.yrange),[-2, 2, -5, 3])

    def test_Craterplotset_summary(self):
        cp = cst.Craterplot(cratercount=self.cc,range=[.2,.7],type='d-fit')
        cps = cst.Craterplotset(cf=self.cf,pf=self.pf,craterplot=[cp])

        def get_cps_summary():
            with patch('sys.stdout', new_callable=io.StringIO) as m:
                cps.create_summary_table()
                table=m.getvalue().split('\n')
                return table[1][:145] # trim off source file

        # the following results are not fundamental, but verified against CraterstatsII (see demo plots):
        # name,area,binning,d_min,d_max,method,resurf,n,n_event,age,age-,age+,a0,a0-,a0+,N(1)
        self.assertEqual(get_cps_summary(),
            'Pickering                  3036.6 pseudo-log   0.2   0.7     d-fit      0   313.0       313   0.668   0.613   0.722 -3.488 -3.525 -3.453 3.25e-04')
        cp.UpdateSettings(binning='10/decade')
        self.assertEqual(get_cps_summary(),
            'Pickering                  3036.6  10/decade   0.2  0.79     d-fit      0   313.0       313   0.645   0.593   0.697 -3.503 -3.539 -3.469 3.14e-04')
        cp.UpdateSettings(type='c-fit', resurf=1, binning='pseudo-log')
        self.assertEqual(get_cps_summary(),
            'Pickering                  3036.6 pseudo-log   0.2   0.7     c-fit      1   313.0       313   0.691   0.653   0.729 -3.472 -3.497 -3.449 3.37e-04')
        cp.UpdateSettings(type='poisson', range=[.22,.43], resurf=0)
        self.assertEqual(get_cps_summary(),
            'Pickering                  3036.6 pseudo-log  0.22  0.43   poisson      0     223       223   0.684    0.64   0.731 -3.477 -3.506 -3.448 3.34e-04')
        cp.UpdateSettings(type='poisson', range=[.22, .43], binning='none')
        self.assertEqual(get_cps_summary(),
            'Pickering                  3036.6       none  0.22  0.43   poisson      0     223       223   0.684    0.64   0.731 -3.477 -3.506 -3.448 3.34e-04')

max_stars_repo_name: AndrewAnnex/craterstats

max_issues_repo_name: AndrewAnnex/craterstats

max_forks_repo_name: AndrewAnnex/craterstats


text:
import numpy as np

# parameters ranges

min_aoa = 7
max_aoa = 25
min_Re = 20
max_Re = 2000

naca4_max_camber_min = 0
naca4_max_camber_max = 9.5
naca4_max_camber_position_min = 0
naca4_max_camber_position_max = 90
naca4_thickness_min = 1
naca4_thickness_max = 40

L = 1 # chord line 1 m
nu = 1e-05 # air cinematic viscosity m²/s

def make_aoa():
    return np.random.uniform(low=0, high=30, size=1)[0]
    if np.random.uniform(low=0, high=1, size=1)[0] < 0.75:
        return abs(np.random.normal(7.2, 4**0.5, size=1)[0])
    else:
        # bimodal
        if np.random.uniform(low=0, high=1, size=1)[0] < 0.33:
            return abs(np.random.normal(2, 2**0.5, size=1)[0])
        elif np.random.uniform(low=0, high=1, size=1)[0] > 0.66:
            return abs(np.random.normal(22, 4**0.5, size=1)[0])
        else:
            return np.random.normal((min_aoa+max_aoa)/2, 4.5**0.5, size=1)[0]

def make_Re():
    #return 1000
    # bimodal
    #if np.random.uniform(low=0, high=1, size=1)[0] > 0.75:
    #    return np.random.uniform(low=min_Re, high=min_Re*10, size=1)[0]
    #else:
    return np.random.uniform(low=min_Re, high=max_Re, size=1)[0]

def make_naca4():
    return [
        4, #np.random.uniform(low=naca4_max_camber_min, high=naca4_max_camber_max, size=1)[0],
        40, #np.random.uniform(low=naca4_max_camber_position_min, high=naca4_max_camber_position_max, size=1)[0],
        18 #np.random.uniform(low=naca4_thickness_min, high=naca4_thickness_max, size=1)[0],
    ]

def Re_to_U(Re, _L=L, _nu=nu):
    u = Re*_nu/_L
    print("Freestream vel (m/s): " + str(u))
    return Re*_nu/_L

if __name__ == '__main__':
    for i in range(0,30):
        aoa = make_aoa()
        print(aoa)


text:
import toml
import argparse
import numpy as np

from scipy.stats import entropy

from pom import POM
from sample_script import get_points_covered_by_lidar_config

def evaluate(map, pom_params, lidar_params, config):
        points = get_points_covered_by_lidar_config(
            pom_params, lidar_params, config,  lidar_params['lidar_nos']
        )

        H_entropy = 0.0 # -plogp-(1-p)log(1-p)

        total_entropy = 0.0
        for x in map:
            for xy in x:
                for xyz in xy:
                    if xyz == 0.0 or xyz == 1:
                        continue
                    total_entropy += entropy([1 - xyz, xyz])

        for point in range(len(points)):
            p = map[points[point]]
            if p == 0.0 or p == 1:
                continue
            H_entropy += entropy([1-p,p])

        return H_entropy, total_entropy, total_entropy - H_entropy

if __name__ == "__main__":
    parser = argparse.ArgumentParser()

    parser.add_argument('-p','--params', type=str, default="multihyper.toml", help="Params")
    parser.add_argument('-c','--configuration', type=str, default="config.toml", help="Configuration")

    args = parser.parse_args()
    params = toml.load(args.params)
    configs = toml.load(args.configuration)['config']

    pom_car, num_valid_frames_car = POM(
        random=True, pom_params=params["pom"], lidar_params=params["lidar"]
    ).create_data_from_logs(
        "./routes/square/vehicle"
    )
    print(111)
    pom_car = pom_car.astype(float) / num_valid_frames_car

    pom_ped, num_valid_frames_ped = POM(
        random=True, pom_params=params["pom"], lidar_params=params["lidar"]
    ).create_data_from_logs(
        "./routes/square/pedestrian"
    )

    pom_ped = pom_ped.astype(float) / num_valid_frames_ped

    pom_cyc, num_valid_frames_cyc = POM(
        random=True, pom_params=params["pom"], lidar_params=params["lidar"]
    ).create_data_from_logs(
        "./routes/square/cyclist"
    )

    pom_cyc = pom_cyc.astype(float) / num_valid_frames_cyc
    type = ['square', 'center', 'line', 'pyramid', 'trapezoid', 'line_roll', 'pyramid_roll',
            'pyramid_pitch']
    pom_list = [('car', pom_car), ('ped', pom_ped), ('cyc', pom_cyc)]
    for key, config in configs.items():
        for pom in pom_list:
            H_entropy, total_entropy, IG = evaluate(pom[1], params['pom'], params['lidar'], config)
            print(
                f"Key {type[int(key)]}, {pom[0]}: H_entropy {H_entropy}, total_entropy {total_entropy}, IG {IG}")


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Oct 24 09:53:59 2019

@authors: Dan Wendling

Last modified: 2020-01-11

------------------------------------------------------
 ** Semantic Search Analysis: Generate suggestions **
------------------------------------------------------

This script: For items that were not matched, derive "ListToCheck" - the highest
frequency unmatched terms, that you have time to work on, and then generate
options from MetaMapLite, CSpell and FuzzyWuzzy. Now that the safe bets have
been tagged and removed from view, let's run more liberal guessing and use
human review to resolve SOME trickier tagging.

RECOMMENDATION: Run the processes for MetaMapLite and CSpell, put the output
files in place, then run this whole script at once.

Consider how much work will benefit you. Solving 75-80% of your total search
volume, might be the limit to what you can accomplish in a reasonable amount
of time. Terms with multiple requests are more likely to be real things, such
as your product or service terminology. Items searched only once or twice may
not be classify-able. Real examples: "halloween genetics"; "yoga nutrition". ??

INPUTS:
    - data/interim/UnmatchedAfterMetathesaurus.xlsx (or current file of unmatched)
    - data/interim/result_mm.txt
    - data/interim/result_cspell.txt
    - data/matchFiles/SemanticNetworkReference.xlsx

OUTPUTS:
    - data/interim/03_suggestions.xlsx

----------------
SCRIPT CONTENTS
----------------

1. Start-up / What to put into place, where
2. Build suggestions from MetaMap Lite API, suggestionsFromMM
3. Build suggestions from CSpell spelling checker, suggestionsFromCSpell
4. Build suggestions from FuzzyWuzzy, suggestionsFromFuzzyLoose
5. Unite all suggestions into one Excel file for manual work
6. Update multiple resources, re-run Phase 1, etc.
"""

#%%
# ============================================
# 1. Start-up / What to put into place, where
# ============================================

import pandas as pd
import numpy as np
import os
# For MetaMap
import sys
import argparse
import requests # http://python-requests.org/
# For fuzzy matching
from fuzzywuzzy import fuzz, process

from pathlib import *
# To be used with str(Path.home())

# Set working directory and directories for read/write
home_folder = str(Path.home()) # os.path.expanduser('~')
os.chdir(home_folder + '/Projects/classifysearches')

dataMatchFiles = 'data/matchFiles/' # Permanent helper files; both reading and writing required
dataInterim = 'data/interim/' # Save to disk as desired, to re-start easily

# Open file of unmatched terms
unmatchedTerms = pd.read_excel(dataInterim + 'UnmatchedAfterMetathesaurus.xlsx')
unmatchedTerms.columns
'''
'AdjustedQueryTerm', 'TotalSearchFreq'
'''

'''
Eyeball the df and determine the appropriate cutoff; focus on solving the most
frequently searched items THAT YOU HAVE TIME FOR. Terms searched more than 5
times, the top 500 unresolved per month - something along those lines.
RECOMMENDATION: Start with the top 100 terms by frequency and run them all the
way through to updating the match files; see you can improve the earlier
processing.
'''

# Pilot site, month analyzed: 1148 untagged terms were searched 5 times or more.
ListToCheck = unmatchedTerms.iloc[0:1148].reset_index()
# ListToCheck = unmatchedTerms.loc[(unmatchedTerms['TotalSearchFreq'] >= 5)].reset_index()
ListToCheck = ListToCheck.rename(columns={'index': 'tempID'})
ListToCheck.drop(['TotalSearchFreq'], axis=1, inplace=True)

# Send to file
ListToCheck.to_csv(r'data/interim/SuggestionsNeeded.txt', header=False, index=False, sep='|')

#%%
# ===============================================================
# 2. Build suggestions from MetaMap Lite API, suggestionsFromMM
# ===============================================================
'''
The version of this for license holders may be better, but requires a
license agrement and local installation. The code below is for an https-only
API. (Switch later if you need an extra boost.)

RECOMMENDATION: To see what MetaMapLite is doing, go to
https://ii-public1.nlm.nih.gov/metamaplite/rest and run a search for smallpox.
The result will be:

C0037354: smallpox: Smallpox : [dsyn]: [LCH, NCI_NICHD, MTH, CHV, CSP, MSH, MEDLINEPLUS, LCH_NW, NCI, ICD9CM, NDFRT, DXP]: 0 : 8 : 0.0
C0037355: smallpox: Smallpox Vaccine : [imft, phsu]: [MTH, LCH, CHV, CSP, MSH, LCH_NW, RXNORM, NDFRT, HL7V2.5, HL7V3.0, VANDF]: 0 : 8 : 0.0
C0037356: smallpox: Smallpox Viruses : [virs]: [MTH, CHV, LNC, NCI_CDISC, CSP, MSH, NCI, NCBI]: 0 : 8 : 0.0

If the best match for a search term is the vaccine entry, what's important to
you will be:
    C0037355 - the CUI
    Smallpox Vaccine - the Preferred term
    [imft, phsu] - the Semantic Type abbreviations

In this example, you have classified your term as both Immunologic Factor (imft)
and Pharmacologic Substance (phsu).

The MetaMap home page is https://metamap.nlm.nih.gov

The most important change from NLM's MetaMap script is this line, which allows
you to link dataframe terms with MetaMap output (sldiwi):
    '--docformat', default='sldiwi',
'''

# -----------------------------------
# PROCESS THE FILE OUTSIDE OF PYTHON
# -----------------------------------

'''
Go to the terminal / command line and run mmlrestclient.py, like this:

cd ~/Projects/classifysearches/data/interim

python mmlrestclient.py https://ii-public1.nlm.nih.gov/metamaplite/rest/annotate SuggestionsNeeded.txt --output result_mm.txt
'''

# ------------------------------------
# Merge source and result into one df
# ------------------------------------

'''
89|MMI|0.46|Calliphoridae|C0322495|[euka]|"blowfly"-text-0-"blowfly"-NN-0|0/7|
C0037354: smallpox: Smallpox : [dsyn]: [LCH, NCI_NICHD, MTH, CHV, CSP, MSH, MEDLINEPLUS, LCH_NW, NCI, ICD9CM, NDFRT, DXP]: 0 : 8 : 0.0

Described at https://github.com/NCBI-Codeathons/Use-UMLS-and-Python-to-classify-website-visitor-queries-into-measurable-categories/wiki/4.-Understanding-MetaMap-output
'''

mmCols = ['tempID', 'Score', 'SuggestedTerm', 'CUI', 'SemTypeList']
result_mm = pd.read_csv(dataInterim + 'result_mm.txt', sep='|', usecols=[0,2,3,4,5], names=mmCols)

# Join
suggestionsFromMM = pd.merge(ListToCheck, result_mm, how='left', on=['tempID'])
suggestionsFromMM.columns
'''
'tempID', 'AdjustedQueryTerm', 'Score', 'SuggestedTerm', 'CUI',
       'SemTypeList'
'''

# Drop
suggestionsFromMM.drop(['tempID'], axis=1, inplace=True)
suggestionsFromMM = suggestionsFromMM.dropna(subset=['SuggestedTerm'])

# Add source col
suggestionsFromMM['Source'] = 'MetaMap'

# Make SemTypeList content look like SemanticNetworkReference.SemanticTypeAbr
suggestionsFromMM['SemTypeList'] = suggestionsFromMM['SemTypeList'].str.replace(', ', '|')
suggestionsFromMM['SemTypeList'] = suggestionsFromMM['SemTypeList'].str.replace('\[', '')
suggestionsFromMM['SemTypeList'] = suggestionsFromMM['SemTypeList'].str.replace('\]', '')

# Bring in Sem Type info
SemTypes = pd.read_excel(dataMatchFiles + 'SemanticNetworkReference.xlsx')
# Reduce to join cols
SemTypes = SemTypes[['SemanticTypeAbr', 'SemanticType']]
# dropna
SemTypes = SemTypes.dropna(subset=['SemanticTypeAbr'])

# Join
suggestionsFromMM2 = pd.merge(suggestionsFromMM, SemTypes, how='left', left_on='SemTypeList', right_on='SemanticTypeAbr')
suggestionsFromMM2.columns
'''
'AdjustedQueryTerm', 'Score', 'SuggestedTerm', 'CUI', 'SemTypeList',
       'Source', 'SemanticTypeAbr', 'SemanticType'
'''

# Some combo types are not in SemanticNetworkReference
# Pull out nulls
fixSem = suggestionsFromMM2[suggestionsFromMM2['SemanticType'].isnull()]
# Dupe off col
fixSem['SemanticType'] = fixSem['SemTypeList']

# Move abbreviation into index
SemTypes.set_index('SemanticTypeAbr', inplace=True)

# Find and replace within fixSem using SemTypes df
for key, value in SemTypes.iterrows():
    currKey = key
    currST = value['SemanticType']
    fixSem['SemanticType'] = fixSem['SemanticType'].str.replace(currKey, currST)

#  drop na from main df
suggestionsFromMM2 = suggestionsFromMM2.dropna(subset=['SemanticType'])

# append the new content
suggestionsFromMM2 = suggestionsFromMM2.append([fixSem])

# drop abbrev cols
suggestionsFromMM2.drop(['SemanticTypeAbr', 'SemTypeList'], axis=1, inplace=True)
suggestionsFromMM2 = suggestionsFromMM2.reset_index(drop=True)

del [[result_mm, suggestionsFromMM, fixSem]]

#%%
# =========================================================================
# 3. Build suggestions from CSpell spelling checker, suggestionsFromCSpell
# =========================================================================
'''
See wiki for installation, help, etc.:
    https://github.com/NCBI-Codeathons/Use-UMLS-and-Python-to-classify-website-visitor-queries-into-measurable-categories/wiki/3.-Installing-and-running-CSpell

Install the Java application CSpell, "a distributable spell checker for consumer
language," from https://lsg3.nlm.nih.gov/LexSysGroup/Projects/cSpell/current/web/index.html

Specifically, https://lsg3.nlm.nih.gov/LexSysGroup/Projects/cSpell/current/web/download.html
No data yet on file size limits; 10,000 rows can be processed in around 12
minutes on a normal workstation, which uses these options:

cspell -I:cspell_infile.txt -si -o:result_cspell.txt

A 30,000-row file was processed in around 1 hour. RECOMMENDATION: Limit size
to what fits into your available work time. Here we use the same highest-freq
unmatched terms, as above.

# reduce to one col
cspell_infile = ListToCheck[['AdjustedQueryTerm']]

# Write out
cspell_infile.to_csv(dataInterim + 'cspell_infile.txt', encoding='utf-8', index=False)
'''

# ---------------------------------------
# NOW PROCESS THE FILE OUTSIDE OF PYTHON
# ---------------------------------------
'''
Run from terminal command line; location similar to:

cd C:\tools\cSpell2018\bin

cspell -i:SuggestionsNeeded.txt -si -o:result_cspell.txt
'''

# ------------------------------------
# Merge source and result into one df
# ------------------------------------

cspellCols = ['AdjustedQueryTerm', 'SuggestedTerm']
result_cspell = pd.read_csv(dataInterim + 'result_cspell.txt', sep='|', index_col=False, usecols=[1,3], names=cspellCols)  # skiprows=7,
result_cspell.columns
'''
'AdjustedQueryTerm', 'SuggestedTerm'
'''

# Drop rows where input and output are the same
# https://stackoverflow.com/questions/43951558/remove-rows-that-two-columns-have-the-same-values-by-pandas
suggestionsFromCSpell = result_cspell[result_cspell['AdjustedQueryTerm'] != result_cspell['SuggestedTerm']]

# Add source col
suggestionsFromCSpell['Source'] = 'CSpell'

del [[result_cspell]]

#%%
# ====================================================================
# 3. Build suggestions from fuzzy matching, suggestionsFromFuzzyLoose
# ====================================================================
'''
Over time you can change the parameters to match your time and desired level
of effort. You can reduce the list, change the type of match (full phrase or
any word), and change the score, to change the number of candidates to
match how much time you want to spend making manual assignments.

When starting with a new site you should probably spend a good deal of time
here, to make connections the other steps can't make. Decisions you make
here will provide training data that the machine learning component can use.

Some configuration options are described at
https://www.neudesic.com/blog/fuzzywuzzy-using-python/.
See for example fuzz.ratio (conservative) vs. fuzz.partial_ratio (medium) vs.
fuzz.token_set_ratio (any single word in the phrases, very liberal). The more
liberal you get here, the more you will see multiple-concept searches, which
you don't need to see at this point. This is not a good time to solve those.

5,000 in ~25 minutes; ~10,000 in ~50 minutes (at score_cutoff=85)

Python's FuzzyWuzzy was written for single inputs to a web form; here, however,
we use it to compare one dataframe column to another dataframe's column.
Takes extra lines of code to match the tokenized function output back
to the original untokenized term, which is necessary for this work.

# Quick test, if you want - punctuation difference
fuzz.ratio('Testing FuzzyWuzzy', 'Testing FuzzyWuzzy!!')

FuzzyWuzzyResults - What the results of this function mean:
('hippocratic oath', 100, 2987)
('Best match string from dataset_2' (PastMatches), 'Score of best match', 'Index of best match string in PastMatches')

Re-start:
uniqueUnassignedAfterStep2 = pd.read_excel(dataInterim + 'uniqueUnassignedAfterStep2.xlsx')
'''

# -----
# Prep
# -----

# Bring in
SiteSpecificMatches = pd.read_excel(dataMatchFiles + 'SiteSpecificMatches.xlsx')
SiteSpecificMatches.columns
'''
'AdjustedQueryTerm', 'PreferredTerm', 'SemanticType'
'''

# Bring in
PastMatches = pd.read_excel(dataMatchFiles + 'PastMatches.xlsx')
PastMatches.columns
'''
'SemanticType', 'AdjustedQueryTerm', 'PreferredTerm', 'ui'
'''

# Append
TermsInUseFull = SiteSpecificMatches.append([PastMatches])

# Reduce
TermsInUseList = TermsInUseFull[['AdjustedQueryTerm']]

# Dedupe
TermsInUseList = TermsInUseList.drop_duplicates(subset=['AdjustedQueryTerm'], keep='first')

del [[SiteSpecificMatches, PastMatches]] # free memory

# ------------------------------------
# Fuzzy match using fuzz.ratio
# ------------------------------------

# Compare ListToCheck to TermsInUseList
def fuzzy_match(x, choices, scorer, cutoff):
    return process.extractOne(
        x, choices=choices, scorer=scorer, score_cutoff=cutoff
    )

# Create series FuzzyWuzzyResults
FuzzyWuzzyProcResult1 = ListToCheck.loc[:, 'AdjustedQueryTerm'].apply(
        fuzzy_match,
    args=( TermsInUseList.loc[:, 'AdjustedQueryTerm'],
            fuzz.ratio,
            75 # Items must have this score or higher to appear in the results
        )
)

# Convert FuzzyWuzzyResults Series to df
FuzzyWuzzyProcResult2 = pd.DataFrame(FuzzyWuzzyProcResult1)

# Move Index (IDs) into 'FuzzyIndex' col because Index values will be discarded
FuzzyWuzzyProcResult2 = FuzzyWuzzyProcResult2.reset_index()
FuzzyWuzzyProcResult2 = FuzzyWuzzyProcResult2.rename(columns={'index': 'FuzzyIndex'})

# Remove nulls
FuzzyWuzzyProcResult2 = FuzzyWuzzyProcResult2[FuzzyWuzzyProcResult2.AdjustedQueryTerm.notnull() == True] # remove nulls

# Move tuple output into 3 cols
FuzzyWuzzyProcResult2[['SuggestedTerm', 'Score', 'TermsInUseListIndex']] = FuzzyWuzzyProcResult2['AdjustedQueryTerm'].apply(pd.Series)
FuzzyWuzzyProcResult2.drop(['AdjustedQueryTerm'], axis=1, inplace=True) # drop tuples

# Merge result to the orig source list cols
FuzzyWuzzyProcResult3 = pd.merge(FuzzyWuzzyProcResult2, ListToCheck, how='left', left_index=True, right_index=True)
FuzzyWuzzyProcResult3.columns
'''
'FuzzyIndex', 'SuggestedTerm', 'Score', 'TermsInUseListIndex',
       'tempID', 'AdjustedQueryTerm'
'''

# Reduce/re-order
suggestionsFromFuzzyWuzzy = FuzzyWuzzyProcResult3[['AdjustedQueryTerm', 'SuggestedTerm', 'Score']]

# Add source info
suggestionsFromFuzzyWuzzy['Source'] = 'FuzzyWuzzy'

# In case the suggestion is correct, get PreferredTerm and Sem Type
suggestionsFromFuzzyWuzzy2 = pd.merge(suggestionsFromFuzzyWuzzy, TermsInUseFull, how='left', left_on=['SuggestedTerm'], right_on='AdjustedQueryTerm')
suggestionsFromFuzzyWuzzy2.columns
'''
'AdjustedQueryTerm_x', 'SuggestedTerm', 'Score', 'Source',
       'AdjustedQueryTerm_y', 'PreferredTerm', 'SemanticType', 'ui'
'''

# drop _y and rename _x
suggestionsFromFuzzyWuzzy2.drop(['AdjustedQueryTerm_y'], axis=1, inplace=True)
suggestionsFromFuzzyWuzzy2.rename(columns={'AdjustedQueryTerm_x': 'AdjustedQueryTerm'}, inplace=True)

# del [[]]

#%%
# =============================================================
# 5. Unite all suggestions into one Excel file for manual work
# =============================================================
'''
Append into one table and prep for manual work.
'''

# Append into one df
suggestions = suggestionsFromMM2.append([suggestionsFromCSpell, suggestionsFromFuzzyWuzzy2])
suggestions.columns
'''
'AdjustedQueryTerm', 'CUI', 'Score', 'SemanticType', 'Source',
       'SuggestedTerm'
'''

# Add frequency from unmatchedTerms
unmatchedTerms = pd.read_excel(dataInterim + 'UnmatchedAfterMetathesaurus.xlsx')
unmatchedTerms.columns
'''
'AdjustedQueryTerm', 'TotalSearchFreq'
'''

# Join
suggestions2 = pd.merge(suggestions, unmatchedTerms, how='left', left_on=['AdjustedQueryTerm'], right_on='AdjustedQueryTerm')
suggestions2.columns
'''
'AdjustedQueryTerm', 'SuggestedTerm', 'SemanticType', 'CUI', 'Score',
       'Source', 'TotalSearchFreq'
'''

# Adjust
suggestions2['UseYN'] = ''
suggestions2 = suggestions2[['UseYN', 'TotalSearchFreq', 'AdjustedQueryTerm', 'SuggestedTerm', 'SemanticType', 'CUI', 'Source', 'Score']]
suggestions2 = suggestions2.sort_values(by=['TotalSearchFreq', 'AdjustedQueryTerm', 'SemanticType'], ascending=[False,True,True])
suggestions2 = suggestions2.reset_index(drop=True)

del [[ListToCheck, SemTypes, suggestions, suggestionsFromCSpell, suggestionsFromFuzzyWuzzy,
      suggestionsFromFuzzyWuzzy2, suggestionsFromMM2, unmatchedTerms]]
# SemanticNetworkReference, TermsInUse,

#%%
# =============================================================
# 6. Update multiple resources, re-run Phase 1, etc.
# =============================================================
'''
Some terms should be added to the match files; do this work and re-run Phase 1.

For now we have to do this manually.

A wireframe for a suggestion-selector user interface is in the GitHub wiki.

Do what you have time for.
'''

# Write out
writer = pd.ExcelWriter(dataInterim + '03_suggestions.xlsx')
suggestions2.to_excel(writer,'Suggestions', index=False)
# df2.to_excel(writer,'Sheet2')
writer.save()


text:
import numpy as np
import torch
import torch.distributed as dist

def tensor(x, device):
    if isinstance(x, torch.Tensor):
        return x.to(device)

    x = np.asarray(x, dtype=np.float)
    x = torch.tensor(x, device=device, dtype=torch.float32)
    return x

def input_preprocessing(x, device):
    x = tensor(x, device)
    x = x.float()
    x /= 255.0
    return x

def to_np(t):
    return t.cpu().detach().numpy()

def random_seed(seed=None):
    np.random.seed(seed)
    torch.manual_seed(np.random.randint(int(1e6)))

def restore_model(model, save_path):
    checkpoint = torch.load(save_path)
    model.network.load_state_dict(checkpoint["model_state_dict"])
    model.optimizer.load_state_dict(checkpoint["optimizer_state_dict"])
    update = checkpoint["update"]
    return update

def sync_initial_weights(model):
    for param in model.parameters():
        dist.broadcast(param.data, src=0)

def sync_gradients(model):
    for param in model.parameters():
        dist.all_reduce(param.grad.data, op=dist.ReduceOp.SUM)

def cleanup():
    dist.destroy_process_group()

def sync_values(tensor_sum_values, tensor_nb_values):
    dist.reduce(tensor_sum_values, dst=0)
    dist.reduce(tensor_nb_values, dst=0)
    return tensor_sum_values / tensor_nb_values

def range_tensor(t, device):
    return torch.arange(t).long().to(device)

def zeros(shape, dtype):
    """Attempt to return torch tensor of zeros, or if numpy dtype provided,
    return numpy array or zeros."""
    try:
        return torch.zeros(shape, dtype=dtype)
    except TypeError:
        return np.zeros(shape, dtype=dtype)


text:
#!/usr/bin/env python
#
import unittest
from femagtools import moo
import numpy as np
import matplotlib as mp
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d as mpl
import matplotlib.cm as cm

class FesProblem(moo.Problem):
    def __init__(self):
        super(FesProblem, self).__init__(4, 0, 2)

    def objfun(self, x):
        D = len(x)
        return (sum([abs(x[i]-np.exp(((i+1.)/D) ** 2)/3) ** 0.5
                    for i in range(D)]),
                sum([(x[i]-0.5*np.cos(10*(i+1)*np.pi/D)-0.5) ** 2
                    for i in range(D)]))

def surface():
    "show a nice 3d plot just for demo purposes"
    prob = FesProblem()
    steps = (15, 15)
    lower = (0.0, 0.0)
    upper = (1.0, 1.0)
    domain = [np.linspace(l, u, s)
              for s, l, u in zip(steps, lower, upper)]

    X = domain[0]
    Y = domain[1]

    L = [len(d) for d in domain]
    s = []
    e = 1
    for d in domain:
        L = L[1:]
        LS = np.prod(L)
        s.append(np.repeat(d.tolist()*LS, e))
        if len(L) > 0:
            e = e*L[0]

    z = np.array([prob.objfun(x) for x in np.array(s).T]).T[1]
    z = np.reshape(z, steps)

    fig = plt.figure()
    ax = mpl.Axes3D(fig)  # fig.add_subplot(111,projection='3d')
    ax.plot_surface(X, Y, z, rstride=1,
                    cstride=1, cmap=cm.jet, linewidth=0)
    plt.show()

class AlgorithmTest(unittest.TestCase):
    def test_evolve(self):
        prob = FesProblem()
        pop = moo.Population(prob, 20)
        pop.eval()

        F = np.array([ind.cur_f for ind in pop.individuals]).T
        #plt.scatter(F[0], F[1], c="b", s=50, alpha=0.5)
        F = np.array([ind.cur_f
                      for ind in pop.individuals if ind.rank == 0]).T
        #plt.scatter(F[0], F[1], c="g", s=100, alpha=0.5)

        algo = moo.Nsga2()
        for i in range(25):
            newpop = algo.evolve(pop)
            newpop.eval()
            pop.merge(newpop)

        #F = np.array([ind.cur_f for ind in pop.individuals]).T
        #plt.scatter(F[0], F[1], c="r", s=200, alpha=0.5)

        #plt.xlabel("$f^{(1)}$")
        #plt.ylabel("$f^{(2)}$")
        #plt.savefig('moo.png')

        # TODO: make these random tests reproducible
        #np.testing.assert_almost_equal(pop.compute_norm_dist(), 1.0 )
        #self.assertEqual(pop.dom_count, [0]*pop.size())
        #self.assertEqual(pop.pareto_rank, [0]*pop.size() )
        #self.assertEqual(pop.best_idx(), [0, 1, 3, 4, 2, 6, 5, 7] )

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

max_stars_repo_name: dapu/femagtools

max_issues_repo_name: dapu/femagtools

max_forks_repo_name: dapu/femagtools


text:
# Copyright 2017 The TensorFlow 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.
# ==============================================================================
"""Tests for cyclegan.train."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

from absl import flags
import numpy as np
import tensorflow as tf

import train

FLAGS = flags.FLAGS
mock = tf.test.mock
tfgan = tf.contrib.gan

def _test_generator(input_images):
  """Simple generator function."""
  return input_images * tf.get_variable('dummy_g', initializer=2.0)

def _test_discriminator(image_batch, unused_conditioning=None):
  """Simple discriminator function."""
  return tf.contrib.layers.flatten(
      image_batch * tf.get_variable('dummy_d', initializer=2.0))

train.networks.generator = _test_generator
train.networks.discriminator = _test_discriminator

class TrainTest(tf.test.TestCase):

  @mock.patch.object(tfgan, 'eval', autospec=True)
  def test_define_model(self, mock_eval):
    FLAGS.batch_size = 2
    images_shape = [FLAGS.batch_size, 4, 4, 3]
    images_x_np = np.zeros(shape=images_shape)
    images_y_np = np.zeros(shape=images_shape)
    images_x = tf.constant(images_x_np, dtype=tf.float32)
    images_y = tf.constant(images_y_np, dtype=tf.float32)

    cyclegan_model = train._define_model(images_x, images_y)
    self.assertIsInstance(cyclegan_model, tfgan.CycleGANModel)
    self.assertShapeEqual(images_x_np, cyclegan_model.reconstructed_x)
    self.assertShapeEqual(images_y_np, cyclegan_model.reconstructed_y)

  @mock.patch.object(train.networks, 'generator', autospec=True)
  @mock.patch.object(train.networks, 'discriminator', autospec=True)
  @mock.patch.object(
      tf.train, 'get_or_create_global_step', autospec=True)
  def test_get_lr(self, mock_get_or_create_global_step,
                  unused_mock_discriminator, unused_mock_generator):
    FLAGS.max_number_of_steps = 10
    base_lr = 0.01
    with self.test_session(use_gpu=True) as sess:
      mock_get_or_create_global_step.return_value = tf.constant(2)
      lr_step2 = sess.run(train._get_lr(base_lr))
      mock_get_or_create_global_step.return_value = tf.constant(9)
      lr_step9 = sess.run(train._get_lr(base_lr))

    self.assertAlmostEqual(base_lr, lr_step2)
    self.assertAlmostEqual(base_lr * 0.2, lr_step9)

  @mock.patch.object(tf.train, 'AdamOptimizer', autospec=True)
  def test_get_optimizer(self, mock_adam_optimizer):
    gen_lr, dis_lr = 0.1, 0.01
    train._get_optimizer(gen_lr=gen_lr, dis_lr=dis_lr)
    mock_adam_optimizer.assert_has_calls([
        mock.call(gen_lr, beta1=mock.ANY, use_locking=True),
        mock.call(dis_lr, beta1=mock.ANY, use_locking=True)
    ])

  @mock.patch.object(tf.summary, 'scalar', autospec=True)
  def test_define_train_ops(self, mock_summary_scalar):
    FLAGS.batch_size = 2
    FLAGS.generator_lr = 0.1
    FLAGS.discriminator_lr = 0.01

    images_shape = [FLAGS.batch_size, 4, 4, 3]
    images_x = tf.zeros(images_shape, dtype=tf.float32)
    images_y = tf.zeros(images_shape, dtype=tf.float32)

    cyclegan_model = train._define_model(images_x, images_y)
    cyclegan_loss = tfgan.cyclegan_loss(
        cyclegan_model, cycle_consistency_loss_weight=10.0)
    train_ops = train._define_train_ops(cyclegan_model, cyclegan_loss)

    self.assertIsInstance(train_ops, tfgan.GANTrainOps)
    mock_summary_scalar.assert_has_calls([
        mock.call('generator_lr', mock.ANY),
        mock.call('discriminator_lr', mock.ANY)
    ])

  @mock.patch.object(tf, 'gfile', autospec=True)
  @mock.patch.object(train, 'data_provider', autospec=True)
  @mock.patch.object(train, '_define_model', autospec=True)
  @mock.patch.object(tfgan, 'cyclegan_loss', autospec=True)
  @mock.patch.object(train, '_define_train_ops', autospec=True)
  @mock.patch.object(tfgan, 'gan_train', autospec=True)
  def test_main(self, mock_gan_train, mock_define_train_ops, mock_cyclegan_loss,
                mock_define_model, mock_data_provider, mock_gfile):
    FLAGS.image_set_x_file_pattern = '/tmp/x/*.jpg'
    FLAGS.image_set_y_file_pattern = '/tmp/y/*.jpg'
    FLAGS.batch_size = 3
    FLAGS.patch_size = 8
    FLAGS.generator_lr = 0.02
    FLAGS.discriminator_lr = 0.3
    FLAGS.train_log_dir = '/tmp/foo'
    FLAGS.master = 'master'
    FLAGS.task = 0
    FLAGS.cycle_consistency_loss_weight = 2.0
    FLAGS.max_number_of_steps = 1

    mock_data_provider.provide_custom_data.return_value = (
        tf.zeros([3, 2, 2, 3], dtype=tf.float32),
        tf.zeros([3, 2, 2, 3], dtype=tf.float32))

    train.main(None)
    mock_data_provider.provide_custom_data.assert_called_once_with(
        ['/tmp/x/*.jpg', '/tmp/y/*.jpg'], batch_size=3, patch_size=8)
    mock_define_model.assert_called_once_with(mock.ANY, mock.ANY)
    mock_cyclegan_loss.assert_called_once_with(
        mock_define_model.return_value,
        cycle_consistency_loss_weight=2.0,
        tensor_pool_fn=mock.ANY)
    mock_define_train_ops.assert_called_once_with(
        mock_define_model.return_value, mock_cyclegan_loss.return_value)
    mock_gan_train.assert_called_once_with(
        mock_define_train_ops.return_value,
        '/tmp/foo',
        get_hooks_fn=mock.ANY,
        hooks=mock.ANY,
        master='master',
        is_chief=True)

if __name__ == '__main__':
  tf.test.main()

max_stars_repo_name: zcdzcdzcd/models

max_issues_repo_name: zcdzcdzcd/models

max_forks_repo_name: zcdzcdzcd/models


text:
import os
import numpy as np
from pwtools import io, parse,common
from pwtools.test.tools import assert_attrs_not_none
from pwtools.constants import Bohr,Ang,Ha,eV
from .testenv import testdir
from pwtools.test.tools import unpack_compressed
pj = common.pj

def test_cp2k_scf():
    attr_lst = parse.Cp2kSCFOutputFile().attr_lst
    for base in ['cp2k.scf.out.print_low', 'cp2k.scf.out.print_medium']:
        fn = 'files/cp2k/scf/%s' %base
        print("testing: %s" %fn)
        print(common.backtick("gunzip %s.gz" %fn))
        st = io.read_cp2k_scf(fn)
        assert_attrs_not_none(st, attr_lst=attr_lst)

def test_cp2k_md():
    attr_lst = parse.Cp2kMDOutputFile().attr_lst
    # This parser and others have get_econst(), but not all, so ATM it's not
    # part of the Trajectory API
    attr_lst.pop(attr_lst.index('econst'))
    for dr in ['files/cp2k/md/npt_f_print_low', 'files/cp2k/md/nvt_print_low']:
        base = os.path.dirname(dr)
        fn = '%s/cp2k.out' %dr
        print("testing: %s" %fn)
        print(common.backtick('tar -C {0} -xzf {1}.tgz'.format(base,dr)))
        tr = io.read_cp2k_md(fn)
        assert_attrs_not_none(tr, attr_lst=attr_lst)
        pp = parse.Cp2kMDOutputFile(fn)
        forces_outfile = pp._get_forces_from_outfile()*Ha/Bohr/eV*Ang
        assert np.allclose(forces_outfile, tr.forces, rtol=1e-3)

def test_cp2k_cell_opt():
    attr_lst = parse.Cp2kRelaxOutputFile().attr_lst
    attr_lst.pop(attr_lst.index('econst'))
    # There is no PROJECT-frc-1.xyz file, but the input file has
    #    &force_eval
    #        &print
    #            &forces
    #            &end forces
    #        &end print
    #    &end force_eval
    # therefore we can parse forces from the outfile.
    none_attrs = ['ekin',
                  'temperature',
                  'timestep',
                  'velocity',
                  ]
    for dr in ['files/cp2k/cell_opt/cell_opt']:
        base = os.path.dirname(dr)
        fn = '%s/cp2k.out' %dr
        print("testing: %s" %fn)
        print(common.backtick('tar -C {0} -xzf {1}.tgz'.format(base,dr)))
        tr = io.read_cp2k_relax(fn)
        assert_attrs_not_none(tr, attr_lst=attr_lst, none_attrs=none_attrs)

def test_cp2k_txt_vs_dcd():
    # Two exactly equal NPT runs, nstep=16, natoms=57. The dcd run uses
    #   motion/print/trajectory format dcd_aligned_cell
    # the other the default xyz format. So we have
    #   PROJECT-pos-1.dcd
    #   PROJECT-pos-1.xyz
    # Since the cell changes and we use dcd_aligned_cell, the coords from the
    # xyz run are NOT the same as the coords in the dcd file, which HAS to be
    # like this. Only coords_frac can be compared, and coords after the xyz
    # run's cell has been aligned to [[x,0,0],[xy,y,0],[xz,yz,z]].
    dir_xyz = unpack_compressed('files/cp2k/dcd/npt_xyz.tgz')
    dir_dcd = unpack_compressed('files/cp2k/dcd/npt_dcd.tgz')
    tr_xyz = io.read_cp2k_md(pj(dir_xyz, 'cp2k.out'))
    tr_dcd = io.read_cp2k_md_dcd(pj(dir_dcd, 'cp2k.out'))
    for name in ['natoms', 'nstep', 'timestep', 'symbols']:
        assert getattr(tr_xyz,name) == getattr(tr_dcd, name)
    assert tr_xyz.timestep == 1.0
    assert tr_xyz.natoms == 57
    assert tr_xyz.nstep == 16

    # coords are 32bit float in dcd files (single prec, so coords_frac can only
    # be accurate to that precision, which ~1e-8). cryst_const is double
    # precision in the dcd file, so atol can be lower
    assert np.allclose(tr_xyz.coords_frac, tr_dcd.coords_frac, rtol=0,
                       atol=1.5e-7)
    assert np.allclose(tr_xyz.cryst_const, tr_dcd.cryst_const, rtol=0,
                       atol=7e-10)
    assert not np.allclose(tr_xyz.coords, tr_dcd.coords, rtol=0,
                           atol=5e-6)
    assert not np.allclose(tr_xyz.cell, tr_dcd.cell, rtol=0,
                           atol=7e-10)

    # align xyz cell, now cell and coords are the same
    tr_xyz.coords = None; tr_xyz.cell=None; tr_xyz.set_all()

    assert np.allclose(tr_xyz.coords_frac, tr_dcd.coords_frac, rtol=0,
                       atol=1.5e-7)
    assert np.allclose(tr_xyz.cryst_const, tr_dcd.cryst_const, rtol=0,
                       atol=7e-10)
    assert np.allclose(tr_xyz.coords, tr_dcd.coords, rtol=0,
                       atol=5e-6)
    assert np.allclose(tr_xyz.cell, tr_dcd.cell, rtol=0,
                       atol=7e-10)

max_stars_repo_name: elcorto/pwtools

max_issues_repo_name: elcorto/pwtools

max_forks_repo_name: elcorto/pwtools