text:
import click
import pandas as pd
import numpy as np
from collections import namedtuple

from scipy.interpolate import LinearNDInterpolator
from scipy.stats import gaussian_kde
from sklearn.mixture import GaussianMixture
from scipy.stats import rankdata

from pyprophet.stats import pnorm, pemp, pi0est #, lfdr, qvalue
# from pyprophet.stats import stat_metrics
from pyprophet.stats import find_nearest_matches
from pyprophet.stats import posterior_chromatogram_hypotheses_fast

from pyprophet.ipf import compute_model_fdr
from pyprophet.optimized import count_num_positives

DensityModel = namedtuple('DensityModel', ('model', 'evaluate'))

class ErrorStatisticsCalculator:
    def __init__(self, scored_table,
                 density_estimator='gmm',
                 grid_size=256,
                 parametric=False,
                 pfdr=False,
                 tric_chromprob=False,
                 **kwargs):
        self.grid_size = grid_size
        self.parametric = parametric
        self.pi0_args = {
            (k if k == 'pi0_method' \
             else ('lambda_' if k == 'pi0_lambda' else k[len('pi0_'):])
            ): v
            for k, v in kwargs.items()
            if k.startswith('pi0_')
        }
        self.lfdr_args = {
            ('trunc' if k == 'lfdr_truncate' \
             else ('transf' if k == 'lfdr_transformation' else k[len('lfdr_'):])
            ): v
            for k, v in kwargs.items()
            if k.startswith('lfdr_')
        }
        self.pfdr = pfdr
        self.tric_chromprob = tric_chromprob

        self.scored_table = scored_table

        if density_estimator == 'kde':
            def density_model_func(scores, decoy_peptide=0, decoy_glycan=0):
                model = gaussian_kde(scores.T)
                return DensityModel(
                    model=model,
                    evaluate=lambda scores: model.evaluate(scores.T)
                )
        elif density_estimator == 'gmm':
            def density_model_func(scores, decoy_peptide=0, decoy_glycan=0):
                if decoy_peptide and decoy_glycan:
                    n_components = 1
                elif not decoy_peptide and not decoy_glycan:
                    n_components = 4
                else:
                    n_components = 2
                model = GaussianMixture(n_components=n_components)
                model.fit(scores)
                return DensityModel(
                    model=model,
                    evaluate=lambda scores: np.exp(model.score_samples(scores))
                )
        else:
            raise ValueError('invalid density_estimator: ' + str(density_estimator))
        self.density_model_func = density_model_func

    def get_top_scores(self, decoy_peptide, decoy_glycan,
                       score=['d_score_peptide', 'd_score_glycan', 'd_score_combined']):
        row = self.scored_table.apply(lambda x: True, axis=1)
        if 'peak_group_rank' in self.scored_table.columns:
            row &= self.scored_table['peak_group_rank'] == 1
        if decoy_peptide:
            row &= self.scored_table['decoy_peptide'] == 1
        elif decoy_peptide is not None:
            row &= self.scored_table['decoy_peptide'] == 0
        if decoy_glycan:
            row &= self.scored_table['decoy_glycan'] == 1
        elif decoy_glycan is not None:
            row &= self.scored_table['decoy_glycan'] == 0

        scores = self.scored_table.loc[row, score]
        scores = scores.dropna()
        return scores.values

    def calculate_partial_stat(self, part):
        pvalue = pnorm if self.parametric else pemp

        if part == 'both':
            score_part = 'd_score_combined'
            scores_target = self.get_top_scores(
                decoy_peptide=0,
                decoy_glycan=0,
                score=score_part
            )
            scores_decoy = self.get_top_scores(
                decoy_peptide=1,
                decoy_glycan=1,
                score=score_part
            )
        elif part == 'peptide':
            score_part = 'd_score_' + part
            scores_target = self.get_top_scores(
                decoy_peptide=0,
                decoy_glycan=None,
                score=score_part
            )
            scores_decoy = self.get_top_scores(
                decoy_peptide=1,
                decoy_glycan=None,
                score=score_part
            )
        elif part == 'glycan':
            score_part = 'd_score_' + part
            scores_target = self.get_top_scores(
                decoy_peptide=None,
                decoy_glycan=0,
                score=score_part
            )
            scores_decoy = self.get_top_scores(
                decoy_peptide=None,
                decoy_glycan=1,
                score=score_part
            )
        else:
            raise click.ClickException("Unspecified scoring part selected.")

        p_value = pvalue(scores_target, scores_decoy)
        pi0 = pi0est(p_value, **self.pi0_args)
        # pep = lfdr(p_value, pi0=pi0['pi0'], **self.lfdr_args)

        stat = pd.DataFrame(scores_target, columns=[score_part])
        stat['p_value_' + part] = p_value
        # stat['pep_' + part] = pep

        if not hasattr(self, 'stats'):
            self.stats = {}
        self.stats[part] = stat

        if not hasattr(self, 'pi0'):
            self.pi0 = {}
        self.pi0[part] = pi0

        return stat, pi0

    def fit_density_models(self):
        decoy_dict = {
            'target': (0, 0),
            'decoy_peptide': (1, 0),
            'decoy_glycan': (0, 1),
            'decoy_both': (1, 1)
        }

        density_models = {}

        for decoy_type, decoy_args in decoy_dict.items():
            scores = self.get_top_scores(
                decoy_peptide=decoy_args[0],
                decoy_glycan=decoy_args[1],
                score=['d_score_peptide', 'd_score_glycan']
            )

            model = self.density_model_func(
                scores,
                decoy_peptide=decoy_args[0],
                decoy_glycan=decoy_args[1]
            )

            density_models[decoy_type] = model

        self.density_models = density_models

        return density_models

    def create_score_grid(self):
        def get_nonoutlier_range(x, lower=True, upper=False):
            q3, q1 = np.percentile(x, [75, 25])
            iqr = q3 - q1
            min_ = np.min(x)
            max_ = np.max(x)
            if lower:
                min_ = np.max([min_, q1 - 1.5 * iqr])
            if upper:
                max_ = np.min([max_, q3 + 1.5 * iqr])
            return (min_, max_)

        def get_grid_cutoffs(total_range, nonoutlier_range, num_cutoffs):
            margin = (nonoutlier_range[1] - nonoutlier_range[0]) * 0.05
            if total_range[0] < nonoutlier_range[0] and \
                total_range[1] > nonoutlier_range[0]:
                cutoffs = np.concatenate((
                    [total_range[0] - margin],
                    np.linspace(
                        nonoutlier_range[0] - margin,
                        nonoutlier_range[1] + margin,
                        num_cutoffs - 2
                    ),
                    [total_range[1] + margin]
                ))
            elif total_range[0] < nonoutlier_range[0]:
                cutoffs = np.concatenate((
                    [total_range[0] - margin],
                    np.linspace(
                        nonoutlier_range[0] - margin,
                        nonoutlier_range[1] + margin,
                        num_cutoffs - 1
                    )
                ))
            elif total_range[1] > nonoutlier_range[0]:
                cutoffs = np.concatenate((
                    np.linspace(
                        nonoutlier_range[0] - margin,
                        nonoutlier_range[1] + margin,
                        num_cutoffs - 1
                    ),
                    [total_range[1] + margin]
                ))
            else:
                cutoffs = np.linspace(
                    nonoutlier_range[0] - margin,
                    nonoutlier_range[1] + margin,
                    num_cutoffs - 1
                )
            return cutoffs

        decoy_dict = {
            'target': (0, 0),
            'decoy_peptide': (1, 0),
            'decoy_glycan': (0, 1),
            'decoy_both': (1, 1)
        }
        scores = {
            decoy_type: self.get_top_scores(
                decoy_peptide=decoy_args[0],
                decoy_glycan=decoy_args[1],
                score=['d_score_peptide', 'd_score_glycan']
            )
            for decoy_type, decoy_args in decoy_dict.items()
        }
        nonoutlier_ranges = np.array([
            (get_nonoutlier_range(s[:, 0]), get_nonoutlier_range(s[:, 1]))
            for decoy_type, s in scores.items()
        ])

        x_cutoffs = get_grid_cutoffs(
            nonoutlier_range=(
                np.min(nonoutlier_ranges[:, 0, 0]),
                np.max(nonoutlier_ranges[:, 0, 1])
            ),
            total_range=(
                self.scored_table['d_score_peptide'].min(),
                self.scored_table['d_score_peptide'].max()
            ),
            num_cutoffs=self.grid_size[0] \
                if isinstance(self.grid_size, tuple) or \
                    isinstance(self.grid_size, list) \
                else self.grid_size
        )
        y_cutoffs = get_grid_cutoffs(
            nonoutlier_range=(
                np.min(nonoutlier_ranges[:, 1, 0]),
                np.max(nonoutlier_ranges[:, 1, 1])
            ),
            total_range=(
                self.scored_table['d_score_glycan'].min(),
                self.scored_table['d_score_glycan'].max()
            ),
            num_cutoffs=self.grid_size[1] \
                if isinstance(self.grid_size, tuple) or \
                    isinstance(self.grid_size, list) \
                else self.grid_size
        )
        return np.meshgrid(x_cutoffs, y_cutoffs)

    def calculate_total_stat(self):
        pi00 = self.pi0['both']['pi0']
        pi01 = self.pi0['peptide']['pi0'] - pi00
        pi10 = self.pi0['glycan']['pi0'] - pi00
        pi11 = 1 - pi01 - pi10 - pi00

        if pi01 < 0:
            # raise ValueError('peptide pi0 < both pi0')
            pi01 = 0
        if pi10 < 0:
            # raise ValueError('glycan pi0 < both pi0')
            pi10 = 0
        if pi11 < 0:
            # raise ValueError('peptide pi0 + glycan pi0 - both pi0 > 1')
            pi11 = 0

        if not hasattr(self, 'density_models'):
            self.fit_density_models()

        X, Y = self.create_score_grid()
        scores = np.column_stack((X.ravel(), Y.ravel()))
        stat = pd.DataFrame(
            scores,
            columns=['d_score_peptide', 'd_score_glycan']
        )
        for decoy_type, model in self.density_models.items():
            stat['density_' + decoy_type] = model.evaluate(scores)

        stat['density_peptide_null_glycan_nonnull'] = np.maximum(
            (stat['density_decoy_peptide'].values - \
             (pi00 + pi10) * stat['density_decoy_both'].values) \
            / (pi01 + pi11), 0
        )
        stat['density_peptide_nonnull_glycan_null'] = np.maximum(
            (stat['density_decoy_glycan'].values - \
             (pi00 + pi01) * stat['density_decoy_both'].values) \
            / (pi10 + pi11), 0
        )
        stat['density_nonnull'] = np.maximum(
            (stat['density_target'].values - \
             pi00 * stat['density_decoy_both'].values - \
             pi01 * stat['density_peptide_null_glycan_nonnull'].values - \
             pi10 * stat['density_peptide_nonnull_glycan_null'].values) \
            / pi11, 0
        )

        pep_peptide = np.divide(
            pi00 * stat['density_decoy_both'].values + \
                pi01 * stat['density_peptide_null_glycan_nonnull'].values,
            stat['density_target'].values,
            out=np.zeros_like(stat['density_target'].values),
            where=stat['density_target'].values != 0
        )
        pep_glycan = np.divide(
            pi00 * stat['density_decoy_both'].values + \
                pi10 * stat['density_peptide_nonnull_glycan_null'].values,
            stat['density_target'].values,
            out=np.zeros_like(stat['density_target'].values),
            where=stat['density_target'].values != 0
        )
        pep_both = np.divide(
            pi00 * stat['density_decoy_both'].values,
            stat['density_target'].values,
            out=np.zeros_like(stat['density_target'].values),
            where=stat['density_target'].values != 0
        )
        pep_total = pep_peptide + pep_glycan - pep_both

        if self.lfdr_args.get('trunc', True):
            pep_peptide = np.minimum(pep_peptide, 1)
            pep_glycan = np.minimum(pep_glycan, 1)
            pep_both = np.minimum(pep_both, 1)
            pep_total = np.minimum(pep_total, 1)

        def monotonize(values, i0, j0, i_ascending=False, j_ascending=False):
            if i_ascending:
                for j in range(values.shape[1]):
                    for i in range(i0, 0, -1):
                        if values[i - 1, j] > values[i, j]:
                            values[i - 1, j] = values[i, j]
                    for i in range(i0, values.shape[0] - 1):
                        if values[i + 1, j] < values[i, j]:
                            values[i + 1, j] = values[i, j]
            elif i_ascending is not None:
                for j in range(values.shape[1]):
                    for i in range(i0, 0, -1):
                        if values[i - 1, j] < values[i, j]:
                            values[i - 1, j] = values[i, j]
                    for i in range(i0, values.shape[0] - 1):
                        if values[i + 1, j] > values[i, j]:
                            values[i + 1, j] = values[i, j]
            if j_ascending:
                for i in range(values.shape[0]):
                    for j in range(j0, 0, -1):
                        if values[i, j - 1] > values[i, j]:
                            values[i, j - 1] = values[i, j]
                    for i in range(j0, values.shape[1] - 1):
                        if values[i, j + 1] < values[i, j]:
                            values[i, j + 1] = values[i, j]
            elif j_ascending is not None:
                for i in range(values.shape[0]):
                    for j in range(j0, 0, -1):
                        if values[i, j - 1] < values[i, j]:
                            values[i, j - 1] = values[i, j]
                    for j in range(j0, values.shape[1] - 1):
                        if values[i, j + 1] > values[i, j]:
                            values[i, j + 1] = values[i, j]
            return values

        if self.lfdr_args.get('monotone', True):
            i0, j0 = np.unravel_index(
                np.argmax(stat['density_nonnull'].values),
                X.shape
            )
            monotonize(pep_peptide.reshape(X.shape), i0=i0, j0=j0)
            monotonize(pep_glycan.reshape(X.shape), i0=i0, j0=j0)
            monotonize(pep_both.reshape(X.shape), i0=i0, j0=j0)
            monotonize(pep_total.reshape(X.shape), i0=i0, j0=j0)

        stat['pep_both'] = pep_both
        stat['pep_peptide'] = pep_peptide
        stat['pep_glycan'] = pep_glycan
        stat['pep'] = pep_total

        def qvalue(pep, density, X, Y):
            dX = np.diff(X, axis=1)
            dX = np.column_stack((dX, np.min(dX, axis=1)))
            dY = np.diff(Y, axis=0)
            dY = np.row_stack((dY, np.min(dY, axis=0)))
            dS = dX * dY
            ds = dS.ravel()

            pep = pep.ravel()
            density = density.ravel()
            order = np.argsort(pep)
            frac_fp = np.cumsum(pep[order] * density[order] * ds[order])
            frac_positive = np.cumsum(density[order] * ds[order])
            fdr_ordered = np.divide(
                frac_fp,
                frac_positive,
                out=np.zeros_like(pep),
                where=frac_positive != 0
            )
            ranks = rankdata(pep, method='max')
            fdr = fdr_ordered[ranks - 1]
            return fdr

        stat['q_value'] = qvalue(
            stat['pep'].values,
            stat['density_target'].values,
            X, Y
        )

        if not hasattr(self, 'stats'):
            self.stats = {}
        self.stats['total'] = stat

        return stat

    def lookup_stat_table(self, scored_table, stat, by, value=None):
        if value is None:
            value = [
                col for col in stat.columns
                if col not in scored_table.columns
            ]
        if isinstance(value, str):
            value = [value]

        row_na = scored_table[by].isna()
        idx = find_nearest_matches(
            np.float32(stat[by].values),
            np.float32(scored_table[by].loc[~row_na].values)
        )

        for col in value:
            scored_table[col] = scored_table[by].map(lambda x: np.nan)
            scored_table.loc[~row_na, col] = \
                stat[col].iloc[idx].values

#        stat = stat.drop_duplicates(subset=by).sort_values(by=by)
#        for col in value:
#            f = np.interp(
#                scored_table[by].loc[~row_na].values,
#                stat[by].values, stat[col].values,
#            )
#            scored_table[col] = scored_table[by].map(lambda x: np.nan)
#            scored_table.loc[~row_na, col] = f

        return scored_table

    def lookup_stat_grid(self, scored_table, stat,
                         by_x='d_score_peptide', by_y='d_score_glycan',
                         value=None):
        if value is None:
            value = [
                col for col in stat.columns
                if col not in scored_table.columns
            ]
        if isinstance(value, str):
            value = [value]

        row_na = scored_table[by_x].isna() | \
            scored_table[by_y].isna()
        for col in value:
            ip = LinearNDInterpolator(
                points=stat[[by_x, by_y]].values,
                values=stat[col].values
            )
            f = ip(
                scored_table.loc[~row_na, [by_x, by_y]].values
            )
            scored_table[col] = scored_table[by_x] \
                .map(lambda x: np.nan)
            scored_table.loc[~row_na, col] = f
        return scored_table

    def lookup_error_stat(self, scored_table, stats):
        for part in ['peptide', 'glycan', 'both']:
            if part == 'both':
                score_part = 'd_score_combined'
            else:
                score_part = 'd_score_' + part
            stat = stats.get(part, None)
            if stat is not None:
                scored_table = self.lookup_stat_table(
                    scored_table, stat,
                    by=score_part
                )

        stat = stats.get('total', None)
        if stat is not None:
            scored_table = self.lookup_stat_grid(
                scored_table, stat,
                value=['pep_peptide', 'pep_glycan', 'pep_both', 'pep']
            )

        return scored_table

    def qvalue(self, scored_table):
        target = (scored_table['decoy_peptide'] == 0) & \
            (scored_table['decoy_glycan'] == 0)
        if 'peak_group_rank' in scored_table.columns:
            target &= scored_table['peak_group_rank'] == 1

        qvalue = scored_table['pep'].map(lambda x: 0.0)
        qvalue.loc[~target] = self.lookup_stat_grid(
            scored_table.loc[~target].copy(), self.stats['total'],
            value=['q_value']
        )['q_value']

        qvalue.loc[target] = compute_model_fdr(scored_table.loc[target, 'pep'])

        scored_table['q_value'] = qvalue
        return scored_table

    def chromatogram_probabilities(self, scored_table):
        pi0 = self.pi0['peptide']['pi0'] + self.pi0['glycan']['pi0'] - \
            self.pi0['both']['pi0']

        texp = namedtuple('DummyExperiment', ['df'])(
            df = scored_table[['group_id', 'precursor_id', 'pep']] \
                .rename(columns={
                    'precursor_id': 'tg_num_id'
                })
        )
        allhypothesis, h0 = posterior_chromatogram_hypotheses_fast(
            texp, pi0
        )
        scored_table['h_score'] = allhypothesis
        scored_table['h0_score'] = h0
        return scored_table

    def calculate_stat_metrics(self):
        pi0_both = self.pi0['both']['pi0']
        pi0_peptide = self.pi0['peptide']['pi0']
        pi0_glycan = self.pi0['glycan']['pi0']
        pi0_total = pi0_peptide + pi0_glycan - pi0_both

#        decoy_dict = {
#            'target': (0, 0),
#            'decoy_peptide': (1, 0),
#            'decoy_glycan': (0, 1),
#            'decoy_both': (1, 1)
#        }
#        scores = {
#            decoy_type: -self.get_top_scores(
#                decoy_peptide=decoy_args[0],
#                decoy_glycan=decoy_args[1],
#                score='pep'
#            )
#            for decoy_type, decoy_args in decoy_dict.items()
#        }
#
#        pvalue = pnorm if self.parametric else pemp
#        fpr_both = pvalue(scores['target'], scores['decoy_both'])
#        fpr_peptide = pvalue(scores['target'], scores['decoy_peptide'])
#        fpr_glycan = pvalue(scores['target'], scores['decoy_glycan'])
#        fpr = (fpr_peptide * pi0_peptide + fpr_glycan * pi0_glycan - \
#            fpr_both * pi0_both) / pi0_total
#
#        stat = pd.DataFrame(
#            np.column_stack((-scores['target'], fpr)),
#            columns=['pep', 'p_value']
#        )
#        stat['q_value'] = compute_model_fdr(stat['pep'].values)
#        stat.sort_values(by='pep', ascending=False, inplace=True)
#        stat.reset_index(drop=True, inplace=True)
#        metrics = stat_metrics(stat['p_value'].values, pi0=pi0_total, pfdr=self.pfdr)
#        stat = pd.concat(
#            (stat, metrics),
#            axis=1
#        )

        stat = pd.DataFrame(
            self.get_top_scores(
                decoy_peptide=0,
                decoy_glycan=0,
                score=['pep', 'q_value']
            ),
            columns=['pep', 'q_value']
        )
        stat.sort_values(by='pep', ascending=False, inplace=True)
        stat.reset_index(drop=True, inplace=True)

        num_total = len(stat['pep'].values)
        num_positives = count_num_positives(stat['pep'].values)
        num_negatives = num_total - num_positives
        num_null = pi0_total * num_total

        fp = num_positives * stat['q_value'].values
        tp = num_positives - fp
        fpr = np.minimum(fp / num_null, 1.0)
        tn = num_null * (1.0 - fpr)
        fn = num_negatives - num_null * (1.0 - fpr)

        fdr = np.divide(fp, num_positives, out=np.zeros_like(fp), where=num_positives!=0)
        fnr = np.divide(fn, num_negatives, out=np.zeros_like(fn), where=num_negatives!=0)

        if self.pfdr:
            fdr /= (1.0 - (1.0 - fpr) ** num_total)
            fdr[fpr == 0] = 1.0 / num_total

            fnr /= 1.0 - fpr ** num_total
            fnr[fpr == 0] = 1.0 / num_total

        sens = tp / (num_total - num_null)
        sens[sens < 0.0] = 0.0
        sens[sens > 1.0] = 1.0

        fdr[fdr < 0.0] = 0.0
        fdr[fdr > 1.0] = 1.0
        fdr[num_positives == 0] = 0.0

        fnr[fnr < 0.0] = 0.0
        fnr[fnr > 1.0] = 1.0
        fnr[num_positives == 0] = 0.0

        svalues = pd.Series(sens)[::-1].cummax()[::-1]
        metrics = pd.DataFrame({
            'tp': tp, 'fp': fp, 'tn': tn, 'fn': fn,
            'fpr': fpr, 'fdr': fdr, 'fnr': fnr,
            'svalue': svalues
        })
        stat = pd.concat(
            (stat, metrics),
            axis=1
        )

        self.stat_metrics = stat
        return stat

    def final_error_stat(self):
        stat = self.lookup_stat_table(
            self.stats['total'].copy(),
            self.stat_metrics,
            by='pep'
        )
        return stat

    def summary_error_stat(self,
                           qvalues=[0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5]):
        qvalues = np.array(qvalues)
        idx = find_nearest_matches(
            np.float32(self.stat_metrics['q_value'].values),
            np.float32(qvalues)
        )

        stat_sub = self.stat_metrics.iloc[idx].copy()
        for i_sub, (i0, i1) in enumerate(zip(idx, idx[1:])):
            if i1 == i0:
                stat_sub.iloc[i_sub + 1, :] = None
        stat_sub['q_value'] = qvalues
        stat_sub.reset_index(inplace=True, drop=True)

        summary_stat = stat_sub[[
            'q_value', 'pep',
            'svalue',
            'fdr', 'fnr', 'fpr',
            'tp', 'tn', 'fp', 'fn'
        ]]
        return summary_stat

    def error_statistics(self):
        for part in ['peptide', 'glycan', 'both']:
            self.calculate_partial_stat(part)

        self.calculate_total_stat()

        self.scored_table = self.lookup_error_stat(self.scored_table, self.stats)
        self.scored_table = self.qvalue(self.scored_table)

        if self.tric_chromprob:
            self.scored_table = \
                self.chromatogram_probabilities(self.scored_table)

        self.calculate_stat_metrics()
        final_statistics =  self.final_error_stat()
        summary_statistics = self.summary_error_stat()

        result = {
            'scored_table': self.scored_table,
            'final_statistics': final_statistics,
            'summary_statistics': summary_statistics
        }
        pi0 = self.pi0

        return result, pi0

max_stars_repo_name: lmsac/GproDIA

max_issues_repo_name: lmsac/GproDIA

max_forks_repo_name: lmsac/GproDIA


text:
import fileinput
import numpy as np
import scipy

import cv2
import os.path

def main():
    for file in fileinput.input():
        # Get filename.
        filename = fileinput.filename()

        # Open current file as a grayscale image.
        img = cv2.imread(filename, 0)
        print 'Opened ' + filename

        # Split filename into prefix and extension.
        (name_pref, name_ext) = os.path.splitext(filename)

        # Bilateral smoothing filter.
        img = cv2.bilateralFilter(img,9,75,75)

        # Write image to file.
        newname = name_pref + '_smooth=bilat' + name_ext
        cv2.imwrite(newname, img)

        fileinput.nextfile()

    fileinput.close()

if __name__ == '__main__':
    main()

ext: py

max_stars_repo_name: stcastle/shell-detection

max_issues_repo_name: SamTCastle/shell-detection

max_forks_repo_name: SamTCastle/shell-detection

reason: import numpy,import scipy


text:
import numpy as np

from common.activations import sigmoid

from .decision_tree import decision_tree, traverse
from .regression import logistic_regression

def binary_classifier(X, y, threshold=0.5, **kwargs):
    """Trains a binary classifier from the dataset provided using logistic regression.

    In order to classify new data, use:
        y_hat = model(X)

    Args:
        X (ndarray[m, n]):              ground features
        y (ndarray[m, 1]):              ground labels
        threshold (float, optional):    default classification threshold, can be
                                        overridden if passed as argument to the model
                                        (defaults to 0.5)
        **kwargs (...):                 logistic regression arguments, see documentation

    Returns:
        model (X -> y_hat):             callable binary classifier
    """
    classes = np.unique(y)
    W, b = logistic_regression(X, y, **kwargs)
    def model(X_pred, threshold=threshold):
        indices = (sigmoid(np.dot(X_pred, W) + b) >= threshold).astype(int)
        return classes[indices].reshape(-1, 1)
    return model

def multinomial_classifier(X, y, **kwargs):
    """Trains a multi-class classifier from the dataset
    provided using a set of logistic regression classifiers.

    In order to classify new data, use:
        y_hat = model(X)

    Args:
        X (ndarray[m, n]):              ground features
        y (ndarray[m, 1]):              ground labels
        **kwargs (...):                 logistic regression arguments, see documentation

    Returns:
        model (X -> y_hat):             callable multi-class classifier
    """
    classes = np.unique(y)

    # Build logistic regression models
    lr_params = []
    for cls in classes:
        y_model = (y == cls).astype(int).reshape(-1, 1)
        W, b = logistic_regression(X, y_model, **kwargs)
        lr_params.append((W, b))

    def model(X_pred):
        """Evaluates each logistic regression model, picking the label
        supplied by the most confident model for each input.
        """
        activations = [sigmoid(np.dot(X_pred, W) + b) for W, b in lr_params]
        indices = np.argmax(np.hstack(activations), axis=1)
        return classes[indices].reshape(-1, 1)
    return model

def knn_classifier(X, y, default_k=3):
    """Builds a K-nearest neighbour classifier from the dataset provided.

    In order to classify new data, use:
        y_hat = model(X)

    Args:
        X (ndarray[m, n]):      ground features
        y (ndarray[m, 1]):      ground labels
        k (int, optional):      number of neighbours to consider (defaults to 3)

    Returns:
        model (X -> y_hat):     callable knn classifier
    """
    classes = np.unique(y)

    def model(X_pred, k=default_k):
        # Euclidean distances from every point to every labeled point
        distances = np.linalg.norm(X - X_pred[:, np.newaxis, :], axis=-1)

        # Indices of k nearest neighbours (np.argpartition yields indices
        # s.t. first k elements along axis are ordered)
        indices = np.argpartition(distances, k, axis=-1)[..., :k]

        # Classes of k nearest neighbours
        k_nearest = np.take(y, indices)

        # Find class index of most frequent label among k nearest neighbours
        mode_indices = np.apply_along_axis(lambda x: np.bincount(x).argmax(), 1, k_nearest)
        return classes[mode_indices].reshape(-1, 1)
    return model

def decision_tree_classifier(X, y, attrs=None):
    """Builds a decision tree classifier from the dataset provided.

    In order to classify new data, use:
        y_hat = model(X)

    Args:
        X (ndarray[m, n]):          ground features
        y (ndarray[m, 1]):          ground labels
        attrs (tuple, optional):    tuple of attributes to branch on (defaults to all)

    Returns:
        model (X -> y_hat):         callable decision tree classifier
    """
    tree = decision_tree(X, y, attrs=attrs)
    def model(X_pred):
        return np.apply_along_axis(traverse, axis=0, arr=X_pred, args=(tree, ))
    return model

def svm(X, y):
    # Delegates convex problem optimization to LIBSVM via Sklearn.
    pass

max_stars_repo_name: patrikkj/algorithms

max_issues_repo_name: patrikkj/algorithms

max_forks_repo_name: patrikkj/algorithms

reason: import numpy


text:
#  Copyright 2021 Amazon.com, Inc. or its affiliates. All Rights Reserved.
#
#  Licensed under the Apache License, Version 2.0 (the "License"). You may not use this file except in compliance
#  with the License. A copy of the License is located at
#
#  http://aws.amazon.com/apache2.0/
#
#  or in the "license" file accompanying this file. This file is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES
#  OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions
#  and limitations under the License.
import json
import logging
import multiprocessing as mp
import os
from itertools import chain, repeat

import numpy as np
import scipy.sparse as smat
from pecos.utils import smat_util
from pecos.core import clib
from sklearn.preprocessing import normalize as sk_normalize

LOGGER = logging.getLogger(__name__)

class ClusterChain(object):
    """Utility class for representing a hierarchical clustering as a list of CSC matrices.

    Attributes:
        chain (list): A list of sparse matrices that form a hierarchical clustering chain.
    """

    def __init__(self, chain):
        """Initialization

        Args:
            chain (list): A list of sparse matrices that form a valid hierarchical clustering chain.
        """

        if isinstance(chain, type(self)):
            chain = chain.chain
        elif isinstance(chain, smat.spmatrix):
            chain = [chain.tocsc()]

        # check the validity of the clustering chain
        assert isinstance(chain, list), "clustering chain shall be a list of CSC matrices"
        for i in range(len(chain) - 1):
            assert (
                chain[i].shape[0] == chain[i + 1].shape[1]
            ), "matrices do not form a valid hierarchical clustering chain"

        self.chain = chain

    def __delitem__(self, key):
        del self.chain[key]

    def __getitem__(self, key):
        return self.chain[key]

    def __setitem__(self, key, val):
        self.chain[key] = val

    def __iter__(self):
        return iter(self.chain)

    def __len__(self):
        return len(self.chain)

    def __eq__(self, other):
        """
        Returns:
            True if `ClusterChain`s are of same length and their matrices have the same shapes and number of nonzeros, False otherwise.
        """

        if len(self) != len(other):
            return False

        for C_A, C_B in zip(self, other):
            if C_A.shape != C_B.shape or (C_A != C_B).nnz != 0:
                return False
        return True

    def save(self, folder):
        """Save to disk.

        Args:
            folder (str): Folder to save to.
        """

        os.makedirs(folder, exist_ok=True)
        with open(os.path.join(folder, "config.json"), "w", encoding="utf-8") as fout:
            fout.write(json.dumps({"len": len(self)}))

        for i, C in enumerate(self):
            smat_util.save_matrix(os.path.join(folder, f"C{i}.npz"), C)

    @classmethod
    def load(cls, path_to_cluster):
        """Load from disk.

        Args:
            path_to_cluster (str): Folder where `ClusterChain` was saved to using `ClusterChain.save`.

        Returns:
            ClusterChain: The loaded object.
        """

        if os.path.isfile(path_to_cluster):
            C = smat_util.load_matrix(path_to_cluster)
            return cls.from_partial_chain(C)

        config_path = os.path.join(path_to_cluster, "config.json")
        if not os.path.exists(config_path):
            raise ValueError(f"Cluster config file, {config_path}, does not exist")

        with open(config_path, "r", encoding="utf-8") as fin:
            config = json.loads(fin.read())
            length = config.get("len", None)
            if length is None:
                raise ValueError(
                    f'Cluster config file, {config_path}, does not have "len" parameter'
                )

        chain = []
        for i in range(length):
            chain.append(
                smat_util.load_matrix(os.path.join(path_to_cluster, f"C{i}.npz"))
                .tocsc()
                .astype(np.float32)
            )

        return cls(chain)

    @classmethod
    def from_partial_chain(cls, C, min_codes=None, nr_splits=16):
        """Creates the clustering matrices necessary for a valid hierarchical clustering chain.

        Except for the bottom level of the hierarchy chain provided, all other levels are filled with dummy CSC matrices (all-one).

        Args:
            C (scipy.sparse.base.spmatrix or list): Partial chain for the bottom level of the hierarchy.
                If sparse matrix, this arg should be the matrix representing the bottom level of the hierarchy.
                If list, this arg should be a list of sparse matrices representing the bottom levels of the hierarchy.
            min_codes (int): The number of direct child nodes that the top level of the hierarchy should have.
            nr_splits (int): The out-degree of the clustering matrices that this method will create.

        Returns:
            ClusterChain: The completed hierarchical clustering chain.
        """
        if min_codes is None:
            min_codes = nr_splits

        if isinstance(C, smat.spmatrix):
            cluster_chain = [C.tocsc()]
        else:
            assert isinstance(C, (cls, list, tuple))
            cluster_chain = C
        cur_C = cluster_chain[0]
        if min_codes is None or min_codes <= 1:
            min_codes = cur_C.shape[1]

        # where cur_C.shape == (nr_labels, nr_codes)
        while cur_C.shape[1] > min_codes:
            nr_codes = cur_C.shape[1]
            cur_codes = np.arange(nr_codes)
            new_codes = cur_codes // nr_splits
            shape = (len(cur_codes), new_codes.max() + 1)
            new_C = smat.csc_matrix(
                (np.ones_like(cur_codes), (cur_codes, new_codes)), shape=shape, dtype=np.float32
            )
            cluster_chain = [new_C] + cluster_chain
            cur_C = new_C

        if cur_C.shape[1] > 1:
            new_C = smat.csc_matrix(np.ones((cur_C.shape[1], 1), dtype=np.float32))
            cluster_chain = [new_C] + cluster_chain
        return cls(cluster_chain)

    def matrix_chain_dimension_check(self, M_dict):
        """Check dimension of matrix chain provided by dictionary with keys being number of layers above leaf elements.

        Args:
            M_dict (dict): dictionary of partial matrix chains to check.

        Returns:
            (nr_insts, nr_labels)
        """
        # get/check the dimensions
        assert isinstance(M_dict, dict)
        nr_labels = self.chain[-1].shape[0]
        assert set(M_dict.keys()) <= set(range(len(self) + 1)), "M_dict got invalid key"
        nr_insts = [v.shape[0] for k, v in M_dict.items() if v is not None]
        assert nr_insts.count(nr_insts[0]) == len(nr_insts), "M_dict first dim do not match"
        nr_insts = nr_insts[0]

        # cannot simply use if M_dict.get(0, None) here
        # since truth value of a [MATRIX/ARRAY] is ambiguous
        if M_dict.get(0, None) is not None:
            assert M_dict[0].shape[1] == self.chain[-1].shape[0]
        for i in range(1, len(self) + 1):
            if M_dict.get(i, None) is not None:
                assert (
                    M_dict[i].shape[1] == self.chain[-i].shape[1]
                ), f"{i}: {M_dict[i].shape}!={self.chain[-i].shape}"

        return nr_insts, nr_labels

    def generate_matching_chain(self, M_dict):
        """Generate a chain of instance to cluster matching matrix for user supplied negative (usn) from partial matching chain.

        Args:
            M_dict (dict): dictionary of partial matching chains, with keys being number of layers above leaf elements.
                M_dict[i].shape[0] == nr_inst, for all i.
                M_dict[0].shape[1] == self.chain[-1].shape[0],
                M_dict[i].shape[1] == self.chain[-i].shape[1], for i >= 1
                M_dict.keys() \\subset range(len(self.chain)+1)

        Returns:
            matching_chain: list of csc matrices for user supplied negatives
        """
        matching_chain = [None] * (len(self) + 1)
        # if nothing is given, return a chain of None
        if M_dict is None or all(M_dict[x] is None for x in M_dict):
            return matching_chain

        nr_insts, nr_labels = self.matrix_chain_dimension_check(M_dict)

        # construct matching chain from incomplete chain
        if M_dict.get(0, None) is not None:
            matching_chain[0] = smat_util.binarized(M_dict[0])
        else:
            matching_chain[0] = smat.csc_matrix((nr_insts, nr_labels), dtype=np.float32)
        for i in range(1, len(self) + 1):
            matching_chain[i] = clib.sparse_matmul(matching_chain[i - 1], self.chain[-i])
            if M_dict.get(i, None) is not None:
                matching_chain[i] += smat_util.binarized(M_dict[i])
            matching_chain[i] = matching_chain[i].tocsc().sorted_indices()
        matching_chain.reverse()

        return matching_chain[:-1]

    def generate_relevance_chain(self, R_dict, norm_type=None, induce=True):
        """Generate a chain of instance to cluster relevance matrix for cost sensitive learning from partial relevance chain.

        Args:
            R_dict (dict): dictionary of partial relevance chains, with keys being number of layers above leaf elements.
                R_dict[i].shape[0] == nr_inst, for all i.
                R_dict[0].shape[1] == self.chain[-1].shape[0],
                R_dict[i].shape[1] == self.chain[-i].shape[1], for i >= 1
                R_dict.keys() \\subset range(len(self.chain)+1)
            norm_type (str, optional): row wise normalziation of resulting relevance matrices. Defatult None to ignore.
                Options: ‘l1’, ‘l2’, ‘max’, 'no-norm', None
            induce (bool, optional): whether to induce missing relevance matrix by label aggregation. Default True

        Returns:
            relevance_chain: list of csc matrices for relevance
        """

        relevance_chain = [None] * (len(self) + 1)
        # if nothing is given, return a chain of None
        if R_dict is None or all(R_dict[x] is None for x in R_dict):
            return relevance_chain

        self.matrix_chain_dimension_check(R_dict)

        # construct relevance chain from incomplete chain
        relevance_chain[0] = R_dict.get(0, None)
        for i in range(1, len(self) + 1):
            if R_dict.get(i, None) is not None:
                relevance_chain[i] = R_dict[i]
            elif relevance_chain[i - 1] is not None and induce:
                relevance_chain[i] = clib.sparse_matmul(relevance_chain[i - 1], self.chain[-i])
            else:
                relevance_chain[i] = None
        relevance_chain.reverse()

        if norm_type not in [None, "no-norm"]:
            relevance_chain = [
                sk_normalize(rr.tocsr(), norm=norm_type) if rr is not None else None
                for rr in relevance_chain
            ]

        return relevance_chain[1:]

def hierarchical_kmeans(
    feat_mat,
    max_leaf_size=100,
    imbalanced_ratio=0,
    imbalanced_depth=100,
    spherical=True,
    seed=0,
    kmeans_max_iter=20,
    threads=-1,
):
    """Python implementation of hierarchical 2-means.

    Args:
        feat_mat (numpy.ndarray or scipy.sparse.csr.csr_matrix): Matrix of label features.
        max_leaf_size (int, optional): The maximum size of each leaf node of the tree. Default is `100`.
        imbalanced_ratio (float, optional): Indicates how relaxed the balancedness constraint of 2-means can be.
            Specifically, if an iteration of 2-means is clustering `L` labels, the size of the output 2 clusters will be within approx `imbalanced_ratio * 2 * L` of each other.
            Value between `0.0` and `0.5` (inclusive). Default is `0.0`.
        imbalanced_depth (int, optional): Maximum depth of imbalanced clustering. After depth `imbalanced_depth` is reached, balanced clustering will be used. Default is `100`.
        spherical (bool, optional): True will l2-normalize the centroids of k-means after each iteration. Default is `True`.
        seed (int, optional): Random seed. Default is `0`.
        kmeans_max_iter (int, optional): Maximum number of iterations for each k-means problem. Default is `20`.
        threads (int, optional): Number of threads to use. `-1` denotes all CPUs. Default is `-1`.

    Returns:
        list: List of CSC matrices representing the generated hierarchical cluster chain.
    """

    global run_kmeans

    def run_kmeans(cluster, c1, c2, min_size, kmeans_max_iter, spherical=True):
        indexer = kmeans(feat_mat_global[cluster], c1, c2, min_size, kmeans_max_iter, spherical)
        return cluster[indexer], cluster[~indexer]

    global kmeans

    def kmeans(feat_mat, c1=-1, c2=-1, min_size=50, kmeans_max_iter=20, spherical=True):
        if c1 == -1:
            c1, c2 = np.random.randint(feat_mat.shape[0]), np.random.randint(1, feat_mat.shape[0])
        c1, c2 = feat_mat[c1], feat_mat[(c1 + c2) % feat_mat.shape[0]]
        old_indexer = np.ones(feat_mat.shape[0]) * -1

        for _ in range(kmeans_max_iter):
            scores = np.squeeze(np.asarray(feat_mat.multiply(c1 - c2).sum(1)))
            indexer = scores >= 0
            if indexer.sum() < min_size:
                indexer = np.zeros(feat_mat.shape[0], dtype=np.bool)
                indexer[np.argpartition(-scores, min_size)[:min_size]] = True
            elif (~indexer).sum() < min_size:
                indexer = np.zeros(feat_mat.shape[0], dtype=np.bool)
                indexer[np.argpartition(scores, min_size)[min_size:]] = True
            if np.array_equal(indexer, old_indexer):
                break
            old_indexer = indexer
            c1 = feat_mat[indexer].sum(0)
            c2 = feat_mat[~indexer].sum(0)
            if spherical:
                c1 = sk_normalize(c1)
                c2 = sk_normalize(c2)
        return indexer

    global feat_mat_global
    feat_mat_global = feat_mat
    random = np.random.RandomState(seed)
    cluster_chain = []
    clusters_big, clusters_small = [], []
    if feat_mat.shape[0] > max_leaf_size:
        clusters_big.append(np.arange(feat_mat.shape[0]))
    else:
        clusters_small.append(np.arange(feat_mat.shape[0]))
    threads = threads if threads > 0 else mp.cpu_count()

    while len(clusters_big) > 0:
        # Do balanced clustering beyond imbalanced_depth to ensure reasonably timely termination
        if len(cluster_chain) >= imbalanced_depth:
            imbalanced_ratio = 0

        num_parent_clusters = len(clusters_big) + len(clusters_small)
        new_clusters_big = []
        cols_big, cols_small = [], [x + len(clusters_big) for x in range(len(clusters_small))]
        seeds = [(random.randint(s), random.randint(1, s)) for s in map(len, clusters_big)]
        min_sizes = [int(s * (0.5 - imbalanced_ratio)) for s in map(len, clusters_big)]
        zip_args = zip(
            clusters_big,
            *map(list, zip(*seeds)),
            min_sizes,
            repeat(kmeans_max_iter),
            repeat(spherical),
        )
        local_threads = min(threads, len(clusters_big))
        with mp.Pool(local_threads) as p:
            for col, child_clusters in enumerate(p.starmap(run_kmeans, zip_args)):
                for cluster in child_clusters:
                    if len(cluster) > max_leaf_size:
                        new_clusters_big.append(cluster)
                        cols_big.append(col)
                    elif len(cluster) > 0:
                        clusters_small.append(cluster)
                        cols_small.append(col)

        cols = cols_big + cols_small

        cluster_chain.append(
            smat.csc_matrix(
                (np.ones(len(cols)), (range(len(cols)), cols)),
                shape=(len(new_clusters_big + clusters_small), num_parent_clusters),
                dtype=np.float32,
            )
        )

        clusters_big = new_clusters_big

    C = []
    for col, cluster in enumerate(chain(clusters_big, clusters_small)):
        for row in cluster:
            C.append((row, col))
    cluster_chain.append(
        smat.csc_matrix(
            (np.ones(feat_mat.shape[0]), list(map(list, zip(*C)))),
            shape=(feat_mat.shape[0], len(clusters_big) + len(clusters_small)),
            dtype=np.float32,
        )
    )
    return cluster_chain

max_stars_repo_name: anukaal/pecos

max_issues_repo_name: anukaal/pecos

max_forks_repo_name: anukaal/pecos