text:
# ============================================================================ #
#                                                                              #
#    This is part of the "GrainSizeTools Script"                               #
#    A Python script for characterizing grain size from thin sections          #
#    and paleopiezometry estimates based on grain size.                        #
#                                                                              #
#    Copyright (c) 2014-present   Marco A. Lopez-Sanchez                       #
#                                                                              #
#    Licensed under the Apache License, Version 2.0 (the "License");           #
#    you may not use this file except in compliance with the License.          #
#    You may obtain a copy of the License at                                   #
#                                                                              #
#        http://www.apache.org/licenses/LICENSE-2.0                            #
#                                                                              #
#    Unless required by applicable law or agreed to in writing, software       #
#    distributed under the License is distributed on an "AS IS" BASIS,         #
#    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  #
#    See the License for the specific language governing permissions and       #
#    limitations under the "License".                                          #
#                                                                              #
#    Version 3.0.2                                                             #
#    For details see: http://marcoalopez.github.io/GrainSizeTools/             #
#    download at https://github.com/marcoalopez/GrainSizeTools/releases        #
#                                                                              #
# ============================================================================ #

# Imports
from scipy.stats import bayes_mvs, gaussian_kde, iqr, t, norm
import numpy as np

# ============================================================================ #
# AVERAGES                                                                     #
# ============================================================================ #

def amean(pop, ci=0.95, method='ASTM'):
    """ Returns the arithmetic mean, the Bessel corrected SD,
    and the confidence interval based on the chosen method.

    Parameters
    ----------
    pop : array-like
        the population

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    method : string
        the method to estimate the confidence interval, either
        'ASTM': central limit theorem based (ASTM default)
        'GCI': generalized confidence interval method
        'mCox': modified Cox method

    Assumptions
    -----------
    - arithmetic mean is optimal for normal-like distributions
    - CLT confidence interval is optimized for normal distributions
    - GCI and mCox methods are optimized for lognormal distributions

    Call functions
    --------------
    - CLT_ci
    - GCI_ci
    - mCox_ci

    Returns
    -------
    the arithmetic mean,
    the SD (Bessel corrected),
    the confidence interval (scalar or tuple)
    the confidence interval length (float),
    """

    n = len(pop)
    mean, std = np.mean(pop), np.std(pop, ddof=1)  # SD using n-1 degrees of freedom (Bessel corrected)

    # confidence interval
    if method == 'ASTM':
        conf_int, length = CLT_ci(mean, std, n, ci)
        return mean, std, conf_int, length

    elif method == 'GCI':
        ci_limis, length = GCI_ci(pop, ci)
        return mean, std, ci_limis, length

    elif method == 'mCox':
        ci_limis, length = mCox_ci(pop, ci)
        return mean, std, ci_limis, length

    else:
        raise Exception("ci methods must be 'CLT', 'GCI', or 'mCox'")

def gmean(pop, ci=0.95, method='CLT'):
    """ Returns the geometric mean, the multiplicative (geometric) SD,
    and the confidence interval.

    Parameters
    ----------
    pop : array-like
        the population

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    method : string
       the method to estimate the confidence interval, either
        'CLT': Central limit theorem based
        'bayes': Bayesian based

    Assumptions
    -----------
    - geometric mean is optimal for lognormal-like distributions
    - the multiplicative SD is a measure of the lognormal shape
    - The bayes method is sometimes slighly superior to CLT for
    very small (< 100) sample sizes

    Call functions
    --------------
    - CLT2_ci
    - bayesian_ci

    Returns
    -------
    the geometric mean,
    the multiplicative SD (MSD),
    the confidence interval (tuple),
    the confidence interval length (float),
    """

    # compute statistics of the log-transformed data
    mean_log, n = np.mean(np.log(pop)), len(pop)
    std_log = np.std(np.log(pop), ddof=1)  # Bessel corrected SD (n-1 degrees of freedom)

    # compute the back-transformed values (gmean and mSD in linear scale)
    gmean = np.exp(mean_log)
    mSD = np.exp(std_log)

    # confidence intervals of the back-transformed values
    if method == 'CLT':
        ci_limis, length = CLT2_ci(mean_log, std_log, n, ci)
        return gmean, mSD, ci_limis, length

    elif method == 'bayes':
        ci_limis, length = bayesian_ci(pop, ci)
        return gmean, mSD, ci_limis, length

    else:
        raise Exception("CI methods must be 'CLT' or 'bayes'")

def median(pop, ci=0.95):
    """ Returns the median, the interquartile length, and the confidence
    intervals for the median based on th rule-of-thumb method of Hollander
    and Wolfe (1999).

    Parameters
    ----------
    pop : array-like
        the population

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Assumptions
    -----------
    - median is optimal for both normal and lognormal-like distributions.
    It behaves better than means when data contamination is expected
    (i.e. more robust).
    - the interquertile length/range is a measure of the spread of
    the distribution

    Reference
    ---------
    Hollander and Wolfe (1999) Nonparametric Statistical Methods.
    3rd ed. John Wiley, New York. 787 pp.

    Call functions
    --------------
    - iqr from Scipy

    Returns
    -------
    the median (float),
    the interquartile range,
    the confidence interval (tuple),
    the confidence length (float)
    """
    pop, n = np.sort(pop), len(pop)
    median, iqr_range = np.median(pop), iqr(pop)

    # compute confidence intervals
    ci_limits, length = median_ci(pop, n, ci)

    return median, iqr_range, ci_limits, length

def freq_peak(pop, bandwidth='silverman', max_precision=0.05):
    """ Returns the peak of the frequency ("mode") of a continuous
    distribution based on the Gaussian kernel density estimator. It
    uses Scipy's gaussian kde method.

    Parameters
    ----------
    pop : array_like
        the diameters of the grains

    bandwidth : string, positive scalar or callable
        the method to estimate the bandwidth or a scalar directly defining the
        bandwidth. Methods can be 'silverman' or 'scott'.

    max_precision : positive scalar, default is 0.05
        the maximum precision expected for the "peak" estimator.

    Call functions
    --------------
    - gen_xgrid
    - gaussian_kde from scipy

    Returns
    -------
    the x and y values to contruct the kde (tuple),
    the mode or peak grain size,
    the density value of the peak,
    the bandwidth
    """

    # check bandwidth and estimate Gaussian kernel density function
    if isinstance(bandwidth, (int, float)):
        bw = bandwidth / np.std(pop, ddof=1)
        kde = gaussian_kde(pop, bw_method=bw)

    elif isinstance(bandwidth, str):
        kde = gaussian_kde(pop, bw_method=bandwidth)
        bw = round(kde.covariance_factor() * pop.std(ddof=1), 2)

    else:
        raise ValueError("bandwidth must be integer, float, or plug-in methods 'silverman' or 'scott'")

    # locate and get the frequency peak
    xgrid = gen_xgrid(pop.min(), pop.max(), max_precision)
    densities = kde(xgrid)
    y_max, peak_grain_size = np.max(densities), xgrid[np.argmax(densities)]

    return (xgrid, densities), peak_grain_size, y_max, bw

# ============================================================================ #
# CONFIDENCE INTERVAL METHODS                                                  #
# ============================================================================ #

def CLT_ci(amean, std, n, ci=0.95):
    """ Estimate the error margin for the arithmetic mean based
    on the central limit theorem and the t-statistics. This is
    the method described in the ASTM norm E112-12.

    Parameters
    ----------
    amean : scalar, float
        the arithmetic mean of the population

    std : scalar, float
        the standard deviation of the population

    n : scalar, positive int
        the sample size

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Reference
    ---------
    ASTM-E112-12 (1996) Standard test methods for determining
    average grain size.

    Call
    ----
    calc_t

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """
    t_score = critical_t(confidence=ci, sample_size=n)
    err = t_score * std / np.sqrt(n)

    lower, upper = amean - err, amean + err
    interval = upper - lower

    return (lower, upper), interval

def CLT2_ci(mean_log, std_log, n, ci=0.95):
    """ Returns the error margin for the geometric mean based
    on the central limit theorem and the t-statistics.

    Parameters
    ----------
    mean_log : scalar, float
        the arithmetic mean of the log-transformed data

    std_log : scalar, float
        the standard deviation of the log-transformed data

    n : scalar, positive int
        the sample size

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Reference
    ---------
    ASTM-E112-12 (1996) Standard test methods for determining
    average grain size.

    Call
    ----
    calc_t

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """
    t_score = critical_t(ci, n)

    upper_ci = np.exp(mean_log + t_score * (std_log / np.sqrt(n)))
    lower_ci = np.exp(mean_log - t_score * (std_log / np.sqrt(n)))
    interval = upper_ci - lower_ci

    return (lower_ci, upper_ci), interval

def mCox_ci(data, ci=0.95):
    """ Returns the error margin for the arithmetic mean using the modified
    Cox method. This is a method optimized from lognormal populations. The
    method implemented below uses the Bessel corrected SD as it produces
    safer/robust results for small sample sizes

    Parameters
    ----------
    data : array_like
        the dataset

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Reference
    ---------
    Armstrong (1992) https://doi.org/10.1080/15298669291360003
    Lopez-Sanchez (2020) TODO

    Call
    ----
    calc_t

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """

    n = len(data)
    t = critical_t(confidence=ci, sample_size=n)
    data = np.log(data)
    mean_log, std_log = np.mean(data), np.std(data, ddof=1)

    lower = np.exp(mean_log + 0.5
                   * std_log**2 - t
                   * (std_log / np.sqrt(n))
                   * np.sqrt(1 + (std_log**2 * n) / (2 * (n + 1))))
    upper = np.exp(mean_log + 0.5
                   * std_log**2 + t
                   * (std_log / np.sqrt(n))
                   * np.sqrt(1 + (std_log**2 * n) / (2 * (n + 1))))
    interval = upper - lower

    return (lower, upper), interval

def GCI_ci(data, ci=0.95, runs=10000):
    """ Returns the confidence interval for the arithmetic mean using the
    generalized confidence interval (GCI) method of Krishnamoorthy and Mathew
    (2003). This is a Monte Carlo method optimized for lognormal populations.

    Parameters
    ----------
    data : array_like
        the dataset

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    runs : integer, default=10000
        the number of (Monte Carlo) iterations to generate z and u**2 values

    Reference
    ---------
    Krishnamoorthy and Mathew (2003) https://doi.org/10.1016/S0378-3758(02)00153-2

    Assumptions
    -----------
    - The population follows a lognormal distribution

    Call
    ----
    GCI_equation

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """

    # estimate the log-transformed population y = ln(x) and the degrees of freedom
    data = np.log(data)
    mu_log, var_log, n = np.mean(data), np.var(data), len(data)
    ddof = n - 1
    alpha = 1 - ci

    # Generate random values from the normal N(0,1) distribution
    z_array = np.random.normal(loc=0, scale=1.0, size=runs)

    # Generate random values from (non-central) chi-square distribution
    # with n-1 degrees of freedom
    u2_array = np.random.noncentral_chisquare(df=ddof, nonc=0, size=runs)
    u_array = np.sqrt(u2_array)

    # Compute the T values and sort them
    T_array = GCI_equation(mu_log, var_log, z_array, u_array, n)
    T_array = np.sort(T_array)

    # Estimate confidence limits
    lower = np.percentile(T_array, 100 * (alpha / 2))
    upper = np.percentile(T_array, 100 * (1 - (alpha / 2)))
    interval = upper - lower

    return (lower, upper), interval

def GCI_equation(mu_log, var_log, z, u, n):
    """ Generalized confidence interval (GCI) equation.

    Parameters
    ----------
    mu_log : integer, float
        the mean of the log-transformed population
    var_log : integer, float
        the variance of the log-transformed population
    z : array-like
        random values of the normal N(0,1) distribution
    u : array-like
        random values of the chi-square distribution with n-1 degrees
        of freedom
    n : integer, float
        size of the dataset

    Returns
    -------
    scalar or array-like
    """

    # estimate the second and third terms of the equation
    second_term = (z / (u / np.sqrt(n - 1))) * (np.sqrt(var_log) / np.sqrt(n))
    third_term = 0.5 * var_log / (u**2 / (n - 1))

    return np.exp(mu_log - second_term + third_term)

def bayesian_ci(data, ci=0.95):
    """ Use a bayesian approach to estimate the confidence intervals
    of the geometric mean. It uses the scipy bayes_msv routine over
    the log-transformed data and then estimate the back-transformed
    values.

    Parameters
    ----------
    data : array_like
        the dataset

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Reference
    ---------
    Oliphant (2006) https://scholarsarchive.byu.edu/facpub/278

    Assumptions
    -----------
    - The population follows a lognormal distribution

    Call
    ----
    bayes_mvs from scipy.stats module

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """

    data = np.log(data)
    mu_log, var_log, std_log = bayes_mvs(data, alpha=ci)
    mu, (lower_log, upper_log) = mu_log
    lower, upper = np.exp(lower_log), np.exp(upper_log)
    interval = upper - lower

    return (lower, upper), interval

def median_ci(pop, n, ci=0.95):
    """ Estimate the approximate ci 95% error margins for the median
    using a rule of thumb based on Hollander and Wolfe (1999).

    Parameters
    ----------
    pop : numpy array
        a sorted dataset

    n : scalar, positive int
        the sample size

    ci : float, scalar between 0 and 1
        the confidence interval, default = 0.95

    Reference
    ---------
    Hollander and Wolfe (1999) Nonparametric Statistical Methods.
    3rd ed. John Wiley, New York. 787 pp.

    Call
    ----
    norm.ppf from scipy

    Returns
    -------
    the lower and upper confidence intervals (tuple)
    the interval length (scalar)
    """

    z_score = norm.ppf(1 - (1 - ci) / 2)  # two-tailed z score

    id_upper = 1 + (n / 2) + (z_score * np.sqrt(n)) / 2
    id_lower = (n / 2) - (z_score * np.sqrt(n)) / 2

    if int(np.ceil(id_upper)) >= len(pop):
        upper_ci, lower_ci = pop[-1], pop[int(np.floor(id_lower))]
    else:
        upper_ci, lower_ci = pop[int(np.ceil(id_upper))], pop[int(np.floor(id_lower))]
    interval = upper_ci - lower_ci

    return (lower_ci, upper_ci), interval

# ============================================================================ #
# AUXILIARY FUNCTIONS                                                          #
# ============================================================================ #

def critical_t(confidence, sample_size):
    """Returns the (two-tailed) critical value of t-distribution

    Parameters
    ----------
    confidence : float, scalar between 0 and 1
        the level of confidence. E.g. 0.95 -> 95%

    sample_size : scalar, int
        the sample size

    Assumptions
    -----------
    - the population is symmetric
    """

    # recalculate confidence for the two-tailed t-distribution
    confidence = confidence + ((1 - confidence) / 2)

    return t.ppf(confidence, sample_size)

def gen_xgrid(start, stop, precision):
    """ Returns a mesh of values (i.e. discretize the
    sample space) with a fixed range and desired precision.

    Parameters
    ----------
    start : scalar
        the starting value of the sequence
    stop : scalar
        the end value of the sequence
    precision : scalar, int or float
        the desired precision (density) of the mesh
    """

    rango = stop - start

    # num = range / precision; as long as range > precision
    if rango < precision:
        raise ValueError('The precision must be smaller than the range of grain sizes')
    else:
        n = int(round(rango / precision, 0))

    return np.linspace(start, stop, num=n)

if __name__ == '__main__':
    pass
else:
    print('module averages imported')

max_stars_repo_name: jwestraadt/GrainSizeTools

max_issues_repo_name: jwestraadt/GrainSizeTools

max_forks_repo_name: jwestraadt/GrainSizeTools

reason: import numpy,from scipy


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

"""Utility functions to run the algorithm."""

import logging
from os import listdir
from os.path import isfile, join
import time
import traceback
from typing import Dict, List, Optional, Set, Union
from urllib import parse, request

import bio2bel_phewascatalog
from igraph import Graph
import mygene
import networkx as nx
import pandas as pd
from opentargets import OpenTargetsClient
from pybel import BELGraph, from_path, union
from pybel.dsl import BaseEntity, CentralDogma
from pybel_tools.assembler.reified_graph.assembler import reify_bel_graph

from guiltytargets.gat2vec import gat2vec_paths
from guiltytargets.ppi_network_annotation.model.gene import Gene
from guiltytargets.ppi_network_annotation.parsers import (
    parse_disease_ids, parse_disease_associations
)
from .reified_network import AttributeFileGenerator, LabeledNetworkGenerator, NetworkNx

logger = logging.getLogger(__name__)

def get_base_dir(basedir, path):
    return join(basedir, path)

def get_local_bel_file(basedir, path):
    dir_ = get_base_dir(basedir, path)
    return join(dir_, path + '.bel')

def get_struct_file(basedir, path):
    dir_ = get_base_dir(basedir, path)
    return join(dir_, path + '_graph.adjlist')

def get_attr_file(basedir, path):
    dir_ = get_base_dir(basedir, path)
    return join(dir_, path + '_na.adjlist')

def get_labels_file(basedir, path):
    dir_ = get_base_dir(basedir, path)
    return join(dir_, 'labels_maped.txt')

def generate_bel_network(
        bel_graph_path: str,
        dge_list: List[Gene],
        max_adj_p: float,
        max_log2_fold_change: float,
        min_log2_fold_change: float,
        current_disease_ids_path: Optional[str] = None,
        disease_associations_path: Optional[str] = None,
        collapse: str = None
):
    """Generate the protein-protein interaction network.

    :param bel_graph_path:
    :param dge_list:
    :param max_adj_p:
    :param max_log2_fold_change:
    :param min_log2_fold_change:
    :param current_disease_ids_path:
    :param disease_associations_path:
    :param collapse:
    :return: Interaction network with information on differential expression.
    """
    # Load and reify the graph
    bel_graph = bel_graph_loader(bel_graph_path)
    reified = reify_bel_graph(bel_graph, collapse=collapse)

    if disease_associations_path is not None:
        if current_disease_ids_path:
            current_disease_ids = parse_disease_ids(current_disease_ids_path)
        else:
            current_disease_ids = set()
        disease_associations = parse_disease_associations(disease_associations_path,
                                                          current_disease_ids)
    else:
        disease_associations = None

    # Build an undirected weighted graph with the remaining interactions based on Entrez gene IDs
    network = NetworkNx(
        reified,
        max_adj_p=max_adj_p,
        max_l2fc=max_log2_fold_change,
        min_l2fc=min_log2_fold_change,
    )
    network.set_up_network(dge_list, disease_associations=disease_associations)

    return network

def write_gat2vec_input_files(
        network,
        targets: List[str],
        home_dir: str,
        assoc_score: Optional[Dict] = None
):
    """Write the input files for gat2vec tool.

    :param network: Network object with attributes overlayed on it.
    :param targets:
    :param home_dir:
    :param assoc_score:
    """
    network.write_adj_list(gat2vec_paths.get_adjlist_path(home_dir, "graph"))

    attribute_network = AttributeFileGenerator(network)
    attribute_network.write_attribute_adj_list(gat2vec_paths.get_adjlist_path(home_dir, "na"))

    labeled_network = LabeledNetworkGenerator(network)
    labeled_network.write_index_labels(
        targets,
        gat2vec_paths.get_labels_path(home_dir),
        sample_scores=assoc_score
    )

def add_disease_attribute(graph: nx.Graph, att_mappings: Dict):
    """Add the phenotypes to the Base Entities as attributes."""
    for node in graph:
        if isinstance(node, CentralDogma) and node.name in att_mappings:
            graph.nodes[node]['phenotypes'] = [phtype for _, phtype in att_mappings[node.name]]

def add_dge_attribute(graph: nx.Graph, gene_dict: Dict):
    """ """
    max_adjp = 0.05
    min_l2fc, max_l2fc = -0.5, 0.5
    dge = {
        g.entrez_id or g.symbol: g.log2_fold_change
        for g
        in gene_dict
        if g.padj < max_adjp
    }

    for node in graph:
        if isinstance(node, CentralDogma) and node['name'] in dge.keys():
            if dge[node['name']] > max_l2fc:
                graph.nodes[node]['deg'] = True
                graph.nodes[node]['upreg'] = True
            elif dge[node['name']] < min_l2fc:
                graph.nodes[node]['deg'] = True
                graph.nodes[node]['downreg'] = True

def get_significantly_differentiated(gene_list: List[Gene], max_adjp: float):
    """Returns a dictionary only with significantly differentially expressed genes from the gene list."""
    max_adjp = max_adjp or 0.05

    dge = {g.entrez_id or g.symbol: g.log2_fold_change
           for g in gene_list
           if g.padj < max_adjp}

    return {k: v for k, v in dge.items() if k}

def write_adj_file_attribute(graph, filepath: str, att_mappings: Dict, pred_mapping: Optional[Dict] = None):
    """Write an adjacency file from the attribute graph."""
    if pred_mapping is None:
        pred_mapping = dict()
    with open(filepath, 'w') as file:
        for node in graph.nodes:
            line = f"{node}"
            if 'phenotypes' in graph.nodes[node]:  # There are diseases in the node (gene)
                line += f" {' '.join(str(att_mappings[phe]) for phe in graph.nodes[node]['phenotypes'])}"
            if 'label' in graph.nodes[node] and graph.nodes[node]['label'] in pred_mapping:  # predicate node
                line += f" {pred_mapping[graph.nodes[node]['label']]}"
            print(line, file=file)

# TODO Renamed from download_for_disease. Moved to assembler. Delete.
def generate_targets_file(disease_id, outpath, anno_type: str = 'entrezgene') -> None:
    """Creates a disease list

    :param disease_id: EFO code from the disease.
    :param outpath:
    :param anno_type: `entrezgene` for Entrez Id or `symbol` for Gene symbol.
    :return:
    """
    ot = OpenTargetsClient()
    assoc = ot.get_associations_for_disease(
        disease_id,
        fields=['association_scoredatatypes', 'target.id']
    ).filter(
        datatype='known_drug'
    )
    ensembl_list = [a['target']['id'] for a in assoc]

    mg = mygene.MyGeneInfo()
    id_mappings = mg.getgenes(ensembl_list, fields=anno_type)

    with open(outpath, 'w+') as outfile:
        for mapping in id_mappings:
            if anno_type in mapping.keys():
                outfile.write(mapping[anno_type])
                outfile.write('\n')

# TODO Was renamed from get_association scores. Moved to input. Delete.
def generate_disease_gene_association_file(disease_id, outpath, anno_type: str = 'entrezgene'):
    """Obtain the association scores from the specified disease that are
    stored in the OpenTargets database.

    :param disease_id: The EFO code to the disease.
    :param outpath: The path to the file to be created.
    :param anno_type: `entrezgene` for Entrez Id or `symbol` for Gene symbol.
    :return:
    """
    ot = OpenTargetsClient()
    assoc = ot.get_associations_for_disease(
        disease_id,
        fields=['association_scoreoverall', 'target.id']
    )
    assoc_simple = [
        {
            'id': a['target']['id'],
            'score': a['association_score']['overall']
        }
        for a in assoc
    ]
    ensembl_list = [a['id'] for a in assoc_simple]

    # Obtain the symbols for the genes associated to disease_id
    id_mappings = get_converter_to_entrez(ensembl_list)

    # Get the symbols and the scores
    ensembl_list = [
        (id_mappings[a['id']], a['score'])
        for a in assoc_simple
        if a['id'] in id_mappings
    ]

    with open(outpath, 'w+') as outfile:
        for symbol, score in ensembl_list:
            print(f'{symbol}\t{score}', file=outfile)

# TODO Moved to parsers. Delete.
def get_converter_to_entrez(
        query_list: List[str]
) -> Dict[str, str]:
    """Get a converter from a gene field to Entrez id.

    :param query_list: List of genes to query.
    :return: a dictionary with the query element as keys and entrez id as values.
    """
    mg = mygene.MyGeneInfo()
    gene_query = mg.querymany(query_list, scopes='symbol,entrezgene,ensembl.gene', species=9606, returnall=True)
    gene_query = [x for x in gene_query['out'] if 'entrezgene' in x]
    id_mappings = {
        g['query']: g['entrezgene']
        for g in gene_query
        if 'entrezgene' in g
    }

    logger.debug(f'get_converter_to_entrez - {len(query_list) - len(id_mappings)} not found.')
    logger.debug(f'get_converter_to_entrez - {len(id_mappings)} mappings')

    return id_mappings

# TODO moved this to parsers (GuiltyTargets). Delete.
def parse_gene_list_reified(path: str, graph: nx.Graph) -> list:
    """Parse a list of genes and return them if they are in the network.

    :param str path: The path of input file.
    :param Graph network: The graph with genes as nodes.
    :return list: A list of genes, all of which are in the network.
    """
    # read the file
    df = pd.read_csv(
        path,
        names=['gene'],
        dtype={'gene': str}
    )
    genes = set(df['gene'])

    # get those genes which are in the network
    symbols = [
        node['name']
        for node
        in graph.nodes
        if isinstance(node, CentralDogma) and 'name' in node
    ]

    return list(genes.intersection(symbols))

# TODO moved this to parsers (GuiltyTargets). Delete.
def parse_gene_list(path: str, network) -> list:
    """Parse a list of genes and return them if they are in the network.

    :param str path: The path of input file.
    :param Graph network: The graph with genes as nodes.
    :return list: A list of genes, all of which are in the network.
    """
    # read the file
    df = pd.read_csv(
        path,
        names=['gene'],
        dtype={'gene': str}
    )
    genes = set(df['gene'])

    # get those genes which are in the network
    return network.find_genes(genes)

def write_labels_file(outfile: str, graph: nx.Graph, target_list: List[str]):
    """Write the file with the nodes' ids and the labels (known target vs not).

    :param outfile: Path to the output file.
    :param graph: Graph (with the nodes already converted to integers).
    :param target_list: List with the known targets symbols.
    :return:
    """
    with open(outfile, 'w') as file:
        for node in graph.nodes:
            if (
                isinstance(graph.nodes[node]['bel_node'], BaseEntity) and
                'name' in graph.nodes[node]['bel_node'] and
                graph.nodes[node]['bel_node']['name'] in target_list
            ):
                print(f'{node}\t1', file=file)
            else:
                print(f'{node}\t0', file=file)

# TODO Moved to parsers. Delete.
def generate_phewas_file(out_path: str, anno_type: str = 'symbol'):
    """Create a file with gene-disease association from PheWAS data.

    :param out_path: Path where the file will be created.
    :param anno_type: `entrezgene` for Entrez Id or `symbol` for Gene symbol.
    :return:
    """
    phewas_manager = bio2bel_phewascatalog.Manager()
    pw_dict = phewas_manager.to_dict()
    if anno_type == 'entrezgene':
        to_entrez = get_converter_to_entrez(list(pw_dict.keys()))
        pw_dict = {
            entrez: pw_dict[symbol]
            for symbol, entrez
            in to_entrez.items()
        }
    with open(out_path, 'w') as file:
        for target, diseaselist in pw_dict.items():
            for disease in diseaselist:
                print(f'{target} {disease[1]}', file=file)

# TODO delete this only after notebook BEL2Gat2Vec(DGE) becomes obsolete
def get_phenotype_list(rbg: nx.Graph) -> List:
    """ Get the list of phenotypes per gene in the graph.

    :param start_ind:
    :return:
    """
    pw_dict = bio2bel_phewascatalog.Manager().to_dict()
    start_ind = len(rbg.nodes) + 3
    phenotypes_list = set([phe for _list in pw_dict.values()
                           for odds, phe in _list])
    enum_phenotypes = enumerate(phenotypes_list, start=start_ind)
    phenot_hash = {phenot: num for num, phenot in enum_phenotypes}

    annotated_graph = rbg.copy()
    add_disease_attribute(annotated_graph, pw_dict)

# TODO Moved to Parsers. Delete.
def bel_graph_loader(from_dir: str) -> BELGraph:
    """Obtains a combined BELGraph from all the BEL documents in one folder.

    :param from_dir: The folder with the BEL documents.
    :return: A corresponding BEL Graph.
    """
    logger.info("Loading BEL Graph.")
    files = [
        join(from_dir, file)
        for file
        in listdir(from_dir)
        if isfile(join(from_dir, file))
    ]
    bel_files = [file for file in files if file[-4:].lower() == '.bel']
    bel_graphs = [from_path(file) for file in bel_files]
    return union(bel_graphs)

# TODO moved to converters. Delete.
def uniprot_id_to_entrez_id_converter(uniprot_id_list: Union[List[str], Set[str]]) -> Dict:
    """Converts a list of gene identifiers, from Uniprot to Entrez, using the Uniprot API.

    :param uniprot_id_list: The list of identifiers.
    :return: A converter dictionary from Uniprot id to Entrez id.
    """
    url = 'https://www.uniprot.org/uploadlists/'
    params = {
        'from': 'ACC+ID',
        'to': 'P_ENTREZGENEID',
        'format': 'tab',
        'query': ' '.join(uniprot_id_list)
    }
    data = parse.urlencode(params)
    data = data.encode('utf-8')
    req = request.Request(url, data)
    with request.urlopen(req) as f:
        response = f.read()
    conv = {
        line.split('\t')[0]: line.split('\t')[1]
        for line
        in response.decode('utf-8').split('\n')
        if line
    }
    return conv

# to be used for link prediction
def merge_by_attribute(g1: Graph, g2: Graph, v_attr_name) -> Graph:
    """Merge two graphs by one of the attributes (`name` or `symbol`).

    :param g1: Graph 1.
    :param g2: Graph 2.
    :param v_attr_name: The name of the attribute to be used to merge.
    :return: The merged graph.
    """
    assert len(set(g1.vs[v_attr_name])) == len(g1.vs[v_attr_name]), "Merging attribute must be unique"
    assert len(set(g2.vs[v_attr_name])) == len(g2.vs[v_attr_name]), "Merging attribute must be unique"

    v_attr_1 = g1.vs[v_attr_name]
    v_attr_2 = g2.vs[v_attr_name]

    edge_list_by_attribute_1 = set([
        tuple(sorted([v_attr_1[u], v_attr_1[v]]) + [g1.es[(u, v)]['edge_type']])
        for u, v
        in g1.get_edgelist()
    ])
    edge_list_by_attribute_2 = set([
        tuple(sorted([v_attr_2[u], v_attr_2[v]]) + [g1.es[(u, v)]['edge_type']])
        for u, v
        in g2.get_edgelist()
    ])

    # Edges should not coincide (different edge types)
    edge_list_by_attribute_merged = edge_list_by_attribute_1.union(edge_list_by_attribute_2)

    v_attr_merged = sorted(list(set(v_attr_2).union(set(v_attr_1))))

    attribute_to_ind = {v_attr_merged: i for i, v_attr_merged in enumerate(v_attr_merged)}
    edge_list_merged = {
        (attribute_to_ind[i], attribute_to_ind[j]): e_type
        for i, j, e_type
        in edge_list_by_attribute_merged
    }

    graph_merged = Graph(list(edge_list_merged.keys()))
    for edge, e_type in edge_list_merged.items():
        graph_merged.es[edge]['edge_type'] = e_type
    graph_merged.vs[v_attr_name] = v_attr_merged

    # Include attributes that are in both g1 and g2. If different attribute values are present in a vertex,
    # then the one of g1 is used
    for attr_name_other in set(g2.vs.attributes()).intersection(set(g1.vs.attributes())).difference([v_attr_name]):
        attr_other = dict()
        for v in g2.vs():
            attr_other[attribute_to_ind[v[v_attr_name]]] = v[attr_name_other]

        for v in g1.vs():
            attr_other[attribute_to_ind[v[v_attr_name]]] = v[attr_name_other]
        graph_merged.vs[attr_name_other] = [attr_other[i] for i in range(graph_merged.vcount())]

    return graph_merged

def timed_main_run(main_function, logger_obj=logger):
    """Times the run of a function.

    :param main_function: The function to be run.
    :param logger_obj: The logger that will produce the result.
    :return:
    """
    start_time = time.time()
    logger_obj.info('Starting...')
    try:
        main_function()
    except Exception as e:
        logger_obj.error(type(e))
        logger_obj.error(traceback.format_exc())
    finally:
        logger_obj.info(f"Total time: {time.time() - start_time}")

ext: py

max_stars_repo_name: GuiltyTargets/phewas

max_issues_repo_name: GuiltyTargets/phewas

max_forks_repo_name: GuiltyTargets/phewas

reason: import networkx


text:
import email
import logging
import mailbox
import os
import unittest

import networkx as nx
import pandas as pd
from testfixtures import LogCapture

import bigbang.archive as archive
import bigbang.mailman as mailman
import bigbang.parse as parse
import bigbang.process as process
import bigbang.utils as utils
from bigbang import repo_loader
from config.config import CONFIG

test_txt = ""
TEMP_DIR = os.path.join(CONFIG.test_data_path, "tmp")

class TestGit(unittest.TestCase):
    def test_git_dependancy(self):
        repo_loader.get_repo(
            "https://github.com/sbenthall/bigbang.git",
            in_type="remote",
        )

class TestMailman(unittest.TestCase):
    def setUp(self):
        try:
            os.mkdir(TEMP_DIR)
        except OSError:  # Python 2.7-specific, alas; FileExistsError in py3
            pass  # temporary directory already exists, that's cool

    def tearDown(self):
        # remove all files in the temporary files directory, as cleanup
        temp_files = os.listdir(TEMP_DIR)
        for f in temp_files:
            os.remove(os.path.join(TEMP_DIR, f))

    def test_provenance(self):
        test_list_name = "test-list-name"
        test_list_url = "https://example.com/test-list-url/"
        test_notes = "Test notes."
        mailman.populate_provenance(
            TEMP_DIR,
            list_name=test_list_name,
            list_url=test_list_url,
            notes=test_notes,
        )

        self.assertTrue(
            os.path.exists(
                os.path.join(TEMP_DIR, mailman.PROVENANCE_FILENAME)
            ),
            msg="provenance file should have been created",
        )

        provenance = mailman.access_provenance(TEMP_DIR)
        self.assertTrue(
            provenance is not None, "provenance should be something"
        )
        self.assertTrue(
            provenance["list"]["list_name"] == test_list_name,
            "list name should be in the provenance",
        )
        self.assertTrue(
            provenance["list"]["list_url"] == test_list_url,
            "list url should be in the provenance",
        )
        self.assertTrue(
            provenance["notes"] == test_notes,
            "notes should be in the provenance",
        )

        provenance["notes"] = "modified provenance"
        mailman.update_provenance(TEMP_DIR, provenance)
        provenance_next = mailman.access_provenance(TEMP_DIR)
        self.assertTrue(
            provenance_next["notes"] == "modified provenance",
            "confirm modified provenance",
        )

    def test_valid_urls(self):
        test_urls_path = os.path.join(
            CONFIG.test_data_path, "urls-test-file.txt"
        )
        with LogCapture() as _l:
            urls = mailman.urls_to_collect(test_urls_path)
            self.assertTrue(
                "#ignored" not in urls, msg="failed to ignore a comment line"
            )
            self.assertTrue(
                "http://www.example.com/1" in urls,
                msg="failed to find valid url",
            )

            self.assertTrue(
                "http://www.example.com/2/" in urls,
                msg="failed to find valid url, whitespace strip issue",
            )
            self.assertTrue(
                "https://www.example.com/3/" in urls,
                msg="failed to find valid url, whitespace strip issue",
            )
            self.assertTrue(
                "invalid.com" not in urls, msg="accepted invalid url"
            )
            self.assertTrue(
                len(list(_l.actual())) == 2, msg="wrong number of log entries"
            )
            for (fromwhere, level, msg) in _l.actual():
                self.assertTrue(
                    level == "WARNING",
                    msg="logged something that wasn't a warning",
                )
                self.assertTrue(
                    len(urls) == 3, msg="wrong number of urls parsed from file"
                )

    def test_empty_list_compute_activity_issue_246(self):
        test_df_csv_path = os.path.join(
            CONFIG.test_data_path, "empty-archive-df.csv"
        )
        df = pd.read_csv(test_df_csv_path)

        with self.assertRaises(mailman.MissingDataException):
            empty_archive = archive.Archive(df)
            empty_archive.get_activity()

    def test_mailman_normalizer(self):
        browse_url = "https://mailarchive.ietf.org/arch/browse/ietf/"
        search_url = (
            "https://mailarchive.ietf.org/arch/search/?email_list=ietf"
        )
        random_url = "http://example.com"

        better_url = "https://www.ietf.org/mail-archive/text/ietf/"

        self.assertTrue(
            mailman.normalize_archives_url(browse_url) == better_url,
            msg="failed to normalize",
        )
        self.assertTrue(
            mailman.normalize_archives_url(search_url) == better_url,
            msg="failed to normalize",
        )
        self.assertTrue(
            mailman.normalize_archives_url(random_url) == random_url,
            msg="should not have changed other url",
        )

    def test_mailman_list_name(self):
        ietf_archive_url = "https://www.ietf.org/mail-archive/text/ietf/"
        w3c_archive_url = (
            "https://lists.w3.org/Archives/Public/public-privacy/"
        )
        random_url = "http://example.com"

        self.assertTrue(
            mailman.get_list_name(ietf_archive_url) == "ietf",
            msg="failed to grab ietf list name",
        )
        self.assertTrue(
            mailman.get_list_name(w3c_archive_url) == "public-privacy",
            msg="failed to grab w3c list name",
        )
        self.assertTrue(
            mailman.get_list_name(random_url) == random_url,
            msg="should not have changed other url",
        )

    def test_activity_summary(self):
        list_url = (
            "https://lists.w3.org/Archives/Public/test-activity-summary/"
        )
        activity_frame = mailman.open_activity_summary(
            list_url, archive_dir=CONFIG.test_data_path
        )

        self.assertTrue(
            str(type(activity_frame))
            == "<class 'pandas.core.frame.DataFrame'>",
            msg="not a DataFrame?",
        )
        self.assertTrue(
            len(activity_frame.columns) == 1,
            msg="activity summary should have one column",
        )

class TestArchive(unittest.TestCase):
    def test_mailman_chain(self):
        name = "bigbang-dev-test.txt"

        # archive loaded from mbox
        arx = archive.Archive(name, archive_dir="tests/data", mbox=True)

        arx.save("test.csv")

        # archive loaded from stored csv
        arx2 = archive.load("test.csv")

        print(arx.data.dtypes)
        print(arx.data.shape)

        self.assertTrue(
            arx.data.shape == arx2.data.shape,
            msg="Original and restored archives are different shapes",
        )

        self.assertTrue(
            (arx2.data.index == arx.data.index).all(),
            msg="Original and restored archives have nonidentical indices",
        )

        self.assertTrue(
            [t.get_num_messages() for t in arx.get_threads()] == [3, 1, 2],
            msg="Thread message count in mbox archive is off",
        )
        self.assertTrue(
            [t.get_num_messages() for t in arx2.get_threads()] == [3, 1, 2],
            msg="Thread message count in restored archive is off",
        )

        # smoke test entity resolution
        arx2.resolve_entities()

        os.remove("test.csv")

    def test_clean_message(self):
        name = "2001-November.txt"

        arx = archive.Archive(name, archive_dir="tests/data", mbox=True)

        body = arx.data["Body"]["<E165uMn-0002IJ-00@spock.physics.mcgill.ca>"]

        self.assertTrue(
            "But seemingly it is even stranger than this." in body,
            msg="Selected wrong message",
        )

        self.assertTrue(
            "Is it a problem of lapack3.0 of of" in body,
            msg="Quoted text is not in uncleaned message",
        )

        self.assertTrue(
            "Is it a problem of lapack3.0 of of"
            not in utils.clean_message(body),
            msg="Quoted text is in cleaned message",
        )

    def test_email_entity_resolution(self):
        name = "2001-November.txt"

        arx = archive.Archive(name, archive_dir="tests/data", mbox=True)

        e = process.resolve_sender_entities(arx.get_activity(resolved=False))

        utils.repartition_dataframe(arx.get_activity(), e)

        self.assertTrue(True, msg="email entity resolution crashed")

class TestMailParsing(unittest.TestCase):
    def test_split_references(self):
        refs = " <ye1y9ljtxwk.fsf@orange30.ex.ac.uk>\n\t<055701c16727$b57fed90$8fd6afcf@pixi.com>"
        split = parse.split_references(refs)
        self.assertTrue(len(split) == 2, msg=split)

    def test_name_email_parsing(self):
        from_header = "John_Doe (?) <John.Doe@example.com>"
        (raw_name, raw_email) = email.utils.parseaddr(from_header)
        normalized_email = parse.normalize_email_address(raw_email)
        self.assertTrue(
            normalized_email == "john.doe@example.com",
            msg="normalized email case incorrect",
        )

        clean_name = parse.clean_name(raw_name)
        self.assertTrue(clean_name == "John Doe", msg="name not fully cleaned")

        empty_name = parse.clean_name(" ")
        self.assertTrue(
            empty_name is None, msg="empty name not cleaned to None"
        )

        tokenized_name = parse.tokenize_name(clean_name)
        self.assertTrue(
            tokenized_name == "doe john",
            msg="name not properly normalized and tokenized",
        )

        parse.tokenize_name(str("   "))
        self.assertTrue(
            empty_name is None, msg="empty name not tokenized to None"
        )

class TestMailProcessing(unittest.TestCase):
    def test_from_header_distance(self):
        a = "Fernando.Perez at colorado.edu (Fernando.Perez at colorado.edu)"
        b = "Fernando.Perez at colorado.edu (Fernando.Perez@colorado.edu)"

        self.assertTrue(
            process.from_header_distance(a, b) == 0,
            msg="from_header_distance computing incorrect value",
        )

        a = ""
        b = ""

class TestUtils(unittest.TestCase):
    def test_labeled_blockmodel(self):
        g = nx.DiGraph()

        g.add_edge(0, 1)
        g.add_edge(0, 2)
        g.add_edge(0, 3)
        g.add_edge(0, 4)

        p = {"B": [1, 2, 3, 4], "A": [0]}

        bg = utils.labeled_blockmodel(g, p)

        self.assertTrue(
            list(bg.edges(data=True))[0][2]["weight"] == 4.0,
            msg="Incorrect edge weight in labeled blockmodel",
        )

        self.assertTrue(
            list(bg.edges()) == [("A", "B")],
            msg="Incorrected edges in labeled blockmodel",
        )

max_stars_repo_name: RiccardoNanni/bigbang

max_issues_repo_name: RiccardoNanni/bigbang

max_forks_repo_name: RiccardoNanni/bigbang


text:
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
#   See COPYING file distributed along with the PyMVPA package for the
#   copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
'''Tests for filter mappers'''

from mvpa2.testing import *
skip_if_no_external('scipy')

import numpy as np

from mvpa2.datasets import Dataset, vstack
from mvpa2.mappers.filters import FFTResampleMapper, iir_filter

def test_resample():
    time = np.linspace(0, 2*np.pi, 100)
    ds = Dataset(np.vstack((np.sin(time), np.cos(time))).T,
                 sa = {'time': time,
                       'section': np.repeat(range(10), 10)})
    assert_equal(ds.shape, (100, 2))

    # downsample
    num = 10
    rm = FFTResampleMapper(num, window=('gauss', 50),
                           position_attr='time',
                           attr_strategy='sample')
    mds = rm.forward(ds)
    assert_equal(mds.shape, (num, ds.nfeatures))
    # didn't change the orig
    assert_equal(len(ds), 100)

    # check position-based resampling
    ds_partial = ds[0::10]
    mds_partial = rm.forward(ds_partial)
    # despite different input sampling should yield the same output timepoints
    assert_array_almost_equal(mds.sa.time, mds_partial.sa.time)
    # exclude the first points to prevent edge effects, but the data should be
    # very similar too
    assert_array_almost_equal(mds.samples[2:], mds_partial.samples[2:], decimal=2)
    # simple sample of sa's should give meaningful stuff
    assert_array_equal(mds.sa.section, range(10))

    # and now for a dataset with chunks
    cds = vstack([ds.copy(), ds.copy()])
    cds.sa['chunks'] = np.repeat([0,1], len(ds))
    rm = FFTResampleMapper(num, attr_strategy='sample', chunks_attr='chunks',
                           window=('gauss', 50))
    mcds = rm.forward(cds)
    assert_equal(mcds.shape, (20, 2))
    assert_array_equal(mcds.sa.section, np.tile(range(10),2))
    # each individual chunks should be identical to previous dataset
    assert_array_almost_equal(mds.samples, mcds.samples[:10])
    assert_array_almost_equal(mds.samples, mcds.samples[10:])

def test_iirfilter():
    # dataset with one feature from two waves
    t = np.linspace(0, 1.0, 2001)
    xlow = np.sin(2 * np.pi * 5 * t)
    xhigh = np.sin(2 * np.pi * 250 * t)
    x = xlow + xhigh
    ds = Dataset(x, sa={'sid': np.arange(len(x))}, fa={'fid':['theone']})

    # butterworth filter with a cutoff between the waves
    from scipy import signal
    b, a = signal.butter(8, 0.125)
    mds = iir_filter(ds, b, a, padlen=150)
    # check we get just the slow wave out (compensate for edge artifacts)
    assert_false(np.sum(np.abs(mds.samples[100:-100,0] - xlow[100:-100]) > 0.001))
    assert_equal(len(ds.sa), len(mds.sa))
    assert_equal(len(ds.fa), len(mds.fa))
    assert_array_equal(ds.fa.fid, mds.fa.fid)
    assert_array_equal(ds.sa.sid, mds.sa.sid)

max_stars_repo_name: nno/PyMVPA

max_issues_repo_name: nno/PyMVPA

max_forks_repo_name: nno/PyMVPA


text:
#!/usr/bin/env python3
import argparse
import librosa
from tqdm import tqdm
import io
import logging
from pathlib import Path
import pandas as pd
import numpy as np
import soundfile as sf
from pypeln import process as pr
import gzip
import h5py

parser = argparse.ArgumentParser()
parser.add_argument('input_csv')
parser.add_argument('-o', '--output', type=str, required=True)
parser.add_argument('-c', type=int, default=4)
parser.add_argument('-sr', type=int, default=22050)
parser.add_argument('-col',
                    default='filename',
                    type=str,
                    help='Column to search for audio files')
parser.add_argument('-cmn', default=False, action='store_true')
parser.add_argument('-cvn', default=False, action='store_true')
parser.add_argument('-winlen',
                    default=40,
                    type=float,
                    help='FFT duration in ms')
parser.add_argument('-hoplen',
                    default=20,
                    type=float,
                    help='hop duration in ms')

parser.add_argument('-n_mels', default=64, type=int)
ARGS = parser.parse_args()

DF = pd.read_csv(ARGS.input_csv, sep='\t',
                 usecols=[0])  # only read first cols, allows to have messy csv

MEL_ARGS = {
    'n_mels': ARGS.n_mels,
    'n_fft': 2048,
    'hop_length': int(ARGS.sr * ARGS.hoplen / 1000),
    'win_length': int(ARGS.sr * ARGS.winlen / 1000)
}

EPS = np.spacing(1)

def extract_feature(fname):
    """extract_feature
    Extracts a log mel spectrogram feature from a filename, currently supports two filetypes:

    1. Wave
    2. Gzipped wave

    :param fname: filepath to the file to extract
    """
    ext = Path(fname).suffix
    try:
        if ext == '.gz':
            with gzip.open(fname, 'rb') as gzipped_wav:
                y, sr = sf.read(io.BytesIO(gzipped_wav.read()),
                                dtype='float32')
                # Multiple channels, reduce
                if y.ndim == 2:
                    y = y.mean(1)
                y = librosa.resample(y, sr, ARGS.sr)
        elif ext in ('.wav', '.flac'):
            y, sr = sf.read(fname, dtype='float32')
            if y.ndim > 1:
                y = y.mean(1)
            y = librosa.resample(y, sr, ARGS.sr)
    except Exception as e:
        # Exception usually happens because some data has 6 channels , which librosa cant handle
        logging.error(e)
        logging.error(fname)
        raise
    lms_feature = np.log(librosa.feature.melspectrogram(y, **MEL_ARGS) + EPS).T
    return fname, lms_feature

with h5py.File(ARGS.output, 'w') as store:
    for fname, feat in tqdm(pr.map(extract_feature,
                                   DF[ARGS.col].unique(),
                                   workers=ARGS.c,
                                   maxsize=4),
                            total=len(DF[ARGS.col].unique())):
        basename = Path(fname).name
        store[basename] = feat

max_stars_repo_name: Diego-II/Datadriven-GPVAD

max_issues_repo_name: Diego-II/Datadriven-GPVAD

max_forks_repo_name: Diego-II/Datadriven-GPVAD

reason: import numpy


text:
import mxnet as mx
from mxnet import gluon
import numpy as np
from scipy.sparse import coo_matrix

def peak_signal_to_noise_ratio(max_value=1.):
    def compute_psnr(predictions, labels):
        mean_squared_error = gluon.loss.L2Loss(weight=2.)(predictions, labels)
        with mean_squared_error.context:
            max_value_array = mx.nd.array([max_value**2])
        return 10 * (
            mx.nd.log10(max_value_array) - mx.nd.log10(mean_squared_error))

    return compute_psnr

def softmax_cross_entropy_segmentation():
    def compute_cross_entropy_segmentation(predictions, labels):
        with labels.context:
            weight = (labels >= 0)
        loss = gluon.loss.SoftmaxCrossEntropyLoss(axis=1)(predictions, labels,
                                                          weight)
        return loss * weight.size / mx.nd.sum(weight)

    return compute_cross_entropy_segmentation

def average_end_point_error():
    def compute_aepe(predictions, labels):
        labels, masks = labels
        epe = mx.nd.sqrt(mx.nd.sum((labels - predictions)**2, axis=1))
        epe = epe * masks[:, 0]
        aepe = mx.nd.sum(
            epe, axis=(1, 2)) / mx.nd.sum(
                masks[:, 0], axis=(1, 2))
        return mx.nd.mean(aepe)

    return compute_aepe

class LossMetric(object):
    """Customized evaluation metric for per instance loss functions."""

    def __init__(self, loss_function):
        self.loss_function = loss_function
        self.reset()

    def reset(self):
        """Resets the internal evaluation result to initial state."""
        self.number_samples = 0
        self.sum_loss = 0.0

    def update(self, predictions, labels):
        """Updates the loss metric."""
        current_loss = self.loss_function(predictions, labels)
        self.sum_loss += mx.nd.sum(current_loss).asscalar()
        self.number_samples += current_loss.shape[0]

    def get(self):
        """Gets the current evaluation result."""
        if self.number_samples == 0:
            return float('nan')
        return self.sum_loss / self.number_samples

class MiouMetric(object):
    """Mean intersection over union metric for semantic segmentation."""

    def __init__(self, num_classes):
        self.num_classes = num_classes
        self.reset()

    def reset(self):
        """Resets the internal evaluation result to initial state."""
        self.confusion_matrix = 0

    def update(self, predictions, labels):
        """Updates the loss metric."""
        if predictions.shape[1] > 1:
            predictions = mx.nd.argmax(predictions, axis=1)
        else:
            predictions = mx.nd.clip(
                mx.nd.round(predictions), 0, self.num_classes - 1)
        labels = labels.reshape([-1]).asnumpy()
        predictions = predictions.reshape([-1]).asnumpy()
        valid_labels = (labels >= 0)
        self.confusion_matrix += _compute_confusion_matrix(
            predictions[valid_labels], labels[valid_labels], self.num_classes)

    def get(self):
        """Gets the current evaluation result."""
        sum_rows = np.sum(self.confusion_matrix, 0)
        sum_colums = np.sum(self.confusion_matrix, 1)
        diagonal_entries = np.diag(self.confusion_matrix)
        denominator = sum_rows + sum_colums - diagonal_entries

        valid_classes = (denominator != 0)
        num_valid_classes = np.sum(valid_classes)
        denominator += (1 - valid_classes)
        iou = diagonal_entries / denominator
        if num_valid_classes == 0:
            return float('nan')
        return np.sum(iou) / num_valid_classes

def _compute_confusion_matrix(predictions, labels, num_classes):
    if np.min(labels) < 0 or np.max(labels) >= num_classes:
        raise Exception("Labels out of bound.")

    if np.min(predictions) < 0 or np.max(predictions) >= num_classes:
        raise Exception("Predictions out of bound.")

    # Idea borrowed from tensorflow implementation
    values = np.ones(predictions.shape)
    confusion_matrix = coo_matrix((values, (labels, predictions)),
                                  shape=(num_classes, num_classes)).toarray()
    return confusion_matrix


text:
# -*- coding: utf-8 -*-
"""
Created on Sat Jan 11 11:20:02 2014

Example code snippend
@author: Gouthaman Balaraman
"""

from pandas.io.data import DataReader
from datetime import date
import numpy as np
import pandas as pd

# Grab time series data for 5-year history for the stock (here AAPL)
# and for S&P-500 Index
sdate           = date(2008,12,31)
edate           = date(2013,12,31)
df              = DataReader('WFM','yahoo',sdate,edate)
dfb             = DataReader('^GSPC','yahoo',sdate,edate)

# create a time-series of monthly data points
rts             = df.resample('M',how='last')
rbts            = dfb.resample('M',how='last')
dfsm            = pd.DataFrame({'s_adjclose' : rts['Adj Close'],
                        'b_adjclose' : rbts['Adj Close']},
                        index=rts.index)

# compute returns
dfsm[['s_returns','b_returns']] = dfsm[['s_adjclose','b_adjclose']]/ \
    dfsm[['s_adjclose','b_adjclose']].shift(1) -1
dfsm            = dfsm.dropna()

covmat          = np.cov(dfsm["s_returns"],dfsm["b_returns"])

# calculate measures now
beta            = covmat[0,1]/covmat[1,1]
alpha           = np.mean(dfsm["s_returns"])-beta*np.mean(dfsm["b_returns"])

# r_squared = 1. - SS_res/SS_tot
ypred           = alpha + beta * dfsm["b_returns"]
SS_res          = np.sum(np.power(ypred-dfsm["s_returns"],2))
SS_tot          = covmat[0,0]*(len(dfsm)-1) # SS_tot is sample_variance*(n-1)
r_squared       = 1. - SS_res/SS_tot
# 5- year volatiity and 1-year momentum
volatility      = np.sqrt(covmat[0,0])
momentum        = np.prod(1+dfsm["s_returns"].tail(12).values) -1

# annualize the numbers
prd              = 12. # used monthly returns; 12 periods to annualize
alpha           = alpha*prd
volatility      = volatility*np.sqrt(prd)

print beta,alpha, r_squared, volatility, momentum

max_stars_repo_name: gouthambs/karuth-source

max_issues_repo_name: gouthambs/karuth-source

max_forks_repo_name: gouthambs/karuth-source


text:
'''
Control the UVBot from a bluetooth client
using pyBluez with python 3
Yao Li: yaoli90@illinois.edu
'''

import bluetooth
import numpy as np
import matplotlib.pyplot as plt
import math

serverMACAddress = 'B8:27:EB:C4:D0:94'
port = 1
s = bluetooth.BluetoothSocket(bluetooth.RFCOMM)
s.connect((serverMACAddress, port))

size = 1000.
MAP_SIZE_PIXELS = 300
MAP_SIZE_METERS = 10

# Receive data from server
# If the data is too large to receive in one time, receive data with length
# 'size' every time.
def receive_data(data, server, size = 1000.):
    for i in range(math.ceil(MAP_SIZE_PIXELS*MAP_SIZE_PIXELS/size)):
        new_data = server.recv(int(size))
        new_data_bytes = bytes(new_data)
        data += new_data_bytes
    return data

if __name__ == "__main__":

    '''
    Commands:

    'm': Display may
    'q': Quit from the client
    'qa' Quit from the client and stop the robot

    '''
    while 1:
        text = input()
        text = bytes(text,'utf-8')

        # unknown commands
        if text not in [b'q', b'qa', b'st', b'sp', b'm',b'fs',b'pp']:
            print('unknown command')
            continue

        if text in [b'st', b'sp']:
            s.send(text)

        if text in [b'fs', b'pp']:
            s.send(text)
        '''
        Quit
        - Command 'q': quit from the client
        - Command 'qa': quit from the client and stop the robot
        '''
        if text in [b'q', b'qa']:
            s.send(text)
            out_file = open("output.txt","wb")
            out_file.write(map)
            out_file.close()
            break

        '''
        Display Map
        - Type command 'm' to display the current map.
        - The map updates once per second. Close the figure to stop the map display
          and go back to command input level.
        '''
        if text == b'm':
            fig = plt.figure()
            FIGURE_OPEN = True
            while FIGURE_OPEN:
                text = b'm'
                s.send(text)

                map = b''
                pos = b''

                pos = receive_data(pos, s) # receive robot position from server
                map = receive_data(map, s) # receive map from server

                # convert bytes to image
                pos_image = np.fromstring(pos, dtype=np.uint8)
                pos_image = pos_image.reshape((MAP_SIZE_PIXELS,MAP_SIZE_PIXELS))

                map_image = np.fromstring(map, dtype=np.uint8)
                map_image = map_image.reshape((MAP_SIZE_PIXELS,MAP_SIZE_PIXELS))

                # construct an image: pos in red, map in blue, green is blank
                # blue: free space in the map
                # white: robot trajectory
                # green: obstacle
                blank = 255*np.zeros((MAP_SIZE_PIXELS, MAP_SIZE_PIXELS), dtype=np.uint8)
                image = np.stack((pos_image, blank, map_image), axis=-1)

                plt.clf()
                plt.imshow(image)
                plt.pause(1)
                # check whether the figure is still open
                FIGURE_OPEN = plt.fignum_exists(fig.number)

    s.close()

max_stars_repo_name: wangfanxin/UVBot

max_issues_repo_name: wangfanxin/UVBot

max_forks_repo_name: wangfanxin/UVBot


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import spdiags
from scipy.sparse import linalg as lg

# Quadratic function
class Quad(object):
    def __init__(self, Q, p):
        self.Q = Q
        self.p = p

    def func(self, x):
        r_0 = self.p - np.dot(self.Q, x)
        return r_0

# Solver
class MINRES(object):

    def __init__(self, A, b, x_0, TOL, MAXIT):
        self.func = Quad(A, b)
        self.TOL = TOL
        self.MAXIT = MAXIT
        self.Q = self.func.Q
        self.x = x_0
        self.r_vec = [self.func.func(self.x)]
        self.p_vec = [self.r_vec[-1]]
        self.Qr_vec = [np.dot(self.Q, self.r_vec[-1])]
        self.res_vec = [np.linalg.norm(self.r_vec[-1])]
        self.Qp = np.dot(self.Q, self.p_vec[-1])
        self.k = 1

    def calc(self, ):
        while self.res_vec[-1] > self.TOL and self.k < self.MAXIT:
            alpha = np.divide(np.dot(self.r_vec[-1].T, self.Qr_vec[-1]), np.dot(self.Qp.T, self.Qp))
            self.x += alpha * self.p_vec[-1]
            self.r_vec.append(self.r_vec[-1] - (alpha * self.Qp))
            self.Qr_vec.append(np.dot(self.Q, self.r_vec[-1]))
            self.res_vec.append(np.linalg.norm(self.r_vec[-1]))
            beta = np.divide(np.dot(self.r_vec[-1].T, self.Qr_vec[-1]), np.dot(self.r_vec[-2].T, self.Qr_vec[-2]))
            self.p_vec.append(self.r_vec[-1] + np.dot(beta, self.p_vec[-1]))
            self.Qp = np.dot(self.Q, self.p_vec[-1])
            self.k += 1
        return self.res_vec, self.k, self.x

# Generate Matrix
class gen_mat(object):

    def __init__(self, n):
        self.n = n

    def gen_diag(self, ):
        data = np.array([- np.ones(self.n), 2 * np.ones(self.n), - np.ones(self.n)])
        diags = np.array([-1, 0, 1])
        M = spdiags(data, diags, self.n, self.n).toarray()
        return M

def main(n):
    mat = gen_mat(n)
    run = MINRES(mat.gen_diag(), -np.ones(n), np.zeros(n), 1e-5, 2000)
    sci = lg.minres(mat.gen_diag(), -np.ones(n), np.zeros(n), tol=1e-5, maxiter=2000)
    result = run.calc()
    print('result :', result[2])
    print('result_sci :', sci[0])
    plt.semilogy(range(0, result[1]), np.log(result[0]), label='MINRES')
    plt.legend()
    plt.ylabel('Residuals')
    plt.xlabel('Iterations')
    plt.show()

if __name__ == '__main__':
    main(25)


text:
import math
import numpy as np
import os
from . import BaseSequenceGenerator
from keras.utils import Sequence
from keras import backend as K
from keras.preprocessing.sequence import pad_sequences

class TxtFileSequenceGenerator(BaseSequenceGenerator):
    """Generate sequence of text files from the filesystem to the Keras."""

    def __init__(self, X, y, padding_size=0, padding_position='post', text_sep=' ', preprocessors=[], *args, **kwargs):
        """Class constructor.

        Arguments:
            X {list} -- Sample data ids.
            y {list} -- Encoded data labels.

        Keyword Arguments:
            preprocessors {list} -- An objects list to preprocess the samples data. (default: {[]})
        """

        # Initialize superclass parameters first
        super(TxtFileSequenceGenerator, self).__init__(X, y, *args, **kwargs)

        self.preprocessors = preprocessors
        self.padding_size = padding_size
        self.padding_position = padding_position
        self.text_sep = text_sep

    def process(self, id):
        """Read and process a csv file.

        Arguments:
            id {srt} -- Text file id.

        Returns:
            [numpy.ndarray] -- Text file content.
        """

        with open(id) as f:
            content = " ".join(line.rstrip() for line in f)

        # check to see if our preprocessors are not Empty
        if len(self.preprocessors) > 0:

            # loop over the preprocessors and apply each to the content file
            for p in self.preprocessors:
                content = p.preprocess(content)

        content = content.split(self.text_sep)

        if self.padding_size > 0:
            content = pad_sequences([content], maxlen=self.padding_size, padding='post')[0]
            # TODO - informar com um warning que texto foi truncado se self.padding_size < len(content)

        return np.array(content)

max_stars_repo_name: kleysonr/snsdl

max_issues_repo_name: kleysonr/snsdl

max_forks_repo_name: kleysonr/snsdl


text:
import unittest
from slam.FactorGraphSimulator import G2oToroPoseGraphReader
import numpy as np
from slam.Variables import Variable

class TestCase(unittest.TestCase):
    def test_g2o(self) -> None:
        # constructor
        file_path ="/home/chad/Research/optimalTransport/CouplingSLAM/example/slam/MIT_g2o_dataset/trimmed_MIT.g2o"
        pg = G2oToroPoseGraphReader(file_path)
        self.assertEqual(pg.file_path, file_path)
        self.assertEqual(pg.file_type, "g2o")
        self.assertEqual(pg.node_head, "VERTEX_SE2")
        self.assertEqual(pg.factor_head, "EDGE_SE2")

        self.assertEqual(pg.node_list[0], Variable("0", 3))
        self.assertEqual(pg.node_list[-1], Variable("10", 3))

        first_factor = pg.factor_list[0]
        last_factor = pg.factor_list[-1]
        self.assertEqual(first_factor.vars[0], Variable("0", 3))
        self.assertEqual(first_factor.vars[1], Variable("1", 3))
        self.assertIsNone(
            np.testing.assert_almost_equal(first_factor.observation,
                                           np.array([2.039345, 0.003006, 0.014452])))
        info_mat = np.array([[1.778126, 0.026853,0.000000],
                             [0.026853, 3.846788,0.000000],
                             [0.000000, 0.000000,388.684289]])
        self.assertIsNone(
            np.testing.assert_almost_equal(first_factor.noise_cov,
                                           np.linalg.inv(info_mat)))
        self.assertEqual(last_factor.vars[0], Variable("9", 3))
        self.assertEqual(last_factor.vars[1], Variable("4", 3))
        obs = np.array([-7.500000, -8.000000, 1.570796])
        self.assertIsNone(
            np.testing.assert_almost_equal(last_factor.observation,obs))
        info_mat = np.array([[1.008417, -0.820651,0.000000],
                             [-0.820651, 0.902417,0.000000],
                             [0.000000, 0.000000,60.523586]])
        self.assertIsNone(
            np.testing.assert_almost_equal(last_factor.noise_cov,
                                           np.linalg.inv(info_mat)))

    def test_toro(self) -> None:
        # constructor
        file_path ="/home/chad/Research/optimalTransport/CouplingSLAM/example/slam/CSAIL_toro_dataset/CSAIL_P_toro.graph"
        pg = G2oToroPoseGraphReader(file_path)
        self.assertEqual(pg.file_path, file_path)
        self.assertEqual(pg.file_type, "graph")
        self.assertEqual(pg.node_head, "VERTEX2")
        self.assertEqual(pg.factor_head, "EDGE2")

        self.assertEqual(pg.node_list[0], Variable("0", 3))
        self.assertEqual(pg.node_list[-1], Variable("1044", 3))

        first_factor = pg.factor_list[0]
        last_factor = pg.factor_list[-1]
        self.assertEqual(first_factor.vars[0], Variable("0", 3))
        self.assertEqual(first_factor.vars[1], Variable("1", 3))
        obs = np.array([0.082760, 0.003050, 0.284020])
        self.assertIsNone(
            np.testing.assert_almost_equal(first_factor.observation,obs))

        info_mat = np.array([[123.488234, -2144.807885,0.000000],
                             [-2144.807885, 58242.575768,0.000000],
                             [0.000000, 0.000000,6065.357771]])
        self.assertIsNone(
            np.testing.assert_almost_equal(first_factor.noise_cov,
                                           np.linalg.inv(info_mat)))

        self.assertEqual(last_factor.vars[0], Variable("915", 3))
        self.assertEqual(last_factor.vars[1], Variable("1008", 3))
        obs = np.array([0.462450, 0.580490, 1.824460])
        self.assertIsNone(
            np.testing.assert_almost_equal(last_factor.observation,obs))
        info_mat = np.array([[461.494377, -332.244726,0.000000],
                             [-332.244726, 309.128718,0.000000],
                             [0.000000, 0.000000,1253.513867]])
        self.assertIsNone(
            np.testing.assert_almost_equal(last_factor.noise_cov,
                                           np.linalg.inv(info_mat)))

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

max_stars_repo_name: MarineRoboticsGroup/NF-iSAM

max_issues_repo_name: MarineRoboticsGroup/NF-iSAM

max_forks_repo_name: MarineRoboticsGroup/NF-iSAM


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

def plot2D(array, cmap='viridis'):
    plt.figure(0)
    plt.imshow(np.repeat(array, repeats=5, axis=0), cmap=cmap)
    plt.show()

def plot3D(array, n_cols=3, cmap='viridis'):
    plt.figure(0)
    n, h, w = array.shape
    n_cols = min(n, n_cols)
    n_rows = int(np.ceil(n/n_cols))
    for i in range(n_rows):
        for j in range(n_cols):
            if i*n_cols+j < n:
                plt.subplot2grid((n_rows, n_cols), (i, j))
                plt.imshow(array[i*n_cols+j].reshape((h, w)), cmap=cmap)
    plt.show()

def plot(array, n_cols=3, cmap='viridis'):

    if len(array.shape) == 2:
        plot2D(array, cmap=cmap)
    elif len(array.shape) == 3:
        plot3D(array, n_cols, cmap=cmap)

max_stars_repo_name: arakhmat/unveiler

max_issues_repo_name: arakhmat/unveiler

max_forks_repo_name: arakhmat/unveiler


text:
import sys
import matplotlib
import matplotlib.pyplot as plt
import scipy
from PIL import ImageTk, Image
import Tkinter as tk
import scipy.ndimage as ndimage
import scipy.ndimage.filters as filters
import numpy as np
from load_scan import load_scan, get_pixels_hu
from Lvl_Set_Alg import lvlset

coords = []

def onclick(event):
    # ix, iy = event.xdata, event.ydata
    print 'x = %d, y = %d'%(
        event.x, event.y)
    global coords
    coords.append((event.x, event.y))
    return coords

def main():
    if len(sys.argv) > 1:
        data_dir = sys.argv[1]
    # data_dir = '/hpc/bsha219/lung/Data/Human_PE_Study_HRC/ST12/TLC/Raw/DICOMS'
    patient_scans = load_scan('/hpc/bsha219/lung/Data/Human_PE_Study_HRC/ST12/TLC/Raw/DICOMS')
    patient_images = get_pixels_hu(patient_scans)
    global artery_mask
    artery_mask = []
    for k in range(115, 215):
        global coords
        coords = []
        img = patient_images[k]
        root = tk.Tk()
        im = Image.fromarray(img)
        photo = ImageTk.PhotoImage(im)
        canvas = tk.Canvas(root, width=512, height=512)
        canvas.pack()
        canvas.create_image(256, 256, image=photo)
        canvas.bind("<Button-1>", onclick)
        root.mainloop()
        for i in range(img.shape[0]):
            for j in range(img.shape[0]):
                if img[i][j] > 0:
                    img[i][j] *= 12
        # blurred_f = ndimage.gaussian_filter(img, sigma=3)
        filter_blurred_f = filters.median_filter(img, 3)
        filter_blurred_f2 = ndimage.gaussian_filter(filter_blurred_f, sigma=1)

        alpha = 12
        sharpened = filter_blurred_f + alpha * (filter_blurred_f - filter_blurred_f2)
        seed = np.zeros(img.shape)
        for jj in range(len(coords)):
            seed[int(coords[jj][1]), int(coords[jj][0])] = 1
        seg, phi, its = lvlset(sharpened, seed, max_its=1000, display=True, alpha=0.4)
        segment = seg * 1.0
        # scipy.misc.imsave('/hpc/bsha219/Python/Behdad/Out_arterial_mask/LungMask%.4d.jpg' % k, segment) # binary files
        # matplotlib.image.imsave('/hpc/bsha219/Python/Behdad/Out_arterial_mask/' + str(k) + '.png', segment)  # RGB files
        segment = np.uint8(segment)
        # to save them as a stack of masks in a directory
        binary_im = Image.fromarray(segment)
        binary_im.save('/hpc/bsha219/Python/Behdad/Lung_masks/LungMask%.4d.jpg' % k, quality=100)

if __name__ == '__main__':
    main()

max_stars_repo_name: behdadebsh/LungPhD

max_issues_repo_name: behdadebsh/LungPhD

max_forks_repo_name: behdadebsh/LungPhD

reason: import numpy,import scipy


text:
import hashlib
from random import shuffle

import time

import sys
from sklearn import model_selection
from sklearn.metrics import auc

from sklearn.metrics import roc_curve, log_loss
import numpy as np
from qml.config import *
from qml.helpers import get_engine, load, save

import pandas as pd

class QCV:
    def __init__(self, qmodels):
        _, self.conn = get_engine()
        self.qm = qmodels
        self.nsplits = 5

        self.X_train_cache = {}
        self.Y_train_cache = {}

    def cross_val(self, model_id, data_id, force=False):
        start_time = time.time()

        if not force:
            rows = self.conn.execute(
                "select cv_score, cv_time from qml_results where data_id={} and model_id={}".format(data_id, model_id)
            ).fetchone()
            if rows and rows['cv_score'] is not None:
                return rows['cv_score']

        if data_id > 0:
            if data_id in self.X_train_cache:
                X = self.X_train_cache[data_id]
            else:
                X = pd.read_csv(QML_TRAIN_X_FILE_MASK.format(data_id), index_col=QML_INDEX_COL)
                self.X_train_cache[data_id] = X

            if data_id in self.Y_train_cache:
                Y = self.Y_train_cache[data_id]
            else:
                Y = pd.read_csv(QML_TRAIN_Y_FILE_MASK.format(data_id), index_col=QML_INDEX_COL)
                self.Y_train_cache[data_id] = Y
        else:
            #todo
            if 1 in self.Y_train_cache:
                Y = self.Y_train_cache[1]
            else:
                Y = pd.read_csv(QML_TRAIN_Y_FILE_MASK.format(1), index_col=QML_INDEX_COL)
                self.Y_train_cache[1] = Y

        splits = load(self._get_splits_filename(self.nsplits, Y.index))
        if splits is None:
            splits = []
            kf = model_selection.KFold(n_splits=self.nsplits, shuffle=True)
            for ids1, ids2 in kf.split(Y):
                splits.append([sorted(Y.index[ids1]), sorted(Y.index[ids2])])
            save(splits, self._get_splits_filename(self.nsplits, Y.index))

        scores = []
        i = 0
        for train_indexes, test_indexes in splits:
            i += 1
            if data_id>0:
                # only one data set used, so pass data
                X_train = X.loc[train_indexes]
                Y_train = Y.loc[train_indexes][QML_RES_COL]
                X_test = X.loc[test_indexes]
                Y_test = Y.loc[test_indexes][QML_RES_COL]
                res = self.qm.qpredict(model_id, data_id, data=(X_train, Y_train, X_test), force=force)
            else:
                # many datasets can be used, so pass only ids
                X_train = train_indexes
                X_test = test_indexes
                Y_test = Y.loc[test_indexes][QML_RES_COL]
                res = self.qm.qpredict(model_id, data_id, ids=[X_train, X_test], force=force)

            # fpr, tpr, thresholds = roc_curve(Y_test, res)
            # score = auc(fpr, tpr)
            score = log_loss(Y_test, res.astype(np.float64), eps=1e-14)
            print(i, score)
            sys.stdout.flush()
            scores.append(score)

        total_score = sum(scores) / len(scores)
        total_time = time.time() - start_time

        self.conn.execute(
            """insert into
                    qml_results
                set
                    cv_score={0}, cv_time={1}, data_id={2}, model_id={3}
                on duplicate key update
                    cv_score = values(cv_score),
                    cv_time = values(cv_time)
            """.format(
                round(total_score, 6), round(total_time, 1), data_id, model_id
            )
        )

        return total_score

    def _get_splits_filename(self, nsplists, train_ids):
        """ all datasets with the same ids have the same splits"""
        ids_hash = hashlib.md5(
            '_'.join([str(i) for i in sorted(train_ids)]).encode('utf-8')
        ).hexdigest()
        filename = QML_DATA_DIR + 'cv_splits/splits_{}_trainlen{}__h_{}.csv'.format(
            nsplists, len(train_ids), ids_hash
        )
        return filename


text:
#!./python27-gcc482/bin/python
# coding: utf-8
"""
BAIDU CLOUD action
"""

import os
import sys
import pickle
import json
import time
import functools

import numpy as np

from utils.preprocess import get_images
from utils.config_utils import parse_config, print_configs
import mfcc.feature_extractor as mfcc_extractor

import models.pptsm_infer as image_model
import models.audio_infer as audio_model
import models.bmn_infer as prop_model
import models.lstm_infer as classify_model

import logger
logger = logger.Logger()

def record_time_info(func):
    """decorator func to log cost time for func
    """
    @functools.wraps(func)
    def timer(*args):
        """log cost time for func
        """
        logger.info("function [{}] processing ...".format(func.__name__))
        start_time = time.time()
        retval = func(*args)
        cost_time = round(time.time() - start_time, 5)
        logger.info("function [{}] run time: {:.2f} min".format(func.__name__, cost_time / 60))
        return retval
    return timer

class ActionDetection(object):
    """ModelPredict"""
    def __init__(self, cfg_file="configs/configs.yaml"):
        cfg = parse_config(cfg_file)
        self.configs = cfg
        print_configs(self.configs, "Infer")

        name = 'COMMON'
        self.DEBUG          = cfg[name]['DEBUG']
        self.BMN_ONLY       = cfg[name]['BMN_ONLY']
        self.LSTM_ONLY      = cfg[name]['LSTM_ONLY']
        self.PCM_ONLY       = cfg[name]['PCM_ONLY']
        if self.LSTM_ONLY:
            self.prop_dict = {}
            for dataset in ['EuroCup2016']:
                prop_json = '/home/work/datasets/{}/feature_bmn/prop.json'.format(dataset)
                json_data = json.load(open(prop_json, 'r'))
                for item in json_data:
                    basename = prop_json.replace('feature_bmn/prop.json', 'mp4')
                    basename = basename + '/' + item['video_name'] + '.mp4'
                    self.prop_dict[basename] = item['bmn_results']

    @record_time_info
    def load_model(self):
        """
        load_model
        """
        if not self.DEBUG:
            self.image_model = image_model.InferModel(self.configs)
            if not self.PCM_ONLY:
                self.audio_model = audio_model.InferModel(self.configs)

        if not self.LSTM_ONLY:
            self.prop_model = prop_model.InferModel(self.configs)

        if not self.BMN_ONLY:
            self.classify_model = classify_model.InferModel(self.configs)

        logger.info("==> Action Detection prepared.")

    @record_time_info
    def infer(self, imgs_path, pcm_path, fps=5):
        """
        extract_feature
        """
        self.imgs_path = imgs_path
        self.pcm_path = pcm_path
        self.configs['COMMON']['fps'] = fps

        logger.info("==> input video {}".format(os.path.basename(self.imgs_path)))

        # step 1: extract feature
        video_features = self.extract_feature()

        # step2: get proposal
        bmn_results = self.extract_proposal(video_features)

        # step3: classify
        material = {'feature': video_features, 'proposal': bmn_results}
        action_results = self.video_classify(material)

        return bmn_results, action_results

    @record_time_info
    def video_classify(self, material):
        """video classify"""
        if self.BMN_ONLY:
            return []
        action_results = self.classify_model.predict(self.configs, material=material)
        logger.info('action shape {}'.format(np.array(action_results).shape))
        return action_results

    @record_time_info
    def extract_proposal(self, video_features):
        """extract proposal"""
        if self.LSTM_ONLY:
            basename = self.imgs_path.replace('frames', 'mp4') + '.mp4'
            bmn_results = self.prop_dict[basename]
            return bmn_results
        bmn_results = self.prop_model.predict(self.configs, material=video_features)
        logger.info('proposal shape {}'.format(np.array(bmn_results).shape))
        return bmn_results

    @record_time_info
    def extract_feature(self):
        """extract feature"""
        if not self.DEBUG:
            image_path_list = get_images(self.imgs_path)
            self.configs['PPTSM']['frame_list'] = image_path_list
            self.configs['AUDIO']['pcm_file'] = self.pcm_path
            image_features = self.image_model.predict(self.configs)
            if self.PCM_ONLY:
                sample_rate = self.configs['AUDIO']['sample_rate']
                pcm_features = mfcc_extractor.extract_pcm(self.pcm_path, sample_rate)
                audio_features = []
            else:
                audio_features, pcm_features = self.audio_model.predict(self.configs)

            np_image_features = np.array(image_features, dtype=np.float32)
            np_audio_features = np.array(audio_features, dtype=np.float32)
            np_pcm_features = np.array(pcm_features, dtype=np.float32)

            video_features = {'image_feature': np_image_features,
                              'audio_feature': np_audio_features,
                              'pcm_feature': np_pcm_features}
        else:
            feature_path = self.imgs_path.replace("frames", "features") + '.pkl'
            video_features = pickle.load(open(feature_path, 'rb'))

        logger.info("feature shape {} {} {}".format(video_features['image_feature'].shape,
                                                    video_features['audio_feature'].shape,
                                                    video_features['pcm_feature'].shape))

        return video_features

if __name__ == '__main__':

    model_predict = ActionDetection(cfg_file="../configs/configs.yaml")
    model_predict.load_model()

    imgs_path = "/home/work/datasets/EuroCup2016/frames/1be705a8f67648da8ec4b4296fa80895"
    pcm_path = "/home/work/datasets/EuroCup2016/pcm/1be705a8f67648da8ec4b4296fa80895.pcm"

    bmn_results, action_results = model_predict.infer(imgs_path, pcm_path)
    results = {'bmn_results': bmn_results, 'action_results': action_results}

    with open('results.json', 'w', encoding='utf-8') as f:
       data = json.dumps(results, indent=4, ensure_ascii=False)
       f.write(data)