text:
"""FactorAnalysis class for Metran in Pastas."""

from logging import getLogger

import numpy as np
import scipy.optimize as scopt
from pastas.utils import initialize_logger

logger = getLogger(__name__)
initialize_logger(logger)

class FactorAnalysis:
    """Class to perform a factor analysis for the Pastas Metran model.

    Parameters
    ----------
    maxfactors : int, optional.
        maximum number of factors to select. The default is None.

    Examples
    --------
    A minimal working example of the FactorAnalysis class is shown below:

    >>> fa = FactorAnalysis()
    >>> factors = fa.solve(oseries)
    """

    def __init__(self, maxfactors=None):
        self.maxfactors = maxfactors

    def get_eigval_weight(self):
        """Method to get the relative weight of each eigenvalue.

        Returns
        -------
        numpy.ndarray
            All eigenvalues as a fraction of the sum of eigenvalues.
        """
        return self.eigval / np.sum(self.eigval)

    def solve(self, oseries):
        """Method to perform factor analysis.

        Factor analysis is based on the minres algorithm.
        The number of eigenvalues is determined by MAP test.
        If more than one eigenvalue is used,
        the factors are rotated using orthogonal rotation.

        Parameters
        ----------
        oseries : pandas.DataFrame
            Object containing the time series. The
            series can be non-equidistant.

        Raises
        ------
        Exception
            If no proper factors can be derived from the series.

        Returns
        -------
        factors : numpy.ndarray
            Factor loadings.
        """
        correlation = self._get_correlations(oseries)
        self.eigval, eigvec = self._get_eigval(correlation)

        # Velicer's MAP test
        try:
            nfactors, _ = self._maptest(correlation,
                                        eigvec, self.eigval)
            msg = "Number of factors according to Velicer\'s MAP test: " \
                  + f"{nfactors}"
            logger.info(msg)
            if nfactors == 0:
                nfactors = sum(self.eigval > 1)
                msg = "Number of factors according to Kaiser criterion: " \
                      + f"{nfactors}"
                logger.info(msg)
            if self.maxfactors is not None:
                nfactors = min(nfactors, self.maxfactors)
        except Exception:
            nfactors = 0
        factors = self._minres(correlation, nfactors)

        if ((nfactors > 0) and (factors is not None)
                and (np.count_nonzero(factors) > 0)):
            # factors is not None and does not contain nonzero elements
            if nfactors > 1:
                # perform varimax rotation
                comm = np.zeros(factors.shape[0])
                for i in range(factors.shape[0]):
                    for j in range(nfactors):
                        comm[i] = comm[i] + factors[i, j] ** 2
                    factors[i, :] = factors[i, :] / np.sqrt(comm[i])
                factors = self._rotate(factors[:, :nfactors])
                for i in range(factors.shape[0]):
                    factors[i, :] = factors[i, :] * np.sqrt(comm[i])

            # swap sign if dominant sign is negative
            for j in range(factors.shape[1]):
                facsign = 0
                for i in range(factors.shape[0]):
                    facsign = facsign + factors[i, j]
                if facsign < 0:
                    for i in range(factors.shape[0]):
                        if np.sign(factors[i, j]) != 0:
                            factors[i, j] = -1. * factors[i, j]

            self.factors = np.atleast_2d(factors[:, :nfactors])

            self.fep = 100 * np.sum(self.get_eigval_weight()[:nfactors])

        else:
            msg = "No proper common factors could be derived from series."
            logger.warning(msg)
            self.factors = None

        return self.factors

    @staticmethod
    def _rotate(phi, gamma=1, maxiter=20, tol=1e-6):
        """Internal method to rotate factor loadings.

        Uses varimax, quartimax, equamax, or parsimax rotation.

        Parameters
        ----------
        phi : numpy.ndarray
            Eigenvectors to be rotated
        gamma : float, optional
            Coefficient for rotation. The default is 1.
            Varimax: gamma = 1.
            Quartimax: gamma = 0.
            Equamax: gamma = nfac/2.
            Parsimax: gamma = nvar(nfac - 1)/(nvar + nfac - 2).
        maxiter : integer, optional
            Maximum number of iterations. The default is 20.
        tol : float, optional
            Stop criterion. The default is 1e-6.

        Returns
        -------
        phi_rot : 2-dimensional array
            rotated eigenvectors

        References
        ----------
        Kaiser, H.F. (1958): The varimax criterion for analytic rotation in
        factor analysis. Psychometrika 23: 187–200.
        """
        p, k = phi.shape
        R = np.eye(k)
        d = 0
        for _ in range(maxiter):
            d_old = d
            Lambda = np.dot(phi, R)
            u, s, vh = np.linalg.svd(
                np.dot(phi.T, np.asarray(Lambda) ** 3 - (gamma / p)
                       * np.dot(Lambda, np.diag(
                           np.diag(np.dot(Lambda.T, Lambda))))))
            R = np.dot(u, vh)
            d = np.sum(s)
            if (d_old != 0) and (d / d_old < 1 + tol):
                break

        phi_rot = np.dot(phi, R)
        return phi_rot

    def _minres(self, s, nf, covar=False):
        """Internal method for estimating factor loadings.

        Uses the minimum residuals (minres) algorithm.

        Parameters
        ----------
        s : numpy.ndarray
            Correlation matrix
        nf : integer
            Number of factors
        covar : boolean
            True if S is covar

        Returns
        -------
        loadings : numpy.ndarray
            Estimated factor loadings
        """
        sorg = np.copy(s)
        try:
            ssmc = 1 - 1 / np.diag(np.linalg.inv(s))
            if (not(covar) and np.sum(ssmc) == nf) and (nf > 1):
                start = 0.5 * np.ones(nf, dtype=float)
            else:
                start = np.diag(s) - ssmc
        except:
            return

        bounds = list()
        for _ in range(len(start)):
            bounds.append((0.005, 1))

        res = scopt.minimize(self._minresfun, start, method='L-BFGS-B',
                             jac=self._minresgrad, bounds=bounds,
                             args=(s, nf))

        loadings = self._get_loadings(res.x, sorg, nf)

        return loadings

    @staticmethod
    def _maptest(cov, eigvec, eigval):
        """Internal method to run Velicer's MAP test.

        Determines the number of factors to be used. This method includes
        two variations of the MAP test: the orginal and the revised MAP test.

        Parameters
        ----------
        cov : numpy.ndarray
            Covariance matrix.
        eigvec : numpy.ndarray
            Matrix with columns eigenvectors associated with eigenvalues.
        eigval : numpy.ndarray
            Vector with eigenvalues in descending order.

        Returns
        -------
        nfacts : integer
            Number factors according to MAP test.
        nfacts4 : integer
            Number factors according to revised MAP test.

        References
        ----------
        The original MAP test:

        Velicer, W. F. (1976). Determining the number of components
        from the matrix of partial correlations. Psychometrika, 41, 321-327.

        The revised (2000) MAP test i.e., with the partial correlations
        raised to the 4rth power (rather than squared):

        Velicer, W. F., Eaton, C. A., and Fava, J. L. (2000). Construct
        explication through factor or component analysis: A review and
        evaluation of alternative procedures for determining the number
        of factors or components. Pp. 41-71 in R. D. Goffin and
        E. Helmes, eds., Problems and solutions in human assessment.
        Boston: Kluwer.
        """

        nvars = len(eigval)
        fm = np.array([np.arange(nvars, dtype=float),
                       np.arange(nvars, dtype=float)]).T
        np.put(fm, [0, 1], ((np.sum(np.sum(np.square(cov))) - nvars)
                            / (nvars * (nvars - 1))))
        fm4 = np.copy(fm)
        np.put(fm4, [0, 1],
               ((np.sum(np.sum(np.square(np.square(cov)))) - nvars)
                / (nvars * (nvars - 1))))

        for m in range(nvars - 1):
            biga = np.atleast_2d(eigvec[:, :m + 1])
            partcov = cov - np.dot(biga, biga.T)
            # exit function with nfacts=1 if diag partcov contains negatives
            if np.amin(np.diag(partcov)) < 0:
                return 1, 1
            d = np.diag((1 / np.sqrt(np.diag(partcov))))
            pr = np.dot(d, np.dot(partcov, d))
            np.put(fm, [m + 1, 1], ((np.sum(np.sum(np.square(pr))) - nvars)
                                    / (nvars * (nvars - 1))))
            np.put(fm4, [m + 1, 1], ((np.sum(np.sum(np.square(np.square(pr))))
                                      - nvars) / (nvars * (nvars - 1))))

        minfm = fm[0, 1]
        nfacts = 0
        minfm4 = fm4[0, 1]
        nfacts4 = 0
        for s in range(nvars):
            fm[s, 0] = s
            fm4[s, 0] = s
            if fm[s, 1] < minfm:
                minfm = fm[s, 1]
                nfacts = s
            if fm4[s, 1] < minfm4:
                minfm4 = fm4[s, 1]
                nfacts4 = s
        return nfacts, nfacts4

    @staticmethod
    def _minresfun(psi, s, nf):
        """Function to be minimized in minimum residuals (minres) algorithm.

        Parameters
        ----------
        psi : array
            Vector to be adjusted during optimization
        s : array
            Correlation matrix
        nf : integer
            Number of factors

        Returns
        -------
        obj : array
            objective function defined as sum of residuals
        """
        s2 = np.copy(s)
        np.fill_diagonal(s2, 1 - psi)
        eigval, eigvec = np.linalg.eigh(s2)
        eigval[eigval < np.MachAr().eps] = 100 * np.MachAr().eps
        if nf > 1:
            loadings = np.atleast_2d(np.dot(eigvec[:, :nf],
                                            np.diag(np.sqrt(eigval[:nf]))))
        else:
            loadings = eigvec[:, 0] * np.sqrt(eigval[0])
        model = np.dot(loadings, loadings.T)
        residual = np.square(s2 - model)
        np.fill_diagonal(residual, 0)

        return np.sum(residual)

    def _minresgrad(self, psi, s, nf):
        """Internal method to calculate jacobian of function.

        Jacobian to be minimized in minimum residuals (minres) algorithm.

        Parameters
        ----------
        psi : array
            Vector to be adjusted during optimization.
        s : array
            Correlation matrix.
        nf : integer
            Number of factors.

        Returns
        -------
        jac : array
            Jacobian of minresfun.
        """

        load = self._get_loadings(psi, s, nf)
        g = np.dot(load, load.T) + np.diag(psi) - s
        jac = np.diag(g) / np.square(psi)

        return jac

    @staticmethod
    def _get_loadings(psi, s, nf):
        """Internal method to estimate matrix of factor loadings.

        Based on minimum residuals (minres) algorithm.

        Parameters
        ----------
        psi : numpy.ndarray
            Communality estimate.
        s : numpy.ndarray
            Correlation matrix.
        nf : integer
            Number of factors.

        Returns
        -------
        load : npumy.ndarray
            Estimated factor loadings.
        """
        sc = np.diag(1 / np.sqrt(psi))
        sstar = np.dot(sc, np.dot(s, sc))
        eigval, eigvec = np.linalg.eig(sstar)
        L = eigvec[:, :nf]
        load = np.dot(L, np.diag(np.sqrt(np.maximum(
            np.subtract(eigval[:nf], 1), 0))))
        load = np.dot(np.diag(np.sqrt(psi)), load)
        return load

    @staticmethod
    def _get_correlations(oseries):
        """Internal method to calculate correlations for multivariate series.

        Parameters
        ----------
        oseries : pandas.DataFrame
            Multivariate series

        Returns
        -------
        corr : numpy.ndarray
            Correlation matrix
        """
        corr = np.array(oseries.corr())
        return corr

    @staticmethod
    def _get_eigval(correlation):
        """Internal method to get eigenvalues and eigenvectors.

        Get eigenvalues and eigenvectors based on correlation matrix.

        Parameters
        ----------
        correlation : numpy.ndarray
            Correlation matrix for which eigenvalues
            and eigenvectors need to be derived.

        Raises
        ------
        Exception
            If method results in complex eigenvalues and eigenvectors.

        Returns
        -------
        eigval : numpy.ndarray
            Vector with eigenvalues.
        eigvec : numpy.ndarray
            Matrix with eigenvectors.
        """
        # perform eigenvalue decomposition
        eigval, eigvec = np.linalg.eig(correlation)
        if isinstance(eigval[0], np.complex128):
            msg = "Serial correlation matrix has " + \
                  "complex eigenvalues and eigenvectors. " + \
                  "Factors cannot be estimated for these series."
            logger.error(msg)
            raise Exception(msg)
        # sort eigenvalues and eigenvectors
        evals_order = np.argsort(-eigval)
        eigval = eigval[evals_order]
        eigval[eigval < 0] = 0.
        eigvec = eigvec[:, evals_order]
        eigvec = np.atleast_2d(np.dot(eigvec, np.sqrt(np.diag(eigval))))
        return eigval, eigvec

max_stars_repo_name: pastas/metran

max_issues_repo_name: pastas/metran

max_forks_repo_name: pastas/metran

reason: import numpy,import scipy


text:
#! /usr/bin/env python

"""
Module with utilities to estimate the spectral type and gravity of an object
based on spectral indices.

.. [GOR03]
   | Gorlova et al. 2003
   | **Gravity Indicators in the Near-Infrared Spectra of Brown Dwarfs**
   | *The Astrophysical Journal, Volume 593, Issue 1, pp. 1074-1092*
   | `https://arxiv.org/abs/astro-ph/0305147
     <https://arxiv.org/abs/astro-ph/0305147>`_

.. [SLE04]
   | Slesnick et al. 2004
   | **The Spectroscopically Determined Substellar Mass Function of the Orion
     Nebula Cluster**
   | *The Astrophysical Journal, Volume 610, Issue 1, pp. 1045-1063*
   | `https://arxiv.org/abs/astro-ph/0404292
     <https://arxiv.org/abs/astro-ph/0404292>`_

.. [ALL07]
   | Allers et al. 2007
   | **Characterizing Young Brown Dwarfs Using Low-Resolution Near-Infrared
     Spectra**
   | *The Astrophysical Journal, Volume 657, Issue 1, pp. 511-520*
   | `https://arxiv.org/abs/astro-ph/0611408
     <https://arxiv.org/abs/astro-ph/0611408>`_
"""

__author__ = 'V. Christiaens'
__all__ = ['spectral_idx',
           'sp_idx_to_spt',
           'sp_idx_to_gravity',
           'digit_to_spt',
           'spt_to_digit']

import numpy as np
from .utils_spec import find_nearest

def spectral_idx(lbda, spec, band='H2O-1.5', spec_err=None, verbose=False):
    """
    Computes a spectral index. Implemented so far:
        - the Na 1.1 mu index [ALL07]
        - the H2O 1.3 mu index [GOR03]
        - the H2O 1.5 mu index [ALL07]
        - the H2O 2 index [SLE04]
        - the CO 2.3 index [GOR03].

    Parameters
    ----------
    lbda : numpy ndarray
        1d numpy array containing the wavelengths of the spectrum in microns.
    spec : numpy ndarray
        1d numpy array containing the measured flux (arbitrary units accepted).
    band: str, optional {'H2O-1.5', 'H2O-1.3', 'H2O-2', 'Na-1.1', 'CO-2.3'}
        Name of the band where the spectral index is defined (spectral feature
        + wavelength in mu)
    spec_err: numpy ndarray, optional
        1d numpy array containing the uncertainties on the measured flux
        (arbitrary units accepted). If provided the uncertainty on the spectral
        index will also be returned.
    verbose: bool, optional
        Whther to print more information.

    Returns
    -------
    index : float
        Value of the spectral index
    index_err: float
       [if spec_err is provided] Uncertainty on the spectral index.
    """

    if band == 'Na-1.1':
        lbda_num0 = 1.15
        lbda_num1 = 1.16
        lbda_den0 = 1.134
        lbda_den1 = 1.144
    elif band == 'H2O-1.3':
        lbda_num0 = 1.334
        lbda_num1 = 1.338
        lbda_den0 = 1.32
        lbda_den1 = 1.324
    elif band == 'H2O-1.5':
        lbda_num0 = 1.55
        lbda_num1 = 1.56
        lbda_den0 = 1.492
        lbda_den1 = 1.502
    elif band == 'H2O-2':
        lbda_num0 = 2.035
        lbda_num1 = 2.045
        lbda_den0 = 2.145
        lbda_den1 = 2.155
    elif band == 'CO-2.3':
        lbda_num0 = 2.296
        lbda_num1 = 2.298
        lbda_den0 = 2.286
        lbda_den1 = 2.288
    else:
        raise ValueError("Band not recognised")

    idx_num0 = find_nearest(lbda, lbda_num0)
    idx_num1 = find_nearest(lbda, lbda_num1)
    # set floor/ceil constraints if the indices are the same
    if idx_num0 == idx_num1:
        idx_num0 = find_nearest(lbda, lbda_num0, constraint='floor')
        idx_num1 = find_nearest(lbda, lbda_num0, constraint='ceil')
    idx_den0 = find_nearest(lbda, lbda_den0)
    idx_den1 = find_nearest(lbda, lbda_den1)
    # set floor/ceil constraints if the indices are the same
    if idx_den0 == idx_den1:
        idx_den0 = find_nearest(lbda, lbda_den0, constraint='floor')
        idx_den1 = find_nearest(lbda, lbda_den1, constraint='ceil')

    raise_flag = False
    # Check if enough spectral resolution
    err_msg = "Indices overlap for num and den: not enough spectral resolution"
    if idx_den0 == idx_num1:
        raise_flag = True
        if abs(lbda[idx_den0]-lbda_den0)<abs(lbda[idx_num1]-lbda_num1):
            idx_num1 -= 1
            if idx_num0 == idx_num1+1:
                raise ValueError(err_msg)
        else:
            idx_den0 += 1
            if idx_den0 == idx_den1+1:
                raise ValueError(err_msg)
    if idx_den1 == idx_num0:
        raise_flag = True
        if abs(lbda[idx_den1]-lbda_den1)<abs(lbda[idx_num0]-lbda_num0):
            idx_num0 += 1
            if idx_num0 == idx_num1+1:
                raise ValueError(err_msg)
        else:
            idx_den1 -= 1
            if idx_den0 == idx_den1+1:
                raise ValueError(err_msg)

    if verbose:
        if raise_flag:
            msg = "*WARNING*: spectral resolution is barely sufficient for "
            msg += "non-overlapping indices - results to be taken with caution!"
            print(msg)
        msg = "Closest (non-overlapping) wavelengths to {:.3f}, {:.3f}, {:.3f} and"
        msg+= " {:.3f} mu correspond to {:.3f}, {:.3f}, {:.3f} and {:.3f} mu resp."
        print(msg.format(lbda_num0, lbda_num1, lbda_den0, lbda_den1,
                         lbda[idx_num0], lbda[idx_num1], lbda[idx_den0],
                         lbda[idx_den1]))

    num = np.mean(spec[idx_num0:idx_num1+1])
    den = np.mean(spec[idx_den0:idx_den1+1])
    index = num/den

    if spec_err is None:
        return index
    else:
        num_err = np.sum(np.power(spec_err[idx_num0:idx_num1+1]/((idx_num1+1-idx_num0)*den),2))
        den_err = np.sum(np.power(spec_err[idx_num0:idx_num1+1]*(idx_den1+1-idx_den0)*num/(den**2),2))
        index_err = np.sqrt(num_err**2+den_err**2)
        return index, index_err

def sp_idx_to_spt(idx, name='H2O-1.5', idx_err=None, young=True):
    """
    Estimates a spectral type from a spectral index. Implemented so far:
        - the H2O 1.3 mu index [GOR03]
        - the H2O 1.5 mu index [ALL07]
        - the H2O-2 index [SLE04]

    Note on scale of SpT: 0 = M0, 10 = L0.

    Parameters
    ----------
    idx : float
        Value of spectral index
    name: str, optional {'H2O-1.3', 'H2O-1.5, 'H2O-2'}
        The name of the spectral index.
    idx_err: float, optional
        Uncertainty on the spectral index value
    young: bool, opt
        Whether the object is likely young (only used for 'H2O-1.3' index, which
        is gravity dependent)

    Returns
    -------
    spt : float
        Value of the spectral type
    spt_err: float
       [if idx_err is provided] Uncertainty on the spectral type.

    """

    if name == 'H2O-1.3':
        if young:
            a = 10.99
            da = 6.06
            b = 17.45
            db = 7.68
        else:
            a = 21.23
            da = 7.82
            b = 29.11
            db = 9.56

        spt = a - b*idx
        tot_err = da**2 # err related to a
        tot_err += np.power(db*idx,2)  # err related to b
        if idx_err is not None:
            tot_err += np.power(idx_err*b, 2)
        spt_err = np.sqrt(tot_err)

    elif name == 'H2O-1.5':
        a = 0.75
        da = 0.03
        b = 0.044
        db = 0.004

        spt = (idx-a)/b
        tot_err = np.power(da/b,2) # err related to a
        tot_err += np.power(db*(idx-a)/b**2,2)  # err related to b
        if idx_err is not None:
            tot_err += np.power(idx_err/b, 2)
        spt_err = np.sqrt(tot_err)

    elif name == 'H2O-2':
        a = 34.13
        da = 1.19
        b = 27.10
        db = 1.20
        spt = a - b*idx
        tot_err = da**2 # err related to a
        tot_err += np.power(db*idx,2)  # err related to b
        if idx_err is not None:
            tot_err += np.power(idx_err*b, 2)
        spt_err = np.sqrt(tot_err)

    else:
        raise ValueError("not implemented for this line")

    return spt, spt_err

def sp_idx_to_gravity(idx, name='Na-1.1'):
    """
    Provides a qualitative estimate of the gravity/youth based on a
    gravity-sensitive spectral index. Implemented so far:
        - the Na-1.1 index [ALL07]
        - the CO-2.3 index [GOR03]

    Parameters
    ----------
    idx : float
        Value of spectral index
    name: str, optional {'Na-1.1', 'CO-2.3'}
        The name of the spectral index.

    """

    if name == 'Na-1.1':
        if idx[0] < 1.4:
            msg = "The object's gravity is consistent with being very young "
            msg+= "(Na index less than 1.4; Allers et al. 2007))"
        else:
            msg= "The Na index does not suggest a very low gravity (i.e. it is"
            msg+= " unlikely to be very young; Allers et al. 2007))"
    elif name == 'CO-2.3':
        if idx[0] < 0.81:
            msg = "The CO-2.3 index is consistent with either a red giant or "
            msg += " a field L-dwarf with log(g) > 4.5 (CO-2.3 index less than"
            msg += " 0.88; Gorlova et al. 2003)"
        elif idx[0] < 0.88:
            msg = "The object's gravity is likely higher than 4.5, with a SpT "
            msg += "earlier than L0 (CO-2.3 index less than 0.88; Gorlova et "
            msg += "al. 2003)"
        else:
            msg= "The object is likely young, with a log(g) lower than 4.5 "
            msg+= "(CO-2.3 index larger than 0.88; Gorlova et al. 2003)"

    print(msg)

def spt_to_digit(spt, convention='splat'):
    """
    Converts a string representing spectral type into an integer index.

    Parameters
    ----------
    spt : str
        String representing the spectral index
    convention: str, optional {'splat', 'Allers+07'}
        Which convention to use to convert digit into spectral type.
        Convention from splat: K0 = 0, M0=10, L0=20, T0=30, Y9 = 49
        Convention from [ALL07]: M0 = 0, L0 = 10, ...

    Returns
    -------
    idx : float or int
        Index value of the spectral type
    """

    spt_str = spt.split()[0]
    sub = 0
    if convention == 'Allers+07':
        sub=10

    if 'K' in spt_str:
        idx = 0-sub
        loc = spt_str.find('K')
        if convention == 'Allers+07':
            msg="WARNING: Allers+07 convention doesn't handle SpT earlier than M"
            print(msg)
    elif 'M' in spt_str:
        idx = 10-sub
        loc = spt_str.find('M')
    elif 'L' in spt_str:
        idx = 20-sub
        loc = spt_str.find('L')
    elif 'T' in spt_str:
        idx = 30-sub
        loc = spt_str.find('T')
    elif 'Y' in spt_str:
        idx = 40-sub
        loc = spt_str.find('Y')
    else:
        idx = np.nan
        print("WARNING: Spectral type not between K0 and Y9. Idx set to NaN.")

    if len(spt)>1:
        if len(spt[loc+1:]) > 2:
            if spt[loc+2] == ".":
                idx+=float(spt[loc+1:loc+4])
        elif spt[loc+1].isdigit():
            idx+=float(spt[loc+1])
    else:
        idx+=5 # in case just "M" or "L" is given, let's assume it's M5 or L5

    return idx

def digit_to_spt(idx, convention='splat'):
    """
    Converts an integer index into spectral type.

    Parameters
    ----------
    idx : float or int
        Index value of the spectral type
    convention: str, optional {'splat', 'Allers+07'}
        Which convention to use to convert digit into spectral type.
        Convention from splat: K0 = 0, M0=10, L0=20, T0=30, Y9 = 49
        Convention from Allers+07: M0 = 0, L0 = 10, ...

    Returns
    -------
    spt : str
        String representing the spectral index
    """

    if convention == 'splat':
        if idx < 10:
            spt = 'K'
        elif idx < 20:
            spt = 'M'
        elif idx < 30:
            spt = 'L'
        elif idx < 40:
            spt = 'T'
        elif idx < 50:
            spt = 'Y'
        else:
            raise ValueError("Index not between 0 (K0) and 49 (Y9)")
    elif convention == 'Allers07':
        if idx < 10:
            spt = 'M'
        elif idx < 20:
            spt = 'L'
        else:
            raise ValueError("Index not between 0 (M0) and 20 (L9)")

    if not (idx%10)%1:
       spt+="{:.0f}".format(idx%10)
    else:
        spt+="{:.1f}".format(idx%10)

    return spt

ext: py

reason: import numpy


text:
""" Utility functions for extracting STIS spectral data into lightcurves

"""

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

import os
import numpy as np
import scipy
from scipy.interpolate import interp1d
from datetime import datetime
from numba import jit
import astropy
from astropy.io import fits as fits

from .utils import expand_refname, enlarge
from .cos import extract_index, calc_npixels
from .version import version as  __version__

__all__ = ['extract',
           'stis_corrtag',
           'map_image',
           'epsilon',
           'dqinit']

#-------------------------------------------------------------------------------

def extract(filename, **kwargs):
    """ Extract lightcurve from STIS dataset

    This is the main driver of the lightcuve extracion, and definitely
    needs some better documentation.

    Parameters
    ----------
    filename : str
        name of FITS file to extract from
    **kwargs : dict
        arbitrary keyword arguements for tailored extraction

    Kwarg parameters
    ----------------
    verbosity : int, default=0
        Verbosity level for print output
    step : int, default=1
        timestep in seconds for output Lightcurve
    wlim : tuple

    Returns
    -------
    data, meta : Astropy table, dict
        Table with extracted data and dictionary of metadata pairs
    """

    verbosity = kwargs.get('verbosity', 0)
    step = kwargs.get('step', 1)
    wlim = kwargs.get('wlim', (2, 10000))
    xlim = kwargs.get('xlim', (0, 2048))
    ylim = kwargs.get('ylim', (0, 2048))
    filter_airglow = kwargs.get('filter_airglow', True)

    hdu = fits.open(filename)
    if fits.getval(filename, 'OBSTYPE') == 'IMAGING':
        print("Imaging observation found, resetting limits.")
        xlim = (0, 2048)
        ylim = (0, 2048)
        wlim = (-1, 1)

    input_headers = {'a':{}}

    for i, ext in enumerate(hdu):
        try:
            input_headers['a'][i] = ext.header._cards
        except AttributeError:
            pass

    meta = {'source': filename,
            'instrument' : 'STIS',
            'headers': input_headers,
            'stepsize': step,
            'wlim': wlim,
            'xlim': xlim,
            'ylim': ylim}

    if verbosity:
        print('#-----------------------------------------#')
        print('Running LightCurve extraction for STIS Data')
        print('#-----------------------------------------#')
        print()
        print('Extracting from: {}'.format(filename))
        print('With arguments:')
        print('step : {}'.format(step))
        print('wlim : {}'.format(wlim))
        print('xlim : {}'.format(xlim))
        print('ylim : {}'.format(ylim))
        print('filter_airglow : {}'.format(filter_airglow))
        print()

    SECOND_PER_MJD = 1.15741e-5

    time = hdu[1].data['time']
    #time = np.array([round(val, 3) for val in hdu[1].data['time']]).astype(np.float64)

    if not len(time):
        end = 0
        start = 0
    else:
        end = min(time.max(), hdu[1].header['EXPTIME'])
        start = time.min()

    all_steps = np.arange(start, end+step, step)

    if all_steps[-1] > end:
        truncate = True
    else:
        truncate = False

    ### need to fix sdqflags
    index = extract_index(hdu,
                          xlim[0], xlim[1],
                          ylim[0], ylim[1],
                          wlim[0], wlim[1],
                          hdu[1].header['SDQFLAGS'],
                          filter_airglow=filter_airglow)

    if verbosity:
        print("#{} events".format(len(index)))

    gross = np.histogram(hdu['events'].data['time'][index],
                         all_steps,
                         weights=hdu['events'].data['epsilon'][index])[0]

    n_pixels = calc_npixels(hdu, index, xlim)

    response_array = get_fluxes(hdu, index).mean()
    weights = hdu['events'].data['epsilon'][index] / step / response_array
    print("WARNING: The flux is a lie")
    flux =  np.histogram(hdu['events'].data['time'][index], all_steps, weights=weights)[0] / n_pixels

    print("WARNING: The background is a lie")
    background = np.zeros(gross.shape)
    mjd = hdu[1].header['EXPSTART'] + np.array(all_steps[:-1]) * SECOND_PER_MJD
    bins = np.ones(len(gross)) * step
    times = all_steps[:-1]

    if truncate:
        if verbosity:
            print('Truncating the last event bin')

        gross = gross[:-1]
        flux = flux[:-1]
        background = background[:-1]
        mjd = mjd[:-1]
        bins = bins[:-1]
        times = times[:-1]

    data = {'dataset': np.ones(times.shape),
            'times': times,
            'mjd': mjd,
            'bins': bins,
            'gross': gross,
            'background': background,
            'flux': flux}

    if verbosity:
        print('Finished extraction for {}'.format(filename))
        print()

    return data, meta

#-------------------------------------------------------------------------------

def stis_corrtag(tagfile, clean=True):
    """Create a COS-like corrtag file for STIS data

    Parameters
    ----------
    tagfile, str
        input STIS time-tag data file

    """
    print("Creating STIS corrtag for {}".format(tagfile))
    if '_tag.fits' in tagfile:
        x1d_filename = tagfile.replace('_tag.fits', '_x1d.fits')
    elif '_corrtag.fits' in tagfile:
        x1d_filename = tagfile.replace('_corrtag.fits', '_x1d.fits')

    with fits.open(tagfile, 'readonly') as hdu:
        n_events = len(hdu[1].data['TIME'])

        time_data = hdu[1].data['TIME']

        if 'DETAXIS1' in hdu[1].data.names:
            #-- Y data is same in raw and corrected coordinates
            rawx_data = hdu[1].data['DETAXIS1']
            rawy_data = hdu[1].data['AXIS2']
            xcorr_data = hdu[1].data['AXIS1']
            ycorr_data = hdu[1].data['AXIS2']
        else:
            rawx_data = hdu[1].data['RAWX']
            rawy_data = hdu[1].data['RAWY']
            xcorr_data = hdu[1].data['XCORR'].astype(np.int32)
            ycorr_data = hdu[1].data['YCORR'].astype(np.int32)

        #-- copy over the primary
        header0 = hdu[0].header.copy()
        header1 = hdu[1].header.copy()

    eps_data = epsilon(tagfile)
    dq_data = dqinit(tagfile)
    if header0['OBSTYPE'] == 'IMAGING':
        wave_data = np.zeros(n_events)

    elif not os.path.exists(x1d_filename):
        print("Could not find associated extracted spectrum {}".format(x1d_filename))
        wave_data = np.ones(n_events) * hdu[0].header['CENTRWV']

    elif header0['OPT_ELEM'].startswith('E'):
        #-- Put in average wavelength for now
        with fits.open(x1d_filename) as x1d:
            #--initialize everything to 0 wavelength
            wave_data = np.zeros(n_events)

            offset = header1['shifta2']
            for order in x1d[1].data:
                spec_center = order['a2center'] + offset
                spec_height = order['extrsize']

                index = np.where((ycorr_data < spec_center+spec_height) &
                                 (ycorr_data > spec_center-spec_height))[0]

                int_pix = integerize_pixels(xcorr_data[index])
                wave_data[index] = order['wavelength'][int_pix]

    else:
        #-- Grab wavelengths from the x1d file
        int_pix = integerize_pixels(xcorr_data)
        with fits.open(x1d_filename) as x1d:
            wave_data = x1d[1].data['wavelength'][0][int_pix]

    #-- Writeout corrtag file
    hdu_out = fits.HDUList(fits.PrimaryHDU())

    hdu_out[0].header['GEN_DATE'] = (str(datetime.now()), 'Creation Date')
    hdu_out[0].header['LC_VER'] = (__version__, 'lightcurve version used')
    hdu_out[0].header['AP_VER'] = (astropy.__version__, 'Astropy version used')
    hdu_out[0].header['NP_VER'] = (np.__version__, 'Numpy version used')
    hdu_out[0].header['SP_VER'] = (scipy.__version__, 'Scipy version used')

    hdu_out[0].header.extend(header0, end=True)

    time_col = fits.Column('time', 'D', 'second', array=time_data)
    rawx_col = fits.Column('rawx', 'D', 'pixel', array=rawx_data)
    rawy_col = fits.Column('rawy', 'D', 'pixel', array=rawy_data)
    xcorr_col = fits.Column('xcorr', 'D', 'pixel', array=xcorr_data)
    ycorr_col = fits.Column('ycorr', 'D', 'pixel', array=ycorr_data)
    epsilon_col = fits.Column('epsilon', 'D', 'factor', array=eps_data)
    wave_col = fits.Column('wavelength', 'D', 'angstrom', array=wave_data)
    dq_col = fits.Column('dq', 'I', 'DQ', array=dq_data)

    tab = fits.new_table([time_col,
                            rawx_col,
                            rawy_col,
                            xcorr_col,
                            ycorr_col,
                            epsilon_col,
                            wave_col,
                            dq_col])
    hdu_out.append(tab)

    hdu_out[1].header.extend(header1, end=True)

    hdu_out.writeto(tagfile.replace('_tag.fits', '_corrtag.fits'), clobber=True)
    if clean and '_tag.fits' in tagfile:
        print("Removing input tagfile")
        os.remove(tagfile)

#-------------------------------------------------------------------------------

def integerize_pixels(xcoords):

    int_pix = np.round(xcoords).astype(np.int32)
    int_pix = np.where(int_pix < 0, 0, int_pix)
    int_pix = np.where(int_pix > 2047, 2047, int_pix)
    int_pix //= 2

    return int_pix

#-------------------------------------------------------------------------------

@jit
def map_image(image, xcoords, ycoords, default=0):
    n_coord = len(xcoords)
    out_vals = np.zeros(n_coord)

    for i in range(n_coord):
        x = xcoords[i]
        y = ycoords[i]

        if x < 0 or x >= 2048:
            val = default
        elif y < 0 or y >= 2048:
            val = default
        else:
            val = image[y, x]

        out_vals[i] = val

    return out_vals

#-------------------------------------------------------------------------------

def epsilon(tagfile):
    """Compute the total epsilon factor for each event

    Compute the flatfield correction from the P-flat and L-flat reference files
    (PFLTFILE and LFLTFILE respectively).

    Parameters
    ----------
    tagfile, str
        input STIS time-tag data file

    Returns
    -------
    epsilon, np.ndarray
        array of epsilons

    """

    print("Calculating Epsilon")

    with fits.open(tagfile) as hdu:

        epsilon_out = np.ones(hdu[1].data['time'].shape)

        #-- Flatfield correction
        for ref_flat in ['PFLTFILE', 'LFLTFILE']:
            reffile = expand_refname(hdu[0].header[ref_flat])
            print('FLATFIELD CORRECTION {}: {}'.format(ref_flat, reffile))

            if not os.path.exists(reffile):
                print("{} not found, correction not performed".format(reffile))
                return np.ones(len(hdu[1].data))

            with fits.open(reffile) as image_hdu:
                image = image_hdu[1].data

                if not image.shape == (2048, 2048):
                    x_factor = 2048 // image.shape[1]
                    y_factor = 2048 // image.shape[0]

                    print('Enlarging by {},{}'.format(x_factor, y_factor))
                    image = enlarge(image, x_factor, y_factor)

                #--indexing is 1 off
                if 'AXIS1' in hdu[1].data.names:
                    epsilon_out *= map_image(image,
                                             hdu[1].data['AXIS1'] - 1,
                                             hdu[1].data['AXIS2'] - 1)
                else:
                    epsilon_out *= map_image(image,
                                             hdu[1].data['XCORR'].astype(np.int32) - 1,
                                             hdu[1].data['YCORR'].astype(np.int32) - 1)

    return epsilon_out

#-------------------------------------------------------------------------------

def dqinit(tagfile):
    """Compute the data quality information for each pixel from the BPIXTAB.

    Parameters
    ----------
    tagfile, str
        input STIS time-tag data file

    Returns
    -------
    dq, np.ndarray
        array of bitwise dq flags

    """

    print("Calculating DQ")

    with fits.open(tagfile) as hdu:

        reffile = expand_refname(hdu[0].header['BPIXTAB'])
        print('BPIXTAB used {}'.format(reffile))

        if not os.path.exists(reffile):
            print("{} not found, correction not performed".format(reffile))
            return np.zeros(len(hdu[1].data))

        with fits.open(reffile) as bpix:
            #-- Mama bpix regions are in lo-res pixels
            dq_im = np.zeros((1024, 1024)).astype(np.int32)

            #-- indexing is 1 off
            for line in bpix[1].data:
                lx = line['PIX1'] - 1
                ly = line['PIX2'] - 1
                dpix = line['LENGTH']
                axis = line['AXIS']
                flag = line['VALUE']

                if axis == 1:
                    #-- Along X
                    dq_im[ly, lx:lx + dpix] |= flag
                elif axis == 2:
                    #-- Along Y
                    dq_im[ly:ly+dpix, lx] |= flag

            #-- Needs to be expanded into Hi-res
            dq_im = enlarge(dq_im, 2, 2)

        #-- Map to the events
        #-- indexing is 1 off
        if 'AXIS1' in hdu[1].data.names:
            dq_out = map_image(dq_im,
                               hdu[1].data['AXIS1'] - 1,
                               hdu[1].data['AXIS2'] - 1)
        else:
            dq_out = map_image(dq_im,
                               hdu[1].data['XCORR'].astype(np.int32) - 1,
                               hdu[1].data['YCORR'].astype(np.int32) - 1)

    return dq_out

#-------------------------------------------------------------------------------

def crazy_histogram2d(x, y, bins=(2048, 2048)):
    """  Bin events to a 2d image.

    This is faster than the numpy version.

    http://stackoverflow.com/questions/8805601/efficiently-create-2d-histograms-from-large-datasets

    Parameters
    ----------
    x : array-like
        x-coordinates to bin
    y : array-like
        y-coordinates to bin
    bins : tuple, optional
        bins of output image

    Returns
    -------
    grid, x_bins, y_bins : (np.ndarray, np.ndarray, np.ndarray)
        binned image, x bins, and y bins
    """

    try:
        nx, ny = bins
    except TypeError:
        nx = ny = bins

    xmin, xmax = x.min(), x.max()
    ymin, ymax = y.min(), y.max()
    dx = (xmax - xmin) / (nx - 1.0)
    dy = (ymax - ymin) / (ny - 1.0)

    weights = np.ones(x.size)

    # Basically, this is just doing what np.digitize does with one less copy
    xyi = np.vstack((x,y)).T
    xyi = np.floor(xyi, xyi).T

    # Now, we'll exploit a sparse coo_matrix to build the 2D histogram...
    grid = scipy.sparse.coo_matrix((weights, xyi), shape=(nx, ny)).toarray().T

    return grid, np.linspace(xmin, xmax, nx), np.linspace(ymin, ymax, ny)

#-------------------------------------------------------------------------------

def get_fluxes(hdu, index):
    """ Return response curve and wavelength range for specified mode
    """

    fluxtab = hdu[0].header.get('FLUXTAB', '')

    fluxfile = expand_refname(fluxtab)

    if not len(index):
        return np.ones(hdu['events'].data['time'].shape)[index]

    if (not fluxfile) or (not os.path.exists(fluxfile)):
        print(' WARNING: Fluxfile not available %s,' % fluxfile )
        print(' using unity flux calibration instead.')
        return np.ones(hdu['events'].data['time'].shape)[index]

    flux_hdu = fits.open(fluxfile)
    setting_index = np.where((flux_hdu[1].data['SPORDER'] == hdu[0].header['segment']) &
                             (flux_hdu[1].data['OPT_ELEM'] == hdu[0].header['opt_elem']) &
                             (flux_hdu[1].data['CENWAVE'] == hdu[0].header['cenwave']))[0]

    if len(setting_index) == 0:
        print('No row in fluxtab found for this dataset, no FLUXCAL performed')
        return np.ones(hdu['events'].data['time'].shape)[index]
    elif len(setting_index) > 1:
        raise ValueError('Too many rows found: {}'.format(len(setting_index)))

    resp_wave = flux_hdu[1].data[setting_index]['WAVELENGTH'].flatten()
    response = flux_hdu[1].data[setting_index]['SENSITIVITY'].flatten()

    data_max = hdu['events'].data['wavelength'][index].max()
    data_min = hdu['events'].data['wavelength'][index].min()

    if data_max > resp_wave.max():
        print("Expanding minumum response curve by {}".format(data_max - resp_wave.max()))
        resp_wave[np.argmax(resp_wave)] = data_max

    if data_min < resp_wave.min():
        print("Expanding minimum response curve by {}".format(data_min - resp_wave.min()))
        resp_wave[np.argmin(resp_wave)] = data_min

    interp_func = interp1d(resp_wave, response)
    all_resp = interp_func(hdu['events'].data['wavelength'][index])

    return all_resp

#-------------------------------------------------------------------------------

max_stars_repo_name: codacy-badger/lightcurve

max_issues_repo_name: codacy-badger/lightcurve

max_forks_repo_name: codacy-badger/lightcurve

reason: import numpy,import scipy,from scipy,from numba,import astropy,from astropy


text:
##########################################################################
# Testing of display and capture & storage thread combined.
# Scan for camera
# Aquire 14 images
# Convert to b/w
# Save hdf5 files
##########################################################################
# Results
#
#
##########################################################################

# System
import logging, time, platform
from datetime import datetime

# Matrix Algebra
import numpy as np
from numba import vectorize

# OpenCV
import cv2

configs = {
    'camera_res'      : (1280, 720 ),   # width & height
    'exposure'        : -6,             # -1,0 = auto, 1...max=frame interval,
    'autoexposure'    : 1.0,            # depends on camera: 0.25 or 0.75(auto), -1,0,1
    'fps'             : 30,             # 15, 30, 40, 90, 120, 180
    'fourcc'          : -1,             # n.a.
    'buffersize'      : -1,             # n.a.
    'output_res'      : (-1, -1),       # Output resolution, -1,-1 no change
    'flip'            : 0,              # 0=norotation
                                        # 1=ccw90deg
                                        # 2=rotation180
                                        # 3=cw90
                                        # 4=horizontal
                                        # 5=upright diagonal flip
                                        # 6=vertical
                                        # 7=uperleft diagonal flip
    'displayfps'       : 30             # frame rate for display server
}

if configs['displayfps'] >= configs['fps']:
    display_interval = 0
else:
    display_interval = 1.0/configs['displayfps']

res = configs['camera_res']
height = res[1]
width = res[0]
measure_time = 5.0 # average measurements over 5 secs
camera_index = 0 # default camera starts at 0 by operating system

# Processing
data_cube      = np.zeros((14, height, width), 'uint8')
background     = np.zeros((height, width), 'uint8')                         # where we keep bg
flatfield      = np.cast['uint16'](2**8.*np.random.random((height, width))) # flatfield correction scaled so that 255=100%
inten         = np.zeros(14, 'uint16')                                # helper to find background image
data_cube_corr = np.zeros((14, height, width), 'uint16')                    # resulting data cube on CPU
# Numpy Vectorized
@vectorize(['uint16(uint8, uint16, uint8)'], nopython=True, fastmath=True)
def vector_np(data_cube, background, flatfield):
    return np.multiply(np.subtract(data_cube, background), flatfield) # 16bit multiplication

# Setting up logging
logging.basicConfig(level=logging.DEBUG) # options are: DEBUG, INFO, ERROR, WARNING
logger = logging.getLogger("Main")

# Setting up storage
from camera.streamer.h5storageserver import h5Server
now = datetime.now()
filename = now.strftime("%Y%m%d%H%M%S") + ".hdf5"
hdf5 = h5Server("C:\\temp\\" + filename)
logger.log(logging.INFO, "Starting Storage Server")
hdf5.start()

# Create camera interface
logger.log(logging.INFO, "Starting Capture")
# Computer OS and platform dependent
plat = platform.system()
if plat == 'Linux' and platform.machine() == "aarch64": # this is jetson nano for me
    from camera.capture.nanocapture import nanoCapture
    camera = nanoCapture(configs, camera_index)
else:
    from camera.capture.cv2capture_process import cv2Capture
    camera = cv2Capture(configs, camera_index)
logger.log(logging.INFO, "Getting Images")
camera.start()

# Display
window_name    = 'Main'
font           = cv2.FONT_HERSHEY_SIMPLEX
textLocation0  = (10,height-40)
textLocation1  = (10,height-20)
fontScale      = 1
fontColor      = (255,255,255)
lineType       = 2
cv2.namedWindow(window_name, cv2.WINDOW_AUTOSIZE) # or WINDOW_NORMAL

# Initialize Variables
frame_idx              = 0  # index to create data cube out of individual frames
num_cubes_stored       = 0  # keep track of data cubes sent to storage
num_cubes_generated    = 0  # keep track of data cubes generated
last_time              = time.perf_counter() # keep track of time to dispay performance
last_display           = time.perf_counter() # keeo track of time to display images
num_frames_received    = 0  # keep track of how many captured frames reach the main program
num_frames_displayed   = 0  # keep track of how many frames are displayed
measured_dps           = 0  # computed in main thread, number of frames displayed per second

# Main Loop
stop =  False
while(not stop):
    current_time = time.perf_counter()

    # Camera
    (frame_time, frame) = camera.capture.get(block=True, timeout=None)
    data_cube[frame_idx,:,:] = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
    frame_idx += 1
    while not camera.log.empty():
        (level, msg)=camera.log.get_nowait()
        logger.log(level, "Status:{}".format(msg))

    # When we have a complete dataset:
    if frame_idx >= 14: # 0...13 is populated
        frame_idx = 0
        num_cubes_generated += 1

        # HDF5
        try:
            hdf5.queue.put_nowait((frame_time, data_cube))
            num_cubes_stored += 1
        except:
            logger.log(logging.DEBUG, "Status:Storage Queue is full!")

        # Background and Field Correction
        # Where is my background?
        _ = np.sum(data_cube[:,::64,::64], axis=(1,2), out = inten)
        frame_idx_bg = np.argmin(inten)  # minimum intensity is in this frame
        background = data_cube[frame_idx_bg, :, :]
        # Correct the data
        vector_np(data_cube, flatfield, background, out = data_cube_corr)

    while not hdf5.log.empty():
        (level, msg)=hdf5.log.get_nowait()
        logger.log(level, "Status:{}".format(msg))

    # Display performance in main loop
    if current_time - last_time >= measure_time:
        # how many data cubes did we create
        measured_cps_generated = num_cubes_generated/measure_time
        logger.log(logging.DEBUG, "Status:captured cubes generated per second:{}".format(measured_cps_generated))
        num_cubes_generated = 0
        # how many data cubes did we send to storage
        measured_cps_stored = num_cubes_stored/measure_time
        logger.log(logging.DEBUG, "Status:cubes sent to storage per second:{}".format(measured_cps_stored))
        num_cubes_sent = 0
        # how many frames did we display
        measured_dps = num_frames_displayed/measure_time
        logger.log(logging.DEBUG, "Status:frames displayed per second:{}".format(measured_dps))
        num_frames_displayed = 0
        last_time = current_time

    if (current_time - last_display) >= display_interval:
        display_frame = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
        cv2.putText(display_frame,"Capture FPS:{} [Hz]".format(camera.measured_fps), textLocation0, font, fontScale, fontColor, lineType)
        cv2.imshow(window_name, display_frame)
        if cv2.waitKey(1) & 0xFF == ord('q'): stop = True
        last_display = current_time
        num_frames_displayed += 1

# Cleanup
camera.stop()
hdf5.stop()
cv2.destroyAllWindows()

reason: import numpy,from numba


text:
"""a module for basic I/O operations
"""
__author__ = "Reed Essick (reed.essick@gmail.com)"

#-------------------------------------------------

import os
import shutil

import numpy as np

from . import utils

#-------------------------------------------------

ALPHANUMERIC = 'a b c d e f g h i j k l m n o p q r s t u v w x y z A B C D E F G H I J K L M N O P Q R S T U V W X Y Z 0 1 2 3 4 5 6 7 8 9'.split()
DEFAULT_MAX_NUM_SAMPLES = np.infty

#-------------------------------------------------
# basic utilities for manipulating existing sapmles
#-------------------------------------------------

def column2logcolumn(name):
    return 'log(%s)'%name

def check_columns(present, required, logcolumns=[]):
    required = [column2logcolumn(column) if column in logcolumns else column for column in required]
    for column in required:
        assert column in present, 'required column=%s is missing!'%column

def load(inpath, columns=[], logcolumns=[], max_num_samples=DEFAULT_MAX_NUM_SAMPLES):
    data = []

    ### make sure we load only a single copy of each column; preserve the order as much as possible
    unique = []
    for col in columns:
        if col not in unique:
            unique.append(col)
    columns = unique

    with open(inpath, 'r') as obj:
        cols = [_.strip() for _ in obj.readline().strip().split(',')]
        if columns:
            check_columns(cols, columns)
        else:
            columns = cols

        inds = [cols.index(col) for col in columns]

        count = 0
        for line in obj:
            if line[0]=='#':
                continue
            if count >= max_num_samples:
                break
            fields = line.strip().split(',')
            ans = [] ### downselect to what we actually want and cast to float
            for ind in inds:
                try: ### try casting everything to a float
                    ans.append(float(fields[ind]))
                except ValueError: ### could not cast to a float
                    ans.append(fields[ind])
            data.append(ans)
            count += 1

    data = np.array(data) ### cast as an array

    cols = [] ### figure out column names and map to logs as requested
    for i, col in enumerate(columns):
        if col in logcolumns:
            data[:,i] = np.log(data[:,i])
            cols.append(column2logcolumn(col))
        else:
            cols.append(col)

    return data, cols #, dict((strind, col) for col, strind in strmap.items())

def load_weights(*args, **kwargs):
    """loads and returns weights from multiple columns via  delegation to load_logweights
    normalizes the weights while it's at it
    """
    normalize = kwargs.pop('normalize', True)
    return utils.exp_weights(load_logweights(*args, **kwargs), normalize=normalize)

def load_logweights(inpath, weight_columns, logweightcolumns=[], invweightcolumns=[], max_num_samples=DEFAULT_MAX_NUM_SAMPLES):
    """loads and returns logweights from multiple columns
    """
    data, columns = load(inpath, columns=weight_columns, max_num_samples=max_num_samples) ### load the raw data

    for i, column in enumerate(columns): ### iterate through columns, transforming as necessary
        if column in logweightcolumns:
            if column in invweightcolumns:
                data[:,i] *= -1

        else:
            if column in invweightcolumns:
                data[:,i] = 1./data[:,i]

            data[:,i] = np.log(data[:,i])

    # multiply weights across all samples, which is the same as adding the logs
    return np.sum(data, axis=1)

def write(path, data, cols, footer='', tmpdir='.', tmpdigits=6, verbose=False):
    """a standardized writing function that does an atomic write
    """
    path = os.path.abspath(path) ### make sure have the full path

    tmppath = os.path.join(tmpdir, "."+os.path.basename(path)+'-'+''.join(np.random.choice(ALPHANUMERIC, size=tmpdigits))) ### save to a temporary file
    if verbose:
        print('writing to temporary file location: '+tmppath)
    np.savetxt(tmppath, data, header=','.join(cols), comments='', delimiter=',')

    if verbose:
        print('moving to final location: '+path)
    shutil.move(tmppath, path) ### move to the permanent location

max_stars_repo_name: isaaclegred/universality

max_issues_repo_name: isaaclegred/universality

max_forks_repo_name: isaaclegred/universality


text:
import pandas as pd
import numpy as np
import preprocess as pre
from stocks import tickers

cols = ["Date", "Open", "High", "Low", "Close", "Volume", "Name"]
names = tickers

def detrend():
    df = pre.preprocess()

    # detrends each stock time series using the difference method
    for n in names:
        df[n] = df[n].diff()

    return df

if __name__ == "__main__":
    detrend()


text:
# -*- coding: utf-8 -*-

"""
@author: Guo Jingyao
"""

import numpy as np
import scipy
from ..recommender import Recommender
from ...utils.common import sigmoid
from ...utils.common import scale
from ...exception import ScoreException

class SoRec(Recommender):
    """Social recommendation using Probabilistic Matrix Factorization.

    Parameters
    ----------
    k: int, optional, default: 5
        The dimension of the latent factors.
    max_iter: int, optional, default: 100
        Maximum number of iterations or the number of epochs for SGD.
    learning_rate: float, optional, default: 0.001
        The learning rate for SGD_RMSProp.
    gamma: float, optional, default: 0.9
        The weight for previous/current gradient in RMSProp.
    lamda: float, optional, default: 0.001
        The regularization parameter.
    lamda_C: float, optional, default: 10
        The parameter balances the information from the user-item rating matrix and the user social network.
    name: string, optional, default: 'SOREC'
        The name of the recommender model.
    variant: {"linear","non_linear"}, optional, default: 'non_linear'
        Pmf variant. If 'non_linear', the Gaussian mean is the output of a Sigmoid function.\
        If 'linear' the Gaussian mean is the output of the identity function.
    trainable: boolean, optional, default: True
        When False, the model is not trained and Cornac assumes that the model already \
        pre-trained (U, V and Z are not None).
    verbose: boolean, optional, default: False
        When True, some running logs are displayed.
    init_params: dictionary, optional, default: {'U':None,'V':None}
        List of initial parameters, e.g., init_params = {'U':U, 'V':V,'Z':Z}. \
        U: a ndarray of shape (n_users,l), containing the user latent factors. \
        V: a ndarray of shape (n_items,l), containing the item latent factors. \
        Z: a ndarray of shape (n_users,l), containing the social network latent factors. \

    References
    ----------
    * H. Ma, H. Yang, M. R. Lyu, and I. King. SoRec:Social recommendation using probabilistic matrix factorization. \
     CIKM ’08, pages 931–940, Napa Valley, USA, 2008.
    """

    def __init__(self, name="SoRec", k=5, max_iter=100, learning_rate=0.001, lamda_C=10, lamda=0.001, gamma=0.9,
                 trainable=True, verbose=False, init_params={'U': None, 'V': None, 'Z': None}):
        Recommender.__init__(self, name=name, trainable=trainable, verbose=verbose)
        self.k = k
        self.init_params = init_params
        self.max_iter = max_iter
        self.learning_rate = learning_rate
        self.lamda_C = lamda_C
        self.lamda = lamda
        self.gamma = gamma

        self.ll = np.full(max_iter, 0)
        self.eps = 0.000000001
        self.U = init_params['U']  # matrix of user factors
        self.V = init_params['V']  # matrix of item factors
        self.Z = init_params['V']  # matrix of social network factors

        if self.U is None:
            print("random initialize user factors")
        elif self.U.shape[1] != self.k:
            raise ValueError('initial parameters U dimension error')

        if self.V is None:
            print("random initialize item factors")
        elif self.V.shape[1] != self.k:
            raise ValueError('initial parameters V dimension error')

        if self.Z is None:
            print("random initialize social factors")
        elif self.Z.shape[1] != self.k:
            raise ValueError('initial parameters Z dimension error')

    def fit(self, train_set):
        """Fit the model to observations.
        Parameters
        ----------
        train_set: object of type TrainSet, required
            An object contraining the user-item preference in csr scipy sparse format,\
            as well as some useful attributes such as mappings to the original user/item ids.\
            Please refer to the class TrainSet in the "data" module for details.
        """
        import math
        from cornac.models.sorec import sorec

        Recommender.fit(self, train_set)

        if self.trainable:
            # user-item interactions
            (rat_uid, rat_iid, rat_val) = train_set.uir_tuple
            # user-user social network
            map_uid = train_set.uid_list
            social_net = train_set.user_graph.get_train_triplet(map_uid, map_uid)

            social_raw = scipy.sparse.csc_matrix((social_net[:, 2], (social_net[:, 0], social_net[:, 1])),
                                                 shape=(len(map_uid), len(map_uid)))
            outdegree = np.array(social_raw.sum(axis=1)).flatten()
            indegree = np.array(social_raw.sum(axis=0)).flatten()
            weighted_social = []
            for ui, uk, cik in social_net:
                i_out = outdegree[int(ui)]
                k_in = indegree[int(uk)]
                cik_weighted = math.sqrt(k_in / (k_in + i_out)) * cik
                weighted_social.append(cik_weighted)

            (net_uid, net_jid, net_val) = (social_net[:, 0], social_net[:, 1], weighted_social)

            if [self.train_set.min_rating, self.train_set.max_rating] != [0, 1]:
                if self.train_set.min_rating == self.train_set.max_rating:
                    rat_val = scale(rat_val, 0., 1., 0., self.train_set.max_rating)
                else:
                    rat_val = scale(rat_val, 0., 1., self.train_set.min_rating, self.train_set.max_rating)

            rat_val = np.array(rat_val, dtype='float32')
            rat_uid = np.array(rat_uid, dtype='int32')
            rat_iid = np.array(rat_iid, dtype='int32')

            net_val = np.array(net_val, dtype='float32')
            net_uid = np.array(net_uid, dtype='int32')
            net_jid = np.array(net_jid, dtype='int32')

            if self.verbose:
                print('Learning...')

            res = sorec.sorec(rat_uid, rat_iid, rat_val, net_uid, net_jid, net_val, k=self.k,
                              n_users=train_set.num_users,
                              n_items=train_set.num_items, n_ratings=len(rat_val), n_edges=len(net_val),
                              n_epochs=self.max_iter, lamda_C=self.lamda_C,
                              lamda=self.lamda, learning_rate=self.learning_rate, gamma=self.gamma,
                              init_params=self.init_params, verbose=self.verbose)

            self.U = np.asarray(res['U'])
            self.V = np.asarray(res['V'])
            self.Z = np.asarray(res['Z'])

            if self.verbose:
                print('Learning completed')
        elif self.verbose:
            print('%s is trained already (trainable = False)' % self.name)

    def score(self, user_id, item_id=None):
        """Predict the scores/ratings of a user for an item.
        Parameters
        ----------
        user_id: int, required
            The index of the user for whom to perform score prediction.
        item_id: int, optional, default: None
            The index of the item for that to perform score prediction.
            If None, scores for all known items will be returned.
        Returns
        -------
        res : A scalar or a Numpy array
            Relative scores that the user gives to the item or to all known items
        """
        if item_id is None:
            if self.train_set.is_unk_user(user_id):
                raise ScoreException("Can't make score prediction for (user_id=%d)" % user_id)

            known_item_scores = self.V.dot(self.U[user_id, :])
            return known_item_scores
        else:
            if self.train_set.is_unk_user(user_id) or self.train_set.is_unk_item(item_id):
                raise ScoreException("Can't make score prediction for (user_id=%d, item_id=%d)" % (user_id, item_id))

            user_pred = self.V[item_id, :].dot(self.U[user_id, :])
            user_pred = sigmoid(user_pred)
            if self.train_set.min_rating == self.train_set.max_rating:
                user_pred = scale(user_pred, 0., self.train_set.max_rating, 0., 1.)
            else:
                user_pred = scale(user_pred, self.train_set.min_rating, self.train_set.max_rating, 0., 1.)

            return user_pred

max_stars_repo_name: GuoJingyao/cornac

max_issues_repo_name: GuoJingyao/cornac

max_forks_repo_name: GuoJingyao/cornac


text:
# day 11: cellular automata, whee
import numpy as np
import pprint

# custom pretty-print wrapper for convenience)
def pprint_c(array):
    pprint.pprint(["".join(seq) for seq in array])

def check_neighbors(inarray, loc):
    # loc is a tuple (row_i, col_j)

    # define limits in each direction
    Nlim = -1
    Wlim = -1
    Slim, Elim = inarray.shape

    # construct neighbors
    north = (loc[0], loc[1] + 1)
    east = (loc[0] + 1, loc[1])
    south = (loc[0], loc[1] - 1)
    west = (loc[0] - 1, loc[1])
    ne = (loc[0] + 1, loc[1] + 1)
    se = (loc[0] + 1, loc[1] - 1)
    sw = (loc[0] - 1, loc[1] - 1)
    nw = (loc[0] - 1, loc[1] + 1)

    adjFilled = 0
    for nbr in (north, east, south, west, ne, se, sw, nw):
        if nbr[0] < Slim and nbr[0] > Nlim:
            if nbr[1] < Elim and nbr[1] > Wlim:
                # check each neighbor now
                if inarray[nbr] == "#":
                    adjFilled += 1
    return adjFilled

def check_visible(inarray, inloc):
    # loc is a tuple (row_i, col_j)
    def step_loc(inloc, dirstr):
        # subfunction to step in a direction based on a string, why not
        # string can be composed of the characters NESW. I won't pass it
        # nonsensical stuff like NS or NESW or whatever, so I won't write it
        # to check the input logic very carefully
        i, j = inloc
        if "N" in dirstr:
            ii = i - 1
        elif "S" in dirstr:
            ii = i + 1
        else:
            ii = i

        if "W" in dirstr:
            jj = j - 1
        elif "E" in dirstr:
            jj = j + 1
        else:
            jj = j

        return (ii, jj)

    # define limits in each direction
    Nlim = -1
    Wlim = -1
    Slim, Elim = inarray.shape

    visibleFilled = 0
    for dirstr in ["N", "NE", "E", "SE", "S", "SW", "W", "NW"]:
        loc = step_loc(inloc, dirstr)
        while (loc[0] < Slim and loc[0] > Nlim) and (loc[1] < Elim and loc[1] > Wlim):
            if inarray[loc] == "#":
                visibleFilled += 1
                break
            elif inarray[loc] == "L":
                break
            loc = step_loc(loc, dirstr)

    return visibleFilled

def main():
    inraw = np.loadtxt("input/day11.txt", dtype="unicode_")
    inarray = np.array([[a for a in inraw[i]] for i in range(len(inraw))])
    Nrows, Ncols = inarray.shape
    # print(M, N)
    # pprint_c(inarray)
    # test = check_neighbors(inarray, (0, 0))
    # print(test)

    # make two arrays, one for step M, and one for step M+1 (N)
    arrayM = inarray.copy()
    arrayN = inarray.copy()
    iterationNum = 0
    while True:
        for loc, spot in np.ndenumerate(arrayM):
            if spot == "L":
                if check_neighbors(arrayM, loc) == 0:
                    arrayN[loc] = "#"
            elif spot == "#":
                if check_neighbors(arrayM, loc) >= 4:
                    arrayN[loc] = "L"

        iterationNum += 1
        if (arrayM == arrayN).all():
            break
        else:
            arrayM = arrayN.copy()
    # pprint_c(arrayN)
    # print(iterationNum)
    print(f"Part 1: {np.count_nonzero(arrayN == '#')}")

    # part 2: annoying sightline stuff instead of simple adjacency. Repeat the
    # same basic loop as before though.

    # make two arrays, one for step M, and one for step M+1 (N)
    arrayM = inarray.copy()
    arrayN = inarray.copy()
    iterationNum = 0
    while True:
        for loc, spot in np.ndenumerate(arrayM):
            if spot == "L":
                if check_visible(arrayM, loc) == 0:
                    arrayN[loc] = "#"
            elif spot == "#":
                if check_visible(arrayM, loc) >= 5:
                    arrayN[loc] = "L"

        # pprint_c(arrayN[:, :])
        # print()
        iterationNum += 1
        if (arrayM == arrayN).all():
            break
        else:
            arrayM = arrayN.copy()
    print(f"Part 2: {np.count_nonzero(arrayN == '#')}")

if __name__ == '__main__':
    main()


text:
import numpy as np
from torch.utils import data
from torch.utils.data.sampler import WeightedRandomSampler
from data.data_loader import SourceDataset, ReferenceDataset

def make_weighted_sampler(labels):
    class_counts = np.bincount(labels)
    class_weights = 1. / class_counts
    weights = class_weights[labels]
    return WeightedRandomSampler(weights, len(weights))

def create_data_loader(opt, which='source', type=None):
    print('Preparing %s dataset during %s phase...' % (which, type))
    if which == 'source':
        dataset = SourceDataset(opt, type)
    elif which == 'reference':
        dataset = ReferenceDataset(opt, type)
    else:
        raise NotImplementedError

    if type == 'train':
        sampler = make_weighted_sampler(dataset.targets)
        return data.DataLoader(dataset=dataset,
                            batch_size=opt.batch_size,
                            sampler=sampler,
                            num_workers=opt.num_workers,
                            pin_memory=True,
                            drop_last=True)
    elif type is not None:
        return data.DataLoader(dataset=dataset, batch_size=opt.batch_size)
    else:
        raise NotImplementedError('Please specify dataset type!')

max_stars_repo_name: joverwey/Diverse-Motion-Stylization

max_issues_repo_name: joverwey/Diverse-Motion-Stylization

max_forks_repo_name: joverwey/Diverse-Motion-Stylization


text:
# -*- coding: utf-8 -*-
'''
Created on Fri Dec 14 11:35:28 2018

@author:

    Visa Suomi
    Turku University Hospital
    December 2018

@description:

    This model is used to predict radiation dose from pre-treatment patient
    parameters

'''

#%% clear variables

%reset -f
%clear

#%% import necessary libraries

import xgboost as xgb
import pandas as pd
import numpy as np
import sklearn as sk
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, mean_squared_error
from sklearn.utils.class_weight import compute_sample_weight
import scipy as sp
import time
import os

from save_load_variables import save_load_variables
from plot_regression_performance import plot_regression_performance
from analyse_statistics import analyse_statistics
from analyse_correlation import analyse_correlation
from analyse_feature_correlation import analyse_feature_correlation
from plot_feature_importance import plot_feature_importance

#%% define logging and data display format

pd.options.display.max_rows = 10
pd.options.display.float_format = '{:.1f}'.format

#%% read data

df = pd.read_csv(r'C:\Users\visa\Documents\TYKS\Machine learning\Radiation dose\radiation-dose\radiation_dose_data.csv', sep = ',')
df_orig = df.copy()

#%% check for duplicates

duplicates = any(df.duplicated())

#%% create synthetic features

df['BSA'] = 0.007184 * df['paino'].pow(0.425) * df['pituus'].pow(0.725)

#%% calculate nan percent for each label

nan_percent = pd.DataFrame(df.isnull().mean() * 100, columns = ['% of NaN'])

# drop nan values

#df = df.dropna()
#df = df.dropna(subset = ['paino'])
#df = df.dropna(subset = ['AHA_cto'])
#df = df.dropna(subset = ['AHA_b1'])
#df = df.dropna(subset = ['Patient_sex'])
#df = df.dropna(subset = ['FN2BA'])
#df = df.dropna(subset = ['I20.81_I21.01_I21.11_or_I21.41'])
#df = df.dropna(subset = ['add_stent_2_tai_yli'])
#df = df.dropna(subset = ['n_tmp_3'])
#df = df.dropna(subset = ['sten_post_100'])
#df = df.dropna(subset = ['suonia_2_tai_yli'])
#df = df.dropna(subset = ['pituus'])

#%% randomise and divive data for cross-validation

# split data

split_ratio = 0.2
training_set, holdout_set = train_test_split(df, test_size = split_ratio)
validation_set, testing_set = train_test_split(holdout_set, test_size = 0.5)
del holdout_set

# obtain sizes

n_training = training_set.shape[0]
n_validation = validation_set.shape[0]
n_testing = testing_set.shape[0]

#%% calculate correlation and standard deviation matrices

std_mat, corr_mat, most_corr = analyse_correlation(training_set, 13, 'Korjattu_DAP_GYcm2')

#%% analyse individual feature correlations

analyse_feature_correlation(training_set, 'paino', 'Korjattu_DAP_GYcm2', False)
analyse_feature_correlation(training_set, 'AHA_cto', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'Patient_sex', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'FN2BA', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'I20.81_I21.01_I21.11_or_I21.41', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'add_stent_2_tai_yli', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'n_tmp_3', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'sten_post_100', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'suonia_2_tai_yli', 'Korjattu_DAP_GYcm2', True)
analyse_feature_correlation(training_set, 'pituus', 'Korjattu_DAP_GYcm2', False)

#%% analyse target

analyse_statistics(training_set[['Korjattu_DAP_GYcm2']])

#%% replace missing values in all datasets

# create dictionary for impute values based only on training data

impute_values = {'BSA': training_set['BSA'].mean(),
                 'paino': training_set['paino'].mean(),
                 'pituus': training_set['pituus'].mean(),
                 'ind_pci_in_stemi': training_set['ind_pci_in_stemi'].mode()[0],
                 'ind_flap_failure': training_set['ind_flap_failure'].mode()[0],
                 'ind_nstemi': training_set['ind_nstemi'].mode()[0],
                 'ind_diag': training_set['ind_diag'].mode()[0],
                 'ind_uap': training_set['ind_uap'].mode()[0],
                 'ind_heart_failure': training_set['ind_heart_failure'].mode()[0],
                 'ind_stemi_other': training_set['ind_stemi_other'].mode()[0],
                 'ind_stable_ap': training_set['ind_stable_ap'].mode()[0],
                 'ind_arrhythmia_settl': training_set['ind_arrhythmia_settl'].mode()[0],
                 'suonia_2_tai_yli': training_set['suonia_2_tai_yli'].mode()[0],
                 'lm_unprotected': training_set['lm_unprotected'].mode()[0],
                 'im': training_set['im'].mode()[0],
                 'lada': training_set['lada'].mode()[0],
                 'ladb': training_set['ladb'].mode()[0],
                 'ladc': training_set['ladc'].mode()[0],
                 'lcxa': training_set['lcxa'].mode()[0],
                 'lcxb': training_set['lcxb'].mode()[0],
                 'lcxc': training_set['lcxc'].mode()[0],
                 'ld1': training_set['ld1'].mode()[0],
                 'ld2': training_set['ld2'].mode()[0],
                 'lita': training_set['lita'].mode()[0],
                 'lm': training_set['lm'].mode()[0],
                 'lom1': training_set['lom1'].mode()[0],
                 'lom2': training_set['lom2'].mode()[0],
                 'lpd': training_set['lpd'].mode()[0],
                 'lpl': training_set['lpl'].mode()[0],
                 'ram_rv': training_set['ram_rv'].mode()[0],
                 'rcaa': training_set['rcaa'].mode()[0],
                 'rcab': training_set['rcab'].mode()[0],
                 'rcac': training_set['rcac'].mode()[0],
                 'rita': training_set['rita'].mode()[0],
                 'rpd': training_set['rpd'].mode()[0],
                 'rpl': training_set['rpl'].mode()[0],
                 'vgrca_ag': training_set['vgrca_ag'].mode()[0],
                 'vglca1_ag': training_set['vglca1_ag'].mode()[0],
                 'vglca2_ag': training_set['vglca2_ag'].mode()[0],
                 'restenosis': training_set['restenosis'].mode()[0],
                 'stent_dimension': training_set['stent_dimension'].mean(),
                 'ball_dimension': training_set['ball_dimension'].mean(),
                 'add_stent_1': 0,
                 'add_stent_2_tai_yli': 0}

# combine datasets for imputing

df = training_set.append([validation_set, testing_set])

# impute data

for key, val in impute_values.items():
    df[key] = df[key].fillna(val)
del key, val

#%% fill in mutually exclusive categorical values

# obtain categorical impute values

sten_post_training = training_set[['sten_post_0', 'sten_post_25', 'sten_post_60',
                                   'sten_post_85', 'sten_post_100']].idxmax(axis = 1)
impute_values['sten_post'] = sten_post_training.mode()[0]

sten_pre_training = training_set[['sten_pre_100', 'sten_pre_85',
                                  'sten_pre_60']].idxmax(axis = 1)
impute_values['sten_pre'] = sten_pre_training.mode()[0]

AHA_training = training_set[['AHA_a', 'AHA_b1', 'AHA_b2',
                             'AHA_c', 'AHA_cto']].idxmax(axis = 1)
impute_values['AHA'] = AHA_training.mode()[0]

del sten_post_training, sten_pre_training, AHA_training

# impute data

sten_post = df[['sten_post_0', 'sten_post_25', 'sten_post_60',
                'sten_post_85', 'sten_post_100']].idxmax(axis = 1)
sten_post = sten_post.fillna(impute_values['sten_post'])
sten_post = pd.get_dummies(sten_post).astype(int)
sten_post = sten_post[['sten_post_0', 'sten_post_25', 'sten_post_60',
                       'sten_post_85', 'sten_post_100']]
df[['sten_post_0', 'sten_post_25', 'sten_post_60', 'sten_post_85',
    'sten_post_100']] = sten_post

sten_pre = df[['sten_pre_100', 'sten_pre_85', 'sten_pre_60']].idxmax(axis = 1)
sten_pre = sten_pre.fillna(impute_values['sten_pre'])
sten_pre = pd.get_dummies(sten_pre).astype(int)
sten_pre = sten_pre[['sten_pre_100', 'sten_pre_85', 'sten_pre_60']]
df[['sten_pre_100', 'sten_pre_85', 'sten_pre_60']] = sten_pre

AHA = df[['AHA_a', 'AHA_b1', 'AHA_b2', 'AHA_c', 'AHA_cto']].idxmax(axis = 1)
AHA = AHA.fillna(impute_values['AHA'])
AHA = pd.get_dummies(AHA).astype(int)
AHA = AHA[['AHA_a', 'AHA_b1', 'AHA_b2', 'AHA_c', 'AHA_cto']]
df[['AHA_a', 'AHA_b1', 'AHA_b2', 'AHA_c', 'AHA_cto']] = AHA

del sten_post, sten_pre, AHA

#%% check for nan values

df.isnull().values.any()

#%% split impute data back to training, validation and testing

training_set = df[:n_training]
validation_set = df[n_training:n_training+n_validation]
testing_set = df[-n_testing:]

#%% define feature and target labels

#feature_labels = ['AHA_b1', 'paino', 'sten_pre_100', 'Age', 'AHA_cto',
#                  'Patient_sex', 'FN1AC', 'FN2BA', 'I20.81_I21.01_I21.11_or_I21.41']

feature_labels = ['BSA', 'paino', 'pituus', 'Patient_sex', 'Age',
                  'I20.81_I21.01_I21.11_or_I21.41', 'I35.0', 'FN1AC', 'FN2BA',
                  'FN2AA', 'TFC00', 'n_tmp_1', 'n_tmp_2', 'n_tmp_3',
                  'ind_pci_in_stemi', 'ind_flap_failure', 'ind_nstemi',
                  'ind_diag', 'ind_uap', 'ind_heart_failure', 'ind_stemi_other',
                  'ind_stable_ap', 'ind_arrhythmia_settl', 'suonia_2_tai_yli',
                  'lm_unprotected', 'Aiempi_ohitusleikkaus', 'im', 'lada',
                  'ladb', 'ladc', 'lcxa', 'lcxb', 'lcxc', 'ld1', 'ld2', 'lita',
                  'lm', 'lom1', 'lom2', 'lpd', 'lpl', 'ram_rv', 'rcaa', 'rcab',
                  'rcac', 'rita', 'rpd', 'rpl', 'vgrca_ag', 'vglca1_ag',
                  'vglca2_ag', 'restenosis', 'stent_dimension', 'ball_dimension',
                  'add_stent_1', 'add_stent_2_tai_yli', 'sten_post_0',
                  'sten_post_25', 'sten_post_60', 'sten_post_85', 'sten_post_100',
                  'sten_pre_100', 'sten_pre_85', 'sten_pre_60', 'AHA_a', 'AHA_b1',
                  'AHA_b2', 'AHA_c', 'AHA_cto', 'IVUS', 'OCT']

target_label = ['Korjattu_DAP_GYcm2']

#%% extract features and targets

training_features = training_set[feature_labels]
validation_features = validation_set[feature_labels]
testing_features = testing_set[feature_labels]

training_targets = training_set[target_label]
validation_targets = validation_set[target_label]
testing_targets = testing_set[target_label]

#%% calculate sample weights

hist, bin_edges = np.histogram(training_targets, bins = 10)
classes = training_targets.apply(lambda x: pd.cut(x, bin_edges, labels = False,
                                                  include_lowest = True)).values
sample_weights = compute_sample_weight('balanced', classes)

#%% scale features

feature_transform = 'log'

if feature_transform == 'z-score':

    t_mean = training_features.mean()
    t_std = training_features.std()

    training_features = (training_features - t_mean) / t_std
    validation_features = (validation_features - t_mean) / t_std
    testing_features = (testing_features - t_mean) / t_std

if feature_transform == 'log':

    training_features = np.log1p(training_features)
    validation_features = np.log1p(validation_features)
    testing_features = np.log1p(testing_features)

if feature_transform == 'box-cox':

    lmbda = 0.15

    training_features = sp.special.boxcox1p(training_features, lmbda)
    validation_features = sp.special.boxcox1p(validation_features, lmbda)
    testing_features = sp.special.boxcox1p(testing_features, lmbda)

#%% scale targets (for skewed data)

target_transform = 'log'

if target_transform == 'log':

    training_targets = np.log1p(training_targets)
    validation_targets = np.log1p(validation_targets)
    testing_targets = np.log1p(testing_targets)

if target_transform == 'box-cox':

    lmbda = 0.15

    training_targets = sp.special.boxcox1p(training_targets, lmbda)
    validation_targets = sp.special.boxcox1p(validation_targets, lmbda)
    testing_targets = sp.special.boxcox1p(testing_targets, lmbda)

#%% build and train xgboost model

param = {
        'objective': 'reg:linear',
        'eta': 0.02,
        'eval_metric': 'rmse',
        'max_depth': 5,
        'min_child_weight': 3,
        'subsample': 0.8,
        'colsample_bytree': 0.8,
        'silent': 1,
        'seed': 123,
        'alpha': 0.0,
        'lambda': 0.01,
        }

trn = xgb.DMatrix(training_features, label = training_targets, weight = sample_weights)
vld = xgb.DMatrix(validation_features, label = validation_targets)

res = xgb.cv(param, trn, nfold = 4, num_boost_round = 2000, early_stopping_rounds = 50,
             show_stdv = True, metrics = {'rmse'}, maximize = False)

min_index = np.argmin(res['test-rmse-mean'])

evals_result = {}

timestr = time.strftime('%Y%m%d-%H%M%S')

model = xgb.train(param, trn, min_index, [(trn, 'training'), (vld,'validation')],
                  evals_result = evals_result, verbose_eval = 10)

#%% evaluate model performance

# make predictions

training_predictions = model.predict(trn)
training_predictions = pd.DataFrame(training_predictions, columns = target_label,
                                    index = training_features.index, dtype = float)

validation_predictions = model.predict(vld)
validation_predictions = pd.DataFrame(validation_predictions, columns = target_label,
                                      index = validation_features.index, dtype = float)

# calculate loss metrics

training_error = np.sqrt(mean_squared_error(training_targets, training_predictions))
validation_error = np.sqrt(mean_squared_error(validation_targets, validation_predictions))

# convert log targets to linear units (for skewed data)

if target_transform == 'log':

    training_targets_lin = np.expm1(training_targets)
    validation_targets_lin = np.expm1(validation_targets)

    training_predictions_lin = np.expm1(training_predictions)
    validation_predictions_lin = np.expm1(validation_predictions)

# convert box-cox targets to linear units (for skewed data)

if target_transform == 'box-cox':

    training_targets_lin = sp.special.inv_boxcox1p(training_targets, lmbda)
    validation_targets_lin = sp.special.inv_boxcox1p(validation_targets, lmbda)

    training_predictions_lin = sp.special.inv_boxcox1p(training_predictions, lmbda)
    validation_predictions_lin = sp.special.inv_boxcox1p(validation_predictions, lmbda)

# plot training performance

if (target_transform == 'log') or (target_transform == 'box-cox'):

    f1 = plot_regression_performance('xgboost', evals_result, training_targets_lin, training_predictions_lin,
                                     validation_targets_lin, validation_predictions_lin)
else:

    f1 = plot_regression_performance('xgboost', evals_result, training_targets, training_predictions,
                                     validation_targets, validation_predictions)

# plot feature importance

f2 = plot_feature_importance(model, training_features)

#%% save model

model_dir = 'XGBoost models\\%s_TE%d_VE%d' % (timestr,
                                              round(training_error),
                                              round(validation_error))

if not os.path.exists(model_dir):
    os.makedirs(model_dir)

f1.savefig(model_dir + '\\' + 'evaluation_metrics.pdf', dpi = 600, format = 'pdf',
                    bbox_inches = 'tight', pad_inches = 0)
f2.savefig(model_dir + '\\' + 'feature_importance.pdf', dpi = 600, format = 'pdf',
                    bbox_inches = 'tight', pad_inches = 0)

variables_to_save = {'param': param,
                     'trn': trn,
                     'vld': vld,
                     'res': res,
                     'min_index': min_index,
                     'evals_result': evals_result,
                     'nan_percent': nan_percent,
                     'duplicates': duplicates,
                     'most_corr': most_corr,
                     'corr_mat': corr_mat,
                     'std_mat': std_mat,
                     'split_ratio': split_ratio,
                     'sample_weights': sample_weights,
                     'feature_transform': feature_transform,
                     'target_transform': target_transform,
                     'timestr': timestr,
                     'model_dir': model_dir,
                     'df': df,
                     'df_orig': df_orig,
                     'impute_values': impute_values,
                     'feature_labels': feature_labels,
                     'target_label': target_label,
                     'n_training': n_training,
                     'n_validation': n_validation,
                     'n_testing': n_testing,
                     'training_set': training_set,
                     'training_features': training_features,
                     'training_targets': training_targets,
                     'validation_set': validation_set,
                     'validation_features': validation_features,
                     'validation_targets': validation_targets,
                     'testing_set': testing_set,
                     'testing_features': testing_features,
                     'testing_targets': testing_targets}

save_load_variables(model_dir, variables_to_save, 'variables', 'save')

model.save_model(model_dir + '\\' + 'xgboost_model.h5')

max_stars_repo_name: vsuomi/radiation-dose

max_forks_repo_name: vsuomi/radiation-dose


text:
import numpy as np
import SBCcode as sbc
import os
import re
from SBCcode.DataHandling.WriteBinary import WriteBinaryNtupleFile as wb
import pdb

recondir = '/bluearc/storage/recon/devel/SBC-15/output'

dd = sbc.read_bin(os.path.join(recondir, 'PMTpulseAnalysis_all.bin'))
ev_float = np.float64(dd['runid'][:, 0] - 20161000) +\
    np.float64(dd['runid'][:, 1]) * 1e-3 +\
    np.float64(dd['ev']) * 1e-6

firsthit = dd['iPMThit'] < 2
lasthit = np.ones(firsthit.shape, dtype=np.bool)
lasthit[:-1] = firsthit[1:]

dd_oneline = {k: dd[k][lasthit] for k in dd.keys()}

wb(os.path.join(recondir, 'PMTpulseAnalysis_all_oneline.bin'), [dd_oneline],
    rowdef=1, initialkeys=['runid', 'ev'], drop_first_dim=True)

max_stars_repo_name: RunzZhang/SBCcode

max_issues_repo_name: RunzZhang/SBCcode

max_forks_repo_name: RunzZhang/SBCcode


text:
"""Test the Sqlite Wrapper with the CITY ontology."""

import os
import unittest2 as unittest
import sqlite3
import numpy as np
from pathlib import Path
from osp.core.session.db.sql_migrate import check_supported_schema_version, \
    detect_current_schema_version, versions
from osp.core.session.db.sql_migrate import SqlMigrate
from osp.wrappers.sqlite import SqliteSession

try:
    from osp.core.namespaces import city
except ImportError:
    from osp.core.ontology import Parser
    from osp.core.ontology.namespace_registry import namespace_registry
    Parser().parse("city")
    city = namespace_registry.city

DB = "test_sqlite_migration.db"

class TestSqliteCity(unittest.TestCase):
    """Test the sqlite wrapper with the city ontology."""

    def tearDown(self):
        """Remove the database file."""
        if os.path.exists(DB):
            os.remove(DB)

    def run_migration(self, schema_version):
        """Load the data from the given schema version + run the migration."""
        filename = f"sqlite_schema_v{schema_version}.sql"
        with open(Path(__file__).parent / filename, encoding="utf-8") as f:
            with sqlite3.connect(DB) as con:
                con.executescript(f.read())

        self.assertRaises(RuntimeError, check_supported_schema_version,
                          SqliteSession(DB))
        self.assertRaises(RuntimeError, city.CityWrapper,
                          session=SqliteSession(DB))
        self.assertEqual(detect_current_schema_version(
            SqliteSession(DB)._get_table_names("")), schema_version
        )

        with SqliteSession(DB) as session:
            m = SqlMigrate(session)
            m.run()
            self.assertEqual(detect_current_schema_version(
                SqliteSession(DB)._get_table_names("")), max(versions.values())
            )
            self.assertTrue(check_supported_schema_version(SqliteSession(DB)))

    def run_test(self):
        """Test whether all data can be loaded correctly."""
        with SqliteSession(DB) as session:
            w = city.CityWrapper(session=session)
            cities = w.get(rel=city.hasPart)
            c = cities[0]
            self.assertEqual(len(cities), 1)
            self.assertTrue(c.is_a(city.City))
            self.assertEqual(c.name, "Freiburg")
            self.assertEqual(c.uid.hex,
                             "affb72ee61754028bd7e39a92ba3bb77")
            self.assertEqual(c.get(rel=city.isPartOf), [w])
            np.testing.assert_equal(c.coordinates, np.array([42, 12]))

            neighborhoods = c.get(oclass=city.Neighborhood)
            n = neighborhoods[0]
            self.assertEqual(len(neighborhoods), 1)
            self.assertEqual(n.uid.hex,
                             "e30e0287f52b49f396b939a85fc9460d")
            self.assertEqual(n.name, "Zähringen")
            self.assertEqual(n.get(rel=city.isPartOf), [c])
            np.testing.assert_equal(n.coordinates, np.array([0, 0]))

            streets = n.get()
            s = streets[0]
            self.assertEqual(len(streets), 1)
            self.assertEqual(s.uid.hex,
                             "25cb6116e9d04ceb81cdd8cfcbead47b")
            self.assertEqual(s.name, "Le street")
            self.assertEqual(s.get(rel=city.isPartOf), [n])
            np.testing.assert_equal(s.coordinates, np.array([1, 98]))

            citizen = c.get(rel=city.hasInhabitant)
            citizen = sorted(citizen, key=lambda p: p.name)
            self.assertEqual([p.name for p in citizen],
                             ["Hans", "Michel", "Peter"])
            for p in citizen:
                self.assertEqual(p.get(rel=city.hasInhabitant.inverse), [c])
                self.assertEqual(p.get(), [])

    def test_migrate_v0(self):
        """Test migration from schema from v0."""
        # load old schema and run migration

        # connect to db and check if
        self.run_migration(0)
        self.run_test()

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

max_stars_repo_name: simphony/osp-core

max_issues_repo_name: simphony/osp-core

max_forks_repo_name: simphony/osp-core


text:
import numpy as np
from PIL import Image
from utils import load_graph_model, get_input_tensors, get_output_tensors
import tensorflow as tf
try:
    from tflite_runtime.interpreter import Interpreter
except:
    from tensorflow.lite.python.interpreter import Interpreter

def load_frozen_graph(frozen_graph_filename):
    # We load the protobuf file from the disk and parse it to retrieve the
    # unserialized graph_def
    with tf.io.gfile.GFile(frozen_graph_filename, "rb") as f:
        graph_def = tf.compat.v1.GraphDef()
        graph_def.ParseFromString(f.read())

    # Then, we import the graph_def into a new Graph and returns it
    with tf.Graph().as_default() as graph:
        # The name var will prefix every op/nodes in your graph
        # Since we load everything in a new graph, this is not needed
        tf.import_graph_def(graph_def, name="")
    return graph

kWidth = 256
kHeight = 144
tf_model_path = "model_v2.pb"
tflite_model_path = "saved_model_144x256/model_float32.tflite"

print("Loading model...", end="")
graph = load_frozen_graph(tf_model_path)
print("done.\n", end="")

# Get input and output tensors
input_tensor_names = get_input_tensors(graph)
print(input_tensor_names)
output_tensor_names = get_output_tensors(graph)
print(output_tensor_names)

image = Image.open('test.jpg')
h = image.size[1]
w = image.size[0]
img = image.resize((kWidth, kHeight))
img = np.asarray(img)
img = img / 255.
img = img.astype(np.float32)
img = img[np.newaxis,:,:,:]

# Tensorflow Lite
interpreter = Interpreter(tflite_model_path, num_threads=4)
interpreter.allocate_tensors()
input_details = interpreter.get_input_details()[0]['index']
output_details = interpreter.get_output_details()[0]['index']

interpreter.set_tensor(input_details, img)
interpreter.invoke()
output = interpreter.get_tensor(output_details)

print(output.shape)
out1 = output[0][:, :, 0]
out2 = output[0][:, :, 1]

out1 = (out1 > 0.5) * 255
out2 = (out2 > 0.5) * 255

print('out1:', out1.shape)
print('out2:', out2.shape)

out1 = Image.fromarray(np.uint8(out1)).resize((w, h))
out2 = Image.fromarray(np.uint8(out2)).resize((w, h))

out1.save('out1.jpg')
out2.save('out2.jpg')

# evaluate the loaded model directly
print("Running inference...", end="")
input_tensor = graph.get_tensor_by_name(input_tensor_names[0])
with tf.compat.v1.Session(graph=graph) as sess:
    results = sess.run(output_tensor_names, feed_dict={
                       input_tensor: img})
print("done. {} outputs received".format(len(results)))

segments = np.squeeze(results[0], 0)
print('segments', segments.shape)
segments1 = segments[:, :, 0]
segments2 = segments[:, :, 1]
segments1 = (segments1 > 0.5) * 255
segments2 = (segments2 > 0.5) * 255

segments1 = Image.fromarray(np.uint8(segments1)).resize((w, h))
segments2 = Image.fromarray(np.uint8(segments2)).resize((w, h))
segments1.save('segments1.jpg')
segments2.save('segments2.jpg')

mask_img = segments2.resize((w, h), Image.LANCZOS).convert("RGB")
mask_img = tf.keras.preprocessing.image.img_to_array(mask_img, dtype=np.uint8)
fg = np.bitwise_and(np.array(image), np.array(mask_img))
fg_img = Image.fromarray(np.uint8(fg)).resize((w, h))
fg_img.save('fg.jpg')


text:
import time
import numpy as np

from qunetsim.components import Host
from qunetsim.components import Network
from qunetsim.objects import Logger
from qunetsim.objects import Qubit

WAIT_TIME = 10

def checksum_sender(host, q_size, receiver_id, checksum_size_per_qubit):
    bit_arr = np.random.randint(2, size=q_size)
    Logger.get_instance().log('Bit array to be sent: ' + str(bit_arr))
    qubits = []
    for i in range(q_size):
        q_tmp = Qubit(host)
        if bit_arr[i] == 1:
            q_tmp.X()
        qubits.append(q_tmp)

    check_qubits = host.add_checksum(qubits, checksum_size_per_qubit)
    checksum_size = int(q_size / checksum_size_per_qubit)
    qubits.append(check_qubits)
    checksum_cnt = 0
    for i in range(q_size + checksum_size):
        if i < q_size:
            q = qubits[i]
        else:
            q = qubits[q_size][checksum_cnt]
            checksum_cnt = checksum_cnt + 1

        host.send_qubit(receiver_id, q, await_ack=True)

def checksum_receiver(host, q_size, sender_id, checksum_size_per_qubit):
    qubits = []
    checksum_size = int(q_size / checksum_size_per_qubit)
    while len(qubits) < (q_size + checksum_size):
        q = host.get_data_qubit(sender_id, wait=WAIT_TIME)
        qubits.append(q)
        Logger.get_instance().log(str(host.host_id) + ': received qubit')

    checksum_qubits = []
    checksum_cnt = 0

    for i in range(len(qubits)):
        if checksum_cnt < checksum_size:
            checksum_qubits.append(qubits[q_size + i]['q'])
            checksum_cnt = checksum_cnt + 1

    checksum_cnt = 1
    for i in range(len(qubits) - checksum_size):
        qubits[i]['q'].cnot(checksum_qubits[checksum_cnt - 1])
        if i == (checksum_cnt * checksum_size_per_qubit - 1):
            checksum_cnt = checksum_cnt + 1

    errors = 0
    for i in range(len(checksum_qubits)):
        if checksum_qubits[i].measure() != 0:
            errors += 1

    print('---------')
    if errors == 0:
        Logger.get_instance().log('No error exist in UDP packet')
    else:
        Logger.get_instance().log('There were errors in the UDP transmission')
    print('---------')

    rec_bits = []
    for i in range(len(qubits) - checksum_size):
        rec_bits.append(qubits[i]['q'].measure())

    if errors == 0:
        print('---------')
        Logger.get_instance().log('Receiver received the classical bits: ' + str(rec_bits))
        print('---------')
        return True

    return

def main():
    network = Network.get_instance()
    nodes = ["Alice", "Bob", "Eve", "Dean"]
    network.start(nodes)
    network.delay = 0.5

    host_alice = Host('Alice')
    host_alice.add_connection('Bob')
    host_alice.max_ack_wait = 30
    host_alice.delay = 0.2
    host_alice.start()

    host_bob = Host('Bob')
    host_bob.max_ack_wait = 30
    host_bob.delay = 0.2
    host_bob.add_connection('Alice')
    host_bob.add_connection('Eve')
    host_bob.start()

    host_eve = Host('Eve')
    host_eve.max_ack_wait = 30
    host_eve.delay = 0.2
    host_eve.add_connection('Bob')
    host_eve.add_connection('Dean')
    host_eve.start()

    host_dean = Host('Dean')
    host_dean.max_ack_wait = 30
    host_dean.delay = 0.2
    host_dean.add_connection('Eve')
    host_dean.start()

    network.add_host(host_alice)
    network.add_host(host_bob)
    network.add_host(host_eve)
    network.add_host(host_dean)

    network.x_error_rate = 0
    network.packet_drop_rate = 0

    q_size = 6
    checksum_per_qubit = 2

    host_alice.run_protocol(checksum_sender, (q_size, host_dean.host_id, checksum_per_qubit))
    host_dean.run_protocol(checksum_receiver, (q_size, host_alice.host_id, checksum_per_qubit))

    start_time = time.time()
    while time.time() - start_time < 150:
        pass

    network.stop(stop_hosts=True)
    exit()

if __name__ == '__main__':
    main()