text:
import numpy as np
from scipy.interpolate import interp1d
from .baseClass import Hysteresis, SimpleCycle, MonotonicCurve

# This is kind of unorganized right now...

# Todo:
    # make a concatenate to Cycle object.

def concatenateHys(*argv):
    """
    This function creates a new hysteresis from the xy data of a series of
    monotonic curves, or xy curves.

    Parameters
    ----------
    *argv : SimpleCycle objects, or XY data
        A number of monotonic cycle objects to be combined into a hysteresis.
        These curves should be

    Returns
    -------
    hysteresis : Hysteresis Object
        The ouput hysteresis object which is a combination of all of the input
        curves

    """

    # # TODO: enhance hysteresis functionality
    # I would expect that the new curve has all the propreties of the old curves.
    # Here that won't be the case, which is akward.

    xyList = [None]*len(argv)

    for ii, vector in enumerate(argv):
        # Try to read the xy data from a monotonic curve object
        try:
            tempxy = vector.xy
        except:
            tempxy = vector

        # for curves after the first curve, we skip the first value
        # I think we want to do this for all curves?
        if ii >= 1:
            xyList[ii] = tempxy[1:,:]
        else:
            xyList[ii] = tempxy

    # Create new hysteresis, then add objects to the list
    xy = np.concatenate(xyList)
    hysteresis = Hysteresis(xy)

    return hysteresis

# =============================================================================
# resample functions
# =============================================================================

def _linInterpolate(x,y, Nsamples):
    """
    A linear interpolation function that takes a target input curve to a
    target sample curve. The sample curve goes between x0 and xN with Nsamples.
    """

    f = interp1d(x, y)

    outputx = np.linspace(x[0],x[-1], Nsamples)
    outputy = f(outputx)
    outputxy = np.column_stack((outputx, outputy))
    return outputxy

def _getNsamples(Targetdx, dxNet):
    if Targetdx >= abs(dxNet/2):
        print('Targetdx is larger than dxNet/2, no intermediate points made for the target dx = ' +str(dxNet))
        Nsamples = 2
    else:
        Nsamples = int(round(abs(dxNet/Targetdx))) + 1
    return Nsamples

# =============================================================================
# Tools for creating a load protocol
# =============================================================================

def _RemoveNeg(x, y, direction):

    difference = np.append(0, np.diff(x))

    condition = np.where(0 <= difference*direction)

    xOut = x[condition]
    yOut = y[condition]

    xy = np.column_stack([xOut, yOut])

    return xy

def removeNegative(Curve):

    """
    Removes intermitent negative values in a simple curve.
    """

    # Get sample parameters, then pass those to the new curve.

    # if the curve is a SimpleCycle
    if isinstance(Curve, SimpleCycle):

        x = Curve.xy[:,0]
        y = Curve.xy[:,1]
        direction = Curve.direction
        Output = SimpleCycle(_RemoveNeg(x, y, direction))

    # if the curve is a hysteresis, we recursively create a series of Cycles
    elif isinstance(Curve, Hysteresis):
        outputCycles = [None]*Curve.NCycles
        for ii, Cycle in enumerate(Curve.Cycles):
            outputCycles[ii] = removeNegative(Cycle)
        Output = concatenateHys(*outputCycles)    # If Curve

    # if the curve is a Monotonic Cycle
    elif isinstance(Curve, MonotonicCurve):

        x = Curve.xy[:,0]
        y = Curve.xy[:,1]
        direction = Curve.direction

        print('Monotonic curves should have no reversals...')
        Output = MonotonicCurve(_RemoveNeg(x, y, direction))

    # if it is a np array
    elif isinstance(Curve, np.ndarray):
        x = Curve[:,0]
        y = Curve[:,1]

        # TODO: Create a standardized get direction function
        if x[0] <= x[-1]:
            direction = 1
        else:
            direction = -1

        Output = _RemoveNeg(x, y, direction)

    return  Output

max_stars_repo_name: cslotboom/hysteresis

max_issues_repo_name: cslotboom/hysteresis

max_forks_repo_name: cslotboom/hysteresis

reason: import numpy,from scipy


text:
"""Mean covariance estimation."""
from copy import deepcopy
import numpy as np

from .base import sqrtm, invsqrtm, logm, expm
from .ajd import ajd_pham
from .distance import distance_riemann
from .geodesic import geodesic_riemann

def _get_sample_weight(sample_weight, data):
    """Get the sample weights.

    If none provided, weights init to 1. otherwise, weights are normalized.
    """
    if sample_weight is None:
        sample_weight = np.ones(data.shape[0])
    if len(sample_weight) != data.shape[0]:
        raise ValueError("len of sample_weight must be equal to len of data.")
    sample_weight /= np.sum(sample_weight)
    return sample_weight

def mean_riemann(covmats, tol=10e-9, maxiter=50, init=None,
                 sample_weight=None):
    r"""Return the mean covariance matrix according to the Riemannian metric.

    The procedure is similar to a gradient descent minimizing the sum of
    Riemannian distances to the mean.

    .. math::
        \mathbf{C} = \arg\min{(\sum_i \delta_R ( \mathbf{C} , \mathbf{C}_i)^2)}

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param tol: the tolerance to stop the gradient descent
    :param maxiter: The maximum number of iteration, default 50
    :param init: A covariance matrix used to initialize the gradient descent. If None the Arithmetic mean is used
    :param sample_weight: the weight of each matrix
    :returns: the mean covariance matrix

    """  # noqa
    # init
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    if init is None:
        C = np.mean(covmats, axis=0)
    else:
        C = init

    k = 0
    nu = 1.0
    tau = np.finfo(np.float64).max
    crit = np.finfo(np.float64).max
    # stop when J<10^-9 or max iteration = 50
    while (crit > tol) and (k < maxiter) and (nu > tol):
        k = k + 1
        C12 = sqrtm(C)
        Cm12 = invsqrtm(C)
        J = np.zeros((n_channels, n_channels))

        for index in range(n_matrices):
            tmp = np.dot(np.dot(Cm12, covmats[index]), Cm12)
            J += sample_weight[index] * logm(tmp)

        crit = np.linalg.norm(J, ord='fro')
        h = nu * crit
        C = np.dot(np.dot(C12, expm(nu * J)), C12)
        if h < tau:
            nu = 0.95 * nu
            tau = h
        else:
            nu = 0.5 * nu

    return C

def mean_logeuclid(covmats, sample_weight=None):
    r"""Return the mean covariance matrix according to the log-Euclidean
    metric.

    .. math::
        \mathbf{C} = \exp{(\frac{1}{N} \sum_i \log{\mathbf{C}_i})}

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    """
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    T = np.zeros((n_channels, n_channels))
    for index in range(n_matrices):
        T += sample_weight[index] * logm(covmats[index])
    C = expm(T)

    return C

def mean_kullback_sym(covmats, sample_weight=None):
    """Return the mean covariance matrix according to KL divergence.

    This mean is the geometric mean between the Arithmetic and the Harmonic
    mean, as shown in [1]_.

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    References
    ----------
    .. [1]  Moakher, Maher, and Philipp G. Batchelor. "Symmetric
        positive-definite matrices: From geometry to applications and
        visualization." In Visualization and Processing of Tensor Fields, pp.
        285-298. Springer Berlin Heidelberg, 2006.
    """
    C_Arithmetic = mean_euclid(covmats, sample_weight)
    C_Harmonic = mean_harmonic(covmats, sample_weight)
    C = geodesic_riemann(C_Arithmetic, C_Harmonic, 0.5)

    return C

def mean_harmonic(covmats, sample_weight=None):
    r"""Return the harmonic mean of a set of covariance matrices.

    .. math::
        \mathbf{C} = \left(\frac{1}{N} \sum_i {\mathbf{C}_i}^{-1}\right)^{-1}

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    """
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    T = np.zeros((n_channels, n_channels))
    for index in range(n_matrices):
        T += sample_weight[index] * np.linalg.inv(covmats[index])
    C = np.linalg.inv(T)

    return C

def mean_logdet(covmats, tol=10e-5, maxiter=50, init=None, sample_weight=None):
    r"""Return the mean covariance matrix according to the logdet metric.

    This is an iterative procedure where the update is:

    .. math::
        \mathbf{C} = \left(\sum_i \left( 0.5 \mathbf{C} + 0.5 \mathbf{C}_i \right)^{-1} \right)^{-1}

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param tol: the tolerance to stop the gradient descent
    :param maxiter: The maximum number of iteration, default 50
    :param init: A covariance matrix used to initialize the iterative procedure. If None the Arithmetic mean is used
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    """  # noqa
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    if init is None:
        C = np.mean(covmats, axis=0)
    else:
        C = init
    k = 0
    crit = np.finfo(np.float64).max
    # stop when J<10^-9 or max iteration = 50
    while (crit > tol) and (k < maxiter):
        k = k + 1

        J = np.zeros((n_channels, n_channels))

        for index, Ci in enumerate(covmats):
            J += sample_weight[index] * np.linalg.inv(0.5 * Ci + 0.5 * C)

        Cnew = np.linalg.inv(J)
        crit = np.linalg.norm(Cnew - C, ord='fro')

        C = Cnew
    return C

def mean_wasserstein(covmats, tol=10e-4, maxiter=50, init=None,
                     sample_weight=None):
    r"""Return the mean covariance matrix according to the Wasserstein metric.

    This is an iterative procedure where the update is [1]_:

    .. math::
        \mathbf{K} = \left(\sum_i \left( \mathbf{K} \mathbf{C}_i \mathbf{K} \right)^{1/2} \right)^{1/2}

    with :math:`\mathbf{K} = \mathbf{C}^{1/2}`.

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param tol: the tolerance to stop the gradient descent
    :param maxiter: The maximum number of iteration, default 50
    :param init: A covariance matrix used to initialize the iterative procedure. If None the Arithmetic mean is used
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    References
    ----------
    .. [1] Barbaresco, F. "Geometric Radar Processing based on Frechet distance:
        Information geometry versus Optimal Transport Theory", Radar Symposium
        (IRS), 2011 Proceedings International.
    """  # noqa
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    if init is None:
        C = np.mean(covmats, axis=0)
    else:
        C = init
    k = 0
    K = sqrtm(C)
    crit = np.finfo(np.float64).max
    # stop when J<10^-9 or max iteration = 50
    while (crit > tol) and (k < maxiter):
        k = k + 1

        J = np.zeros((n_channels, n_channels))

        for index, Ci in enumerate(covmats):
            tmp = np.dot(np.dot(K, Ci), K)
            J += sample_weight[index] * sqrtm(tmp)

        Knew = sqrtm(J)
        crit = np.linalg.norm(Knew - K, ord='fro')
        K = Knew
    if k == maxiter:
        print('Max iter reach')
    C = np.dot(K, K)
    return C

def mean_euclid(covmats, sample_weight=None):
    r"""Return the mean covariance matrix according to the Euclidean metric :

    .. math::
        \mathbf{C} = \frac{1}{N} \sum_i \mathbf{C}_i

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    """
    return np.average(covmats, axis=0, weights=sample_weight)

def mean_ale(covmats, tol=10e-7, maxiter=50, sample_weight=None):
    """Return the mean covariance matrix according using the AJD-based
    log-Euclidean Mean (ALE). See [1].

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param tol: the tolerance to stop the gradient descent
    :param maxiter: The maximum number of iteration, default 50
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    Notes
    -----
    .. versionadded:: 0.2.4

    References
    ----------
    [1] M. Congedo, B. Afsari, A. Barachant, M. Moakher, 'Approximate Joint
    Diagonalization and Geometric Mean of Symmetric Positive Definite
    Matrices', PLoS ONE, 2015

    """
    sample_weight = _get_sample_weight(sample_weight, covmats)
    n_matrices, n_channels, _ = covmats.shape
    crit = np.inf
    k = 0

    # init with AJD
    B, _ = ajd_pham(covmats)
    while (crit > tol) and (k < maxiter):
        k += 1
        J = np.zeros((n_channels, n_channels))

        for index, Ci in enumerate(covmats):
            tmp = logm(np.dot(np.dot(B.T, Ci), B))
            J += sample_weight[index] * tmp

        update = np.diag(np.diag(expm(J)))
        B = np.dot(B, invsqrtm(update))

        crit = distance_riemann(np.eye(n_channels), update)

    A = np.linalg.inv(B)

    J = np.zeros((n_channels, n_channels))
    for index, Ci in enumerate(covmats):
        tmp = logm(np.dot(np.dot(B.T, Ci), B))
        J += sample_weight[index] * tmp

    C = np.dot(np.dot(A.T, expm(J)), A)
    return C

def mean_alm(covmats, tol=1e-14, maxiter=100,
             verbose=False, sample_weight=None):
    r"""Return Ando-Li-Mathias (ALM) mean

    Find the geometric mean recursively [1]_, generalizing from:

    .. math::
        \mathbf{C} = A^{\frac{1}{2}}(A^{-\frac{1}{2}}B^{\frac{1}{2}}A^{-\frac{1}{2}})^{\frac{1}{2}}A^{\frac{1}{2}}

    require a high number of iterations.

    This is the adaptation of the Matlab code proposed by Dario Bini and
    Bruno Iannazzo, http://bezout.dm.unipi.it/software/mmtoolbox/
    Extremely slow, due to the recursive formulation.

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param tol: the tolerance to stop iterations
    :param maxiter: maximum number of iteration, default 100
    :param verbose: indicate when reaching maxiter
    :param sample_weight: the weight of each matrix

    :returns: the mean covariance matrix

    Notes
    -----
    .. versionadded:: 0.2.8

    References
    ----------
    .. [1] T. Ando, C.-K. Li and R. Mathias, "Geometric Means", Linear Algebra
        Appl. 385 (2004), 305-334.
    """  # noqa
    sample_weight = _get_sample_weight(sample_weight, covmats)
    C = covmats
    C_iter = np.zeros_like(C)
    n_matrices = covmats.shape[0]
    if n_matrices == 2:
        alpha = sample_weight[1] / sample_weight[0] / 2
        X = geodesic_riemann(covmats[0], covmats[1], alpha=alpha)
        return X
    else:
        for k in range(maxiter):
            for h in range(n_matrices):
                s = np.mod(np.arange(h, h + n_matrices - 1) + 1, n_matrices)
                C_iter[h] = mean_alm(C[s], sample_weight=sample_weight[s])

            norm_iter = np.linalg.norm(C_iter[0] - C[0], 2)
            norm_c = np.linalg.norm(C[0], 2)
            if (norm_iter / norm_c) < tol:
                break
            C = deepcopy(C_iter)
        else:
            if verbose:
                print('Max number of iterations reached')
        return C_iter.mean(axis=0)

def mean_identity(covmats, sample_weight=None):
    r"""Return the identity matrix corresponding to the covmats sit size

    .. math::
        \mathbf{C} = \mathbf{I}_d

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :returns: the identity matrix of size n_channels

    """
    C = np.eye(covmats.shape[1])
    return C

def mean_covariance(covmats, metric='riemann', sample_weight=None, *args):
    """Return the mean covariance matrix according to the metric

    :param covmats: Covariance matrices, (n_matrices, n_channels, n_channels)
    :param metric: the metric (default 'riemann'), can be : 'riemann',
        'logeuclid', 'euclid', 'logdet', 'identity', 'wasserstein', 'ale',
        'alm', 'harmonic', 'kullback_sym' or a callable function
    :param sample_weight: the weight of each matrix
    :param args: the argument passed to the sub function
    :returns: the mean covariance matrix, (n_channels, n_channels)

    """
    if callable(metric):
        C = metric(covmats, sample_weight=sample_weight, *args)
    else:
        C = mean_methods[metric](covmats, sample_weight=sample_weight, *args)
    return C

mean_methods = {'riemann': mean_riemann,
                'logeuclid': mean_logeuclid,
                'euclid': mean_euclid,
                'identity': mean_identity,
                'logdet': mean_logdet,
                'wasserstein': mean_wasserstein,
                'ale': mean_ale,
                'harmonic': mean_harmonic,
                'kullback_sym': mean_kullback_sym,
                'alm': mean_alm}

def _check_mean_method(method):
    """checks methods """
    if isinstance(method, str):
        if method not in mean_methods.keys():
            raise ValueError('Unknown mean method')
        else:
            method = mean_methods[method]
    elif not hasattr(method, '__call__'):
        raise ValueError('mean method must be a function or a string.')
    return method

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

def _get_mask_from_nan(covmat):
    nan_col = np.all(np.isnan(covmat), axis=0)
    nan_row = np.all(np.isnan(covmat), axis=1)
    if not np.array_equal(nan_col, nan_row):
        raise ValueError('NaN values are not symmetric.')
    nan_inds = np.where(nan_col)
    subcovmat_ = np.delete(covmat, nan_inds, axis=0)
    subcovmat = np.delete(subcovmat_, nan_inds, axis=1)
    if np.any(np.isnan(subcovmat)):
        raise ValueError('NaN values must fill rows and columns.')
    mask = np.delete(np.eye(covmat.shape[0]), nan_inds, axis=1)
    return mask

def _get_masks_from_nan(covmats):
    masks = []
    for i in range(len(covmats)):
        masks.append(_get_mask_from_nan(covmats[i]))
    return masks

def _apply_masks(covmats, masks):
    maskedcovmats = []
    for i in range(len(covmats)):
        maskedcovmats.append(masks[i].T @ covmats[i] @ masks[i])
    return maskedcovmats

def nanmean_riemann(covmats, tol=10e-9, maxiter=50, init=None,
                    sample_weight=None):
    """Return the Riemannian NaN-mean of covariance matrices.

    The Riemannian NaN-mean is the masked Riemannian mean applied to covariance
    matrices corrupted by symmetric NaN values [1]_.

    Parameters
    ----------
    covmats : ndarray, shape (n_matrices, n_channels, n_channels)
        Covariance matrices, potentially corrupted by symmetric NaN values [1]_
    tol : float (default 10e-9)
        The tolerance to stop the gradient descent.
    maxiter : int (default 50)
        The maximum number of iteration.
    init : None | ndarray, shape (n_channels, n_channels) (default None)
        A covariance matrix used to initialize the gradient descent.
        If None, a regularized Euclidean NaN-mean is used.
    sample_weight : None | ndarray, shape (n_matrices) (default None)
        The weight of each matrix.

    Returns
    -------
    C : ndarray, shape (n_channels, n_channels)
        Riemannian NaN-mean.

    Notes
    -----
    .. versionadded:: 0.2.8

    References
    ----------
    .. [1] F. Yger, S. Chevallier, Q. Barthélemy, S. Sra. "Geodesically-convex
        optimization for averaging partially observed covariance matrices",
        ACML 2021.
    """
    n_matrices, n_channels, _ = covmats.shape
    if init is None:
        Cinit = np.nanmean(covmats, axis=0) + 1e-6 * np.eye(n_channels)
    else:
        Cinit = init

    C = maskedmean_riemann(
        np.nan_to_num(covmats),  # avoid nan contamination in matmul
        _get_masks_from_nan(covmats),
        tol=tol,
        maxiter=maxiter,
        init=Cinit,
        sample_weight=sample_weight
    )
    return C

def maskedmean_riemann(covmats, masks, tol=10e-9, maxiter=50, init=None,
                       sample_weight=None):
    """Return the masked Riemannian mean of covariance matrices.

    Given masks defined as semi-orthogonal matrices, the masked Riemannian mean
    of covariance matrices is obtained with a gradient descent minimizing the
    sum of Riemannian distances between masked covariance matrices and the
    masked mean [1]_.

    Parameters
    ----------
    covmats : ndarray, shape (n_matrices, n_channels, n_channels)
        Covariance matrices.
    masks : list of n_matrices ndarray of shape (n_channels, n_channels_i), \
            with different n_channels_i, such that n_channels_i <= n_channels
        Masks, defined as semi-orthogonal matrices. See [1]_.
    tol : float (default 10e-9)
        The tolerance to stop the gradient descent.
    maxiter : int (default 50)
        The maximum number of iteration.
    init : None | ndarray, shape (n_channels, n_channels) (default None)
        A covariance matrix used to initialize the gradient descent.
        If None, the Identity is used.
    sample_weight : None | ndarray, shape (n_matrices) (default None)
        The weight of each matrix.

    Returns
    -------
    C : ndarray, shape (n_channels, n_channels)
        Masked Riemannian mean.

    Notes
    -----
    .. versionadded:: 0.2.8

    References
    ----------
    .. [1] F. Yger, S. Chevallier, Q. Barthélemy, S. Sra. "Geodesically-convex
        optimization for averaging partially observed covariance matrices",
        ACML 2021.
    """
    sample_weight = _get_sample_weight(sample_weight, covmats)
    maskedcovmats = _apply_masks(covmats, masks)
    n_matrices, n_channels, _ = covmats.shape
    if init is None:
        C = np.eye(n_channels)
    else:
        C = init

    k = 0
    nu = 1.0
    tau = np.finfo(np.float64).max
    crit = np.finfo(np.float64).max
    while (crit > tol) and (k < maxiter) and (nu > tol):
        k = k + 1
        maskedC = _apply_masks(np.tile(C, (n_matrices, 1, 1)), masks)
        J = np.zeros((n_channels, n_channels))
        for i in range(n_matrices):
            C12, Cm12 = sqrtm(maskedC[i]), invsqrtm(maskedC[i])
            tmp = C12 @ logm(Cm12 @ maskedcovmats[i] @ Cm12) @ C12
            J += sample_weight[i] * masks[i] @ tmp @ masks[i].T
        C12, Cm12 = sqrtm(C), invsqrtm(C)
        C = C12 @ expm(Cm12 @ (nu * J) @ Cm12) @ C12

        crit = np.linalg.norm(J, ord='fro')
        h = nu * crit
        if h < tau:
            nu = 0.95 * nu
            tau = h
        else:
            nu = 0.5 * nu

    return C

ext: py

max_stars_repo_name: c-cameron/pyRiemann

max_issues_repo_name: c-cameron/pyRiemann

max_forks_repo_name: c-cameron/pyRiemann

reason: import numpy


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

from ..utils import NotFittedError, UnsupportedModelError
import numpy as np
import pandas as pd
import copy

from sklearn.cluster import KMeans
from probatus.utils import shap_helpers

def return_confusion_metric(y_true, y_score, normalize = False):
    """
    Computes a confusion metric as absolute difference between the y_true and y_score.
    If normalize eis set to tru, it will normalize y_score to the maximum value in the array
    Args:
        y_true: (np.ndarray or pd.Series) true targets
        y_score: (np.ndarray or pd.Series) model output
        normalize: boolean, normalize or not to the maximum value

    Returns: (np.ndarray or pd.Series) conflusion metric

    """

    if normalize:
        y_score = y_score/y_score.max()

    return np.abs(y_true - y_score)

class BaseInspector(object):
    def __init__(self, algotype, **kwargs):
        self.fitted = False
        self.algotype = algotype
        # TODO fix compilatiopn issue on  for hdbscan
        # if algotype =='dbscan':
        #     self.clusterer = hdbscan.HDBSCAN(prediction_data=True,**kwargs)
        if algotype =='kmeans':
            self.clusterer = KMeans(**kwargs)
        else:
            raise UnsupportedModelError("The algorithm {} is not supported".format(algotype))

    def __repr__(self):
        repr_ = "{},\n\t{}".format(self.__class__.__name__,self.algotype)
        if self.fitted:
            repr_ += "\n\tTotal clusters {}".format(np.unique(self.clusterer.labels_).shape[0])
        return repr_

    def fit_clusters(self, X):
        """
        Perform the fit of the clusters with the algorithm specified in the constructor
        Args:
            X: input features

        Returns: cluster labels

        """
        self.clusterer.fit(X)
        self.fitted = True

        return self

    def predict_clusters(self,X):
        if not self.fitted:
            raise NotFittedError("Inspector not fitter. Run .fit()")

        labels = None
        if self.algotype == 'kmeans':
            labels =  self.clusterer.predict(X)
        if self.algotype == 'dbscan':
            raise NotImplementedError("Implementation not finished (note the hdbscan package is not imported yet!)")
            #labels, strengths = hdbscan.approximate_predict(self.clusterer, X)

        return labels

    @staticmethod
    def assert_is_dataframe(df):

        if isinstance(df,pd.DataFrame):
            return df

        elif isinstance(df,np.ndarray) and len(df.shape)==2:
            return pd.DataFrame(df)

        else:
            raise NotImplementedError("Sorry, X needs to be a pd.DataFrame for for a 2 dimensional numpy array")

    @staticmethod
    def assert_is_series(series, index=None):

        if isinstance(series, pd.Series):
            return series
        elif isinstance(series, pd.DataFrame) and series.shape[1] == 1:
            return pd.Series(series.values.ravel(), index=series.index)
        elif isinstance(series, np.ndarray) and len(series.shape) == 1 and index is not None:
            return pd.Series(series, index=index)
        else:
            raise TypeError(
                "The object should be a pd.Series, a dataframe with one collumn or a 1 dimensional numpy array")

class InspectorShap(BaseInspector):
    """
    Class to perform inspection of the model prediction based on Shapley values.

    It uses the calculated shapley values for the train model to build clusters in the shap space.
    For each cluster, an average confusion, average predicted probability and observed rate of a single class is
    calculated.
    Every sub cluster can be retrieved with the function slice_cluster to perform deeper analysis.

    The original dataframe indexing is used in slicing the dataframe, ensuring easy filtering

    Args:
            model: (obj) pretrained model (with sklearn-like API)
            algotype: (str) clustering algorithm (supported are kmeans and hdbscan)
            confusion_metric: (str) Confusion metric to use:
                - "proba": it will calculate the confusion metric as the absolute value of the target minus the predicted
                           probability. This provides a continuous measure of confusion, where 0 indicated correct predictions
                           and the closer the number is to 1, the higher the confusion
            normalize_probability: (boolean) if true, it will normalize the probabilities to the max value when computing
                the confusion metric
            cluster_probabilities: (boolean) if true, uses the model prediction as an input for the cluster prediction
            **kwargs: keyword arguments for the clustering algorithm

    """

    def __init__(self, model, algotype='kmeans', confusion_metric = 'proba',
                 normalize_probability=False,cluster_probability = False, **kwargs):

        super().__init__(algotype, **kwargs)
        self.model = model
        self.isinspected = False
        self.hasmultiple_dfs = False
        self.normalize_proba = normalize_probability
        self.cluster_probabilities = cluster_probability
        self.agg_summary_df = None
        self.set_names = None
        self.confusion_metric = confusion_metric
        self.cluster_report = None
        self.y = None
        self.predicted_proba = None
        self.X_shap = None
        self.clusters = None
        self.init_eval_set_report_variables()

        if confusion_metric not in ['proba']:
            #TODO implement the target method
            raise NotImplementedError("confusion metric {} not supported. See docstrings".format(confusion_metric))

    def __repr__(self):
        repr_ = "{},\n\t{}".format(self.__class__.__name__, self.algotype)
        if self.fitted:
            repr_ += "\n\tTotal clusters {}".format(np.unique(self.clusterer.labels_).shape[0])
        return repr_

    def init_eval_set_report_variables(self):
        self.X_shaps = list()
        self.clusters_list = list()
        self.ys = list()
        self.predicted_probas = list()

    def compute_probabilities(self, X):
        """
        Compute the probabilities for the model using the sklearn API
        Args:
            X: Feature set

        Returns: (np.array) probability

        """

        return self.model.predict_proba(X)[:,1]

    def fit_clusters(self, X):
        """
        Perform the fit of the clusters with the algorithm specified in the constructor
        Args:
            X: input features

        """
        X = copy.deepcopy(X)

        if self.cluster_probabilities:
            X['probs'] = self.predicted_proba

        return super().fit_clusters(X)

    def predict_clusters(self,X):
        """
        Predicts the clusters of the dataset X
        Args:
            X: features

        Returns: cluster labels

        """
        X = copy.deepcopy(X)

        if self.cluster_probabilities:
            X['probs'] = self.predicted_proba

        return super().predict_clusters(X)

    def fit(self, X, y=None, eval_set = None, sample_names=None, **shap_kwargs):
        """
        Fits and orchestrates the cluster calculations
        Args:
            X: (pd.DataFrame) with the features set used to train the model
            y: (pd.Series, default=None): targets used to train the model
            eval_set: (list, default=None). list of tuples in the shape (X,y) containing evaluation samples, for example
                a test sample, validation sample etc... X corresponds to the feature set of the sample, y corresponds
                to the targets of the samples
            sample_names: (list of strings, default=None): list of suffixed for the samples. If none, it will be labelled with
                sample_{i}, where i corresponds to the index of the sample.
                List length must match that of eval_set
            **shap_kwargs:  kwargs to pass to the Shapley Tree Explained
        """

        self.set_names = sample_names
        if sample_names is not None:
            # Make sure that the amount of eval sets matches the set names
            assert len(eval_set) == len(sample_names), "set_names must be the same length as eval_set"

        self.y, self.predicted_proba, self.X_shap, self.clusters = self.perform_fit_calc(X=X, y=y,
                                                                                fit_clusters=True,
                                                                                **shap_kwargs)
        if eval_set is not None:
            assert isinstance(eval_set, list), "eval_set needs to be a list"

            self.hasmultiple_dfs = True
            # Reset lists in case inspect run multiple times
            self.init_eval_set_report_variables()

            for X_, y_ in eval_set:
                y_, predicted_proba_, X_shap_, clusters_ = self.perform_fit_calc(X=X_, y=y_, fit_clusters=False,
                                                                                **shap_kwargs)

                self.X_shaps.append(X_shap_)
                self.ys.append(y_)
                self.predicted_probas.append(predicted_proba_)
                self.clusters_list.append(clusters_)

        return self

    def perform_fit_calc(self, X, y, fit_clusters=False, **shap_kwargs):
        """
        Performs cluster calculations for a specific X and y
        Args:
            X: pd.DataFrame with the features set used to train the model
            y: pd.Series (default None): targets used to train the model
            fit_clusters: flag indicating whether clustering algorithm should be trained with computed shap values
            **shap_kwargs:  kwargs to pass to the Shapley Tree Explained
        """

        X = self.assert_is_dataframe(X)
        y = self.assert_is_series(y, index = X.index)

        # Compute probabilities for the input X using model
        predicted_proba = pd.Series(self.compute_probabilities(X), index = y.index,name = 'pred_proba')

        # Compute SHAP values and cluster them
        X_shap = shap_helpers.shap_to_df(self.model, X, **shap_kwargs)
        if fit_clusters:
            self.fit_clusters(X_shap)
        clusters = pd.Series(self.predict_clusters(X_shap), index=y.index, name='cluster_id')
        return y, predicted_proba, X_shap, clusters

    def _compute_report(self):
        """
        Helper function to compute the report of the ispector - performs aggregations per cluster id

        """

        self.summary_df = self.create_summary_df(self.clusters, self.y, self.predicted_proba, normalize=self.normalize_proba)
        self.agg_summary_df = self.aggregate_summary_df(self.summary_df)

        if self.hasmultiple_dfs:

            self.summary_dfs = [
                self.create_summary_df(clust, y, pred_proba, normalize=self.normalize_proba)
                for  clust, y, pred_proba in zip(self.clusters_list, self.ys, self.predicted_probas)
            ]

            self.agg_summary_dfs = [
                self.aggregate_summary_df(df)
                for df in self.summary_dfs
            ]

    def compute(self):
        """
        Calculates a report containing the information per cluster.
        Includes the following:
            - cluster id
            - total number of observations in the cluster
            - total number of target 1 in the cluster
            - target 1 rate (ration of target 1 counts/observations)
            - average predicted probabilitites
            - average confusion

        If multiple eval_sets were passed in the inspect() functions, the output will contain those aggregations as well
        The output names will use the sample names provided in the inspect function. Otherwise they will be labelled by
        the suffix sample_{i}, where i is the index of the sample

        Returns: (pd.DataFrame) with above mentioned aggregations.

        """

        if self.cluster_report is not None:
            return self.cluster_report

        self._compute_report()
        out = copy.deepcopy(self.agg_summary_df)

        if self.hasmultiple_dfs:

            for ix, agg_summary_df in enumerate(self.agg_summary_dfs):
                if self.set_names is None:
                    sample_suffix = "sample_{}".format(ix+1)
                else: sample_suffix = self.set_names[ix]

                out = pd.merge(out, agg_summary_df, how="left", on='cluster_id',
                               suffixes = ('','_{}'.format(sample_suffix)))

        self.cluster_report = out
        return self.cluster_report

    def slice_cluster(self, cluster_id,  summary_df=None, X_shap=None, y=None, predicted_proba=None,
                      complementary = False):
        """
        Slices the input dataframes by the cluster.

        Args:
            cluster_id: (int or list for multiple cluster_id) cluster ids to to slice
            summary_df: Optional parameter - the summary_df on which the masking should be performed.
                if not passed the slicing is performed on summary generated by inspect method on X and y
            X_shap: Optional parameter - the SHAP values generated from on X on which the masking should be performed.
                if not passed the slicing is performed on X_shap generated by inspect method on X and y
            y: Optional parameter - the y on which the masking should be performed.
                if not passed the slicing is performed on y passed to inspect
            predicted_proba: Optional parameter - the predicted_proba on which the masking should be performed.
                if not passed the slicing is performed on predicted_proba generated by inspect method on X and y
            complementary: flag that returns the cluster_id if set to False, otherwise the complementary dataframe (ie
                those with ~mask)

        Returns: tuple: Dataframe of sliced shapley values, series of sliced targets, sliced probabilities

        """

        if self.cluster_report is None:
            self.compute()

        # Check if input specified by user, otherwise use the ones from self
        if summary_df is None:
            summary_df = self.summary_df
        if X_shap is None:
            X_shap = self.X_shap
        if y is None:
            y = self.y
        if predicted_proba is None:
            predicted_proba = self.predicted_proba

        mask = self.get_cluster_mask(summary_df, cluster_id)
        if not complementary:
            return X_shap[mask], y[mask], predicted_proba[mask]
        else:
            return X_shap[~mask], y[~mask], predicted_proba[~mask]

    def slice_cluster_eval_set(self, cluster_id, complementary=False):
        """
        Slices the input dataframes passed in the eval_set in the inspect function by the cluster id.

        Args:
            cluster_id: (int or list for multiple cluster_id) cluster ids to to slice
            complementary: flag that returns the cluster_id if set to False, otherwise the complementary dataframe (ie
                those with ~mask)

        Returns: list of tuplse: each element of the list containst
                    Dataframe of sliced shapley values, series of sliced targets, sliced probabilities

        """

        if not self.hasmultiple_dfs:
            raise NotFittedError("You did not fit the eval set. Please add an eval set when calling inspect()")

        output = []
        for X_shap, y, predicted_proba, summary_df in zip(
                self.X_shaps, self.ys, self.predicted_probas, self.summary_dfs):
            output.append(self.slice_cluster(cluster_id=cluster_id, summary_df=summary_df, X_shap=X_shap, y=y,
                                             predicted_proba=predicted_proba, complementary=complementary))
        return output

    @staticmethod
    def get_cluster_mask(df, cluster_id):
        """
        Returns the mask to filter the cluster id
        Args:
            df: dataframe with 'cluster_id' in it
            cluster_id: int or list of cluster ids to mask

        Returns:

        """

        if not isinstance(cluster_id, list):
            cluster_id = [cluster_id]

        mask = df['cluster_id'].isin(cluster_id)
        return mask

    @staticmethod
    def create_summary_df(cluster,y, probas, normalize=False):
        """
        Creates a summary by concatenating the cluster series, the targets, the probabilities and the measured confusion
        Args:
            cluster: pd.Series od clusters
            y: pd.Series od targets
            probas: pd.Series of predicted probabilities of the model
            normalize: boolean (if the predicted probabilities should be normalized to the max value

        Returns: pd.DataFrame (concatenation of the inputs)

        """

        confusion = return_confusion_metric(y, probas, normalize = normalize).rename("confusion")

        summary = [
            cluster,
            y.rename("target"),
            probas,
            confusion
        ]

        return pd.concat(summary, axis=1)

    @staticmethod
    def aggregate_summary_df(df):
        """
        Performs the aggregations at the cluster_id level needed to generate the report of the inspection
        Args:
            df: input df to aggregate

        Returns: pd.Dataframe with aggregation results
        """

        out = df.groupby("cluster_id").agg(
            total_label_1=pd.NamedAgg(column='target', aggfunc="sum"),
            total_entries=pd.NamedAgg(column='target', aggfunc="count"),
            label_1_rate=pd.NamedAgg(column='target', aggfunc="mean"),
            average_confusion=pd.NamedAgg(column='confusion', aggfunc="mean"),
            average_pred_proba=pd.NamedAgg(column='pred_proba', aggfunc="mean"),
        ).reset_index().rename(columns={"index": "cluster_id"}).sort_values(by='cluster_id')

        return out

    def fit_compute(self, X, y=None, eval_set=None, sample_names=None, **shap_kwargs):
        """
        Fits and orchestrates the cluster calculations and returns the computed report

        Args:
            X: (pd.DataFrame) with the features set used to train the model
            y: (pd.Series, default=None): targets used to train the model
            eval_set: (list, default=None). list of tuples in the shape (X,y) containing evaluation samples, for example
                a test sample, validation sample etc... X corresponds to the feature set of the sample, y corresponds
                to the targets of the samples
            sample_names: (list of strings, default=None): list of suffixed for the samples. If none, it will be labelled with
                sample_{i}, where i corresponds to the index of the sample.
                List length must match that of eval_set
            **shap_kwargs:  kwargs to pass to the Shapley Tree Explained

        Returns:
            (pd.DataFrame) Report with aggregations described in compute() method.
        """
        self.fit(X, y, eval_set, sample_names, **shap_kwargs)
        return self.compute()

max_stars_repo_name: sbjelogr/probatus

max_issues_repo_name: sbjelogr/probatus

max_forks_repo_name: sbjelogr/probatus


text:
"""Define FastAPI object."""
import jax
from fastapi import FastAPI
from fastapi.middleware.cors import CORSMiddleware
from pydantic import BaseSettings

app = FastAPI()

origins = ['http://localhost:3000']

app.add_middleware(
    CORSMiddleware,
    allow_origins=origins,
    allow_methods=["*"],
    allow_headers=["*"],
)

class Settings(BaseSettings):
    """Pydantic model of app settings."""

    device: str = jax.default_backend()
    target: str = 'cuda' if jax.default_backend() == 'gpu' else 'cpu'

    class Config:
        """Pydantic config."""

        allow_mutation = False

SETTINGS = Settings()

@app.get('/settings/gpu')
def check_gpu_enabled() -> bool:
    """Check if gpu is enabled."""
    return SETTINGS.device == 'gpu'

reason: import jax


text:
import asyncio
import logging
import struct
from typing import Optional

import numpy as np
import pytest

from tickit.core.device import DeviceUpdate
from tickit.core.typedefs import SimTime
from tickit.devices.cryostream.cryostream import CryostreamDevice
from tickit.devices.cryostream.states import PhaseIds
from tickit.devices.cryostream.status import ExtendedStatus

# # # # # Cryostream Tests # # # # #

@pytest.fixture
def cryostream() -> CryostreamDevice:
    return CryostreamDevice()

def test_cryostream_constructor():
    CryostreamDevice()

def test_cryostream_update_hold(cryostream):
    starting_temperature = cryostream.gas_temp
    device_update: DeviceUpdate = cryostream.update(time=SimTime(1e9), inputs={})
    device_update.outputs["temperature"] == starting_temperature

@pytest.mark.asyncio
async def test_cryostream_update_cool(cryostream: CryostreamDevice):
    starting_temperature = cryostream.gas_temp
    target_temperature = starting_temperature - 50
    await cryostream.cool(target_temperature)

    time = SimTime(0)
    device_update: DeviceUpdate
    while cryostream.phase_id != PhaseIds.HOLD.value:
        logging.info(f"Running time step: {time}")
        device_update = cryostream.update(time, inputs={})
        time_update: Optional[SimTime] = device_update.call_at

        if time_update is None:
            time = SimTime(int(time) + int(1e9))
            continue
        else:
            time = time_update

    max_diff = 10
    margin_of_error = np.array([+max_diff, -max_diff])
    assert any(
        (device_update.outputs["temperature"] + margin_of_error) == target_temperature
    )

@pytest.mark.asyncio
async def test_cryostream_update_end(cryostream: CryostreamDevice):
    starting_temperature = cryostream.default_temp_shutdown - 100
    cryostream.gas_temp = starting_temperature
    await cryostream.end(cryostream.default_ramp_rate)
    assert cryostream.phase_id == PhaseIds.RAMP.value
    assert cryostream.gas_flow == 5

    time = SimTime(0)
    device_update: DeviceUpdate
    while cryostream.phase_id != PhaseIds.HOLD.value:
        logging.info(f"Running time step: {time}")
        device_update = cryostream.update(time, inputs={})
        time_update: Optional[SimTime] = device_update.call_at

        if time_update is None:
            continue
        else:
            time = time_update

    max_diff = 10
    target = cryostream.default_temp_shutdown
    assert (
        target - max_diff <= device_update.outputs["temperature"] <= target + max_diff
    )

@pytest.mark.asyncio
async def test_cryostream_update_plat(cryostream: CryostreamDevice):
    starting_temperature = cryostream.gas_temp
    await cryostream.plat(5)
    assert cryostream.phase_id == PhaseIds.PLAT.value

    time = SimTime(0)
    time_update: Optional[SimTime] = time
    device_update: DeviceUpdate
    while time_update is not None:
        logging.info(f"Running time step: {time}")
        device_update = cryostream.update(time, inputs={})
        time_update = device_update.call_at

        if time_update is None:
            continue
        else:
            time = time_update

    assert device_update.outputs["temperature"] == starting_temperature

@pytest.mark.asyncio
@pytest.mark.parametrize(
    "tickit_task", ["examples/configs/cryo-tcp.yaml"], indirect=True
)
async def test_cryostream_system(tickit_task):
    reader, writer = await asyncio.open_connection("localhost", 25565)

    async def write(data: bytes):
        writer.write(data)
        await writer.drain()

    async def get_status() -> ExtendedStatus:
        return ExtendedStatus.from_packed(await reader.read(42))

    # Check we start at 300K
    status = await get_status()
    assert status.gas_temp == 30000
    assert status.turbo_mode == 0
    # Start a ramp
    await write(b"\x06\x0b" + struct.pack(">HH", 360, 30030))
    # Check we go to to 300.3K, eventually
    status = await get_status()
    assert status.gas_temp == pytest.approx(30015, rel=5)
    status = await get_status()
    assert status.gas_temp == pytest.approx(30025, rel=5)
    status = await get_status()
    assert status.gas_temp == pytest.approx(30030, rel=5)
    # Turbo on
    await write(b"\x03\x14\x01")
    status = await get_status()
    assert status.gas_temp == pytest.approx(30030, rel=5)
    assert status.turbo_mode == 1
    # Cool to 300K
    await write(b"\x04\x0e" + struct.pack(">H", 30000))
    status = await get_status()
    assert status.gas_temp == pytest.approx(30000, rel=5)

max_stars_repo_name: dls-controls/tickit

max_issues_repo_name: dls-controls/tickit

max_forks_repo_name: dls-controls/tickit


text:
import pygame
from bisect import bisect_left
from PIL import Image
from bresenham import bresenham
import numpy as np
import delegator
import math
from threading import Thread
from subprocess import Popen, PIPE
from pyo import *
from time import time
from collections import namedtuple

class WacomValues:

    def __init__(self):
        "docstring"
        self.tilt_x = 0
        self.tilt_y = 0
        self.pressure = 0
        self.thread = Thread(target=lambda : self.update_wacom_stuff())
        self.thread.start()

    def update_wacom_stuff(self):
        libin = Popen(["unbuffer", "libinput", "debug-events"], stdout=PIPE)
        while True:
            l = libin.stdout.readline().decode("utf-8")
            # tilt 3
            # pressure 4
            try:
                if "TABLET_TOOL" in l:
                    self.tilt_x, self.tilt_y = [min(max((float(num.strip())+60)/120, -1.0), 1.0) for num in l.split("\t")[3].split("tilt:")[1].replace("*","").split("/")]
                    if "pressure" in l:
                        self.pressure = float(l.split("\t")[4].split("pressure:")[1].split(" ")[1].replace("*",""))
            except:
                pass

class Brush:
    def preview(self, points):
        pyo_sine_pan.setPan(points[0].pan)
        for point in points:
            faderidx = int((1-point.y)*current_surface.get_height()*preview_quality_factor)
            if faderidx < len(pyo_sine_faders):
                fader = pyo_sine_faders[faderidx]
                fader.mul = fader_loudness*(float(point.amp))
                fader.play()

class PixelBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        points = [Point(y, 1, x)]
        current_spec_unit.modify_or_insert_line(points, x)

class HarmonicBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        points = [Point(y, 1, x)]
        ratio = 2
        bandwidth = math.log(max_freq/min_freq,2)
        next_y = (bandwidth*y - math.log(ratio,2)) / bandwidth
        while next_y > 0 and ratio < 12:
            points.append(Point(next_y, 1/(ratio/1.5), x))
            ratio+=1
            next_y = (bandwidth*y - math.log(ratio,2)) / bandwidth

        current_spec_unit.modify_or_insert_line(points, x)

class WSampleLoopBrush(Brush):
    def __init__(self):

        loadt = Thread(target=lambda:self.load())
        loadt.start()

    def load(self):
        self.loaded = False
        bpo = (current_surface.get_height()-1)/math.log(max_freq/min_freq,2)
        self.sndTable = SndTable('sample.wav')
        delegator.run("arss sample.wav sample.bmp --min-freq {} --max-freq {} --pps 150 --bpo {}".format(min_freq, max_freq, bpo))
        self.sample = np.asarray(Image.open("sample.bmp"))
        self.idx = 0
        print("loaded")
        self.loaded = True

    def preview(self, offset):
        if pyo_sample_reader.table != self.sndTable:
            pyo_sample_reader.setTable(self.sndTable)
            pyo_sample_reader.setFreq(self.sndTable.getRate())

        # number of octaves
        bandwidth = math.log(max_freq/min_freq,2)
        # number of semitones in the bandwidth
        note_num = bandwidth*12
        transpo = offset * note_num
        pyo_sample_transpo.setTranspo(transpo)
        pyo_sample_fader.mul = wacom.pressure
        pyo_sample_reader.setPhase(self.idx/len(self.sample[0]))
        pyo_sample_pan.setPan(wacom.tilt_x)
        pyo_sample_fader.play()

    def draw(self, x, y):
        if not self.loaded:
            return
        x = x/current_surface.get_width()
        offset = y/current_surface.get_height() - 0.5
        increment = 1/current_surface.get_height()
        points = []
        for i in range(len(self.sample)):
            pointy = i * increment + offset
            if wacom.pressure * (self.sample[i][self.idx][0]) <=3:
                continue
            if pointy >= 1 or pointy < 0:
                continue
            points.append(Point(pointy, wacom.pressure * (self.sample[i][self.idx][0])/255, wacom.tilt_x))
        # advance one point in the self.sample and loop
        self.preview((1 - y/current_surface.get_height())-0.5)
        current_spec_unit.modify_or_insert_line(points, x)
        self.idx+=1
        self.idx%=len(self.sample[0])

class WHarmonicBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        points = [Point(y, wacom.pressure, wacom.tilt_x)]
        ratio = 2
        # frequency bandwidth in number of octaves
        bandwidth = math.log(max_freq/min_freq,2)
        next_y = (bandwidth*y - math.log(ratio,2)) / bandwidth
        while next_y > 0 and ratio < 12:
            multip = 1
            if ratio %2 == 0:
                multip = wacom.tilt_y * 2
            else:
                multip = (1 - wacom.tilt_y) * 2
            amp = wacom.pressure/max(((ratio/1.5) * multip),1)
            points.append(Point(next_y, amp, wacom.tilt_x))
            ratio+=1
            next_y = (bandwidth*y - math.log(ratio,2)) / bandwidth
        self.preview(points)
        current_spec_unit.modify_or_insert_line(points, x)

class WSineBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        points = [Point(y, wacom.pressure, wacom.tilt_x)]
        self.preview(points)
        current_spec_unit.modify_or_insert_line(points, x)

class WNoiseBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        increment = 1/current_surface.get_height()
        points = []
        h = int(wacom.tilt_y*(current_surface.get_height()/3))
        highest = 1 - max(y - (h*increment),0)
        lowest = 1 - min(y + (h*increment),1)
        for i in range(-h,h+1):
            points.append(Point(min(max(y+i*increment, 0), 1 - increment), wacom.pressure, wacom.tilt_x))
        self.preview(highest, lowest, wacom.pressure, wacom.tilt_x)
        current_spec_unit.modify_or_insert_line(points, x)

    def preview(self, highest, lowest, amp, pan):
        # number of octaves
        bandwidth = math.log(max_freq/min_freq,2)
        # number of midi note in the bandwidth
        note_num = bandwidth*12
        # the target midi note is the normalized location in y + the base midi note offset

        highest * bandwidth
        pyo_noise_hp.setFreq(midi_to_freq((note_num * lowest) + freq_to_midi(min_freq)))
        pyo_noise_lp.setFreq(midi_to_freq((note_num * highest) + freq_to_midi(min_freq)))
        pyo_noise_pan.setPan(pan)
        pyo_noise_fader.mul = amp
        pyo_noise_fader.play()

class NoiseBrush(Brush):
    def draw(self, x, y):
        x = x/current_surface.get_width()
        y = y/current_surface.get_height()
        increment = 1/current_surface.get_height()
        points = []
        for i in range(-10,11):
            points.append(Point(y+i*increment, 1, x))
        current_spec_unit.modify_or_insert_line(points, x)

class Point:
    def __init__(self, y, amp, pan):
        self.y = y
        self.amp = amp
        self.pan = pan
        self.dirty = True
        self.gesture_id = current_gesture

    def render(self, x_value):
        if not self.dirty:
            return
        self.dirty = False
        y_value = round(self.y * current_surface.get_height())
        l_amp, r_amp = [self.amp * math.sin(self.pan*math.pi/2), self.amp * math.cos(self.pan*math.pi/2)]
        pixel_color = [l_amp*255,0 , r_amp*255]

        canvas_l[y_value][x_value] += pixel_color[2]
        canvas_l[y_value][x_value] = min(canvas_l[y_value][x_value], 255)
        canvas_r[y_value][x_value] += pixel_color[0]
        canvas_r[y_value][x_value] = min(canvas_r[y_value][x_value], 255)
        try:
            current_surface.fill([canvas_r[y_value][x_value],0,canvas_l[y_value][x_value]], ((x_value,y_value),(1,1)))
        except:
            pass

class Line:
    def __init__(self, x, width):
        self.x = x
        self.width = 1
        self.points = []
        self.points_y = []
        self.dirty = True
        self.gesture_id = current_gesture

    def render(self):
        if not self.dirty:
            return
        self.dirty = False

        x_value = round(self.x * current_surface.get_width())

        for point in self.points:
            point.render(x_value)

    def add_point(self, point):
        self.dirty = True
        idx = bisect_left(self.points_y, point.y)
        adjacent_point = self.points[idx] if idx < len(self.points) else None
        if adjacent_point and adjacent_point.y == point.y:
            adjacent_point.pan = point.pan
            adjacent_point.amp = point.amp
            adjacent_point.dirty = True
        else:
            self.points.insert(idx, point)
            self.points_y.insert(idx, point.y)

    def add_points(self, points):
        for point in points:
            self.add_point(point)

class SpectrumUnit:
    def __init__(self):
        self.lines = []
        # keeping this alows us to not build this list every time we need to bissect
        self.lines_x = []

    def modify_or_insert_line(self, points, x):
        line = None
        # right now, can only add one line at a time
        idx = bisect_left(self.lines_x, x)
        adjacent_line = self.lines[idx] if idx < len(self.lines) else None
        if adjacent_line and adjacent_line.x == x:
            line = adjacent_line
        else:
            line = Line(x, 1/current_surface.get_width())
            self.lines.insert(idx, line)
            self.lines_x.insert(idx, x)
        line.add_points(points)

    def invalidate_canvas(self):
        global canvas_l
        global canvas_r
        canvas_l = [[0 for i in range(w)] for i in range(h)]
        canvas_r = [[0 for i in range(w)] for i in range(h)]
        current_surface.fill([0,0,0])
        for line in self.lines:
            line.dirty = True
            for point in line.points:
                point.dirty = True

    def remove_last_gesture(self):
        gesture_to_remove = gesture_stack.pop()
        for line in self.lines:
            line.points = [point for point in line.points if point.gesture_id != gesture_to_remove]
        self.invalidate_canvas()

h = 1000
w = 1800

max_freq = 20000
min_freq = 27.5

canvas_l = [[0 for i in range(w)] for i in range(h)]
canvas_r = [[0 for i in range(w)] for i in range(h)]

gesture_stack = []
current_spec_unit = SpectrumUnit()
current_surface = None
current_brush = None
wacom = WacomValues()
current_gesture = 0
pyo_server = Server(audio="jack").boot()
pyo_server.start()
def midi_to_freq(note_num):
    return (440 / 32) * (2 ** ((note_num - 9) / 12))

def freq_to_midi(freq):
    return math.log(freq/440.0,2) * 12 + 69

preview_quality_factor = 0.3
bpo = ((preview_quality_factor*h)-1)/math.log(max_freq/min_freq,2)
midi_notes_increment = 12/bpo
fader_loudness = 0.1

pyo_sine_faders = [Fader(fadeout=0.03, dur=0.04, mul=fader_loudness) for i in range(int(h*preview_quality_factor))]
pyo_sines = FastSine(mul=pyo_sine_faders,freq=[midi_to_freq(freq_to_midi(min_freq)+i*midi_notes_increment) for i in range(int(h*preview_quality_factor))])
pyo_sine_pan = Pan(pyo_sines).out()

pyo_noise_fader = Fader(fadeout=0.03, dur=0.04, mul=1)
pyo_noise = Noise(mul=pyo_noise_fader)
pyo_noise_hp = Biquadx(pyo_noise, q=1, type=1)
pyo_noise_lp = Biquadx(pyo_noise_hp, q=1)
pyo_noise_pan = Pan(pyo_noise_lp).out()

pyo_sample_fader = Fader(fadeout=0.03, dur=0.04, mul=1)
pyo_sample_reader = Osc(HarmTable([1]), mul=pyo_sample_fader)
pyo_sample_transpo = Harmonizer(pyo_sample_reader, transpo=0, winsize=0.03)
pyo_sample_pan = Pan(pyo_sample_transpo).out()

last_grid_len = 33
grid_len = 33

def main():

    global last_grid_len
    global grid_len
    pressing = False

    pygame.init()
    pygame.display.set_caption("dessin")
    global current_gesture
    screen = pygame.display.set_mode((w,h))
    global current_surface
    current_surface = screen
    brushes = [WSineBrush(), WHarmonicBrush(), WNoiseBrush(), WSampleLoopBrush()]
    current_brush = brushes[0]
    last_mouse_pos = None
    running = True
    framenum=0
    ctrl= False
    playing_start = False
    playing = False
    last_x = 0
    grid_changed = True
    while running:
        for line in current_spec_unit.lines:
            line.render()

        if grid_changed:
            for i in range(int(current_surface.get_width()/last_grid_len)):
                current_surface.fill([0,40,0,0.1], ((i * last_grid_len,0),(1,current_surface.get_height())), pygame.BLEND_RGBA_SUB)

            for i in range(int(current_surface.get_width()/grid_len)):
                current_surface.fill([0,40,0,0.1], ((i * grid_len,0),(1,current_surface.get_height())), pygame.BLEND_RGBA_ADD)

            last_grid_len = grid_len
            grid_changed = False

        if playing:
            num_seconds = current_surface.get_width()/150
            x = int(((playing-playing_start)/num_seconds) * current_surface.get_width())
            if x >= current_surface.get_width():
                playing = False
                playing_start = False
                last_x=0
                current_surface.fill([0,100,0,0.3], ((x - (x-last_x),0),(x-last_x,current_surface.get_height())), pygame.BLEND_RGBA_SUB)
                continue
            if x-last_x >=1:
                current_surface.fill([0,100,0,0.3], ((x - (x-last_x),0),(x-last_x,current_surface.get_height())), pygame.BLEND_RGBA_ADD)
            playing = time()
            last_x = x

        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                running = False
            if event.type == pygame.MOUSEBUTTONUP:
                pressing = False
                gesture_stack.append(current_gesture)
                current_gesture += 1
                last_mouse_pos = None
                grid_changed = True
            if event.type == pygame.MOUSEBUTTONDOWN:
                pressing = True
                last_mouse_pos = [-1,-1]
            if event.type == pygame.MOUSEMOTION and pressing:

                if last_mouse_pos != [-1,-1]:
                    points = bresenham(last_mouse_pos[0], last_mouse_pos[1], event.pos[0], event.pos[1])
                else:
                    points = [[event.pos[0], event.pos[1]]]
                for p in list(points)[0:-1]:
                    current_brush.draw(p[0], p[1])
                last_mouse_pos = event.pos
            if event.type == pygame.KEYDOWN:
                num = -1
                try:
                    num = int(event.unicode)
                except:
                    pass
                if num > -1:
                    current_brush = brushes[num-1]
                elif event.key == 122 and pygame.key.get_mods() & pygame.KMOD_CTRL:
                    current_spec_unit.remove_last_gesture()

                elif event.unicode == "+":
                    last_grid_len = grid_len
                    grid_len += 1
                    grid_changed = True

                elif event.unicode == "-":
                    last_grid_len = grid_len
                    grid_len -= 1
                    grid_changed = True

                elif event.unicode == " ":
                    if playing:
                        continue
                    il = Image.fromarray(np.asarray(canvas_l, np.dtype('uint8')))
                    il.convert("RGB").save("out_l.bmp")
                    ir = Image.fromarray(np.asarray(canvas_r, np.dtype('uint8')))
                    ir.convert("RGB").save("out_r.bmp")
                    delegator.run("arss out_l.bmp o_l.wav --sine --min-freq 27.5 --max-freq 20000 -r 44000 -f 16 --pps 150")
                    delegator.run("arss out_r.bmp o_r.wav --sine --min-freq 27.5 --max-freq 20000 -r 44000 -f 16 --pps 150")
                    delegator.run("sox o_l.wav o_r.wav --channels 2 --combine merge o.wav")
                    delegator.run("play o.wav", block=False)
                    playing_start = time()
                    playing = time()
        pygame.display.update()

if __name__=="__main__":
    main()

max_stars_repo_name: adelrune/freqsketch

max_issues_repo_name: adelrune/freqsketch

max_forks_repo_name: adelrune/freqsketch


text:
#! /usr/bin/env python
from numpy import *
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst
import os
from log_it import *

#############################################################################
# Created on 3/24/2014 by Gail Schmidt, USGS/EROS LSRD Project
# Created Python script to open and read the input single band ENVI file
# to obtain particular attributes and data from the file.
#
# History:
#   Updated on 3/26/2014 by Gail Schmidt, USGS/EROS
#       Removed metadata reads from the old TIF files that were not being
#       used
############################################################################
class ENVI_Scene:
    """Class to work with the single band ENVI files.
    """

    filename = ""                 # name of input image file
    dX = 0                        # pixel size in X direction
    dY = 0                        # pixel size in Y direction
    NRow = 0                      # number of lines in the scene
    NCol = 0                      # number of samples in the scene
    NoData = 0                    # fill value
    NorthBoundingCoordinate = 0   # northern bounding coord
    SouthBoundingCoordinate = 0   # southern bounding coord
    WestBoundingCoordinate = 0    # western bounding coord
    EastBoundingCoordinate = 0    # eastern bounding coord

    # datasets created by gdal.Open
    dataset = None

    # bands
    band1 = None

    def __init__(self, fname, log_handler=None):
        """Class constructor to open and read metadata of the ENVI file.
        Description: Class constructor verifies the input file exists, then
            opens it, reads the metadata, and establishes pointers to the
            single band.

        History:
          Created on 3/24/2014 by Gail Schmidt, USGS/EROS LSRD Project

        Args:
          fname - name of the input ENVI file to be processed
          log_handler - open log file for logging or None for stdout

        Returns:
            None - error opening or reading the file via GDAL
            Object - successful processing
        """

        # make sure the file exists
        if not os.path.exists(fname):
            msg = 'Input file does not exist: ' + fname
            logIt (msg, log_handler)
            return None

        # store the filename for this class and open it
        self.filename = fname
        self.dataset = gdal.Open (fname)
        if self.dataset is None:
            msg = 'GDAL could not open input file: ' + fname
            logIt (msg, log_handler)
            return None

        # create connections to the band
        self.band1 = self.dataset.GetRasterBand(1)

        # get scene and coordinate information
        gt = self.dataset.GetGeoTransform()
        self.NCol = self.dataset.RasterXSize
        self.dX = gt[1]
        self.WestBoundingCoordinate = gt[0]
        self.EastBoundingCoordinate = self.WestBoundingCoordinate +  \
            (self.NCol * self.dX)

        self.NRow = self.dataset.RasterYSize
        self.dY = gt[5]
        self.NorthBoundingCoordinate = gt[3]
        self.SouthBoundingCoordinate = self.NorthBoundingCoordinate +  \
            (self.NRow * self.dY)

        self.NoData = self.band1.GetNoDataValue()

    def __del__(self):
        """Class constructor to clear the bands.
        Description: class destructor cleans up all the sub dataset and band
            pointers.

        History:
          Created in 2013 by Jodi Riegle and Todd Hawbaker, USGS Rocky Mountain
              Geographic Science Center

        Args: None

        Returns: Nothing
        """

        self.dataset = None
        self.band1 = None

    def xy2ij(self, x, y):
        """Converts projection space to pixel space.
        Description: xy2ij converts projection coordinates (x,y) to pixel
            space (i,j)

        History:
          Created in 2013 by Jodi Riegle and Todd Hawbaker, USGS Rocky Mountain
              Geographic Science Center

        Args:
          x - x projection coordinate
          y - y projection coordinate

        Returns:
            (col, row) - i,j pixel referring to x,y coordinate
        """

        col = int((x - self.WestBoundingCoordinate) / self.dX)
        row = int((y - self.NorthBoundingCoordinate) / self.dY)
        return([col, row])

    def ij2xy(self, col, row):
        """Converts pixel space to projection space.
        Description: ij2xy converts pixel space values (i,j) to projection
            coordinates (x,y)

        History:
          Created in 2013 by Jodi Riegle and Todd Hawbaker, USGS Rocky Mountain
              Geographic Science Center

        Args:
          row - row (line) pixel space coordinate
          col - column (sample) pixel space coordinate

        Returns:
            (x, y) - x,y projection coord referring to pixel i,j
        """

        x = self.WestBoundingCoordinate + (col * self.dX)
        y = self.NorthBoundingCoordinate + (row * self.dY)
        return([x,y])

    def getRowOfBandValues(self, y):
        """Reads the specified line at specified projection y coordinate.
        Description: getRowOfBandValues reads a row of band values from band1
            at the associated row for the projection y coordinate.  The
            projection y coordinate is converted to the actual row/line in the
            current scene, and then that row is read (reading all the samples
            in the line) and returned.

        History:
          Created in 2013 by Jodi Riegle and Todd Hawbaker, USGS Rocky Mountain
              Geographic Science Center
          Updated on 5/1/2013 by Gail Schmidt, USGS/EROS LSRD Project
              Modified to validate the y projection coordinate falls within
              the current scene.

        Args:
          y - y projection coordinate

        Returns:
            None - if pixel location does not fall within the coordinates of
                this scene
            Array - associated line of data y location as band1
        """

        if (y < self.NorthBoundingCoordinate) and  \
           (y > self.SouthBoundingCoordinate):
            # get the i,j pixel (line) for the y projection coordinates
            ij = self.xy2ij(0, y)

            # read the line from band1 in the the scene
            x1 = self.band1.ReadAsArray(0, ij[1], self.NCol, 1)[0,]

            return (x1)
        else:
            return (None)

    def getBandValues(self):
        """Reads the specified band of data.
        Description: getBandValues reads an entire band of data for the
            band 1.

        History:
          Created on 5/1/2013 by Gail Schmidt, USGS/EROS LSRD Project

        Args: None

        Returns:
            Array - associated data for band 1
        """

        x = zeros ((self.NRow, self.NCol))
        return (self.band1.ReadAsArray())

#### end ENVI_Scene class ####

max_stars_repo_name: heyyy/burned-area

max_issues_repo_name: heyyy/burned-area

max_forks_repo_name: heyyy/burned-area

reason: from numpy


text:
from dataclasses import dataclass, field
from typing import List

import numpy
from marshmallow import Schema, fields, post_load

class Distribution:
    """
    Base type for all distribution classes.
    """

    def dumpsPacked(self):
        raise Exception(f"'dumpsPacked' has not been implemented in {str(type(self))}")

class NormalSchema(Schema):
    _type = fields.String(default="IncrementalInference.PackedNormal")
    mu = fields.Float()
    sigma = fields.Float()

    @post_load
    def marshal(self, data, **kwargs):
        return Normal(**data)

    class Meta:
        ordered = True

@dataclass
class Normal(Distribution):
    mu: float
    sigma: float
    _type: str = "IncrementalInference.PackedNormal"

    def dumps(self):
        return NormalSchema().dumps(self)

    def dump(self):
        return NormalSchema().dump(self)

    @staticmethod
    def load(data):
        return NormalSchema().load(data)

class RayleighSchema(Schema):
    _type = fields.String(default="IncrementalInference.PackedRayleigh")
    sigma = fields.Float()

    @post_load
    def marshal(self, data, **kwargs):
        return Normal(**data)

    class Meta:
        ordered = True

@dataclass
class Rayleigh(Distribution):
    sigma: float
    _type: str = "IncrementalInference.PackedRayleigh"

    def dumps(self):
        return RayleighSchema().dumps(self)

    def dump(self):
        return RayleighSchema().dump(self)

    @staticmethod
    def load(data):
        return RayleighSchema().load(data)

class FullNormalSchema(Schema):
    _type = fields.String(default="IncrementalInference.PackedFullNormal")
    mu = fields.List(fields.Float)
    cov = fields.Method("get_cov", "set_cov", required=True)

    def get_cov(self, obj):
        return obj.cov.flatten().tolist()

    def set_cov(self, obj):
        raise Exception("Deserialization not supported yet.")  # obj.cov.reshape(3,3)?

    class Meta:
        ordered = True

    @post_load
    def marshal(self, data, **kwargs):
        return FullNormal(**data)

@dataclass
class FullNormal(Distribution):
    # TODO: Remove the default initializers.
    mu: numpy.ndarray = field(default_factory=lambda: numpy.zeros(3))
    cov: numpy.ndarray = field(default_factory=lambda: numpy.diag([0.1, 0.1, 0.1]))
    _type: str = "IncrementalInference.PackedFullNormal"

    def dumps(self):
        return FullNormalSchema().dumps(self)

    def dump(self):
        return FullNormalSchema().dump(self)

    @staticmethod
    def load(data):
        return FullNormalSchema().load(data)

@dataclass
class Uniform(Distribution):
    a: float
    b: float

    def dumps(self):
        return UniformSchema().dumps(self)

    def dump(self):
        return UniformSchema().dump(self)

    @staticmethod
    def load(data):
        return UniformSchema().load(data)

class UniformSchema(Schema):
    _type = fields.String(default="IncrementalInference.PackedUniform")
    a = fields.Float()
    b = fields.Float()
    PackedSamplableTypeJSON = fields.String(
        default="IncrementalInference.PackedUniform"
    )

    class Meta:
        ordered = True

@dataclass
class Categorical(Distribution):
    p: List[float]

    def dumps(self):
        return CategoricalSchema().dumps(self)

    def dump(self):
        return CategoricalSchema().dump(self)

    @staticmethod
    def load(data):
        return CategoricalSchema().load(data)

class CategoricalSchema(Schema):
    _type = fields.String(default="IncrementalInference.PackedCategorical")
    p = fields.List(fields.Float)

    class Meta:
        ordered = True


text:
import numpy as np
import pandas as pd
import seaborn as sns
import xlsxwriter
import psycopg2
from sqlalchemy import create_engine

engine = create_engine('postgresql://#:############@#####.###:5432/######',echo=False)
n = int(input('колличество строк: '))

df = pd.DataFrame({
				    #'Номер пробы'     : [],
				    #'Номер скважины'  : [],
				    #'гл. *** м.'      : [],
				   	#'кальций-ион'     : [],
				   	#'магний-ион'      : []
					'q1'      : [],
				    'q2'      : [],
				    'q3'      : [],
				   	'q4'      : [],
				   	'q5'      : [],

	}
	)

for i in range(n):
	#k1 = int  (input('Номер пробы: '     ))
	#k2 = int  (input('Номер скважины: '  ))
	#k3 = float(input('гл. *** м.: '      ))
	#k4 = int  (input('кальций-ион: '     ))
	#k5 = int  (input('магний-ион: '      ))
	k1 = int  (input('Номер пробы: '     ))
	k2 = int  (input('Номер скважины: '  ))
	k3 = float(input('гл. *** м.: '      ))
	k4 = int  (input('кальций-ион: '     ))
	k5 = int  (input('магний-ион: '      ))

	df.loc[i+1] = k1, k2, k3, k4, k5

df['q6'] = df['q3'] - df['q4']

print(df)

df.to_sql(name='test', con=engine, if_exists = 'replace', index=False)

writer = pd.ExcelWriter('pandas_simple.xlsx', engine='xlsxwriter')

# Convert the dataframe to an XlsxWriter Excel object.
df.to_excel(writer, sheet_name='Sheet1')

# Close the Pandas Excel writer and output the Excel file.
writer.save()

print(df)

df.to_sql(name='test', con=engine, if_exists = 'replace', index=False)

writer = pd.ExcelWriter('pandas_simple.xlsx', engine='xlsxwriter')

# Convert the dataframe to an XlsxWriter Excel object.
df.to_excel(writer, sheet_name='Sheet1')

# Close the Pandas Excel writer and output the Excel file.
writer.save()

max_stars_repo_name: FoxProklya/BeginPandas

max_issues_repo_name: FoxProklya/BeginPandas

max_forks_repo_name: FoxProklya/BeginPandas


text:
"""
Module containing class representing a 1D distribution that jumps from integers
to integers. Its PMF is represented by: $$f(x,y) = \\begin{cases}\
1-p & y=x \\\\ p/2 & y=x-1\\wedge x\\ne \\text{max} \\\\\
p & y=x-1\\wedge x=\\text{max}\\\\ p/2 & y=x+1\\wedge x\\ne \\text{min}\\\\\
p & y=x+1\\wedge x=\\text{min}\\\\ 0 & \\text{otherwise}\\end{cases}$$

**File**: $DISTPY/distpy/jumping/AdjacencyJumpingDistribution.py
**Author**: Keith Tauscher
**Date**: 3 Jul 2021
"""
import numpy as np
from ..util import int_types, numerical_types
from .JumpingDistribution import JumpingDistribution

log2 = np.log(2)

class AdjacencyJumpingDistribution(JumpingDistribution):
    """
    Class representing a 1D distribution that jumps from integers to integers.
    Its PMF is represented by: $$f(x,y) = \\begin{cases} 1-p & y=x \\\\\
    p/2 & y=x-1\\wedge x\\ne \\text{max} \\\\ p & y=x-1\\wedge x=\\text{max}\
    \\\\ p/2 & y=x+1\\wedge x\\ne \\text{min}\\\\ p &\
    y=x+1\\wedge x=\\text{min}\\\\ 0 & \\text{otherwise}\\end{cases}$$
    """
    def __init__(self, jumping_probability=0.5, minimum=None, maximum=None):
        """
        Initializes an `AdjacencyJumpingDistribution` with the given jumping
        probability (and extrema, if applicable).

        Parameters
        ----------
        jumping_probability : float
            number between 0 and 1 (exclusive) describing the probability with
            which the destination is different from the source.
        minimum : int or None
            - if None, no minimum is used
            - if an int, the minimum integer ever drawn by the distribution
        maximum : int or None
            - if None, no maximum is used
            - if an int, the maximum integer ever drawn by the distribution
        """
        self.jumping_probability = jumping_probability
        self.minimum = minimum
        self.maximum = maximum

    @property
    def jumping_probability(self):
        """
        The probability, \\(0<p<1\\) with which the destination is different
        than the source.
        """
        if not hasattr(self, '_jumping_probability'):
            raise AttributeError("jumping_probability referenced before it " +\
                "was set.")
        return self._jumping_probability

    @jumping_probability.setter
    def jumping_probability(self, value):
        """
        Setter for `AdjacencyJumpingDistribution.jumping_probability`.

        Parameters
        ----------
        value : float
            number greater than 0 and less than 1
        """
        if type(value) in numerical_types:
            if (value > 0) and (value < 1):
                self._jumping_probability = value
            else:
                raise ValueError("jumping_probability, jp, doesn't satisfy " +\
                    "0<jp<1.")
        else:
            raise TypeError("jumping_probability was set to a non-number.")

    @property
    def log_jumping_probability(self):
        """
        The natural logarithm of the jumping probability, \\(\\ln{p}\\).
        """
        if not hasattr(self, '_log_jumping_probability'):
            self._log_jumping_probability = np.log(self.jumping_probability)
        return self._log_jumping_probability

    @property
    def log_of_complement_of_jumping_probability(self):
        """
        The natural logarithm of the complement of the jumping probability,
        \\(\\ln{(1-p)}\\).
        """
        if not hasattr(self, '_log_of_complement_of_jumping_probability'):
            self._log_of_complement_of_jumping_probability =\
                np.log(1 - self.jumping_probability)
        return self._log_of_complement_of_jumping_probability

    @property
    def minimum(self):
        """
        Either None (if this distribution should be able to jump towards
        negative infinity) or the minimum integer this distribution should ever
        draw.
        """
        if not hasattr(self, '_minimum'):
            raise AttributeError("minimum referenced before it was set.")
        return self._minimum

    @minimum.setter
    def minimum(self, value):
        """
        Setter for `AdjacencyJumpingDistribution.minimum`.

        Parameters
        ----------
        value : int or None
            - if None, no minimum is used
            - if an int, the minimum integer ever drawn by the distribution
        """
        if (type(value) is type(None)) or (type(value) in int_types):
            self._minimum = value
        else:
            raise TypeError("minimum was set to a non-int.")

    @property
    def maximum(self):
        """
        Either None (if this distribution should be able to jump towards
        infinity) or the maximum integer this distribution should ever draw.
        """
        if not hasattr(self, '_maximum'):
            raise AttributeError("maximum referenced before it was set.")
        return self._maximum

    @maximum.setter
    def maximum(self, value):
        """
        Setter for `AdjacencyJumpingDistribution.maximum`.

        Parameters
        ----------
        value : int or None
            - if None, no maximum is used
            - if an int, the maximum integer ever drawn by the distribution
        """
        if type(value) is type(None):
            self._maximum = value
        elif type(value) in int_types:
            if (type(self.minimum) is type(None)) or (value > self.minimum):
                self._maximum = value
            else:
                raise ValueError("maximum wasn't greater than minimum.")
        else:
            raise TypeError("maximum was set to a non-int.")

    def draw_single_value(self, source, random=np.random):
        """
        Draws a single value from this distribution.

        Parameters
        ----------
        source : int
            single starting point
        random : numpy.random.RandomState
            the random number generator to use (default: `numpy.random`)

        Returns
        -------
        destination : int
            single int satisfying \\((y-x)\\in\\{-1,0,1\\}\\)
        """
        uniform = random.rand()
        if uniform < self.jumping_probability:
            if source == self.minimum:
                return (source + 1)
            elif source == self.maximum:
                return (source - 1)
            elif uniform < (self.jumping_probability / 2.):
                return (source - 1)
            else:
                return (source + 1)
        else:
            return source

    def draw_shaped_values(self, source, shape, random=np.random):
        """
        Draws arbitrary shape of random values given the source point.

        Parameters
        ----------
        source : int
            a single integer number from which to jump
        shape : tuple
            tuple of ints \\((n_1,n_2,\\ldots,n_k)\\),
            \\(\\prod_{m=1}^kn_m\\) destinations are returned as a
            `numpy.ndarray` of shape \\((n_1,n_2,\\ldots,n_k)\\)
        random : numpy.random.RandomState
            the random number generator to use (default: `numpy.random`)

        Returns
        -------
        destinations : int or numpy.ndarray
            drawn destination(s)
        """
        uniform = random.rand(*shape)
        jumps = np.where(uniform < self.jumping_probability, 1, 0)
        if source == self.minimum:
            pass
        elif source == self.maximum:
            jumps = -jumps
        else:
            jumps[np.where(uniform < self.jumping_probability / 2.)[0]] = -1
        return source + jumps

    def draw(self, source, shape=None, random=np.random):
        """
        Draws random value(s) given the source point.

        Parameters
        ----------
        source : int
            a single integer number from which to jump
        shape : None or int or tuple
            - if None, a single int destination is returned
            - if int \\(n\\), \\(n\\) destinations are returned in a 1D
            `numpy.ndarray` of length \\(n\\)
            - if tuple of ints \\((n_1,n_2,\\ldots,n_k)\\),
            \\(\\prod_{m=1}^kn_m\\) destinations are returned as a
            `numpy.ndarray` of shape \\((n_1,n_2,\\ldots,n_k)\\)
        random : numpy.random.RandomState
            the random number generator to use (default: `numpy.random`)

        Returns
        -------
        destinations : int or numpy.ndarray
            drawn destination(s)
        """
        if type(shape) is type(None):
            return self.draw_single_value(source, random=random)
        if type(shape) in int_types:
            shape = (shape,)
        return self.draw_shaped_values(source, shape, random=random)

    def log_value(self, source, destination):
        """
        Computes the log-PDF of jumping from `source` to `destination`. It must
        be implemented by all subclasses.

        Parameters
        ----------
        source : int
            source point
        destination : int
            destination point

        Returns
        -------
        log_pdf : float
            if the distribution is \\(f(x,y)=\\text{Pr}[y|x]\\), `source` is
            \\(x\\) and `destination` is \\(y\\), then `log_pdf` is given by
            \\(\\ln{f(x,y)}\\)
        """
        displacement = destination - source
        if displacement == 0:
            return self.log_of_complement_of_jumping_probability
        elif displacement == 1:
            return_value = self.log_jumping_probability
            if source != self.minimum:
                return_value -= log2
            return return_value
        elif displacement == -1:
            return_value = self.log_jumping_probability
            if source != self.maximum:
                return_value -= log2
            return return_value
        else:
            return -np.inf

    def log_value_difference(self, source, destination):
        """
        Computes the difference in the log-PDF of jumping from `source` to
        `destination` and of jumping from `destination` to `source`.

        Parameters
        ----------
        source : int
            source point
        destination : int
            destination point

        Returns
        -------
        log_pdf_difference : float
            if the distribution is \\(f(x,y)=\\text{Pr}[y|x]\\), `source` is
            \\(x\\) and `destination` is \\(y\\), then `log_pdf_difference` is
            given by \\(\\ln{f(x,y)}-\\ln{f(y,x)}\\)
        """
        displacement = destination - source
        if displacement == 0:
            return 0.
        elif displacement == 1:
            source_is_minimum = (source == self.minimum)
            destination_is_maximum = (destination == self.maximum)
            if source_is_minimum == destination_is_maximum:
                return 0.
            elif source_is_minimum:
                return log2
            else:
                return -log2
        elif displacement == -1:
            source_is_maximum = (source == self.maximum)
            destination_is_minimum = (destination == self.minimum)
            if source_is_maximum == destination_is_minimum:
                return 0.
            elif source_is_maximum:
                return log2
            else:
                return -log2
        else:
            raise ValueError("source and destination could not connected " +\
                "by only a single jump.")

    @property
    def numparams(self):
        """
        The integer number of parameters described by this distribution. Since
        this distribution is univariate, this property is 1.
        """
        return 1

    def __eq__(self, other):
        """
        Tests for equality between this jumping distribution and other.

        Parameters
        ----------
        other : object
            object with which to check for equality

        Returns
        -------
        result : bool
            True if and only if `other` is an `AdjacencyJumpingDistribution`
            with the same `AdjacencyJumpingDistribution.minimum`,
            `AdjacencyJumpingDistribution.maximum`, and
            `AdjacencyJumpingDistribution.jumping_probability`
        """
        if isinstance(other, AdjacencyJumpingDistribution):
            jumping_probabilities_equal = np.isclose(self.jumping_probability,\
                other.jumping_probability, atol=1e-6)
            minima_equal = (self.minimum == other.minimum)
            maxima_equal = (self.maximum == other.maximum)
            return (jumping_probabilities_equal and (minima_equal and\
                maxima_equal))
        else:
            return False

    @property
    def is_discrete(self):
        """
        Boolean describing whether this `AdjacencyJumpingDistribution`. Since
        it exists, on a grid, this is always True.
        """
        return True

    def fill_hdf5_group(self, group):
        """
        Fills the given hdf5 file group with information about this
        distribution.

        Parameters
        ----------
        group : h5py.Group
            hdf5 file group to fill with information about this distribution
        """
        group.attrs['class'] = 'AdjacencyJumpingDistribution'
        group.attrs['jumping_probability'] = self.jumping_probability
        if type(self.minimum) is not type(None):
            group.attrs['minimum'] = self.minimum
        if type(self.maximum) is not type(None):
            group.attrs['maximum'] = self.maximum

    @staticmethod
    def load_from_hdf5_group(group):
        """
        Loads an `AdjacencyJumpingDistribution` from the given hdf5 file group.

        Parameters
        ----------
        group : h5py.Group
            the same hdf5 file group which
            `AdjacencyJumpingDistribution.fill_hdf5_group` was called on

        Returns
        -------
        loaded : `AdjacencyJumpingDistribution`
            an `AdjacencyJumpingDistribution` object created from the
            information in the given group
        """
        try:
            assert group.attrs['class'] == 'AdjacencyJumpingDistribution'
        except:
            raise ValueError("The given group does not seem to contain an " +\
                "AdjacencyJumpingDistribution.")
        jumping_probability = group.attrs['jumping_probability']
        if 'minimum' in group.attrs:
            minimum = group.attrs['minimum']
        else:
            minimum = None
        if 'maximum' in group.attrs:
            maximum = group.attrs['maximum']
        else:
            maximum = None
        return\
            AdjacencyJumpingDistribution(jumping_probability, minimum, maximum)

max_stars_repo_name: CU-NESS/distpy

max_issues_repo_name: CU-NESS/distpy

max_forks_repo_name: CU-NESS/distpy


text:
from numpy import linalg, zeros, ones, hstack, asarray
import itertools

def basis_vector(n, i):
    """ Return an array like [0, 0, ..., 1, ..., 0, 0]

    >>> from multipolyfit.core import basis_vector
    >>> basis_vector(3, 1)
    array([0, 1, 0])
    >>> basis_vector(5, 4)
    array([0, 0, 0, 0, 1])
    """
    x = zeros(n, dtype=int)
    x[i] = 1
    return x

def as_tall(x):
    """ Turns a row vector into a column vector """
    return x.reshape(x.shape + (1,))

def multipolyfit(xs, y, deg, full=False, model_out=False, powers_out=False):
    """
    Least squares multivariate polynomial fit

    Fit a polynomial like ``y = a**2 + 3a - 2ab + 4b**2 - 1``
    with many covariates a, b, c, ...

    Parameters
    ----------

    xs : array_like, shape (M, k)
         x-coordinates of the k covariates over the M sample points
    y :  array_like, shape(M,)
         y-coordinates of the sample points.
    deg : int
         Degree o fthe fitting polynomial
    model_out : bool (defaults to True)
         If True return a callable function
         If False return an array of coefficients
    powers_out : bool (defaults to False)
         Returns the meaning of each of the coefficients in the form of an
         iterator that gives the powers over the inputs and 1
         For example if xs corresponds to the covariates a,b,c then the array
         [1, 2, 1, 0] corresponds to 1**1 * a**2 * b**1 * c**0

    See Also
    --------
        numpy.polyfit

    """
    y = asarray(y).squeeze()
    rows = y.shape[0]
    xs = asarray(xs)
    num_covariates = xs.shape[1]
    xs = hstack((ones((xs.shape[0], 1), dtype=xs.dtype), xs))

    generators = [basis_vector(num_covariates + 1, i)
                  for i in range(num_covariates + 1)]

    # All combinations of degrees
    powers = map(sum, itertools.combinations_with_replacement(generators, deg))

    # Raise data to specified degree pattern, stack in order
    A = hstack(asarray([as_tall((xs ** p).prod(1)) for p in powers]))

    beta = linalg.lstsq(A, y)[0]

    if model_out:
        return mk_model(beta, powers)

    if powers_out:
        return beta, powers
    return beta

def mk_model(beta, powers):
    """ Create a callable python function out of beta/powers from multipolyfit

    This function is callable from within multipolyfit using the model_out flag
    """

    # Create a function that takes in many x values
    # and returns an approximate y value
    def model(*args):
        num_covariates = len(powers[0]) - 1
        if len(args) != (num_covariates):
            raise ValueError("Expected %d inputs" % num_covariates)
        xs = asarray((1,) + args)
        return sum([coeff * (xs ** p).prod()
                    for p, coeff in zip(powers, beta)])

    return model

def mk_sympy_function(beta, powers):
    from sympy import symbols, Add, Mul, S
    num_covariates = len(powers[0]) - 1
    xs = (S.One,) + symbols('x0:%d' % num_covariates)
    return Add(*[coeff * Mul(*[x ** deg for x, deg in zip(xs, power)])
                 for power, coeff in zip(powers, beta)])

reason: from numpy,from sympy


text:
#need to install "scikit-image" first
import originpro as op
import numpy as np
from skimage.util import invert

fn = op.path('e') + r'Samples\Image Processing and Analysis\Flower.jpg'
iw=op.new_image()
iw.from_file(fn)
iw.lname='original image'
#make a copy to put result
iw2=iw.duplicate()
iw2.lname='inverted'
img = iw2.to_np()
inv = invert(img)
iw2.from_np(inv)
#put into new empty image window
iw3=op.new_image()
#need to first setup the image window as 3 channels and not multi-frames
iw3.setup(3,False)
iw3.lname='invert back'
inv = invert(inv)
iw3.from_np(inv)

max_stars_repo_path: 6. Images/Inverting an Image with skimage.py

max_stars_repo_name: originlab/Python

max_issues_repo_name: originlab/Python

max_forks_repo_name: originlab/Python


text:
import requests
from unidecode import unidecode
import pandas as pd
import numpy as np

# Carregando dataset
arquivo = "santos_rent.csv"
santos = pd.read_csv(arquivo)
arquivo2 = "santos_local.csv"
end = pd.read_csv(arquivo2)

santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Primeiro de Maio','1 de Maio'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('º de Maio','1 de Maio'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Rua Engenheiro Manoel Ferramenta Júnior','Manoel Ferramenta Júnior'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Matias','Mathias'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Doutor Tolentino Filgueiras','Tolentino Filgueiras'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Doutor Galeão Carvalhal','Galeão Carvalhal'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace(' xx',''))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Costa x','Costa'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('J. Menino','José Menino'))
santos['endereco'] = santos['endereco'].apply(lambda x: x.replace('Doutor Cunha','Cunha'))

data = santos.merge(end, on='endereco')
data.head()

data = data.drop(['Unnamed: 0_x','Unnamed: 0_y'] , axis = 1)
data

data.isna().sum()

data = data.dropna(how = 'any', subset=['rooms', 'bathroom'])

data.isna().sum()

import requests
from time import sleep

def get_geo(cep):
    try:
        url = "https://www.cepaberto.com/api/v3/cep?cep="
        # O seu token está visível apenas pra você
        headers = {'Authorization': 'Token token=461f19961c9388e6f3383a8504e9fd25'}
        response = requests.get(url+cep, headers=headers)
        sleep(0.5)
        x = response.json()['latitude']
        y = response.json()['longitude']
        print(x + ',' + y)
        return x + ',' + y
    except:
        print(cep)
        return np.nan

end['CEP'] = end['CEP'].astype(str)
end['CEP'] = end['CEP'].apply(lambda x: x.replace('-',''))

end['coord'] = end['CEP'].apply(lambda x: get_geo(x))

# Novo dataframe contendo os valores da geolocalização divididos em 2 colunas
new = end["coord"].str.split(",", n = 1, expand = True)

# Inserindo os valores de latitude e longitude no dataframe "end"
end["lat"]= new[0]
end["long"]= new[1]

# Retirando a coluna antigas
end.drop(columns =["coord"], inplace = True)

end.to_csv('santos_geo.csv')