text:
import warnings
from dataclasses import dataclass
from typing import Callable, Dict, Union, List

import numpy as np
from scipy.optimize.slsqp import _epsilon  # noqa

from .modelspec import Parameter, ParameterSet
from .util import ReprMixin

@dataclass
class Configuration(ReprMixin):
    """
    Configuration for the ReMeta toolbox

    Parameters
    ----------
    *** Basic definition of the model ***
    skip_meta : bool
        If False, only fit the sensory level.
    meta_noise_type : str
        Whether the model considers noise at readout or report.
        Possible values: 'noisy_report', 'noisy_readout'
    meta_noise_model : str
        Metacognitive noise distribution.
        Possible valus: 'norm', 'gumbel', 'lognorm', 'lognorm_varstd', 'beta', 'beta_std', 'betaprime', 'gamma',
                        'censored_norm', 'censored_gumbel', 'censored_lognorm', 'censored_lognorm_varstd',
                        'censored_betaprime', 'censored_gamma',
                        'truncated_norm', 'truncated_norm_lookup', 'truncated_norm_fit',
                        'truncated_gumbel', 'truncated_gumbel_lookup',
                        'truncated_lognorm', 'truncated_lognorm_varstd'
    meta_link_function : str
        Metacognitive link function. In case of criterion-based link functions {x} refers to the number of criteria.
        Possible values: 'probability_correct', 'identity', 'tanh', 'normcdf', 'erf', 'alg', 'guder', 'linear',
                         'detection_model_linear', 'detection_model_mean', 'detection_model_mode',
                         'detection_model_full', 'detection_model_ideal'
                         '{x}_criteria', '{x}_criteria_linear', '{x}_criteria_linear_tanh', '{x}_criteria_variable'
    detection_model : bool
        Experimental. Define the model as a detection model.
    detection_model_nchannels : int
        Defines the number of sensory channels in case of a detection model.

    *** Enable or disable specific parameters ***
    * Each setting can take the values 0, 1 or 2:
    *    0: Disable parameter.
    *    1: Enable parameter.
    *    2: Enable parameter and fit separate values for the negative and positive stimulus category (in the case of
    *       sensory parameters, indicated by the suffix '_sens') or negative and positive decision values (in the
    *       case of metacognitive parameters, indicated by the suffix '_meta').
    enable_warping_sens : int (default: 0)
        Fit a non-linear transducer prior to sensory processing (the nonlinear function is defined via
        `function_warping_sens`).
    enable_noise_sens : int (default: 1)
        Fit separate sensory noise parameters for both stimulus categories.
    enable_noise_transform_sens : int (default: 0)
        Fit an additional sensory noise parameter for signal-dependent sensory noise (the type of dependency is is
        defined via `function_noise_transform_sens`).
    enable_thresh_sens : int (default: 0)
        Fit a sensory threshold.
    enable_bias_sens : int (default: 1)
        Fit a sensory bias towards one of the stimulus categories.
    enable_noise_meta : int (default: 1)
        Fit a metacognitive noise parameter
    enable_noise_transform_meta : int (default: 0)
        Fit an additional metacognitive noise parameter for signal-dependent metacognitive noise (the type of dependency
        is is defined via `function_noise_transform_meta`). Note: at present, enable_noise_transform_meta=2 leads to
        biased results and is therefore discouraged.
    enable_readout_term_meta : int (default: 0)
        Fit an additive metacognitive bias at readout.
    enable_slope_meta : int (default: 0)
        Fit a slope parameter for the link function. In the case of a criterion-based link function, multiple
        parameters are fitted for each criterion/confidence level.
    enable_confidence_term_meta : int (default: 0)
        Fit an additive metacognitive bias at report.
    enable_scaling_meta : int (default: 0)
        Fit a confidence scaling parameter.
    enable_criteria_meta : int (default: 0)
        Fit criteria for a criterion-based link function. Note that the number of criteria is set via the link function
        argument (see `meta_link_function`)
    enable_levels_meta : int (default: 0)
        Fit confidence levels for a criterion-based link function. If disabled, equidistant confidence levels are
        assumed (e.g. 1/3, 2/3, 1 for three criteria).

    *** Define fitting characteristics of the parameters ***
    * The fitting of each parameter is characzerized as follows:
    *     1) An initial guess.
    *     2) Lower and upper bound.
    *     3) Grid range, i.e. the range of values that are tested during the initial gridsearch search.
    * Sensible default values are provided for all parameters. To tweak those, one can either define an entire
    * ParameterSet, which is a container for a set of parameters, or each parameter individually. Note that the
    * parameters must be either defined as a Parameter instance or as List[Parameter] in case when separate values are
    * fitted for the positive and negative stimulus category/decision value).
    paramset_sens : ParameterSet
        Parameter set for the sensory level.
    paramset_meta : ParameterSet
        Parameter set for the metacognitive level.

    warping_sens : Union[Parameter, List[Parameter]]
        Parameter for the nonlinear transducer.
    noise_sens : Union[Parameter, List[Parameter]]
        Parameter for sensory noise.
    noise_transform_sens : Union[Parameter, List[Parameter]]
        Parameter for multiplicative sensory noise.
    thresh_sens : Union[Parameter, List[Parameter]]
        Parameter for the sensory threshold.
    bias_sens : Union[Parameter, List[Parameter]]
        Parameter for the sensory bias.
    noise_meta : Union[Parameter, List[Parameter]]
        Parameter for metacognitive noise.
    noise_transform_meta : Union[Parameter, List[Parameter]]
        Parameter for multiplicative megacognitive noise.
    readout_term_meta : Union[Parameter, List[Parameter]]
        Parameter for the metacognitive readout term.
    slope_meta : Union[Parameter, List[Parameter]]
        Parameter for the link function slope.
    confidence_term_meta : Union[Parameter, List[Parameter]]
        Parameter for the metacognitive report term.
    scaling_meta : Union[Parameter, List[Parameter]]
        Parameter for confidence scaling.
    criterion{x}_meta : Union[Parameter, List[Parameter]]
        Parameter for the xth confidence criterion (0-based).
    level{x}_meta : Union[Parameter, List[Parameter]]
        Parameter for the xth confidence level (0-based).

    *** Constraints ***
    * Constraints for the parameters are provided in the form of a function, as defined below.
    constraints_sens_callable : Callable
        Constraints for sensory parameters. See method `_constraints_sens_default` in this module for the expected
        format of this Callable. See method `_constraints_sens_default` in this module for the expected
        format of this Callable. The function must return a list of dictionaries, where each
        dictionary specifies a scipy optimize constraint.
    constraints_meta_callable : Callable
        Constraints for metacognitive parameters. The function must take the configuration object as the first argument
        and the stimulus array as the second argument. The function must return a list of dictionaries, where each
        dictionary specifies a scipy optimize constraint.

    *** Methodoligcal aspects of parameter fitting ***
    * Note: this applies to the fitting of metacognitive parameters only.
    gridsearch : bool (default: True)
        If True, perform initial (usually coarse) gridsearch search, based on the gridsearch defined for a Parameter.
    fine_gridsearch : bool (default: False)
        If True, perform an iteratively finer gridsearch search for each parameter.
    grid_multiproc : bool (default: False)
        If True, use all available cores for the gridsearch search. If False, use a single core.
    global_minimization : bool (default: False)
        If True, use a global minimization routine.
    gradient_free : bool (default: None)
        If True, use a gradien-free optimization routine.
    slsqp_epsilon : float (default: None)
        Set parameter epsilon parameter for the slsqp optimization method.

    *** Transformation functions ***
    function_warping_sens: str (default: 'power')
        Can be one of 'power', 'exponential' or 'identity'.
    function_noise_transform_sens: str (default: 'multiplicative')
        Can be one of 'multiplicative', 'power', 'exponential' or 'logarithm'.
    function_noise_transform_meta: str (default: 'multiplicative')
        Can be one of 'multiplicative', 'power', 'exponential' or 'logarithm'.

    *** Preprocessing ***
    normalize_stimuli_by_max : bool (default: True)
        If True, normalize provided stimuli by their maximum value.
    confidence_bounds_error : float
        Set ratings < confidence_bounds_error to 0 and > confidence_bounds_erroror to 1. This might be useful, if
        'finger errors' are assumed for values close to the scale's extremes.

    *** Parameters for the metacognitive likelihood computation ***
    max_dv_deviation : int
        Number of standard deviations around the mean considered for sensory uncertainty.
    nbins_dv : int
        Number of discrete decision values bins that are considered to represent sensory uncertainty.
    binsize_meta : float
        Integration bin size for the computation of the likelihood around empirical confidence values (noisy-report)
        or metacognitive evidence (noisy-readout).

    *** Other ***
    true_params : Dict
        Pass true (known) parameter values. This can be useful for testing to compare the likelihood of true and
        fitted parameters. The likelihood of true parameters is returned (and printed).
    force_settings : bool
        Some setting combinations are known to be incompatible and/or to produce biased fits. If True, fit the model
        nevertheless.
    settings_ignore_warnings : bool (default: False)
        If True, ignore warnings about user-specified settings.
    print_configuration : bool (default: True)
        If True, print the configuration at instatiation of the ReMeta class.
    """

    skip_meta: bool = False
    meta_noise_type: str = 'noisy_report'
    meta_noise_model: str = 'truncated_norm'
    meta_link_function: str = 'probability_correct'
    detection_model: bool = False
    detection_model_nchannels: int = 10

    enable_warping_sens: int = 0
    enable_noise_sens: int = 1
    enable_noise_transform_sens: int = 0
    enable_thresh_sens: int = 0
    enable_bias_sens: int = 1
    enable_noise_meta: int = 1
    enable_noise_transform_meta: int = 0
    enable_readout_term_meta: int = 0
    enable_slope_meta: int = 1
    enable_confidence_term_meta: int = 0
    enable_scaling_meta: int = 0
    enable_criteria_meta: int = 0
    enable_levels_meta: int = 0

    paramset_sens: ParameterSet = None
    paramset_meta: ParameterSet = None

    warping_sens: Union[Parameter, List[Parameter]] = None
    noise_sens: Union[Parameter, List[Parameter]] = None
    noise_transform_sens: Union[Parameter, List[Parameter]] = None
    thresh_sens: Union[Parameter, List[Parameter]] = None
    bias_sens: Union[Parameter, List[Parameter]] = None
    noise_meta: Union[Parameter, List[Parameter]] = None
    noise_transform_meta: Union[Parameter, List[Parameter]] = None
    readout_term_meta: Union[Parameter, List[Parameter]] = None
    slope_meta: Union[Parameter, List[Parameter]] = None
    confidence_term_meta: Union[Parameter, List[Parameter]] = None
    scaling_meta: Union[Parameter, List[Parameter]] = None

    constraints_sens_callable: Callable = None
    constraints_meta_callable: Callable = None

    function_warping_sens: str = 'power'
    function_noise_transform_sens: str = 'multiplicative'
    function_noise_transform_meta: str = 'multiplicative'

    gridsearch: bool = True
    fine_gridsearch: bool = False
    grid_multiproc: bool = False
    global_minimization: bool = False
    gradient_free: bool = None
    slsqp_epsilon: float = None

    normalize_stimuli_by_max: bool = True
    confidence_bounds_error: float = 0

    min_likelihood_meta: float = 1e-10
    min_likelihood_sens: float = 1e-10
    binsize_meta: float = 1e-1
    max_dv_deviation: int = 5
    nbins_dv: int = 101

    true_params: Dict = None
    force_settings: bool = False
    settings_ignore_warnings: bool = False
    print_configuration: bool = False

    noise_sens_default: float = 0.001
    noise_meta_default: float = 0.1
    noise_meta_min: float = 0.001

    _warping_sens_default: Parameter = Parameter(guess=0.1, bounds=(-10, 10), grid_range=np.arange(-10, 11, 5))
    _noise_transform_sens_default: Parameter = Parameter(guess=0, bounds=(0, 10), grid_range=np.arange(0, 1.1, 0.25))
    _noise_sens_default: Parameter = Parameter(guess=0.1, bounds=(1e-3, 100), grid_range=np.arange(0.1, 0.9, 0.25))
    _thresh_sens_default: Parameter = Parameter(guess=0, bounds=(0, 1), grid_range=np.arange(0, 0.41, 0.2))
    _bias_sens_default: Parameter = Parameter(guess=0, bounds=(-1, 1), grid_range=np.arange(-0.1, 0.11, 0.1))
    _noise_meta_default: Parameter = Parameter(guess=0.2, bounds=(1e-5, 50), grid_range=np.arange(0.05, 1, 0.1))
    _noise_transform_meta_default: Parameter = Parameter(guess=0, bounds=(0, 10),
                                                         grid_range=np.arange(0, 1.1, 0.25))
    _readout_term_meta_default: Parameter = Parameter(guess=0, bounds=(-1, 1), grid_range=np.arange(-0.2, 0.21, 0.1))
    _slope_meta_default: Parameter = Parameter(guess=1, bounds=(0.1, 50), grid_range=np.arange(0.2, 1.71, 0.3))
    _confidence_term_meta_default: Parameter = Parameter(guess=0, bounds=(-1, 1), grid_range=np.arange(-0.3, 0.31, 0.15))
    _scaling_meta_default: Parameter = Parameter(guess=1, bounds=(0.1, 10), grid_range=np.arange(0.5, 2.01, 0.3))
    _criterion_meta_default: Parameter = Parameter(guess=0, bounds=(1e-6, 50),
                                                   grid_range=np.exp(np.linspace(0, np.log(2), 8)) - 0.9)
    _level_meta_default: Parameter = Parameter(guess=0, bounds=(1e-6, 1),
                                               grid_range=np.exp(np.linspace(0, np.log(2), 8)) - 0.9)

    def __post_init__(self):

        if self.meta_link_function == 'probability_correct_ideal':
            self.enable_slope_meta = False

        if self.gradient_free is None:
            if '_criteria' in self.meta_link_function or '_transform' in self.meta_noise_model:
                self.gradient_free = True
            else:
                self.gradient_free = False

        if self.slsqp_epsilon is None:
            if self.meta_noise_type == 'noisy_readout':
                self.slsqp_epsilon = 1e-4
            else:
                self.slsqp_epsilon = _epsilon

        if self._scaling_meta_default is None and (self.meta_noise_type == 'noisy_readout'):
            self._scaling_meta_default = Parameter(guess=1, bounds=(0.1, 10), grid_range=np.arange(0.4, 1.01, 0.2))

        self._check_compatibility()

        self._prepare_params_sens()
        self._prepare_params_meta()
        if self.print_configuration:
            self.print()

    def _check_compatibility(self):
        if self.enable_noise_transform_meta:
            text = 'Fitting signal-dependent metacognitive noise parameters leads to biased estimates (for currently ' \
                   'unknown reasons) and is thus discouraged.'
            if self.force_settings:
                warnings.warn(text)
            else:
                raise ValueError(text)
        # if self.enable_noise_transform_meta == 2:
        #     warnings.warn('Fitting separate signal-dependent metacognitive noise parameters for the two stimulus '
        #                   'categories leads to biased estimates (for currently unknown reasons) and is thus '
        #                   'discouraged.')

        if self.detection_model and 'detection_model' not in self.meta_link_function:
            raise ValueError('Detection Models require a compatible metacogntive link function. Choose one of: '
                             'detection_model_linear, detection_model_mean, detection_model_mode, '
                             'detection_model_full, detection_model_ideal')

        if not self.enable_criteria_meta and '_criteria' in self.meta_link_function:
            raise ValueError('A criterion-based link function was set, but confidence criteria were not enabled.')

        if not self.settings_ignore_warnings:
            if self.enable_criteria_meta and '_criteria' not in self.meta_link_function:
                self.enable_criteria_meta = 0
                warnings.warn('Confidence criteria were enabled but no criterion-based link function was set -> '
                              'auto-setting enable_criteria_meta = 0.')
            if self.enable_levels_meta and '_criteria' not in self.meta_link_function:
                self.enable_levels_meta = 0
                warnings.warn('Confidence criteria were enabled but no criterion-based link function was set -> '
                              'auto-setting enable_levels_meta = 0.')
            if self.enable_levels_meta and not self.enable_criteria_meta:
                self.enable_criteria_meta = self.enable_levels_meta
                warnings.warn(f'Confidence leves were enabled, but confidence criteria not -> auto-setting'
                              f'enable_criteria_meta = {self.enable_levels_meta}')

            if self.meta_noise_type == 'noisy_readout':
                if self.enable_confidence_term_meta:
                    warnings.warn('The setting enable_confidence_term_meta has been enabled for a model of type '
                                  'noisy-readout. This setting results in a link function that is not monotonically '
                                  'increasing and likely leads to biased estimates.')
                if self.enable_scaling_meta:
                    warnings.warn('The setting enable_scaling_meta has been enabled for a model of type noisy-readout. '
                                  'Use this only if you have strong reasons to belief that values of the confidence '
                                  'scaling parameter are <= 1, since true values > 1 cannot be recovered for '
                                  'noisy-readout models.')

            if self.enable_criteria_meta and self.enable_slope_meta:
                self.enable_slope_meta = 0
                warnings.warn('Confidence criteria were enabled, which is in conflict with enable_slope_meta > 0 -> '
                              'auto-setting enable_slope_meta = 0.')

    def _prepare_params_sens(self):
        if self.paramset_sens is None:

            param_names_sens = []
            params_sens = ('warping', 'noise', 'noise_transform', 'thresh', 'bias')
            for param in params_sens:
                if getattr(self, f'enable_{param}_sens'):
                    param_names_sens += [f'{param}_sens']
                    if getattr(self, f'{param}_sens') is None:
                        param_default = getattr(self, f'_{param}_sens_default')
                        if getattr(self, f'enable_{param}_sens') == 2:
                            setattr(self, f'{param}_sens', [param_default, param_default])
                        else:
                            setattr(self, f'{param}_sens', param_default)

            parameters = {k: getattr(self, k) for k in param_names_sens}
            self.paramset_sens = ParameterSet(parameters, param_names_sens)

        if self.constraints_sens_callable is None:
            self.constraints_sens_callable = _constraints_sens_default

    def _prepare_params_meta(self):

        if self.paramset_meta is None:

            if self.enable_noise_meta and self.noise_meta is None:
                if self.meta_noise_model == 'beta':
                    self._noise_meta_default.bounds = (1e-5, 0.5)
                    self._noise_meta_default.grid_range = np.arange(0.05, 0.5, 0.05)
                elif self.meta_noise_type == 'noisy_readout':
                    self._noise_meta_default.bounds = (1e-5, 250)

            param_names_meta = []
            params_meta = ('noise', 'noise_transform', 'readout_term', 'slope', 'confidence_term', 'scaling')
            for param in params_meta:
                if getattr(self, f'enable_{param}_meta'):
                    param_names_meta += [f'{param}_meta']
                    if getattr(self, f'{param}_meta') is None:
                        param_default = getattr(self, f'_{param}_meta_default')
                        if getattr(self, f'enable_{param}_meta') == 2:
                            setattr(self, f'{param}_meta', [param_default.copy(), param_default.copy()])
                        else:
                            setattr(self, f'{param}_meta', param_default.copy())

            if self.enable_criteria_meta and '_criteria' in self.meta_link_function:  # noqa
                ncriteria_meta = int(self.meta_link_function.split('_')[0])
                guess_criteria = np.linspace(1 / (ncriteria_meta + 1), 1 - 1 / (ncriteria_meta + 1), ncriteria_meta)

                for i in range(ncriteria_meta):
                    if not hasattr(self, f'criterion{i}_meta') or getattr(self, f'criterion{i}_meta') is None:  # noqa
                        param_names_meta += [f'criterion{i}_meta']
                        param_default = self._criterion_meta_default.copy()
                        param_default.guess = guess_criteria[i]
                        if self.enable_criteria_meta == 2:
                            setattr(self, f'criterion{i}_meta', [param_default.copy(), param_default.copy()])
                        else:
                            setattr(self, f'criterion{i}_meta', param_default)
                    if self.enable_levels_meta and \
                            (not hasattr(self, f'level{i}_meta') or getattr(self, f'level{i}_meta') is None):
                        param_names_meta += [f'level{i}_meta']
                        param_default = self._level_meta_default.copy()
                        param_default.guess = guess_criteria[i]
                        if self.enable_levels_meta == 2:
                            setattr(self, f'level{i}_meta', [param_default.copy(), param_default.copy()])
                        else:
                            setattr(self, f'level{i}_meta', param_default)

            parameters = {k: getattr(self, k) for k in param_names_meta}
            self.paramset_meta = ParameterSet(parameters, param_names_meta)

        if self.constraints_meta_callable is None:
            self.constraints_meta_callable = _constraints_meta_default

    def print(self):
        # print('***********************')
        print(f'{self.__class__.__name__}')
        for k, v in self.__dict__.items():
            if not self.skip_meta or ('meta_' not in k and '_meta' not in k):
                print('\n'.join([f'\t{k}: {v}']))
        # print('***********************')

    def __repr__(self):
        txt = f'{self.__class__.__name__}\n'
        txt += '\n'.join([f'\t{k}: {v}' for k, v in self.__dict__.items()])
        return txt

def _constraints_sens_default(remeta):
    constraints = []
    if remeta.cfg.enable_noise_transform_sens:
        constraints += [{'type': 'ineq', 'fun': lambda x: np.min(remeta._negll_sens(x, return_noise=True)) - 0.001}]  # noqa
    return constraints

def _constraints_meta_default(remeta):
    constraints = []

    # Set constraints for the metacognitive noise slope
    # Note that metacognitive noise is based on confidence in case of the noisy-report model and
    # on metacognitive decision values in case of the noisy-readout model
    if remeta.cfg.enable_noise_transform_meta:
        constraints += [{'type': 'ineq', 'fun': lambda x: np.min(remeta.fun_meta(x, return_noise=True,
                                                                                constraint_mode=True)) - 0.001}]

    if '_criteria' in remeta.cfg.meta_link_function:
        ncrit_meta = int(remeta.cfg.meta_link_function.split('_')[0])
        # in case of a linear_tanh link function, the last parameter criterion corresponds to the slope of the tanh
        ncrit_meta_ = ncrit_meta - 1 if 'linear_tanh' in remeta.cfg.meta_link_function else ncrit_meta
        if remeta.cfg.enable_criteria_meta == 2:
            def fun_criteria(x):
                crit = remeta.fun_meta(x, return_criteria=True, constraint_mode=True)
                return np.sum([(crit[0][i] > (0 if i == 0 else crit[0][i - 1])) for i in range(ncrit_meta_)]) + \
                    np.sum([(crit[1][i] > (0 if i == 0 else crit[1][i - 1])) for i in range(ncrit_meta_)]) - \
                    2 * ncrit_meta_
        else:
            def fun_criteria(x):
                crit = remeta.fun_meta(x, return_criteria=True, constraint_mode=True)
                return np.sum([(crit[i] > (0 if i == 0 else crit[i - 1])) for i in range(ncrit_meta_)]) - ncrit_meta_
        constraints += [{'type': 'eq', 'fun': fun_criteria}]
        if remeta.cfg.enable_levels_meta:
            if remeta.cfg.enable_levels_meta == 2:
                def fun_levels(x):
                    lev = remeta.fun_meta(x, return_levels=True, constraint_mode=True)
                    return np.sum([(lev[0][i] > (0 if i == 0 else lev[0][i - 1])) for i in range(ncrit_meta)]) + \
                        np.sum([(lev[1][i] > (0 if i == 0 else lev[1][i - 1])) for i in range(ncrit_meta)]) - \
                        2 * ncrit_meta
            else:
                def fun_levels(x):
                    lev = remeta.fun_meta(x, return_levels=True, constraint_mode=True)
                    return np.sum([(lev[i] > (0 if i == 0 else lev[i - 1])) for i in range(ncrit_meta)]) - ncrit_meta
            constraints += [{'type': 'eq', 'fun': fun_levels}]
    return constraints

max_stars_repo_name: m-guggenmos/remeta

max_issues_repo_name: m-guggenmos/remeta

max_forks_repo_name: m-guggenmos/remeta

reason: import numpy,from scipy


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

# input parameter
den = 8880.0
cp = 386.0
cond = 398.0
temp_bc = 100.0
temp_init = 0.0
lx = 1.0
nx = 101
tend = 20000.0
dt = 0.1
tout = 100.0

alpha = cond / (den * cp)
dx = lx / (nx - 1)
nt = int(tend / dt)
nout = int(tout / dt)

#initial condition
temp = np.full(nx, temp_init)
time = 0.0

temp_new = np.zeros(nx)

# Boundary condition
temp[0] = temp_bc # Dirichlet @ x=0
temp[nx-1] = temp[nx-2] # Neumann @ x=Lx

# graph data array
ims = []
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)

gx = np.zeros(nx)
for i in range(nx):
    gx[i] = i * dx

# time loop
for n in range(1, nt+1):
    # FTCS
    for i in range(1, nx-1):
        temp_new[i] = temp[i] + dt * alpha * (temp[i+1] - 2.0 * temp[i] + temp[i-1]) / (dx * dx)

    # update
    for i in range(1, nx-1):
        temp[i] = temp_new[i]

    # Boundary condition
    temp[0] = temp_bc # Dirichlet @ x=0
    temp[nx-1] = temp[nx-2] # Neumann @ x=Lx

    time += dt

    if n % nout == 0:
        print('n: {0:7d}, time: {1:8.1f}, temp: {2:10.6f}'.format(n, time, temp[nx-1]))
        im_line = ax.plot(gx, temp, 'b')
        im_time = ax.text(0, 110, 'Time = {0:8.1f} [s]'.format(time))
        ims.append(im_line + [im_time])

# graph plot
ax.set_xlabel('x [m]')
ax.set_ylabel('Temperature [C]')
ax.grid()
anm = animation.ArtistAnimation(fig, ims, interval=50)
anm.save('animation.gif', writer='pillow')
plt.show()

ext: py

reason: import numpy


text:
import numpy as np
import tkinter as tk
from tkinter import ttk
from stellapy.GUI.graph_tools.display_information import display_information

#################################################################
#                  CLASS TO CHOOSE SIMULATIONS
#################################################################

class Simulations:

    def __init__(self, parent, frame):
        ''' This options frame controls which experiments/simulations are plotted. '''

        # Make the parents available
        self.tab = parent
        self.root = parent.root

        # Attributes
        self.experiment = None
        self.simulation = None
        self.experiment_id = "All experiments"
        self.simulation_id = "All simulations"
        self.options_kx = ["0.0"]
        self.options_ky = ["0.0", "100"]
        self.options_experiment = ["First plot data"]
        self.options_simulation = ["First plot data"]
        self.options_simulationsids = []

        # Variables
        self.var_experiment = tk.StringVar(value=self.options_experiment[0])
        self.var_simulation = tk.StringVar(value=self.options_simulation[0])
        self.var_kxmin = tk.StringVar(value=self.options_kx[0])
        self.var_kxmax = tk.StringVar(value=self.options_kx[-1])
        self.var_kymin = tk.StringVar(value=self.options_ky[0])
        self.var_kymax = tk.StringVar(value=self.options_ky[-1])

        # Configure the frame
        tk.Grid.columnconfigure(frame, 0, weight=1)

        # Create subframes
        self.frame1 = ttk.Frame(frame)
        self.frame2 = ttk.Frame(frame)
        self.frame1.grid(row=0, column=0, padx=(0,0), pady=(0,0), stick='NSEW')
        self.frame2.grid(row=2, column=0, padx=(0,0), pady=(0,0), stick='NSEW')

        # Choice between the experiments
        self.lbl_experiment = ttk.Label(self.frame1, text="Experiment: "); width=20
        self.mnu_experiment = ttk.OptionMenu(self.frame1, self.var_experiment, self.options_experiment[0], *self.options_experiment, style='option.TMenubutton')
        self.mnu_experiment["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_experiment.config(width=width)
        self.var_experiment.trace('w', self.change_plottedExperiment) # link function to a change of the dropdown options

        # Choice between the simulations
        self.lbl_simulation = ttk.Label(self.frame1, text="Simulation: ")
        self.mnu_simulation = ttk.OptionMenu(self.frame1, self.var_simulation, self.options_simulation[0], *self.options_simulation, style='option.TMenubutton')
        self.mnu_simulation["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_simulation.config(width=width)
        self.var_simulation.trace('w', self.change_plottedSimulation) # link function to a change of the dropdown options

        # Show the range of the kx and the k values
        self.lbl_kxmin = ttk.Label(self.frame2, text="kx min: ", style='opt_sign.TLabel'); width=4
        self.lbl_kxmax = ttk.Label(self.frame2, text=" kx max: ", style='opt_sign.TLabel')
        self.lbl_kymin = ttk.Label(self.frame2, text="ky min: ", style='opt_sign.TLabel')
        self.lbl_kymax = ttk.Label(self.frame2, text=" ky max: ", style='opt_sign.TLabel')
        self.mnu_kxmin = ttk.OptionMenu(self.frame2, self.var_kxmin, self.options_kx[0], *self.options_kx, style='option.TMenubutton')
        self.mnu_kxmin["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_kxmin.config(width=width)
        self.mnu_kxmax = ttk.OptionMenu(self.frame2, self.var_kxmax, self.options_kx[0], *self.options_kx, style='option.TMenubutton')
        self.mnu_kxmax["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_kxmax.config(width=width)
        self.mnu_kymin = ttk.OptionMenu(self.frame2, self.var_kymin, self.options_ky[0], *self.options_ky, style='option.TMenubutton')
        self.mnu_kymin["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_kymin.config(width=width)
        self.mnu_kymax = ttk.OptionMenu(self.frame2, self.var_kymax, self.options_ky[1], *self.options_ky, style='option.TMenubutton')
        self.mnu_kymax["menu"].config(bg=self.root.color['bbg'], fg=self.root.color['fg'], activebackground=self.root.color['bg'], activeforeground=self.root.color['fg'])
        self.mnu_kymax.config(width=width)
        self.var_kxmin.trace('w', self.change_selectedModes) # link function to a change of the dropdown options
        self.var_kxmax.trace('w', self.change_selectedModes) # link function to a change of the dropdown options
        self.var_kymin.trace('w', self.change_selectedModes) # link function to a change of the dropdown options
        self.var_kymax.trace('w', self.change_selectedModes) # link function to a change of the dropdown options

        # Configure the subframe
        tk.Grid.columnconfigure(self.frame1, 1, weight=1)
        tk.Grid.columnconfigure(self.frame2, 1, weight=1, uniform="2")
        tk.Grid.columnconfigure(self.frame2, 3, weight=1, uniform="2")

        # Place the widgets in the frame
        self.lbl_experiment.grid(row=0, column=0, stick='NSEW', padx=(0,0), pady=(2,2))
        self.mnu_experiment.grid(row=0, column=1, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=2, ipady=1)
        self.lbl_simulation.grid(row=1, column=0, stick='NSEW', padx=(0,0), pady=(2,2))
        self.mnu_simulation.grid(row=1, column=1, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=2, ipady=0)
        self.lbl_kxmin.grid(row=0, column=0, stick='W',    padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.mnu_kxmin.grid(row=0, column=1, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.lbl_kxmax.grid(row=0, column=2, stick='W',    padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.mnu_kxmax.grid(row=0, column=3, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.lbl_kymin.grid(row=1, column=0, stick='W',    padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.mnu_kymin.grid(row=1, column=1, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.lbl_kymax.grid(row=1, column=2, stick='W',    padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        self.mnu_kymax.grid(row=1, column=3, stick='NSEW', padx=(0,0), pady=(2,2), ipadx=1, ipady=0)
        if True: return

#################################################################
#                          METHODS
#################################################################

    def change_plottedExperiment(self, *args):

        # Change the experiment that is plotted
        self.experiment_id = self.var_experiment.get()

        # Get a reference to the experiment
        self.experiment = None
        for experiment in self.root.Research.experiments:
            if experiment.id == self.experiment_id:
                self.experiment = experiment
        if self.experiment==None and self.experiment_id!="All experiments":
            self.experiment = self.root.Research.experiments[0]
            self.experiment_id = self.experiment.id
        if self.experiment==None and self.experiment_id=="All experiments":
            self.experiment = self.root.Research.experiments[0]

        # Reset the graph classes
        self.tab.Graph[0].load_defaults()
        self.tab.Graph[1].load_defaults()
        self.tab.Graph[2].load_defaults()

    #------------------------------------
    def change_plottedSimulation(self, *args):

        # Change the simulation that is plotted
        self.simulation_id = self.var_simulation.get()
        self.simulation_id = self.options_simulationsids[self.options_simulations.index(self.simulation_id)]

        # Get a reference to the simulation
        self.simulation = None
        for simulation in self.experiment.simulations:
            if simulation.id == self.simulation_id:
                self.simulation = simulation
        if self.simulation==None and self.simulation_id!="All simulations":
            self.simulation = self.experiment.simulations[0]
            self.simulation_id = self.simulation.id
        if self.simulation==None and self.simulation_id=="All simulations":
            self.simulation = self.experiment.simulations[0]
        self.replot=True
        return

    #--------------------------------------
    def change_selectedModes(self, *args):
        ''' Tell the graph class which modes to plot.'''

        # Get the kx and k ranges set in the GUI
        kxmin = float(self.var_kxmin.get())
        kxmax = float(self.var_kxmax.get())
        kymin = float(self.var_kymin.get())
        kymax = float(self.var_kymax.get())

        # Make sure we don't scan over both kx and ky
        if kxmin != kxmax and kymin != kymax:
            message = """You can not scan over kx and ky simultaneously\n
            Please make sure one of the ranges is a fixed number."""
            display_information(self.root, "WARNING", message)

        # Tell the graph class which modes to plot.
        self.tab.Graph[0].kx = [ kxmin, kxmax ]
        self.tab.Graph[0].ky = [ kymin, kymax ]
        self.tab.Graph[1].kx = [ kxmin, kxmax ]
        self.tab.Graph[1].ky = [ kymin, kymax ]
        self.tab.Graph[2].kx = [ kxmin, kxmax ]
        self.tab.Graph[2].ky = [ kymin, kymax ]

        # Make sure one of the ranges is a float
        if kxmin == kxmax:
            self.tab.Graph[0].kx = kxmin
            self.tab.Graph[1].kx = kxmin
            self.tab.Graph[2].kx = kxmin
        elif kymin == kymax:
            self.tab.Graph[0].ky = kymin
            self.tab.Graph[1].ky = kymin
            self.tab.Graph[2].ky = kymin
        return

    #-------------------------------------------
    def display_plottedModesAndExperiments(self):
        ''' When the graph is plotted, adjust the menus on the GUI of the modes and experiments. '''

        # Get the experiment plotted by the function
        self.options_experiments = ["All experiments"] + [ e.id for e in self.root.Research.experiments ]
        if self.experiment==None:
            self.experiment = self.root.Research.experiments[0]
            self.var_experiment.set(self.options_experiments[0])
        if self.experiment.id not in self.options_experiments:
            self.experiment = self.root.Research.experiments[0]
            self.var_experiment.set(self.options_experiments[0])
        if self.experiment_id not in self.options_experiments:
            self.experiment = self.root.Research.experiments[0]
            self.var_experiment.set(self.options_experiments[0])
        if self.experiment_id == "All experiments":
            self.experiment = self.root.Research.experiments[0]
        for experiment in self.root.Research.experiments:
            if experiment.id == self.experiment_id:
                self.experiment = experiment
                self.var_experiment.set(self.experiment.id)

        # Reset the options in the menu
        self.mnu_experiment['menu'].delete(0, 'end')
        for experiment in self.options_experiments:
            self.mnu_experiment['menu'].add_command(label=experiment, command=tk._setit(self.var_experiment, experiment))

        # Get the simulation plotted by the function
        self.options_simulationsids = ["All simulations"] + [ s.id for s in self.experiment.simulations ]
        self.options_simulations = ["All simulations"] + [ s.id.split("__")[-1] for s in self.experiment.simulations ]
        self.options_simulations = ["All simulations"] + [ v for v in self.experiment.variedValues ]
        self.options_simulations = [ s.replace('\\', '').replace(',', '').replace('$', '') for s in self.options_simulations ]
        if self.simulation==None:
            self.simulation = self.experiment.simulations[0]
            self.var_simulation.set(self.options_simulations[0])
        if self.simulation.id not in self.options_simulationsids:
            self.simulation = self.experiment.simulations[0]
            self.var_simulation.set(self.options_simulations[0])
        if self.simulation_id not in self.options_simulationsids:
            self.simulation = self.experiment.simulations[0]
            self.var_simulation.set(self.options_simulations[0])
        for simulation in self.experiment.simulations:
            if simulation.id == self.simulation_id:
                self.simulation = simulation
                self.var_simulation.set(self.options_simulations[1+self.experiment.simulations.index(simulation)])

        # Reset the options in the menu
        self.mnu_simulation['menu'].delete(0, 'end')
        for simulation in self.options_simulations:
            self.mnu_simulation['menu'].add_command(label=simulation, command=tk._setit(self.var_simulation, simulation))

        # Update the options for the kx and ky ranges
        self.options_kx = self.experiment.total_kx
        self.options_ky = self.experiment.total_ky
        if self.experiment_id == "All experiments":
            for experiment in self.root.Research.experiments:
                self.options_kx = sorted(list(set(self.options_kx + experiment.total_kx)))
                self.options_ky = sorted(list(set(self.options_ky + experiment.total_ky)))
        self.mnu_kxmin['menu'].delete(0, 'end')
        self.mnu_kxmax['menu'].delete(0, 'end')
        self.mnu_kymin['menu'].delete(0, 'end')
        self.mnu_kymax['menu'].delete(0, 'end')
        for kx in self.options_kx:
            self.mnu_kxmin['menu'].add_command(label=round(kx,2), command=tk._setit(self.var_kxmin, kx))
            self.mnu_kxmax['menu'].add_command(label=round(kx,2), command=tk._setit(self.var_kxmax, kx))
        for ky in self.options_ky:
            self.mnu_kymin['menu'].add_command(label=round(ky,2), command=tk._setit(self.var_kymin, ky))
            self.mnu_kymax['menu'].add_command(label=round(ky,2), command=tk._setit(self.var_kymax, ky))

        # Update the options for the k ranges for the extraction frames!
        for Plot in [self.tab.PlotLinearMap, self.tab.PlotParameterInfluence]:
            mnu_kvalue = Plot.Extraction.mnu_kvalue
            options_k  = Plot.options_k
            var_kvalue = Plot.var_specificKvalue
            if isinstance(self.tab.Graph[0].kx, list): options_k  = self.experiment.total_ky
            if isinstance(self.tab.Graph[0].ky, list): options_k  = self.experiment.total_ky
            mnu_kvalue['menu'].delete(0, 'end')
            for k in options_k:
                mnu_kvalue['menu'].add_command(label=round(k,2), command=tk._setit(var_kvalue, k))

        # Update the GUI variables
        if isinstance(self.tab.Graph[0].kx, list):
            kx_min = round(np.max([self.tab.Graph[0].kx[0], self.options_kx[0]]),2)
            self.var_kxmin.set(kx_min)
            kx_max = round(np.min([self.tab.Graph[0].kx[1], self.options_kx[-1]]),2)
            self.var_kxmax.set(kx_max)
        if isinstance(self.tab.Graph[0].ky, list):
            ky_min = round(np.max([self.tab.Graph[0].ky[0], self.options_ky[0]]),2)
            self.var_kymin.set(ky_min)
            ky_max = round(np.min([self.tab.Graph[0].ky[1], self.options_ky[-1]]),2)
            self.var_kymax.set(ky_max)

        return


text:
import glob
import numpy as np
import matplotlib.pyplot as plt
import Imstats as img

psds = []

path = "../TestImgs/*.jpg"
for fname in glob.glob(path):
	psds.append(img.Imstat(fname).plot_power_spectrum())

print len(psds)

for i in range(0,len(psds)):
	#plt.plot(psds[i])
	plt.semilogy(psds[i])
plt.show()

max_stars_repo_name: giulioungaretti/ArtScope

max_issues_repo_name: giulioungaretti/ArtScope

max_forks_repo_name: giulioungaretti/ArtScope


text:
import argparse
import json
import numpy as np
import random
import itertools
import collections

from lm_eval import models, tasks

def parse_args():
    parser = argparse.ArgumentParser()
    parser.add_argument('--model', required=True)
    parser.add_argument('--model_args', default="")
    parser.add_argument('--tasks', default="all_tasks")
    parser.add_argument('--provide_description', action="store_true")
    parser.add_argument('--num_fewshot', type=int, default=1)
    parser.add_argument('--seed', type=int, default=1234)
    parser.add_argument('--output_path', default=None)
    parser.add_argument('--limit', type=int, default=None)
    return parser.parse_args()

def main():
    args = parse_args()
    random.seed(args.seed)
    np.random.seed(args.seed)

    lm = models.get_model(args.model).create_from_arg_string(args.model_args)
    if args.tasks == "all_tasks":
        task_names = tasks.ALL_TASKS
    else:
        task_names = args.tasks.split(",")
    task_dict = tasks.get_task_dict(task_names)

    # TODO: fall back to test docs
    task_dict_items = [(name, task) for name, task in task_dict.items() if task.has_validation_docs()]

    results = collections.defaultdict(dict)

    requests = collections.defaultdict(list)
    requests_origin = collections.defaultdict(list)

    # if we ever run into issues where the eval tasks don't fit in memory and we can't afford a machine with bigger memory,
    # we can always modify this plumbing to support that, but i didn't want to include it just yet because overengineering is bad
    # (or we could make it write the requests to disk and then read them back out again - probably using an sqlite db because of all the moving parts we have

    # TODO: we need unit tests & sanity checks or something to ensure that the return of `validation_docs` is stable

    docs = {}

    # get lists of each type of requeste
    for task_name, task in task_dict_items:
        for doc_id, doc in enumerate(itertools.islice(task.validation_docs(), 0, args.limit)):
            docs[(task_name, doc_id)] = doc

            ctx = task.fewshot_context(
                doc=doc,
                provide_description=args.provide_description,
                num_fewshot=args.num_fewshot,
            )

            reqs = task.construct_requests(doc, ctx)

            for i, req in enumerate(reqs):
                requests[req.type].append(req)
                # i: index in requests for a single task instance
                # doc_id: unique id that we can get back to a doc using `docs`
                requests_origin[req.type].append((i, task_name, doc, doc_id))

    # all responses for each (task, doc)
    process_res_queue = collections.defaultdict(list)

    # execute each type of request
    for reqtype, reqs in requests.items():
        # TODO: right now, this code runs multiple seperate LM requests for multiple Requests differing
        # only in index. We could implement some kind of caching, but that would be more of a bandaid
        # solution. we could also implement some kind of autogrouping here; they should end up next to each other.

        resps = getattr(lm, reqtype)([req.args for req in reqs])

        resps = [x if req.index is None else x[req.index] for x, req in zip(resps, reqs)]

        for resp, (i, task_name, doc, doc_id) in zip(resps, requests_origin[reqtype]):
            process_res_queue[(task_name, doc_id)].append((i, resp))

    vals = collections.defaultdict(list)

    # unpack results and sort back in order and return control to Task
    for (task_name, doc_id), requests in process_res_queue.items():
        requests.sort(key=lambda x: x[0])
        requests = [x[1] for x in requests]

        task = task_dict[task_name]
        doc = docs[(task_name, doc_id)]

        metrics = task.process_results(doc, requests)
        for metric, value in metrics.items():
            vals[(task_name, metric)].append(value)

    # aggregate results
    for (task_name, metric), items in vals.items():
        task = task_dict[task_name]
        results[task_name][metric] = task.aggregation()[metric](items)

    dumped = json.dumps(results, indent=2)
    print(dumped)
    if args.output_path:
        with open(args.output_path, "w") as f:
            f.write(dumped)

if __name__ == "__main__":
    main()

max_stars_repo_name: VitamintK/lm-evaluation-harness

max_issues_repo_name: VitamintK/lm-evaluation-harness

max_forks_repo_name: VitamintK/lm-evaluation-harness


text:
from enum import Enum
from typing import Union

import numpy as np
import pandas as pd

from experiment_analysis.exceptions import (
    InvalidInputDataframe,
    InvalidParameter,
)

class CappingMethod(Enum):
    STANDARD = "STANDARD"
    WITH_CUPED = "WITH_CUPED"

def apply_capping(
    df: pd.DataFrame,
    metric: str,
    sd_threshold: int = 5,
    method: Union[str, CappingMethod] = CappingMethod.STANDARD,
    metric_raw: str = None,
) -> pd.Series:
    """Calculates metrics capped at threshold of given SDs from mean (upper bound capping only). Intended for non-negative metrics.

    :param df: pd.Dataframe
    :param metric: str Name of dataframe column containing experiment metric to cap
    :param sd_threshold: sd_threshold Number of standard deviations from mean to set as outlier threshold. Default is 5.
    :param method: Union[str, CappingMethod] If "STANDARD" assumed non-CUPED metric. Use with "WITH_CUPED" if want to cap a CUPED metric.
    :param metric_raw: str Name of dataframe column containing metric before CUPED transformation if applying capping to a CUPED metric. Used in "WITH_CUPED" method.
    :return: Capped metric
    :rtype: pd.Series
    """
    user_method: CappingMethod

    if isinstance(method, str):
        user_method = CappingMethod(method.upper())
    else:
        user_method = CappingMethod(method)

    if metric not in df.columns:
        raise InvalidInputDataframe(
            f"input dataframe is missing metric column {metric!r}"
        )

    if not isinstance(sd_threshold, int):
        raise InvalidParameter(f"sd_threshold must be integer")

    if user_method == CappingMethod.STANDARD:
        metric_mean = df.loc[df[metric] > 0][metric].mean()
        metric_stdev = df.loc[df[metric] > 0][metric].std()
    elif user_method == CappingMethod.WITH_CUPED:

        if metric_raw not in df.columns:
            raise InvalidInputDataframe(
                f"input dataframe is missing raw metric column {metric_raw!r}"
            )

        metric_mean = df.loc[df[metric_raw] > 0][metric].mean()
        metric_stdev = df.loc[df[metric_raw] > 0][metric].std()

    upper_bound = metric_mean + sd_threshold * metric_stdev
    capped_metric = pd.Series(
        np.where(df[metric] > upper_bound, upper_bound, df[metric])
    )

    return capped_metric

max_stars_repo_name: Faire/experiment-analysis

max_issues_repo_name: Faire/experiment-analysis

max_forks_repo_name: Faire/experiment-analysis


text:
#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
# by the software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
# Research Corporation, et al.  All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
# license information.
#################################################################################
"""
Tests for solvent condenser unit model.
Authors: Andrew Lee
"""

import pytest
from pyomo.environ import (ConcreteModel,
                           Constraint,
                           Param,
                           TerminationCondition,
                           SolverStatus,
                           units,
                           value)
from pyomo.util.check_units import (assert_units_consistent,
                                    assert_units_equivalent)

from idaes.core import FlowsheetBlock
from idaes.generic_models.properties.core.generic.generic_property import (
        GenericParameterBlock)
from idaes.core.util.model_statistics import (degrees_of_freedom,
                                              number_variables,
                                              number_total_constraints,
                                              number_unused_variables)
from idaes.core.util.testing import initialization_tester
from idaes.core.util import get_solver, scaling as iscale

from idaes.generic_models.unit_models.column_models.solvent_condenser import (
    SolventCondenser)
from idaes.power_generation.carbon_capture.mea_solvent_system.properties.MEA_solvent \
    import configuration as aqueous_mea
from idaes.power_generation.carbon_capture.mea_solvent_system.properties.MEA_vapor \
    import wet_co2

# -----------------------------------------------------------------------------
# Get default solver for testing
solver = get_solver()

# -----------------------------------------------------------------------------
class TestStripperVaporFlow(object):
    @pytest.fixture(scope="class")
    def model(self):
        m = ConcreteModel()
        m.fs = FlowsheetBlock(default={"dynamic": False})

        m.fs.liquid_properties = GenericParameterBlock(default=aqueous_mea)
        m.fs.vapor_properties = GenericParameterBlock(default=wet_co2)

        m.fs.unit = SolventCondenser(default={
            "liquid_property_package": m.fs.liquid_properties,
            "vapor_property_package": m.fs.vapor_properties})

        m.fs.unit.inlet.flow_mol[0].fix(1.1117)
        m.fs.unit.inlet.temperature[0].fix(339.33)
        m.fs.unit.inlet.pressure[0].fix(184360)
        m.fs.unit.inlet.mole_frac_comp[0, "CO2"].fix(0.8817)
        m.fs.unit.inlet.mole_frac_comp[0, "H2O"].fix(0.1183)

        m.fs.unit.reflux.flow_mol[0].fix(0.1083)

        iscale.set_scaling_factor(
            m.fs.unit.vapor_phase.properties_out[0].fug_phase_comp[
                "Vap", "CO2"], 1e-5)
        iscale.set_scaling_factor(
            m.fs.unit.vapor_phase.properties_out[0].fug_phase_comp[
                "Vap", "H2O"], 1e-3)

        iscale.calculate_scaling_factors(m.fs.unit)

        return m

    @pytest.mark.build
    @pytest.mark.unit
    def test_build(self, model):

        assert hasattr(model.fs.unit, "inlet")
        assert len(model.fs.unit.inlet.vars) == 4
        assert hasattr(model.fs.unit.inlet, "flow_mol")
        assert hasattr(model.fs.unit.inlet, "mole_frac_comp")
        assert hasattr(model.fs.unit.inlet, "temperature")
        assert hasattr(model.fs.unit.inlet, "pressure")

        assert hasattr(model.fs.unit, "reflux")
        assert len(model.fs.unit.reflux.vars) == 4
        assert hasattr(model.fs.unit.reflux, "flow_mol")
        assert hasattr(model.fs.unit.reflux, "mole_frac_comp")
        assert hasattr(model.fs.unit.reflux, "temperature")
        assert hasattr(model.fs.unit.reflux, "pressure")

        assert hasattr(model.fs.unit, "vapor_outlet")
        assert len(model.fs.unit.vapor_outlet.vars) == 4
        assert hasattr(model.fs.unit.vapor_outlet, "flow_mol")
        assert hasattr(model.fs.unit.vapor_outlet, "mole_frac_comp")
        assert hasattr(model.fs.unit.vapor_outlet, "temperature")
        assert hasattr(model.fs.unit.vapor_outlet, "pressure")

        assert isinstance(model.fs.unit.unit_material_balance, Constraint)
        assert isinstance(model.fs.unit.unit_enthalpy_balance, Constraint)
        assert isinstance(model.fs.unit.unit_temperature_equality, Constraint)
        assert isinstance(model.fs.unit.unit_pressure_balance, Constraint)
        assert isinstance(model.fs.unit.zero_flow_param, Param)

        assert number_variables(model.fs.unit) == 55
        assert number_total_constraints(model.fs.unit) == 49
        assert number_unused_variables(model.fs.unit) == 0

    @pytest.mark.component
    def test_units(self, model):
        assert_units_consistent(model)
        assert_units_equivalent(model.fs.unit.heat_duty[0], units.W)

    @pytest.mark.unit
    def test_dof(self, model):
        assert degrees_of_freedom(model) == 0

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_initialize(self, model):
        initialization_tester(model)

    # @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, model):
        results = solver.solve(model)

        # Check for optimal solution
        assert results.solver.termination_condition == \
            TerminationCondition.optimal
        assert results.solver.status == SolverStatus.ok

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, model):
        assert (pytest.approx(0.1083, rel=1e-5) ==
                value(model.fs.unit.reflux.flow_mol[0]))
        assert (pytest.approx(0, abs=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'CO2']))
        assert (pytest.approx(0, abs=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'MEA']))
        assert (pytest.approx(1, rel=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'H2O']))
        assert (pytest.approx(184360, rel=1e-5) ==
                value(model.fs.unit.reflux.pressure[0]))
        assert (pytest.approx(303.244, rel=1e-5) ==
                value(model.fs.unit.reflux.temperature[0]))

        assert (pytest.approx(1.0034, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.flow_mol[0]))
        assert (pytest.approx(0.976758, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.mole_frac_comp[0, 'CO2']))
        assert (pytest.approx(0.0232423, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.mole_frac_comp[0, 'H2O']))
        assert (pytest.approx(184360, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.pressure[0]))
        assert (pytest.approx(303.244, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.temperature[0]))

        assert (pytest.approx(-6264.72, rel=1e-5) ==
                value(model.fs.unit.heat_duty[0]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, model):
        assert abs(value(model.fs.unit.inlet.flow_mol[0] -
                         model.fs.unit.reflux.flow_mol[0] -
                         model.fs.unit.vapor_outlet.flow_mol[0])) <= 1e-6

        assert (abs(value(model.fs.unit.inlet.flow_mol[0] *
                          model.fs.unit.inlet.mole_frac_comp[0, "CO2"] -
                          model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "CO2"] -
                          model.fs.unit.vapor_outlet.flow_mol[0] *
                          model.fs.unit.vapor_outlet.mole_frac_comp[0, "CO2"]))
                <= 1e-6)
        assert (abs(value(model.fs.unit.inlet.flow_mol[0] *
                          model.fs.unit.inlet.mole_frac_comp[0, "H2O"] -
                          model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "H2O"] -
                          model.fs.unit.vapor_outlet.flow_mol[0] *
                          model.fs.unit.vapor_outlet.mole_frac_comp[0, "H2O"]))
                <= 1e-6)
        assert (abs(value(model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "MEA"]))
                <= 1e-6)

        assert abs(value(
            model.fs.unit.vapor_phase.properties_in[0]._enthalpy_flow_term[
                "Vap"] -
            model.fs.unit.vapor_phase.properties_out[0]._enthalpy_flow_term[
                "Vap"] -
            model.fs.unit.liquid_phase[0]._enthalpy_flow_term["Liq"] +
            model.fs.unit.heat_duty[0])) <= 1e-6

# -----------------------------------------------------------------------------
class TestStripperHeatDuty(object):
    @pytest.fixture(scope="class")
    def model(self):
        m = ConcreteModel()
        m.fs = FlowsheetBlock(default={"dynamic": False})

        m.fs.liquid_properties = GenericParameterBlock(default=aqueous_mea)
        m.fs.vapor_properties = GenericParameterBlock(default=wet_co2)

        m.fs.unit = SolventCondenser(default={
            "liquid_property_package": m.fs.liquid_properties,
            "vapor_property_package": m.fs.vapor_properties})

        m.fs.unit.inlet.flow_mol[0].fix(1.1117)
        m.fs.unit.inlet.temperature[0].fix(339.33)
        m.fs.unit.inlet.pressure[0].fix(184360)
        m.fs.unit.inlet.mole_frac_comp[0, "CO2"].fix(0.8817)
        m.fs.unit.inlet.mole_frac_comp[0, "H2O"].fix(0.1183)

        m.fs.unit.heat_duty.fix(-6264)

        return m

    @pytest.mark.build
    @pytest.mark.unit
    def test_build(self, model):

        assert hasattr(model.fs.unit, "inlet")
        assert len(model.fs.unit.inlet.vars) == 4
        assert hasattr(model.fs.unit.inlet, "flow_mol")
        assert hasattr(model.fs.unit.inlet, "mole_frac_comp")
        assert hasattr(model.fs.unit.inlet, "temperature")
        assert hasattr(model.fs.unit.inlet, "pressure")

        assert hasattr(model.fs.unit, "reflux")
        assert len(model.fs.unit.reflux.vars) == 4
        assert hasattr(model.fs.unit.reflux, "flow_mol")
        assert hasattr(model.fs.unit.reflux, "mole_frac_comp")
        assert hasattr(model.fs.unit.reflux, "temperature")
        assert hasattr(model.fs.unit.reflux, "pressure")

        assert hasattr(model.fs.unit, "vapor_outlet")
        assert len(model.fs.unit.vapor_outlet.vars) == 4
        assert hasattr(model.fs.unit.vapor_outlet, "flow_mol")
        assert hasattr(model.fs.unit.vapor_outlet, "mole_frac_comp")
        assert hasattr(model.fs.unit.vapor_outlet, "temperature")
        assert hasattr(model.fs.unit.vapor_outlet, "pressure")

        assert isinstance(model.fs.unit.unit_material_balance, Constraint)
        assert isinstance(model.fs.unit.unit_enthalpy_balance, Constraint)
        assert isinstance(model.fs.unit.unit_temperature_equality, Constraint)
        assert isinstance(model.fs.unit.unit_pressure_balance, Constraint)
        assert isinstance(model.fs.unit.zero_flow_param, Param)

        assert number_variables(model.fs.unit) == 55
        assert number_total_constraints(model.fs.unit) == 49
        assert number_unused_variables(model.fs.unit) == 0

    @pytest.mark.component
    def test_units(self, model):
        assert_units_consistent(model)
        assert_units_equivalent(model.fs.unit.heat_duty[0], units.W)

    @pytest.mark.unit
    def test_dof(self, model):
        assert degrees_of_freedom(model) == 0

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_initialize(self, model):
        initialization_tester(model)

    # @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solve(self, model):
        results = solver.solve(model)

        # Check for optimal solution
        assert results.solver.termination_condition == \
            TerminationCondition.optimal
        assert results.solver.status == SolverStatus.ok

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_solution(self, model):
        assert (pytest.approx(0.108291, rel=1e-5) ==
                value(model.fs.unit.reflux.flow_mol[0]))
        assert (pytest.approx(0, abs=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'CO2']))
        assert (pytest.approx(0, abs=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'MEA']))
        assert (pytest.approx(1, rel=1e-3) ==
                value(model.fs.unit.reflux.mole_frac_comp[0, 'H2O']))
        assert (pytest.approx(184360, rel=1e-5) ==
                value(model.fs.unit.reflux.pressure[0]))
        assert (pytest.approx(303.250, rel=1e-5) ==
                value(model.fs.unit.reflux.temperature[0]))

        assert (pytest.approx(1.0034, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.flow_mol[0]))
        assert (pytest.approx(0.976758, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.mole_frac_comp[0, 'CO2']))
        assert (pytest.approx(0.0232509, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.mole_frac_comp[0, 'H2O']))
        assert (pytest.approx(184360, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.pressure[0]))
        assert (pytest.approx(303.250, rel=1e-5) ==
                value(model.fs.unit.vapor_outlet.temperature[0]))

        assert (pytest.approx(-6264, rel=1e-5) ==
                value(model.fs.unit.heat_duty[0]))

    @pytest.mark.solver
    @pytest.mark.skipif(solver is None, reason="Solver not available")
    @pytest.mark.component
    def test_conservation(self, model):
        assert abs(value(model.fs.unit.inlet.flow_mol[0] -
                         model.fs.unit.reflux.flow_mol[0] -
                         model.fs.unit.vapor_outlet.flow_mol[0])) <= 1e-6

        assert (abs(value(model.fs.unit.inlet.flow_mol[0] *
                          model.fs.unit.inlet.mole_frac_comp[0, "CO2"] -
                          model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "CO2"] -
                          model.fs.unit.vapor_outlet.flow_mol[0] *
                          model.fs.unit.vapor_outlet.mole_frac_comp[0, "CO2"]))
                <= 1e-6)
        assert (abs(value(model.fs.unit.inlet.flow_mol[0] *
                          model.fs.unit.inlet.mole_frac_comp[0, "H2O"] -
                          model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "H2O"] -
                          model.fs.unit.vapor_outlet.flow_mol[0] *
                          model.fs.unit.vapor_outlet.mole_frac_comp[0, "H2O"]))
                <= 1e-6)
        assert (abs(value(model.fs.unit.reflux.flow_mol[0] *
                          model.fs.unit.reflux.mole_frac_comp[0, "MEA"]))
                <= 1e-6)

        assert abs(value(
            model.fs.unit.vapor_phase.properties_in[0]._enthalpy_flow_term[
                "Vap"] -
            model.fs.unit.vapor_phase.properties_out[0]._enthalpy_flow_term[
                "Vap"] -
            model.fs.unit.liquid_phase[0]._enthalpy_flow_term["Liq"] +
            model.fs.unit.heat_duty[0])) <= 1e-6

    @pytest.mark.component
    def test_scaling(self, model):
        iscale.set_scaling_factor(
            model.fs.unit.vapor_phase.properties_out[0].fug_phase_comp[
                "Vap", "CO2"], 1e-5)
        iscale.set_scaling_factor(
            model.fs.unit.vapor_phase.properties_out[0].fug_phase_comp[
                "Vap", "H2O"], 1e-3)

        iscale.calculate_scaling_factors(model.fs.unit)

        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_material_balance[0, "CO2"]) == 1
        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_material_balance[0, "H2O"]) == 1
        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_material_balance[0, "MEA"]) == 1e8

        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_phase_equilibrium[0, "CO2"]) == 1e-5
        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_phase_equilibrium[0, "H2O"]) == 1e-3

        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_temperature_equality[0]) == 1e-2

        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_enthalpy_balance[0]) == 1

        assert iscale.get_constraint_transform_applied_scaling_factor(
            model.fs.unit.unit_pressure_balance[0]) == 1e-5

max_stars_repo_name: michaelbynum/idaes-pse

max_issues_repo_name: michaelbynum/idaes-pse

max_forks_repo_name: michaelbynum/idaes-pse

reason: from pyomo


text:
'''
from David Huard's scipy sandbox, also attached to a ticket and
in the matplotlib-user mailinglist  (links ???)

Notes
=====

out of bounds interpolation raises exception and would not be completely
defined ::

>>> scoreatpercentile(x, [0,25,50,100])
Traceback (most recent call last):
...
    raise ValueError("A value in x_new is below the interpolation "
ValueError: A value in x_new is below the interpolation range.
>>> percentileofscore(x, [-50, 50])
Traceback (most recent call last):
...
    raise ValueError("A value in x_new is below the interpolation "
ValueError: A value in x_new is below the interpolation range.

idea
====

histogram and empirical interpolated distribution
-------------------------------------------------

dual constructor
* empirical cdf : cdf on all observations through linear interpolation
* binned cdf : based on histogram
both should work essentially the same, although pdf of empirical has
many spikes, fluctuates a lot
- alternative: binning based on interpolated cdf : example in script
* ppf: quantileatscore based on interpolated cdf
* rvs : generic from ppf
* stats, expectation ? how does integration wrt cdf work - theory?

Problems
* limits, lower and upper bound of support
  does not work or is undefined with empirical cdf and interpolation
* extending bounds ?
  matlab has pareto tails for empirical distribution, breaks linearity

empirical distribution with higher order interpolation
------------------------------------------------------

* should work easily enough with interpolating splines
* not piecewise linear
* can use pareto (or other) tails
* ppf how do I get the inverse function of a higher order spline?
  Chuck: resample and fit spline to inverse function
  this will have an approximation error in the inverse function
* -> does not work: higher order spline does not preserve monotonicity
  see mailing list for response to my question
* pmf from derivative available in spline

-> forget this and use kernel density estimator instead

bootstrap/empirical distribution:
---------------------------------

discrete distribution on real line given observations
what's defined?
* cdf : step function
* pmf : points with equal weight 1/nobs
* rvs : resampling
* ppf : quantileatscore on sample?
* moments : from data ?
* expectation ? sum_{all observations x} [func(x) * pmf(x)]
* similar for discrete distribution on real line
* References : ?
* what's the point? most of it is trivial, just for the record ?

Created on Monday, May 03, 2010, 11:47:03 AM
Author: josef-pktd, parts based on David Huard
License: BSD

'''
import scipy.interpolate as interpolate
import numpy as np

def scoreatpercentile(data, percentile):
    """Return the score at the given percentile of the data.

    Example:
        >>> data = randn(100)
            >>> scoreatpercentile(data, 50)

        will return the median of sample `data`.
    """
    per = np.array(percentile)
    cdf = empiricalcdf(data)
    interpolator = interpolate.interp1d(np.sort(cdf), np.sort(data))
    return interpolator(per/100.)

def percentileofscore(data, score):
    """Return the percentile-position of score relative to data.

    score: Array of scores at which the percentile is computed.

    Return percentiles (0-100).

    Example
            r = randn(50)
        x = linspace(-2,2,100)
        percentileofscore(r,x)

    Raise an error if the score is outside the range of data.
    """
    cdf = empiricalcdf(data)
    interpolator = interpolate.interp1d(np.sort(data), np.sort(cdf))
    return interpolator(score)*100.

def empiricalcdf(data, method='Hazen'):
    """Return the empirical cdf.

    Methods available:
        Hazen:       (i-0.5)/N
            Weibull:     i/(N+1)
        Chegodayev:  (i-.3)/(N+.4)
        Cunnane:     (i-.4)/(N+.2)
        Gringorten:  (i-.44)/(N+.12)
        California:  (i-1)/N

    Where i goes from 1 to N.
    """

    i = np.argsort(np.argsort(data)) + 1.
    N = len(data)
    method = method.lower()
    if method == 'hazen':
        cdf = (i-0.5)/N
    elif method == 'weibull':
        cdf = i/(N+1.)
    elif method == 'california':
        cdf = (i-1.)/N
    elif method == 'chegodayev':
        cdf = (i-.3)/(N+.4)
    elif method == 'cunnane':
        cdf = (i-.4)/(N+.2)
    elif method == 'gringorten':
        cdf = (i-.44)/(N+.12)
    else:
        raise ValueError('Unknown method. Choose among Weibull, Hazen,'
                         'Chegodayev, Cunnane, Gringorten and California.')

    return cdf

class HistDist:
    '''Distribution with piecewise linear cdf, pdf is step function

    can be created from empiricial distribution or from a histogram (not done yet)

    work in progress, not finished

    '''

    def __init__(self, data):
        self.data = np.atleast_1d(data)
        self.binlimit = np.array([self.data.min(), self.data.max()])
        sortind = np.argsort(data)
        self._datasorted = data[sortind]
        self.ranking = np.argsort(sortind)

        cdf = self.empiricalcdf()
        self._empcdfsorted = np.sort(cdf)
        self.cdfintp = interpolate.interp1d(self._datasorted, self._empcdfsorted)
        self.ppfintp = interpolate.interp1d(self._empcdfsorted, self._datasorted)

    def empiricalcdf(self, data=None, method='Hazen'):
        """Return the empirical cdf.

        Methods available:
            Hazen:       (i-0.5)/N
                Weibull:     i/(N+1)
            Chegodayev:  (i-.3)/(N+.4)
            Cunnane:     (i-.4)/(N+.2)
            Gringorten:  (i-.44)/(N+.12)
            California:  (i-1)/N

        Where i goes from 1 to N.
        """

        if data is None:
            data = self.data
            i = self.ranking
        else:
            i = np.argsort(np.argsort(data)) + 1.

        N = len(data)
        method = method.lower()
        if method == 'hazen':
            cdf = (i-0.5)/N
        elif method == 'weibull':
            cdf = i/(N+1.)
        elif method == 'california':
            cdf = (i-1.)/N
        elif method == 'chegodayev':
            cdf = (i-.3)/(N+.4)
        elif method == 'cunnane':
            cdf = (i-.4)/(N+.2)
        elif method == 'gringorten':
            cdf = (i-.44)/(N+.12)
        else:
            raise ValueError('Unknown method. Choose among Weibull, Hazen,'
                             'Chegodayev, Cunnane, Gringorten and California.')

        return cdf

    def cdf_emp(self, score):
        '''
        this is score in dh

        '''
        return self.cdfintp(score)
        #return percentileofscore(self.data, score)

    def ppf_emp(self, quantile):
        '''
        this is score in dh

        '''
        return self.ppfintp(quantile)
        #return scoreatpercentile(self.data, quantile*100)

    #from DHuard http://old.nabble.com/matplotlib-f2903.html
    def optimize_binning(self, method='Freedman'):
        """Find the optimal number of bins and update the bin countaccordingly.
        Available methods : Freedman
                            Scott
        """

        nobs = len(self.data)
        if method=='Freedman':
            IQR = self.ppf_emp(0.75) - self.ppf_emp(0.25) # Interquantile range(75% -25%)
            width = 2* IQR* nobs**(-1./3)

        elif method=='Scott':
            width = 3.49 * np.std(self.data) * nobs**(-1./3)

        self.nbin = (np.ptp(self.binlimit)/width)
        return self.nbin

#changes: josef-pktd
if __name__ == '__main__':
    import matplotlib.pyplot as plt

    nobs = 100
    x = np.random.randn(nobs)

    examples = [2]
    if 1 in examples:
        empiricalcdf(x)
        print(percentileofscore(x, 0.5))
        print(scoreatpercentile(x, 50))
        xsupp = np.linspace(x.min(), x.max())
        pos = percentileofscore(x, xsupp)
        plt.plot(xsupp, pos)
        #perc = np.linspace(2.5, 97.5)
        #plt.plot(scoreatpercentile(x, perc), perc)
        plt.plot(scoreatpercentile(x, pos), pos+1)

        #emp = interpolate.PiecewisePolynomial(np.sort(empiricalcdf(x)), np.sort(x))
        emp=interpolate.InterpolatedUnivariateSpline(np.sort(x),np.sort(empiricalcdf(x)),k=1)
        pdfemp = np.array([emp.derivatives(xi)[1] for xi in xsupp])
        plt.figure()
        plt.plot(xsupp,pdfemp)
        cdf_ongrid = emp(xsupp)
        plt.figure()
        plt.plot(xsupp, cdf_ongrid)

        #get pdf from interpolated cdf on a regular grid
        plt.figure()
        plt.step(xsupp[:-1],np.diff(cdf_ongrid)/np.diff(xsupp))

        #reduce number of bins/steps
        xsupp2 = np.linspace(x.min(), x.max(), 25)
        plt.figure()
        plt.step(xsupp2[:-1],np.diff(emp(xsupp2))/np.diff(xsupp2))

        #pdf using 25 original observations, every (nobs/25)th
        xso = np.sort(x)
        xs = xso[::nobs/25]
        plt.figure()
        plt.step(xs[:-1],np.diff(emp(xs))/np.diff(xs))
        #lower end looks strange

    histd = HistDist(x)
    print(histd.optimize_binning())
    print(histd.cdf_emp(histd.binlimit))
    print(histd.ppf_emp([0.25, 0.5, 0.75]))
    print(histd.cdf_emp([-0.5, -0.25, 0, 0.25, 0.5]))

    xsupp = np.linspace(x.min(), x.max(), 500)
    emp=interpolate.InterpolatedUnivariateSpline(np.sort(x),np.sort(empiricalcdf(x)),k=1)
    #pdfemp = np.array([emp.derivatives(xi)[1] for xi in xsupp])
    #plt.figure()
    #plt.plot(xsupp,pdfemp)
    cdf_ongrid = emp(xsupp)
    plt.figure()
    plt.plot(xsupp, cdf_ongrid)
    ppfintp = interpolate.InterpolatedUnivariateSpline(cdf_ongrid,xsupp,k=3)

    ppfs = ppfintp(cdf_ongrid)
    plt.plot(ppfs, cdf_ongrid)
    #ppfemp=interpolate.InterpolatedUnivariateSpline(np.sort(empiricalcdf(x)),np.sort(x),k=3)
    #Do not use interpolating splines for function approximation
    #with s=0.03 the spline is monotonic at the evaluated values
    ppfemp=interpolate.UnivariateSpline(np.sort(empiricalcdf(x)),np.sort(x),k=3, s=0.03)
    ppfe = ppfemp(cdf_ongrid)
    plt.plot(ppfe, cdf_ongrid)

    print('negative density')
    print('(np.diff(ppfs)).min()', (np.diff(ppfs)).min())
    print('(np.diff(cdf_ongrid)).min()', (np.diff(cdf_ongrid)).min())
    #plt.show()

max_stars_repo_name: aliavni/statsmodels

max_issues_repo_name: aliavni/statsmodels

max_forks_repo_name: aliavni/statsmodels

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
"""
    boxel.config
    ~~~~~~~~~~~~~

    This module provides configurations needed for boxelization. Such as yaml
    palette parsing and palette generation and also dimension calculations.
"""

import math
import os
import yaml
import numpy as np
from PIL import Image
from utils.singleton import singleton

@singleton
class Dimensions:
    """Dimensions calculated upon init of object including the number of pixels
        per boxel, the width in boxels of the image and the dims in boxels

    :param pixel_width:
    """
    def __init__(self, **kwargs):
        self.pixel_width = kwargs['pixel_width']

        # determine boxel size -- rows first
        self.pixels_per_box = math.ceil(
                kwargs['img_size'][0]/float(self.pixel_width))
        self.output_dims = [
                self.pixel_width,
                int(math.ceil(kwargs['img_size'][1]/self.pixels_per_box))]

@singleton
class Palette:
    """A palette generated by a YAML file containing containg an array of RGB
        values to be used as the only colors for boxels. Contains the name of
        the palette, an array of colors, the amount of colors and a image with
        only those colors in it for easy conversion of other images to this
        palette.

    :param file: A filepath to palette.yml these should be in the palette
        directory
    """
    def __init__(self, **kwargs):
        self.name, self.colors, self.palette_img = self._load(kwargs['file'])
        self.num_of_colors = len(self.colors)

    def _load(self, filepath):
        """Loads the YAML file with an array of RGB values and generates an
            array of RGB colors, gets the name of the palette and makes an
            Image with the colors in palette.yml

        :param filepath: A filepath to palette.yml these should be in the
            palette directory
        """
        if filepath:
            with open(filepath, 'r') as stream:
                base = os.path.basename(filepath)
                filename = os.path.splitext(base)[0]
                dct = yaml.load(stream)
                palette = np.asarray(dct[filename], dtype=np.uint8)
                fill = (256*3) - (palette.shape[0] * palette.shape[1])
                palette_img = Image.new('P', (50, 50))
                palette_img.putpalette(
                        # We add zeros since it's expecting 256*3 color values
                        np.pad(
                            np.asarray(palette).flatten(), (0, fill),
                            'constant'))
                return filename, palette, palette_img
        else:
            return None, np.asarray([], dtype=np.uint8)

max_stars_repo_name: strugee/Boxel


text:
import numpy as np
import tensorflow as tf
from utils.sth import sth
from typing import Dict

def make_off_policy_class(mode='share'):
    if mode == 'share':
        from algos.tf2algos.base.share_rnn_cnn_policy import SharedPolicy as BasePolicy
    else:
        from algos.tf2algos.base.policy import Policy as BasePolicy

    class Off_Policy(BasePolicy):
        def __init__(self,
                     s_dim,
                     visual_sources,
                     visual_resolution,
                     a_dim,
                     is_continuous,
                     **kwargs):
            super().__init__(
                s_dim=s_dim,
                visual_sources=visual_sources,
                visual_resolution=visual_resolution,
                a_dim=a_dim,
                is_continuous=is_continuous,
                **kwargs)
            self.buffer_size = int(kwargs.get('buffer_size', 10000))
            self.use_priority = kwargs.get('use_priority', False)
            self.n_step = kwargs.get('n_step', False)
            self.use_isw = bool(kwargs.get('use_isw', False))
            self.train_times_per_step = int(kwargs.get('train_times_per_step', 1))

        def set_buffer(self, buffer):
            '''
            TODO: Annotation
            '''
            self.data = buffer

        def store_data(self, s, visual_s, a, r, s_, visual_s_, done):
            """
            for off-policy training, use this function to store <s, a, r, s_, done> into ReplayBuffer.
            """
            assert isinstance(a, np.ndarray), "store need action type is np.ndarray"
            assert isinstance(r, np.ndarray), "store need reward type is np.ndarray"
            assert isinstance(done, np.ndarray), "store need done type is np.ndarray"
            self._running_average(s)
            self.data.add(
                s,
                visual_s,
                a,
                r[:, np.newaxis],   # 升维
                s_,
                visual_s_,
                done[:, np.newaxis]  # 升维
            )

        def no_op_store(self, s, visual_s, a, r, s_, visual_s_, done):
            assert isinstance(a, np.ndarray), "no_op_store need action type is np.ndarray"
            assert isinstance(r, np.ndarray), "no_op_store need reward type is np.ndarray"
            assert isinstance(done, np.ndarray), "no_op_store need done type is np.ndarray"
            self._running_average(s)
            self.data.add(
                s,
                visual_s,
                a,
                r[:, np.newaxis],
                s_,
                visual_s_,
                done[:, np.newaxis]
            )

        def get_transitions(self, data_name_list=['s', 'visual_s', 'a', 'r', 's_', 'visual_s_', 'done']):
            '''
            TODO: Annotation
            '''
            data = self.data.sample()   # 经验池取数据
            if not self.is_continuous and 'a' in data_name_list:
                a_idx = data_name_list.index('a')
                data[a_idx] = sth.int2one_hot(data[a_idx].astype(np.int32), self.a_dim)
            if 's' in data_name_list:
                s_idx = data_name_list.index('s')
                data[s_idx] = self.normalize_vector_obs(data[s_idx])
            if 's_' in data_name_list:
                s_idx = data_name_list.index('s_')
                data[s_idx] = self.normalize_vector_obs(data[s_idx])

            return dict([
                [n, d] for n, d in zip(data_name_list, list(map(self.data_convert, data)))
            ])

        def get_value_from_dict(self, data_name_list, data_dict):
            '''
            TODO: Annotation
            '''
            return [data_dict.get(n) for n in data_name_list]

        def _learn(self, function_dict: Dict):
            '''
            TODO: Annotation
            '''
            _pre_process = function_dict.get('pre_process_function', lambda *args: args)    # 预处理过程
            _train = function_dict.get('train_function', lambda *args: (None, {}))   # 训练过程
            _update = function_dict.get('update_function', lambda *args: None)  # maybe need update parameters of target networks
            _summary = function_dict.get('summary_dict', {})    # 记录输出到tensorboard的词典
            _sample_data_list = function_dict.get('sample_data_list', ['s', 'visual_s', 'a', 'r', 's_', 'visual_s_', 'done'])  # 需要从经验池提取的经验
            _train_data_list = function_dict.get('train_data_list', ['ss', 'vvss', 'a', 'r', 'done'])  # 需要从经验池提取的经验

            if self.data.is_lg_batch_size:
                # -----------初始化某些变量---------------
                self.intermediate_variable_reset()
                # --------------------------------------

                # --------------------------------------从经验池中获取数据
                data = self.get_transitions(data_name_list=_sample_data_list)  # default: s, visual_s, a, r, s_, visual_s_, done
                # --------------------------------------

                # --------------------------------------如果使用RNN， 就将s和s‘状态进行拼接处理
                if self.use_rnn:
                    data['ss'] = tf.concat([    # [B, T, N], [B, T, N] => [B, T+1, N]
                        data['s'],
                        data['s_'][:, -1:]
                    ], axis=1)
                    data['vvss'] = tf.concat([
                        data['visual_s'],
                        data['visual_s_'][:, -1:]
                    ], axis=1)
                # --------------------------------------如果不使用RNN， 就将s和s‘状态进行堆叠处理
                else:
                    data['ss'] = tf.concat([data['s'], data['s_']], axis=0)  # [B, N] => [2*B, N]
                    data['vvss'] = tf.concat([data['visual_s'], data['visual_s_']], axis=0)
                # --------------------------------------

                # --------------------------------------预处理过程
                data = _pre_process(data)[0]
                # --------------------------------------

                # --------------------------------------好奇心部分
                if self.use_curiosity:
                    # -------------------------------------- 如果使用RNN，那么需要将数据维度从三维转换为二维
                    if self.use_rnn:
                        data['s'] = tf.reshape(data['s'], [-1, data['s'].shape[-1]])    # [B, T, N] => [B*T, N]
                        data['s_'] = tf.reshape(data['s_'], [-1, data['s_'].shape[-1]])
                        data['visual_s'] = tf.reshape(data['visual_s'], [-1, data['visual_s'].shape[-1]])
                        data['visual_s_'] = tf.reshape(data['visual_s_'], [-1, data['visual_s_'].shape[-1]])
                    crsty_r, crsty_loss, crsty_summaries = self.curiosity_model(
                        *self.get_value_from_dict(data_name_list=['s', 'visual_s', 'a', 's_', 'visual_s_'], data_dict=data))
                    data['r'] += crsty_r
                    self.summaries.update(crsty_summaries)
                else:
                    crsty_loss = tf.constant(value=0., dtype=self._data_type)
                # --------------------------------------

                # --------------------------------------优先经验回放部分，获取重要性比例
                if self.use_priority and self.use_isw:
                    _isw = self.data.get_IS_w().reshape(-1, 1)  # [B, ] => [B, 1]
                    _isw = self.data_convert(_isw)
                else:
                    _isw = tf.constant(value=1., dtype=self._data_type)
                # --------------------------------------

                # --------------------------------------获取需要传给train函数的参数
                _training_data = self.get_value_from_dict(data_name_list=_train_data_list, data_dict=data)
                # --------------------------------------

                # --------------------------------------burn in隐状态部分
                if self.use_rnn:
                    cell_state = self.initial_cell_state()
                    if self.burn_in_time_step > 0:
                        _s, _visual_s = self.data.get_burn_in_states()
                        cell_state = self.get_burn_in_feature(_s, _visual_s, cell_state)
                else:
                    cell_state = (None,)
                # --------------------------------------

                # --------------------------------------训练主程序，返回可能用于PER权重更新的TD error，和需要输出tensorboard的信息
                td_error, summaries = _train(_training_data, _isw, crsty_loss, cell_state)
                # --------------------------------------

                # --------------------------------------更新summary
                self.summaries.update(summaries)
                # --------------------------------------

                # --------------------------------------优先经验回放的更新部分
                if self.use_priority:
                    td_error = np.squeeze(td_error.numpy())
                    self.data.update(td_error, self.train_step)
                # --------------------------------------

                # --------------------------------------target网络的更新部分
                _update()
                # --------------------------------------

                # --------------------------------------更新summary
                self.summaries.update(_summary)
                # --------------------------------------

                # --------------------------------------写summary到tensorboard
                self.write_training_summaries(self.global_step, self.summaries)
                # --------------------------------------

    return Off_Policy


text:
# -*- coding: utf-8 -*-
import iris
import numpy as np
import pytest
from dateutil.relativedelta import relativedelta
from iris.time import PartialDateTime

from wildfires.data.datasets import (
    ERA5_DryDayPeriod,
    MonthlyDataset,
    UnexpectedCoordinateError,
    dummy_lat_lon_cube,
)

from .utils import data_availability

@pytest.mark.slow
@data_availability
def test_era_temporal_shifting():
    normal = ERA5_DryDayPeriod()
    shifted = normal.get_temporally_shifted_dataset(months=-3)

    # Since the shifted datasets represents previous months' data, there will not be
    # an overlap at the very beginning of the dataset, but rather at the end.

    # 2-month wide selection window.
    end_date = normal.max_time
    start_date = end_date - relativedelta(months=2)

    normal_cubes = normal.get_monthly_data(
        start=PartialDateTime(year=start_date.year, month=start_date.month),
        end=PartialDateTime(year=end_date.year, month=end_date.month),
    )

    # Shift window by 3 months to take into account the dataset shift.
    shift_end_date = end_date + relativedelta(months=3)
    shift_start_date = start_date + relativedelta(months=3)

    shifted_cubes = shifted.get_monthly_data(
        start=PartialDateTime(year=shift_start_date.year, month=shift_start_date.month),
        end=PartialDateTime(year=shift_end_date.year, month=shift_end_date.month),
    )

    for cube, shifted_cube in zip(normal_cubes, shifted_cubes):
        assert np.all(np.isclose(cube.data, shifted_cube.data))
        if any(hasattr(c, "mask") for c in (cube, shifted_cube)):
            assert all(hasattr(c.data, "mask") for c in (cube, shifted_cube))
            assert np.all(cube.data.mask == shifted_cube.data.mask)

        assert cube.attributes == shifted_cube.attributes

class DummyDataset(MonthlyDataset):
    def __init__(self):
        np.random.seed(1)
        self.cubes = iris.cube.CubeList(
            [dummy_lat_lon_cube(np.random.random((10, 10, 10)), monthly=True)]
        )

@pytest.mark.parametrize("deep", (False, True))
def test_temporal_shifting(deep):
    normal = DummyDataset()
    shifted = normal.get_temporally_shifted_dataset(months=-3, deep=deep)

    # Confirm that the underlying data is identical (not copied) if `deep=False`.
    if deep:
        assert normal.cube.data is not shifted.cube.data
    else:
        assert normal.cube.data is shifted.cube.data

    # Since the shifted datasets represents previous months' data, there will not be
    # an overlap at the very beginning of the dataset, but rather at the end.

    # 2-month wide selection window.
    end_date = normal.max_time
    start_date = end_date - relativedelta(months=2)

    normal_cubes = normal.get_monthly_data(
        start=PartialDateTime(year=start_date.year, month=start_date.month),
        end=PartialDateTime(year=end_date.year, month=end_date.month),
    )

    # Shift window by 3 months to take into account the dataset shift.
    shift_end_date = end_date + relativedelta(months=3)
    shift_start_date = start_date + relativedelta(months=3)

    shifted_cubes = shifted.get_monthly_data(
        start=PartialDateTime(year=shift_start_date.year, month=shift_start_date.month),
        end=PartialDateTime(year=shift_end_date.year, month=shift_end_date.month),
    )

    # Test that the data matches.
    for cube, shifted_cube in zip(normal_cubes, shifted_cubes):
        assert np.all(np.isclose(cube.data, shifted_cube.data))
        if any(hasattr(c, "mask") for c in (cube, shifted_cube)):
            assert all(hasattr(c.data, "mask") for c in (cube, shifted_cube))
            assert np.all(cube.data.mask == shifted_cube.data.mask)

        assert cube.attributes == shifted_cube.attributes

def test_other_temporal_coords():
    """No other temporal coordinates should be permitted prior to shifting.

    The shifting operation may or may not affect these auxiliary coordinates, and
    predicting how they may be affected is outside the scope of this function.

    """
    # Add a month_number coordinate to verify that an Exception should be raised.
    dataset = DummyDataset()
    iris.coord_categorisation.add_month_number(dataset.cube, "time")
    with pytest.raises(UnexpectedCoordinateError):
        # Attempt to get shifted dataset.
        dataset.get_temporally_shifted_dataset(months=-1)

max_stars_repo_name: akuhnregnier/wildfires

max_issues_repo_name: akuhnregnier/wildfires

max_forks_repo_name: akuhnregnier/wildfires