text:
"""trying to develop a function to get the convulsion of the time course"""
# noinspection PyUnresolvedReferences

from __future__ import print_function, division
import scipy.stats
from scipy.stats import gamma
import numpy as np

# hrf function from lecture
def hrf(tr_times):
    """ Return values for HRF at given times """
    # Gamma pdf for the peak
    peak_values = gamma.pdf(tr_times, 6)
    # Gamma pdf for the undershoot
    undershoot_values = gamma.pdf(tr_times, 12)
    # Combine them
    values = peak_values - 0.35 * undershoot_values
    # Scale max to 0.6
    return values / np.max(values) * 0.6

# convolution and deal with overflow
def convolution(neural, hrf_at_trs):
    convolved = np.convolve(neural, hrf_at_trs)
    convolved = convolved[:-(len(hrf_at_trs)-1)]
    return convolved

# get all convolved result and save in a dictionary
def get_all_convolved (dictionary_of_all_neural, hrf_at_trs, path ="."):
    result = {}
    dict_key = dictionary_of_all_neural.keys()
    dict_key = sorted(dict_key)
    for i in dict_key:
        convolved = convolution(dictionary_of_all_neural[i][0], hrf_at_trs)
        result[i] = convolved
        np.savetxt(path + "convolved_%s.txt" % i, convolved)
        print ("convolution for %s: complete" % i)
    return result

# remove outlier in convolved results
def remove_outlier (subid ,convolsion_data_dict, outlier_directory, axis=0):
    """
    Parameters
    ----------
    subid = subject ID
    convolsion_data_dict: the convolsion dictionary
    outlier_directory: outlier dictionary
    axis: in this case, axis=0

    Returns
    -------
    dictionary of convolsion result without outlier
    """
    result = {}
    for i in convolsion_data_dict.keys():
        result[i] = np.delete(convolsion_data_dict[i], outlier_directory[subid + '_' + i[:6]], axis)
    return result

reason: import numpy,import scipy,from scipy


text:
import unittest
import os.path as osp
import numpy as np
from terminaltables import AsciiTable
from tqdm import tqdm

from retinal.config import get_cfg
from retinal.data.fgadr_dataset import FGADRDataset
from retinal.data.data_transform import retinal_lesion

class TestFGADRDataset(unittest.TestCase):
    def test(self):
        data_root = "./data/FGADR-Seg/Seg-set"
        samples_path = osp.join(data_root, "test.txt")
        classes_path = osp.join(data_root, "classes.txt")
        cfg = get_cfg()
        transformer = retinal_lesion(cfg)

        dataset = FGADRDataset(
            data_root, samples_path, classes_path,
            data_transformer=None
        )

        total_area = 0
        class_area = [0 for _ in range(6)]
        for i in tqdm(range(len(dataset))):
            img, target = dataset[i]
            total_area += target.shape[0] * target.shape[1]
            for j in range(6):
                class_area[j] += np.sum(target[:, :, j])

        table_data = [["id"] + ["Class"] + ["Area"] + ["Proportion"]]
        for i in range(6):
            table_data.append(
                [i] + [dataset.classes[i]]
                + [class_area[i]]
                + ["{:.4f}".format(class_area[i] / total_area)]
            )
        table_data.append(
            [""] + ["all"] + [total_area] + [""]
        )
        table = AsciiTable(table_data)
        print(table.table)

        self.assertTrue(True)

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

ext: py

max_stars_repo_name: by-liu/RetinalApp

max_issues_repo_name: by-liu/RetinalApp

max_forks_repo_name: by-liu/RetinalApp

reason: import numpy


text:
import numpy as np
from astropy.io import ascii
import os
from tqdm import tqdm
from astropy.time import Time
import pandas as pd
import warnings
import time
import random
from gatspy import periodic, datasets
import matplotlib.pyplot as plt
from astropy.table import Table
from gatspy import datasets, periodic
import scipy.stats as sci_stat
global data_path, utc_tag

data_path = 'data/plasticc/data/'

# Caution, be careful when ignoring warnings!
import sys
import warnings

if not sys.warnoptions:
    warnings.simplefilter("ignore")

import argparse
import datetime
utc = datetime.datetime.utcnow()
utc_tag = f"{utc.year}_{utc.month}_{utc.day}_{utc.second}"

parser = argparse.ArgumentParser(description="LSP Summarize Arguments")
parser.add_argument('-N', '--NumberInjected', type=int, help='Number of injected light curves (int)')
parser.add_argument('-class', '--classType', type=str, help='Class type. Currently supports rrl, eb, agn, tde')
parser.add_argument('-kmax', '--MaximumFourierComponents', type=int, help='Maximum number of fourier components (>1)')
parser.add_argument('-k_term_base', '--FreqTermsBaseAll', type=int, help='Number of frequency terms to use for the base model common to all bands')
parser.add_argument('-fmin', '--MinSearchPeriod', type=float, help='Minimum search period')
parser.add_argument('-fmax', '--MaxSearchPeriod', type=float, help='Maximum search period')
parser.add_argument('-dur', '--duration', type=int, help='Baseline duration (default 365 days)')
parser.add_argument('-dets', '--detection', type=str, help='Detection type (default all)')

args = parser.parse_args()

def generate_toi_table(data, meta_info, meta_theta_EB, meta_theta_RRL):
    """
    Generate table that contains the light curve ID and transient type.

    Input
    -----
    data: Head data table that contains photometry
    meta_info: Table that contains the meta-data (i.e classification name)
    meta_theta_<TYPE>: Table that contains metadata information (i.e Period)

    """
    id_av_rrl, id_av_eb = [], []

    for uid in tqdm(np.unique(data['object_id'])):
        ww = np.where(meta_theta_EB['object_id'] == uid)
        if np.shape(ww)[-1]==1:
            id_av_eb.append(uid)

    for uid in tqdm(np.unique(data['object_id'])):
        ww = np.where(meta_theta_RRL['object_id'] == uid)
        if np.shape(ww)[-1]==1:
            id_av_rrl.append(uid)

    id_av_rrl, id_av_eb = np.array(id_av_rrl), np.array(id_av_eb)

    _id1 = np.array(['rrl' for _ in range(len(id_av_rrl))])
    _id2 = np.array(['eb' for _ in range(len(id_av_eb))])

    # All ID's and & ID tags
    all_id = np.concatenate([id_av_rrl, id_av_eb])
    _id_all = np.concatenate([_id1, _id2])

    # Final TOI table
    toi_table = Table([all_id, _id_all], names=('obj_id', 'type'))

    return toi_table

# Read ALL PlastiCC data & meta
data = pd.read_csv(data_path + "plasticc_train_lightcurves.csv.gz",
                  compression='gzip',
                  error_bad_lines=False)

meta_info = ascii.read(data_path + "plasticc_train_metadata.csv") # ascii meta since it's smaller
meta_theta_EB = ascii.read(data_path + 'plasticc_modelpar/' + 'plasticc_modelpar_016_EB.csv')
meta_theta_RRL = ascii.read(data_path + 'plasticc_modelpar/' + 'plasticc_modelpar_092_RRL.csv')

# Fetch all toi's
toi_table = generate_toi_table(data, meta_info, meta_theta_EB, meta_theta_RRL)

# Helper functions to read and digest plasticc data
def generate_lc(obj_id, band='all', data_table=data, det=1):
    """Unpack and return PlastiCC data in numpy array format.

    Paremters
    ---------
    obj_id: Object ID
    band: Photometric bandpass filter. 'all' includes ugrizy, or 'ugrizy'
    data_table: Pandas data table containing the light curves
    det: Detection from the image subtraction algorithm. ==1 detection, ==0 not detection (i.e upper limit) or 'all': uses both 0 & 1

    Returns
    -------
    mjd, magnitude, magnitude_error, filter (if band=='all')
    """

    if det==0 or det==1:
        data_table_mod = data_table[data_table['detected_bool']==det]
    elif det=='all':
        data_table_mod = data_table # select both

    # Select light curve based on the ID
    lc = data_table_mod[data_table_mod['object_id']==obj_id]

    lsst_bands = list('ugrizy') # lsst photomeric bands

    lc_array = lc.to_numpy()

    # Capture empty light curve
    assert len(lc_array[:,1])>0, ("Sorry, it seems like your obj_id query was wrong!")

    mjd, flux, flux_err = lc_array[:,1], lc_array[:,3], lc_array[:,4]
    flt = lc_array[:,2].astype(int).astype(str)

    for j in range(6):
        flt[flt==str(j)] = lsst_bands[j]

    if band=='all':
        return mjd, flux, flux_err, flt
    else:
        return mjd[flt==band], flux[flt==band], flux_err[flt==band], flt[flt==band]

def fetch_type(lid, table=toi_table):
    """Fetch the classification type of transient given light curve table and original head TOI table."""
    return table[table['obj_id']==lid]

def fetch_meta_info(lc_id, lc_type):
    """Fetch metadata for transient type.

    Input
    -----
    lc_id: Light curve ID
    lc_type: classification type (i.e rrl, eb)

    Output
    ------
    meta_<type>_table: Table that contains metadata (i.e period and other physical properties)
    """
    if lc_type=='rrl':
        # crossmatch to approprirate table
        xm_ = np.where(meta_theta_RRL['object_id']==lc_id)
        return meta_theta_RRL[xm_]
    elif lc_type=='eb':
        # crossmatch to approprirate table
        xm_ = np.where(meta_theta_EB['object_id']==lc_id)
        return meta_theta_EB[xm_]

# Write a function that will generate N random from each class (equal)
def draw_rand_trans(table, N, class_type='rrl'):
    """Given N this function will draw an equal number of trnasinets.
       Note: It will not draw the same transiennt
    """
    # isolate each unique class
    req_tab = table[table['type']==class_type]

    # Random number generator w/o repeat
    rng = np.random.default_rng()
    rn = rng.choice(len(req_tab), size=N, replace=False)

    return req_tab[rn]

def run_multi_lsp(x, y, err, fts, fmin=0.1, fmax=150, k_band=1, m_base=1, mode='fast', dt_cut=365, k_term_base=0):
    """Run all methods of multiband gatspy Lomb-Scargle Periodogram.

        Input
        ------
        x, y, err, fts: phase, magnitudes/flux, error, filter list
        fmin, fmax: minimum and maximum search period in the Lomb-Scargle
        k (int): Number of Fourier components
        mode (str): LSP method. Currently supports 'fast' and "general"
        dt_cut (int): Maximum baseline time. Default is 1 year worth of photometry.

        Output
        ------
        best_period: Best period found from the highest peak in the LSP
        TODO: False Alarm Probability, TOP N peaks?!
    """

    try:
        # Pre-processing to photometry
        dt = x-x[0] # calculate baseline
        x, y, err, fts = x[dt<=dt_cut], y[dt<=dt_cut], err[dt<=dt_cut], fts[dt<=dt_cut]
        y += -1*min(y) # TODO: PLASTICC light curves can be negative. For now normalize such that they're at least positive
        dt = x-x[0] # evaluate baseline again!
        # Check fmax limit
        if max(dt)<fmax:
            fmax = max(dt)-3
    except:
        return np.nan

    if mode=='fast':
        try:
            model = periodic.LombScargleMultibandFast(fit_period=True,optimizer_kwds={"quiet": True},
                                  Nterms=k)
            model.optimizer.set(period_range=(fmin, fmax))
            model = model.fit(x, y, dy=err, filts=fts)
            return model.best_period
        except:
            return np.nan
    elif mode=='general':
        try:
            model = periodic.LombScargleMultiband(fit_period=True,optimizer_kwds={"quiet": True},
                      Nterms_base=m_base, Nterms_band=k_band)
            model.optimizer.set(period_range=(fmin, fmax))
            model = model.fit(x, y, dy=err, filts=fts)
            return model.best_period
        except:
            return np.nan

def run_single_lsp(x, y, err, fts, band='u', fmin=0.1, fmax=150, k_band=1, mode='fast', dt_cut=365):
    """Run all methods of single-band gatspy Lomb-Scargle Periodogram.

        Input
        ------
        x, y, err, fts: phase, magnitudes/flux, error, filter list
        band (str): Photometric band you want to run LSP. Currently supports 'ugrizy'.
        fmin, fmax: minimum and maximum search period in the Lomb-Scargle
        k (int): Number of Fourier components
        mode (str): LSP method. Currently supports 'fast' and "general"
        dt_cut (int): Maximum baseline time. Default is 1 year worth of photometry.

        Output
        ------
        best_period: Best period found from the highest peak in the LSP
        TODO: False Alarm Probability, TOP N peaks?!
    """

    try:
        # Pre-processing to photometry
        dt = x-x[0] # calculate transient duration
        x, y, err, fts = x[dt<=dt_cut], y[dt<=dt_cut], err[dt<=dt_cut], fts[dt<=dt_cut]
        y += -1*min(y)
        dt = x-x[0] # updated dt

        # isolate photometric band
        x, y, err = x[fts==band], y[fts==band], err[fts==band]

        # Check fmax limit
        if max(dt)<fmax:
            fmax = max(dt)-5

    except:
        return np.nan

    if mode=='fast':
        try:
            model = periodic.LombScargleFast(fit_period=True,optimizer_kwds={"quiet": True},
                                  Nterms=1)
            model.optimizer.set(period_range=(fmin, fmax))
            model = model.fit(x, y, dy=err)
            return model.best_period
        except:
            return np.nan
    elif mode=='general':
        try:
            model = periodic.LombScargle(fit_period=True,optimizer_kwds={"quiet": True}, Nterms=k_band)
            model.optimizer.set(period_range=(fmin, fmax))
            model = model.fit(x, y, dy=err)
            return model.best_period
        except:
            return np.nan

def generate_tags(kmax):
    """Generate titles for master table on LSP analysis"""
    # Create data table
    m_lsp_name_fast_list = []
    m_lsp_name_gen_list = []
    for i in range(kmax):
        m_lsp_name_fast_list.append('multi_lsp_f'+f'{i+1}')
        m_lsp_name_gen_list.append('multi_lsp_g'+f'{i+1}')

    s_lsp_gen_list = []
    for iii in range(kmax):
        for jj, band_name in enumerate(list('ugrizy')):
            s_lsp_gen_list.append('s_lsp_g'+f'{iii+1}'+f'_{band_name}')

    s_lsp_fast_list = []
    for band_name in list('ugrizy'):
        s_lsp_fast_list.append(f's_lsp_f_{band_name}')

    master_names = np.concatenate([['id'], ['ndet'], ['ptrue'], m_lsp_name_fast_list, m_lsp_name_gen_list, s_lsp_gen_list, s_lsp_fast_list])

    return master_names

def calc_all_lsp(N, transient_class='rrl', fmin=0.1, fmax=150, table=toi_table, dur=365, det_type='all', k=1, m=1):
    """

    Will return a master table (ascii.Table) of both General & Fast LSP (single & multi)

    Parameters:
    -------
    Table content will contain the following parameters:
    lc_id: Light curve ID
    N_det: Number of photometric detections (multi: ugrizy; single: band)
    multi_f1...multi_fk: Multiband fast LSP (kth component)
    multi_g1...multi_gk: Multibamd general LSP (kth component)
    ugrizy_g1...urizy_gk: Single general LSP (kth component)
    urizy_f1: Single fast LSP (k=1)
    p_true: True injected period (0 if not-periodic variable)
    """

    # Generalte N class
    toi_special = draw_rand_trans(table,
                                  N=N,
                                  class_type=transient_class)

    _id_unq = toi_special['obj_id'].data # all the ID's of the unique class
    ndet = np.zeros(shape=N)
    ptrue = np.zeros(shape=N)
    M_multi_fast = np.zeros(shape=(N,k_max_comp))
    M_multi_general = np.zeros(shape=(N,k_max_comp))
    M_single_fast = np.zeros(shape=(N, 6)) # ugrizy columns
    M_single_general = np.zeros(shape=(N, k_max_comp, 6))

    i = 0
    for _id in tqdm(_id_unq):

        # Fetch photometry
        t, mag, mag_err, flt = generate_lc(_id,  det=det_type)

        for j in range(k_max_comp):
            # multi-LSP
            m_f_lsp = run_multi_lsp(t, mag, mag_err, flt, mode='fast', fmin=fmin, fmax=fmax, dt_cut=dur, k_band=k,  m_base=m)
            m_g_lsp = run_multi_lsp(t, mag, mag_err, flt, mode='general', fmin=fmin, fmax=fmax, dt_cut=dur, k_band=k,  m_base=m)

            for ii, flt_lst in enumerate(list('ugrizy')):
                # Single-band lsp per band
                X_g_lsp = run_single_lsp(t, mag, mag_err, flt, band=flt_lst, k_band=j, mode='general', fmin=fmin, fmax=fmax, dt_cut=dur) # general
                M_single_general[i, j, ii] = X_g_lsp # append data

            M_multi_fast[i, j]=m_f_lsp
            M_multi_general[i, j]=m_g_lsp

        # fast single_lsp
        for ii, flt_lst in enumerate(list('ugrizy')):
            X_f_lsp = run_single_lsp(t, mag, mag_err, flt, band=flt_lst, mode='fast', fmin=fmin, fmax=fmax, dt_cut=dur) # general
            M_single_fast[i, ii] = X_f_lsp

        if transient_class=='rrl' or transient_class=='eb':
            f_tp = fetch_meta_info(_id, lc_type=transient_class)['PERIOD'].data[0]
        else:
            f_tp = np.nan

        ptrue[i] = f_tp
        ndet[i] = len(t)

        i+=1

    master_names = generate_tags(k_max_comp) # generate names for table
    Table_master = Table(names=master_names)

    for i in range(N):
        # Collapse multi matrix
        M_LSP_F = M_multi_fast[i]
        M_LSP_G = M_multi_general[i]

        # Collapse single matrix
        S_LSP_U = M_single_general[i, :, 0]
        S_LSP_G = M_single_general[i, :, 1]
        S_LSP_R = M_single_general[i, :, 2]
        S_LSP_I = M_single_general[i, :, 3]
        S_LSP_Z = M_single_general[i, :, 4]
        S_LSP_Y = M_single_general[i, :, 5]

        S_LSP_fast = M_single_fast[i]

        master_var_col = np.concatenate([[_id_unq[i]], [ndet[i]], [ptrue[i]],
                                        M_LSP_F, M_LSP_G,
                                        S_LSP_U, S_LSP_G, S_LSP_R, S_LSP_I, S_LSP_Z, S_LSP_Y,
                                        S_LSP_fast])

        Table_master.add_row(master_var_col)

    # Store table with latest date
    Table_master.write(f"data/{transient_class}_master_N{N}_{utc_tag}", format='ascii')
    return Table_master

def main(Ninj, clf_type, Kmax, baseK, Pmin, Pmax, baseline_dur=365, det_type='all'):
    calc_all_lsp(Ninj, transient_class=clf_type, k=Kmax, m=baseK, fmin=Pmin, fmax=Pmax, dur=baseline_dur, det_type=det_type, table=toi_table)

if __name__ == "__main__":
    main(args.NumberInjected, args.classType, args.MaximumFourierComponents, args.FreqTermsBaseAll, args.MinSearchPeriod, args.MaxSearchPeriod, args.duration, args.detection)

max_stars_repo_name: AndyTza/periodos

max_issues_repo_name: AndyTza/periodos

max_forks_repo_name: AndyTza/periodos

reason: import numpy,import scipy,from astropy


text:
import numpy as np
import cv2
import hand_sample
import filter_frame

cap = cv2.VideoCapture(0)
rows = 480
cols = 640
hs = hand_sample.HandSample(rows, cols)

while True:
    ret, frame = cap.read()
    original = frame.copy()
    hs.draw_hand_sample(frame)
    cv2.imshow('frame', frame)

    if cv2.waitKey(1) & 0xFF == ord('q'):
        histogram = hs.get_hand_sample(original)
        break

ff = filter_frame.FilterFrame()

while True:
    ret, frame = cap.read()
    dst = ff.apply_filter(frame, histogram)
    im2, contours, hierarchy = cv2.findContours(dst, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    sortedContours = np.array(sorted(contours, key=cv2.contourArea))
    biggestContour = sortedContours[-1:]
    cnt = sortedContours[-1]
    frame = ff.getFingertips(frame, cnt)
    cv2.drawContours(frame, biggestContour, -1, (0, 255, 0), 3)

    cv2.imshow('frame', dst)
    cv2.imshow('frame', frame)
    if cv2.waitKey(1) & 0xFF == ord('q'):
        break

cap.release()
cv2.destroyAllWindows()


text:
from pathlib import Path
import pandas as pd
import numpy as np
import jinja2

from image_tagging.tagging_from_caption import get_non_wanted_tags

ROOT = Path(__file__).parents[1]

def get_result_df():
    # result_file = ROOT / 'Image_IDs_Marken_Sortiment.tagged.sample.pickle'
    result_file = ROOT / 'Image_IDs_Marken_Sortiment.tagged.pickle'

    df = pd.read_pickle(result_file)
    df = df.assign(image_name=df['BRANDSHOPPROMOTION_IMAGE_ID'])
    df = df.assign(image_url='https://i.otto.de/i/otto/' + df['image_name'] + '?maxW=300&maxH=300&sm=C&fmt=jpg')
    df = df.assign(index=df.index)

    # remove non downloadable image
    df = df[~df['score'].isnull()]
    return df

template = jinja2.Template('''

Image #{{ index }}

                {{ image_name }}

Tags

            {% if tags %}

                    tag
                    word score

                {% for tag in tags %}

                    {{ tag['name'] }}

                                {{ "%.1f%%"|format(tag['score_word']*100|float) }}

                {% endfor %}

            {% endif %}

''')

def render(row):
    tags_dict = []
    non_wanted = get_non_wanted_tags()
    tag_list = row['score']
    tag_list = sorted(tag_list, key=lambda t: t.word_score, reverse=True)
    for tag in tag_list:
        if tag.tag not in non_wanted:
            tags_dict.append({'name': tag.tag, 'score_caption': tag.caption_score_compensated, 'score_word': tag.word_score})
    output = template.render(
        index=row['index'],
        image_name=row['image_name'],
        image_url=row['image_url'],
        tags=tags_dict,
        description_tags=False
    )
    return output

def render_dataframe(df: pd.DataFrame, part: int, of: int):
    print(f'rendering result {part:03d} of {of:03d}')
    if part == 17:
        a=42

    rows = df.apply(render, axis=1)

    contents = '\n\n'.join(rows)

    html = f"""
<!DOCTYPE html>

<title>Page Title</title>

{contents}

"""

    output_file = ROOT / 'image_tagging/results' / f'results_{part:03d}_of_{of:03d}.html'
    output_file.parent.mkdir(parents=True, exist_ok=True)
    output_file.write_text(html)

def visualize():
    df = get_result_df()

    chunksize = 40
    chunk_count = len(df) // chunksize + 1

    for i, chunk in enumerate(np.array_split(df, chunk_count), 1):
        render_dataframe(chunk, i, chunk_count)

def generate_result_csv():
    df = get_result_df()

    def get_result_cols(row):
        non_wanted = get_non_wanted_tags()
        tag_list = row['score']
        tag_list = list(filter(lambda x: x.tag not in non_wanted, tag_list))
        tag_list = sorted(tag_list, key=lambda t: t.word_score, reverse=True)
        tags = '|'.join([t.tag for t in tag_list])
        confidences = '|'.join([str(t.word_score) for t in tag_list])
        return pd.Series({'tags': tags, 'confidences': confidences})

    df = df.apply(get_result_cols, axis=1).join(df)

    df = df[['image_name', 'tags', 'confidences']]

    output_file = ROOT / 'image_tagging/results' / f'results.csv'
    df.to_csv(output_file, index=False)

if __name__ == '__main__':
    # df = visualize()
    generate_result_csv()
    # a = 42


text:
from .mms_read_feeps_sector_masks_csv import mms_read_feeps_sector_masks_csv
from pytplot import get_data, store_data
import numpy as np

def mms_feeps_remove_sun(sensor_eyes, trange, probe='1', datatype='electron', data_units='intensity', data_rate='srvy', level='l2', suffix=''):

    sector_times, spin_sectors = get_data('mms'+probe+'_epd_feeps_' + data_rate + '_' + level + '_' + datatype + '_spinsectnum'+suffix)
    mask_sectors = mms_read_feeps_sector_masks_csv(trange=trange)
    out_vars = []

    top_sensors = [str(eye) for eye in sensor_eyes['top']]
    bot_sensors = [str(eye) for eye in sensor_eyes['bottom']]

    for sensor in top_sensors:
        var_name = 'mms'+str(probe)+'_epd_feeps_'+data_rate+'_'+level+'_'+datatype+'_top_'+data_units+'_sensorid_'+sensor+'_clean'

        top_data_tuple = get_data(var_name)
        if top_data_tuple is None:
            print('skipping variable: ' + var_name)
            continue
        times, top_data, top_energies = top_data_tuple

        if mask_sectors.get('mms'+probe+'imaskt'+sensor) is not None:
            bad_sectors = mask_sectors['mms'+probe+'imaskt'+sensor]

            for bad_sector in bad_sectors:
                this_bad_sector = np.where(spin_sectors == bad_sector)[0]
                if len(this_bad_sector) is not 0:
                    top_data[this_bad_sector] = np.nan

        try:
            store_data(var_name+'_sun_removed'+suffix, data={'x': times, 'y': top_data, 'v': top_energies})
            out_vars.append(var_name+'_sun_removed'+suffix)
        except Warning:
            continue

    if level is not 'sitl':
        for sensor in bot_sensors:
            var_name = 'mms'+str(probe)+'_epd_feeps_'+data_rate+'_'+level+'_'+datatype+'_bottom_'+data_units+'_sensorid_'+sensor+'_clean'

            bot_data_tuple = get_data(var_name)
            if bot_data_tuple is None:
                print('skipping: ' + var_name)
                continue
            times, bot_data, bot_energies = bot_data_tuple

            if mask_sectors.get('mms'+probe+'imaskb'+sensor) is not None:
                bad_sectors = mask_sectors['mms'+probe+'imaskb'+sensor]

                for bad_sector in bad_sectors:
                    this_bad_sector = np.where(spin_sectors == bad_sector)[0]
                    if len(this_bad_sector) is not 0:
                        bot_data[this_bad_sector] = np.nan

            try:
                store_data(var_name+'_sun_removed'+suffix, data={'x': times, 'y': bot_data, 'v': bot_energies})
                out_vars.append(var_name+'_sun_removed'+suffix)
            except Warning:
                continue

max_stars_repo_name: asherp/pyspedas

max_issues_repo_name: asherp/pyspedas

max_forks_repo_name: asherp/pyspedas


text:
from scipy.integrate import ode
import matplotlib.pyplot as plt
from alu_model import *
from models import *
from parameters import *
import numpy as np

def int_to_bool_list(num):
    return [int(bool(num & (1 << n))) for n in range(4)]

def list_to_int(xs):
    v = 0
    f = 1
    for x in xs:
        v += x*f
        f *= 2
    return v

rho_x = 0

params = delta_L, gamma_L_X, n_y, theta_L_X, eta_x, omega_x, m_x, delta_x, rho_x, gamma_x, theta_x, r_X

# simulation parameters
t_end = 1000
N = t_end

sums = []
x_axis = []
y_axis = []
for a in range(4):
    A = int_to_bool_list(a)
    A0 = A[0]
    A1 = A[1]

    sums.insert(0, [])
    x_axis.append(a)
    y_axis.insert(0, a)
    for b in range(4):
        """
        ADDER

        A0, A1, B0, B1, \
        S0_I0, S0_I1, S0_J0, \
        S1_I0, S1_I1, S1_I2, S1_I3, S1_I4, S1_I5, S1_I6, S1_I7, S1_I8, \
        C_I0, C_I1, C_I2, C_I3, C_I4, C_I5, \
        S0_L_A, S0_L_B, S0_L_I0, S0_L_I1, S0_L_J0, \
        S1_L_A0_0, S1_L_A0_1, S1_L_A1_0, S1_L_A1_1, S1_L_B0_0, S1_L_B0_1, S1_L_B1_0, S1_L_B1_1, S1_L_I0, S1_L_I1, S1_L_I2, S1_L_I3, S1_L_I4, S1_L_I5, S1_L_I6, S1_L_I7, S1_L_I8, \
        C_L_A0_0, C_L_A0_1, C_L_B0_0, C_L_B0_1, C_L_I0, C_L_I1, C_L_I2, C_L_I3, C_L_I4, C_L_I5, \
        S0_N_A, S0_N_B, S0_N_I0, S0_N_I1, S0_N_J0, \
        S1_N_A0_0, S1_N_A0_1, S1_N_A1_0, S1_N_A1_1, S1_N_A1_2, S1_N_A1_3, S1_N_B0_0, S1_N_B0_1, S1_N_B1_0, S1_N_B1_1, S1_N_B1_2, S1_N_B1_3, S1_N_I0, S1_N_I1, S1_N_I2, S1_N_I3, S1_N_I4, S1_N_I5, S1_N_I6, S1_N_I7, S1_N_I8,\
        C_N_A0_0, C_N_A0_1, C_N_A1_0, C_N_A1_1, C_N_B0_0, C_N_B0_1, C_N_B1_0, C_N_B1_1, C_N_I0, C_N_I1, C_N_I2, C_N_I3, C_N_I4, C_N_I5, \
        S0, S1, C
        """

        B = int_to_bool_list(b)
        B0 = B[0]
        B1 = B[1]

        Y0 = np.zeros(97)

        Y0[:4] = A0, A1, B0, B1
        # N_X nastavimo vse na 1 - st celic
        Y0[54:-3] = 1

        T = np.linspace(0, t_end, N)

        t1 = t_end
        dt = t_end/N
        T = np.arange(0,t1+dt,dt)
        Y = np.zeros([1+N,len(Y0)])
        Y[0,:] = Y0

        r = ode(adder_model_ODE).set_integrator('zvode', method='bdf')
        r.set_initial_value(Y0, T[0]).set_f_params(params)

        i = 1
        while r.successful() and r.t < t1:
            Y[i,:] = r.integrate(r.t+dt)
            i += 1

        S0 = 1 if Y[:,-3][-1] > 1 else 0
        S1 = 1 if Y[:,-2][-1] > 1 else 0
        C = 1 if Y[:,-1][-1] > 1 else 0

        s = C * 4 + S1 * 2 + S0

        sums[0].append(s)

        print(f' {A1}{A0}')
        print(f' {B1}{B0}')
        print('--------')
        print(f'{C}{S1}{S0}')
        print(f'{a}+{b}={s}')
        print()

    print()

print(x_axis, y_axis, sums)
fig, ax = plt.subplots()
plt.xlabel('A')
plt.ylabel('B')
im = ax.imshow(sums)
ax.set_xticks(np.arange(len(x_axis)))
ax.set_yticks(np.arange(len(y_axis)))
ax.set_xticklabels(x_axis)
ax.set_yticklabels(y_axis)
plt.setp(ax.get_xticklabels(), rotation=45, ha="right",
         rotation_mode="anchor")
for i in range(len(y_axis)):
    for j in range(len(x_axis)):
        text = ax.text(j, i, sums[i][j], ha="center", va="center", color="w")

ax.set_title("Adder Heat map")
fig.tight_layout()
_, labels = plt.yticks()
plt.show()
# plt.savefig('cblb/figs/heatmap_adder.pdf')

# plt.figure()
# plt.plot(T,out)
# plt.show()

reason: import numpy,from scipy


text:
"""VIMS camera model."""

import numpy as np

from .errors import VIMSCameraError

class VIMSCameraAbstract:
    """Abstract VIMS Camera object.

    Parameters
    ----------
    offsets: [int, int]
        Acquisition X and Z offsets.
    swaths: [int, int]
        Acquisition width and length swaths.

    """

    # Boresite location (X, Y)
    BORESITE = np.array([[[32]], [[32]]])

    # Camera pixel size [rad/pixel]
    PIXEL_SIZE = None

    # Pixel scaling factor in (X, Y) direction
    SCALE = (1, 1)

    def __init__(self, offsets, swaths):
        self.offset_x, self.offset_y = offsets
        self.swath_x, self.swath_y = swaths

        self.__grid = None
        self.__pixels = None

    def __str__(self):
        return self.__class__.__name__

    def __repr__(self):
        return '\n - '.join([
            f'<{self}>',
            f'Swath: {self.swath}',
            f'Offset: {self.offset}',
            f'Scaling: {self.scale}',
            f'Pixel size: {self.PIXEL_SIZE * 1e3} mrad/pix',
            f'Boresite: {tuple(self.BORESITE.flatten())}',
        ])

    @property
    def swath(self):
        """Camera swath (x, y)."""
        return self.swath_x, self.swath_y

    @property
    def offset(self):
        """Camera offset (x, y)."""
        return self.offset_x, self.offset_y

    @property
    def scale_x(self):
        """Scaling factor in X direction."""
        return self.SCALE[0]

    @property
    def scale_y(self):
        """Scaling factor in Y direction."""
        return self.SCALE[1]

    @property
    def scale(self):
        """Camera scale (x, y)."""
        return self.scale_x, self.scale_y

    @staticmethod
    def _positions(offset, swath, scale):
        """Scaled initial pixel position (X, Y).

        Pixel position are expressed in the
        camera focal plane.

        """
        if scale == 1:
            start = offset
        elif scale == 2:
            start = offset + (swath // 2) / 2 - 1 / 4
        elif scale == 3:
            start = offset + swath / 3 - 1 / 3
        else:
            raise VIMSCameraError('Scale value must be 1, 2, or 3.')

        stop = start + (swath - 1) / scale
        return np.linspace(start, stop, swath)

    @property
    def _x(self):
        """Scaled sample position on the sensor."""
        return self._positions(self.offset_x, self.swath_x, self.scale_x)

    @property
    def _y(self):
        """Line position on the sensor."""
        return self._positions(self.offset_y, self.swath_y, self.scale_y)

    @property
    def grid(self):
        """Camera pixel grid (X, Y)."""
        if self.__grid is None:
            self.__grid = np.asarray(np.meshgrid(self._x, self._y))
            self.__pixels = None
        return self.__grid

    @property
    def cgrid(self):
        """Camera grid FOV contour.

        Note
        ----
        The corners pixels are set on circle or an ellipse
        based on the shape of the pixel (:py:attr:`scale`).

        """
        x, y = self._x, self._y
        xl, xr, yt, yb = x[0], x[-1], y[0], y[-1]

        ex, ey = .5 / np.array(self.scale)         # Edges
        dx, dy = ex / np.sqrt(2), ey / np.sqrt(2)  # Ellipse corners

        return np.hstack([
            [[xl - dx], [yt - dy]],               # Top-Left corner
            [x, self.swath_x * [yt - ey]],        # Top edge
            [[xr + dx], [yt - dy]],               # Top-Right corner
            [self.swath_y * [xr + ex], y],        # Right edge
            [[xr + dx], [yb + dy]],               # Bottom-Right corner
            [x[::-1], self.swath_x * [yb + ey]],  # Bottom edge
            [[xl - dx], [yb + dy]],               # Bottom-Left corner
            [self.swath_y * [xl - ex], y[::-1]],  # Left edge
            [[xl - dx], [yt - dy]],               # Top-Left corner
        ])

    @property
    def rgrid(self):
        """Camera grid pixel corners.

        Corners order:
            - Top Left
            - Top Right
            - Bottom Right
            - Bottom Left

        Note
        ----
        Corner are defined as diagonal points
        compare to the pixel center.

        """
        x, y = self._x, self._y

        cx, cy = .5 / np.array(self.scale)  # Rectangle corners

        tl = np.meshgrid(x - cx, y - cy)
        tr = np.meshgrid(x + cx, y - cy)
        br = np.meshgrid(x + cx, y + cy)
        bl = np.meshgrid(x - cx, y + cy)

        return np.moveaxis([tl, tr, br, bl], 0, 3)

    @property
    def fgrid(self):
        """Camera grid pixel footprint.

        Footprint order:
            - Top Left
            - Top
            - Top Right
            - Right
            - Bottom Right
            - Bottom
            - Bottom Left
            - Left
            - Top Left

        Note
        ----
        The corners pixels are set on circle or an ellipse
        based on the shape of the pixel (:py:attr:`scale`).

        """
        x, y = self._x, self._y

        ex, ey = .5 / np.array(self.scale)         # Edges
        dx, dy = ex / np.sqrt(2), ey / np.sqrt(2)  # Ellipse corners

        l = np.meshgrid(x - ex, y)
        r = np.meshgrid(x + ex, y)
        t = np.meshgrid(x, y + ey)
        b = np.meshgrid(x, y - ey)

        tl = np.meshgrid(x - dx, y + dy)
        tr = np.meshgrid(x + dx, y + dy)
        bl = np.meshgrid(x - dx, y - dy)
        br = np.meshgrid(x + dx, y - dy)

        return np.moveaxis([tl, t, tr, r, br, b, bl, l, tl], 0, 3)

    @property
    def extent(self):
        """Camera grid extent."""
        return [self._x[0] - .5 / self.scale_x,
                self._x[-1] + .5 / self.scale_x,
                self._y[-1] + .5 / self.scale_y,
                self._y[0] - .5 / self.scale_y]

    def xy2ang(self, x, y):
        """Convert pixel coordinates in camera look vector.

        Parameters
        ----------
        x: float or np.array
            Sample coordinates on the sensor.
        y: float or np.array
            Line coordinates on the sensor.

        Return
        ------
        (float, float, float)
            XYZ normalized pixel vector in Camera frame.

        """
        shape = (3, *np.shape(x))
        if np.ndim(x) == 1:
            x, y = [x], [y]
        elif np.ndim(x) == 3:
            s = (2, 1, np.size(x) // 2)
            x, y = np.reshape(x, s), np.reshape(y, s)

        phi, theta = (np.array([x, y]) - self.BORESITE) * self.PIXEL_SIZE
        return np.reshape([
            np.cos(theta) * np.sin(phi),
            np.sin(theta),
            np.cos(theta) * np.cos(phi),
        ], shape)

    @property
    def pixels(self):
        """Camera pixels orientation in J2000 frame.

        Return
        ------
        array(3, NL, NS)
            Pixel boresights in Camera frame.
        """
        if self.__pixels is None:
            self.__pixels = self.xy2ang(*self.grid)
        return self.__pixels

    @property
    def cpixels(self):
        """Camera contour pixels orientation in J2000 frame.

        Return
        ------
        array(3, N)
            Pixel contours boresights in Camera frame.
        """
        return self.xy2ang(*self.cgrid)

    @property
    def rpixels(self):
        """Camera corner pixels orientation in J2000 frame.

        Return
        ------
        array(3, NL, NS, 4)
            Pixel corners boresights in Camera frame.

        Note
        ----
        Corner are defined as diagonal points
        compare to the pixel center.

        """
        return self.xy2ang(*self.rgrid)

    @property
    def fpixels(self):
        """Camera footprint pixels orientation in J2000 frame.

        Return
        ------
        array(3, NL, NS, 9)
            Pixel footprints boresights in Camera frame.

        """
        return self.xy2ang(*self.fgrid)

    @property
    def pix_res_x(self):
        """Pixel angular resolution in X direction."""
        return 2 * np.tan(self.PIXEL_SIZE / (2 * self.scale_x))

    @property
    def pix_res_y(self):
        """Pixel angular resolution in Y direction."""
        return 2 * np.tan(self.PIXEL_SIZE / (2 * self.scale_y))

    @property
    def pix_res(self):
        """Mean pixel angular resolution."""
        return np.sqrt(self.pix_res_x * self.pix_res_y)

class VIMSCameraVis(VIMSCameraAbstract):
    """VIMS-VIS camera in ``NORMAL`` sampling mode."""

    PIXEL_SIZE = 0.506e-3

class VIMSCameraIr(VIMSCameraAbstract):
    """VIMS-IR camera in ``NORMAL`` sampling mode."""

    PIXEL_SIZE = 0.495e-3

class VIMSCameraVisHR(VIMSCameraVis):
    """VIMS-VIS camera in ``HI-RES`` sampling mode."""

    SCALE = (3, 3)

class VIMSCameraIrHR(VIMSCameraIr):
    """VIMS-IR camera in ``HI-RES`` sampling mode."""

    SCALE = (2, 1)

class VIMSCamera:
    """VIMS Camera object.

    Parameters
    ----------
    channel: str
        Camera channel (``VIS``|``IR``).
    mode: str
        Acquisition sampling mode (``NORMAL``|``HI-RES``).
    offsets: [int, int]
        Acquisition X and Z offsets.
    swaths: [int, int]
        Acquisition width and length swaths.

    """

    def __new__(cls, channel, mode, offsets, swaths):

        if channel not in ['VIS', 'IR']:
            raise VIMSCameraError(f'Unknown channel `{channel}`. '
                                  'Only `VIS` and `IR` are available')

        if mode not in ['NORMAL', 'HI-RES']:
            raise VIMSCameraError(f'Unknown sampling mode `{mode}`. '
                                  'Only `NORMAL` and `HI-RES` are available')

        if mode == 'NORMAL' and channel == 'IR':
            return VIMSCameraIr(offsets, swaths)

        if mode == 'NORMAL' and channel == 'VIS':
            return VIMSCameraVis(offsets, swaths)

        if mode == 'HI-RES' and channel == 'IR':
            return VIMSCameraIrHR(offsets, swaths)

        if mode == 'HI-RES' and channel == 'VIS':
            return VIMSCameraVisHR(offsets, swaths)

        return None

max_stars_repo_name: seignovert/pyvims

max_issues_repo_name: seignovert/pyvims

max_forks_repo_name: seignovert/pyvims


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

from concurrent.futures import ProcessPoolExecutor
from functools import partial
import numpy as np
import os
import librosa
from multiprocessing import cpu_count
import argparse

def build_from_path(in_dir, out_dir, num_workers=1):
    executor = ProcessPoolExecutor(max_workers=num_workers)
    futures = []
    index = 1
    with open(os.path.join(in_dir, 'metadata.csv'), encoding='utf-8') as f:
        for line in f:
            parts = line.strip().split('|')
            wav_path = os.path.join(in_dir, 'wavs', '%s.wav' % parts[0])
            text = parts[2]
            futures.append(executor.submit(
                partial(_process_utterance, out_dir, index, wav_path, text)))
            index += 1
    return [future.result() for future in futures]

def _process_utterance(out_dir, index, wav_path, text):
    # Load the audio to a numpy array:
    wav, sr = librosa.load(wav_path, sr=22050)

    wav = wav / np.abs(wav).max() * 0.999
    out = wav
    constant_values = 0.0
    out_dtype = np.float32
    n_fft = 1024
    hop_length = 256
    reference = 20.0
    min_db = -100

    # Compute a mel-scale spectrogram from the trimmed wav:
    # (N, D)
    mel_spectrogram = librosa.feature.melspectrogram(wav, sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=80, fmin=125, fmax=7600).T

    # mel_spectrogram = np.round(mel_spectrogram, decimals=2)
    mel_spectrogram = 20 * np.log10(np.maximum(1e-4, mel_spectrogram)) - reference
    mel_spectrogram = np.clip((mel_spectrogram - min_db) / (-min_db), 0, 1)

    pad = (out.shape[0] // hop_length + 1) * hop_length - out.shape[0]
    pad_l = pad // 2
    pad_r = pad // 2 + pad % 2

    # zero pad for quantized signal
    out = np.pad(out, (pad_l, pad_r), mode="constant", constant_values=constant_values)
    N = mel_spectrogram.shape[0]
    assert len(out) >= N * hop_length

    # time resolution adjustment
    # ensure length of raw audio is multiple of hop_size so that we can use
    # transposed convolution to upsample
    out = out[:N * hop_length]
    assert len(out) % hop_length == 0

    timesteps = len(out)

    # Write the spectrograms to disk:
    audio_filename = 'ljspeech-audio-%05d.npy' % index
    mel_filename = 'ljspeech-mel-%05d.npy' % index
    np.save(os.path.join(out_dir, audio_filename),
            out.astype(out_dtype), allow_pickle=False)
    np.save(os.path.join(out_dir, mel_filename),
            mel_spectrogram.astype(np.float32), allow_pickle=False)

    # Return a tuple describing this training example:
    return audio_filename, mel_filename, timesteps, text

def preprocess(in_dir, out_dir, num_workers):
    os.makedirs(out_dir, exist_ok=True)
    metadata = build_from_path(in_dir, out_dir, num_workers)
    write_metadata(metadata, out_dir)

def write_metadata(metadata, out_dir):
    with open(os.path.join(out_dir, 'train.txt'), 'w', encoding='utf-8') as f:
        for m in metadata:
            f.write('|'.join([str(x) for x in m]) + '\n')
    frames = sum([m[2] for m in metadata])
    sr = 22050
    hours = frames / sr / 3600
    print('Wrote %d utterances, %d time steps (%.2f hours)' % (len(metadata), frames, hours))
    print('Max input length:  %d' % max(len(m[3]) for m in metadata))
    print('Max output length: %d' % max(m[2] for m in metadata))

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='Preprocessing',
                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('--in_dir', '-i', type=str, default='./', help='In Directory')
    parser.add_argument('--out_dir', '-o', type=str, default='./', help='Out Directory')
    args = parser.parse_args()

    num_workers = cpu_count()
    preprocess(args.in_dir, args.out_dir, num_workers)

max_stars_repo_name: roberthoenig/VQ-VAE-Speech

max_issues_repo_name: roberthoenig/VQ-VAE-Speech

max_forks_repo_name: roberthoenig/VQ-VAE-Speech


text:
from scipy.io import wavfile
import soundfile as sf
import numpy as np
import matplotlib.pyplot as plt
import webrtcvad
import struct
import argparse

def main(opts):
    rate, wav = wavfile.read(opts.in_file)
    vad = webrtcvad.Vad()
    vad.set_mode(opts.vad_mode)
    win_len = int(opts.win_len * rate + .5)
    # y is to plot (wav time resolution), vads contains one point per frame
    y = []
    vads = []
    frames = []
    for beg_i in np.arange(0, len(wav), win_len):
        end_i = min(beg_i + win_len, len(wav))
        frame = wav[beg_i:end_i]
        frames.append(frame)
        if len(frame) < win_len:
            P = win_len - len(frame)
            frame = np.concatenate((frame,
                                    np.zeros((P,), dtype=frame.dtype)),
                                    axis=0)
        xb = struct.pack("%dh" % len(frame), *frame)
        is_speech = vad.is_speech(xb,
                                  sample_rate = rate)
        y_i = 1 if is_speech else 0
        vads.append(y_i)
        y += [y_i] * win_len

    if opts.show:
        plt.subplot(2, 1, 1)
        plt.plot(wav)
        plt.subplot(2, 1, 2)
        plt.plot(y)
        plt.show()

    if opts.trim_sil > 0:
        assert opts.out_file is not None
        # post-process to find silence regions over
        # this value, and delete those pieces of speech
        count0 = 0
        fcandidates = []
        out_samples = []
        max_sil = int(np.ceil(((opts.trim_sil / 1000) * 16000) / win_len))
        frame_len = len(frames[0])
        for idx, (y_i, frame) in enumerate(zip(vads, frames), start=1):
            if y_i == 0:
                count0 += 1
                fcandidates.extend(frame.tolist())
            if y_i == 1 or idx >= len(frames):
                # change detected, process all previous counts
                # and frame candidates, to know whether to
                # discard all but max allowed or to leave them alone
                if count0 >= max_sil:
                    # discard all candidates but last trim_sil ones
                    K = max_sil * frame_len
                    fcandidates = fcandidates[-K:]
                if len(fcandidates) > 0:
                    out_samples += fcandidates[:]
                    fcandidates = []
                frame = frame.tolist()
                out_samples += frame[:]
                count0 = 0
        out_samples = np.array(out_samples, dtype=np.int16)
        R = 100 - ((len(out_samples) / len(wav)) * 100.)
        if opts.verbose:
            print('Input num samples: ', len(wav))
            print('Output num samples: ', len(out_samples))
            print('Discarded {} % of samples'.format(R))
        if opts.show:
            plt.subplot(3,1,1)
            plt.plot(wav)
            plt.subplot(3,1,2)
            plt.plot(y)
            plt.subplot(3,1,3)
            plt.plot(out_samples)
            plt.show()
        sf.write(opts.out_file, out_samples, rate, 'PCM_16')
        if opts.out_log is not None:
            with open(opts.out_log, 'w') as out_log_f:
                out_log_f.write('isamples\tosamples\treduction[%]\n{}\t{}\t{:.2f}\n'
                                ''.format(len(wav), len(out_samples), R))

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('in_file', type=str, default=None,
                        help='Path to wav file to be processed (Def: None).')
    parser.add_argument('--trim_sil', type=float,
                        default=0,
                        help='Silence regions over this value '
                             'will be trimmed. The value 0 has no effect. '
                             'Specify in milisec, please (Def: 0).')
    parser.add_argument('--verbose', action='store_true', default=False)
    parser.add_argument('--out_file', type=str, default=None)
    parser.add_argument('--out_log', type=str, default=None,
                        help='Will write log in this file if specified')
    parser.add_argument('--show', action='store_true', default=False)
    parser.add_argument('--win_len', type=float, default=0.02,
                        help='In seconds (Def: 0.02).')
    parser.add_argument('--vad_mode', type=int, default=3)

    opts = parser.parse_args()
    main(opts)

max_stars_repo_name: ishine/pase

max_issues_repo_name: ishine/pase

max_forks_repo_name: ishine/pase


text:
import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
import matplotlib.image as mpimg
import os
import random
import scipy.sparse
import sklearn.svm as svm

w_size = 256

def read_image(f):
  img = mpimg.imread(f)
  t = img.shape
  factor = w_size / t[0]
  img = cv.resize(img, (int(factor*t[1]), int(factor*t[0])))
  if len(img.shape) == 3:
    img = cv.cvtColor(img, cv.COLOR_BGR2GRAY)
  return img

def show_image(img, label):
  plt.imshow(img, interpolation='nearest')
  plt.title(label)
  plt.show()

class Dataset_old:
  def __init__(self, repo, drop_rate):
    # all training and testing data stored in data_X and data_y
    self.data_X = []
    self.data_y = []
    self.train = []
    self.test = []
    self.labels = []

    folders = os.listdir(repo)
    labels = [i.split(".")[1]  for i in folders]
    for i in range(len(folders)):
      current_folder = repo + folders[i] + '/'
      image_names = os.listdir(current_folder)
      self.labels.append(folders[i])
      for j in image_names:
        if(random.random() > drop_rate):
          image = read_image(current_folder + j)
          self.data_X.append(image)
          self.data_y.append(i)

  def show_image(self, index = None, num = 5):
    if index == None:
      shown = np.random.choice(len(self.data_X), num)
      for i in shown:
        show_image(self.data_X[i], self.data_y[i])
    else:
      show_image(self.data_X[index], self.data_y[index])

  def get_data_X(self):
    return self.data_X

  def seperate_data(self, num_test = 300):
    self.test = list(np.random.choice(len(self.data_X), num_test, replace = False))
    self.train = list(set(range(len(self.data_X)))-set(self.test))

class Dataset:
  def __init__(self, repo, samples = 70, no_clutter = True, feature_function = None, pyramid = 0, verbose = False):
    # all training and testing data stored in data_X for better flexibility
    self.data_X = {}
    self.train_X = []
    self.train_y = []
    self.validation_X = []
    self.validation_y = []
    self.test_X = []
    self.test_y = []
    self.labels = []

    folders = os.listdir(repo)
    labels = [i.split(".")[1]  for i in folders]

    for i in range(len(folders)):
      current_folder = repo + folders[i] + '/'
      image_names = os.listdir(current_folder)
      self.labels.append(folders[i])
      tmp = np.array([])
      if verbose:
        print(folders[i], 'is loaded.')
      # exclude the clutter class
      if i == 256 and no_clutter:
        continue
      t = 0
      for j in image_names:
        image = read_image(current_folder + j)
        if feature_function == None:
          tmp = np.vstack([tmp, image]) if tmp.size else image
        else:
          try:
            added = spatial_pyramid(feature_function,
                              pyramid, image)
          except ValueError as e:
            continue
          tmp = np.vstack([tmp, added]) if tmp.size else added

          t += 1
          if t >= samples:
            break
      self.data_X[i] = tmp

  def show_image(self, index = None, num = 5):
    if index == None:
      # randomly show some images in the test set
      shown = np.random.choice(len(self.test_X), num)
      for i in shown:
        show_image(self.test_X[i], self.test_y[i])
    else:
      show_image(self.test_X[index], self.test_y[index])

  def generate_train_test_samples(self, train_number = 60):
    for i, j in self.data_X.items():
      train_choices = np.random.choice(j.shape[0], train_number, replace = False)
      self.train_X.append(j[train_choices])
      self.train_y.append(np.ones([train_number]) * i)

      validation_choices = list(set(range(j.shape[0])) - set(train_choices))[:5]
      self.validation_X.append(j[validation_choices])
      self.validation_y.append(np.ones([len(validation_choices)]) * i)

      not_choices = list(set(range(j.shape[0])) - set(train_choices) - \
                         set(validation_choices))[:5]
      self.test_X.append(j[not_choices])
      self.test_y.append(np.ones([len(not_choices)]) * i)

  def get_train_X(self):
    return self.train_X

  def train(self, model, scaler, training_number = 15, accuracy = True):
    X_t = self.train_X[0][:training_number, :]
    y_t = self.train_y[0][:training_number]
    for i in range(1, len(self.train_X)):
      X_t = np.vstack((X_t, self.train_X[i][:training_number, :]))
      y_t = np.concatenate((y_t, self.train_y[i][:training_number]))
    if scaler:
      scaler.fit(X_t)
      X_t = scaler.transform(X_t)

    model.fit(X_t, y_t)
    if accuracy:
      return sum(model.predict(X_t) == y_t)/X_t.shape[0]

  def test(self, model, scaler, validation = True):
    if validation:
      X_tmp = self.validation_X
      y_tmp = self.validation_y
    else:
      X_tmp = self.test_X
      y_tmp = self.test_y

    X_t = X_tmp[0][:5, :]
    y_t = y_tmp[0][:5]
    for i in range(1, len(X_tmp)):
      X_t = np.vstack((X_t, X_tmp[i][:5, :]))
      y_t = np.concatenate((y_t, y_tmp[i][:5]))

    if scaler:
      X_t = scaler.transform(X_t)
    return sum(model.predict(X_t) == y_t)/X_t.shape[0]

def spatial_pyramid(feature_function, layer, image):
  ret = feature_function(image)
  for k in range(1, layer):
    t = 2 ** k
    xx = image.shape[0] // t
    yy = image.shape[1] // t
    for i in range(t):
      for j in range(t):
        ret = np.hstack((ret, feature_function(image[i*xx:(i+1)*xx, j*yy:(j+1)*yy])))
  return ret

class HistIntersectionModel:
  def __init__(self, C = 1, shrinking=True, probability=False,
                 tol=1e-3, cache_size=200, class_weight=None,
                 verbose=False, max_iter=-1, decision_function_shape='ovr',
                 break_ties=False,
                 random_state=None, **params):

    self.model = svm.SVC(kernel = 'precomputed', C = C , shrinking = shrinking,
                         probability = probability, tol = tol , cache_size = \
                         cache_size,  class_weight = class_weight,
                         verbose = verbose,  max_iter = max_iter,
                         decision_function_shape = decision_function_shape,
                         break_ties = break_ties,
                         random_state = random_state)
    self.intersected = None

  def fit(self, X, y):
    self.intersected = X
    K = hist_intersection(X, X)
    self.model.fit(K, y)

  def predict(self, X):
    K = hist_intersection(X, self.intersected)
    return self.model.predict(K)

  def get_params(self, deep=True):
    return self.model.get_params(deep)

  def decision_function(self, X):
   K = hist_intersection(X, self.intersected)
   return self.model.decision_function(K)

  def predict_proba(self, X):
    K = hist_intersection(X, self.intersected)
    return self.model.predict_proba(K)

def hist_intersection(X, Y):
  kernel = np.zeros((X.shape[0], Y.shape[0]))

  for i in range(Y.shape[1]):
      c1 = Y[:, i].reshape(-1, 1)
      c2 = Y[:, i].reshape(-1, 1)
      kernel += np.minimum(c1, c2.T)
  return kernel

reason: import numpy,import scipy


text:
import math
from time import time
import numpy as np
import argparse

from time_limit import set_time_limit

parser = argparse.ArgumentParser()
parser.add_argument('-n', type=int, default=1000000)
args = parser.parse_args()

def float_operations(n):
    start = time()
    for i in range(0, n):
        sin_i = math.sin(i)
        cos_i = math.cos(i)
        sqrt_i = math.sqrt(i)
    latency = time() - start
    return latency

@set_time_limit()
def lambda_handler(event, context):
    n = int(event['n'])
    result = float_operations(n)
    # print(result)
    return result

if __name__ == '__main__':
    event = dict()
    event['n'] = args.n

    # print()
    # print("#### test: float_operation ####")
    total = list()
    for i in range(1):
        total.append(lambda_handler(event=event, context=None))
    # print("mean: " + str(np.mean(total)))
    # print("std:  " + str(np.std(total)))


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan  7 10:48:48 2021

@author: rayin
"""

import os, sys
import pandas as pd
from pprint import pprint
import matplotlib.pyplot as plt
import numpy as np
import math
import random
import torch
from collections import Counter
from scipy import interp
from torch import nn

from sklearn import metrics
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import cross_val_predict
from sklearn.model_selection import train_test_split
from imblearn.metrics import geometric_mean_score
from sklearn.metrics import matthews_corrcoef
from sklearn.metrics import roc_curve
from sklearn.metrics import auc
from sklearn.metrics import accuracy_score
from sklearn.metrics import precision_score
from sklearn.metrics import recall_score
from sklearn.metrics import f1_score
from sklearn.metrics import roc_auc_score
from sklearn.metrics import precision_recall_curve
from sklearn.metrics import plot_precision_recall_curve
from sklearn.metrics import average_precision_score
from sklearn.metrics import RocCurveDisplay
from sklearn.metrics import PrecisionRecallDisplay

from sklearn import ensemble
from sklearn.linear_model import LogisticRegression
from sklearn import neighbors

import xgboost as xgb
from xgboost import plot_importance
from matplotlib import pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import StratifiedKFold
from sklearn.model_selection import KFold
from sklearn.impute import SimpleImputer

os.chdir("/Users/rayin/Google Drive/Harvard/5_data/UDN/work")

def is_input_contain_nan(df):
    # for i in range(0, df.shape[0]):
    #     for j in range(0, df.shape[1]):
    #         if str(df.iloc[i][j]).lower() == 'nan':
    #             print('row'%i + 'contains nan at'%j + '\n')
    indices_to_keep = ~df.isin([np.nan, np.inf, -np.inf]).any(1)
    return indices_to_keep

def logistic_regression_cv(data_x, data_y):

    train_x, test_x, train_y, test_y = train_test_split(data_x, data_y, test_size=0.2, random_state=42)
    clf = LogisticRegression()
    kfold = StratifiedKFold(n_splits=5) #shuffle=True

    #training
    train_acc = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='accuracy')
    train_pre = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='precision')
    train_rec = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='recall')
    train_f1 = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='f1')
    train_pred = cross_val_predict(clf, train_x, train_y, cv=kfold)
    train_mcc = matthews_corrcoef(train_y, train_pred)

    #testing
    clf = LogisticRegression().fit(train_x, train_y)
    test_acc = accuracy_score(test_y, clf.predict(test_x))
    test_pre = precision_score(test_y, clf.predict(test_x))
    test_rec = recall_score(test_y, clf.predict(test_x))
    test_f1 = f1_score(test_y, clf.predict(test_x))
    test_mcc = matthews_corrcoef(test_y, clf.predict(test_x))

    return train_acc, train_pre, train_rec, train_f1, train_mcc, test_acc, test_pre, test_rec, test_f1, test_mcc

def randomforest_cv(data_x, data_y):
    #np.random.seed(100)
    train_x, test_x, train_y, test_y = train_test_split(data_x, data_y, test_size=0.1, random_state=42)
    #print(test_y.index.values)
    #cross-validation
    #clf = ensemble.RandomForestClassifier(n_estimators=140)
    clf = ensemble.RandomForestClassifier(n_estimators=340, max_depth=5, max_features=7, min_samples_split=30,
                                  min_samples_leaf=10, random_state=42)
    kfold = StratifiedKFold(n_splits=5) #shuffle=True
    #kfold = KFold(n_splits=5, shuffle=True)
    train_acc = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='accuracy')
    train_pre = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='precision')
    train_rec = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='recall')
    train_f1 = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='f1')
    train_pred = cross_val_predict(clf, train_x, train_y, cv=kfold)
    train_mcc = matthews_corrcoef(train_y, train_pred)
    train_auc = roc_auc_score(train_y, train_pred)

    # print('training set:')
    # print("train_acc: %f" %train_acc.mean() + '\n')
    # print("train_pre: %f" %train_pre.mean() + '\n')
    # print("train_rec: %f" %train_rec.mean() + '\n')
    # print("train_f1 : %f " %train_f1.mean() + '\n')
    # print("train_mcc: %f" %np.mean(train_mcc) + '\n')

    #testing
    clf = ensemble.RandomForestClassifier().fit(train_x, train_y)
    test_acc = accuracy_score(test_y, clf.predict(test_x))
    test_pre = precision_score(test_y, clf.predict(test_x))
    test_rec = recall_score(test_y, clf.predict(test_x))
    test_f1 = f1_score(test_y, clf.predict(test_x))
    test_mcc = matthews_corrcoef(test_y, clf.predict(test_x))
    test_auc = roc_auc_score(test_y, clf.predict(test_x))

    # print('testing set:')
    # print("test_acc: %f" %test_acc.mean() + '\n')
    # print("test_pre: %f" %test_pre.mean() + '\n')
    # print("test_rec: %f" %test_rec.mean() + '\n')
    # print("test_f1 : %f " %test_f1.mean() + '\n')
    # print("test_mcc: %f" %np.mean(test_mcc) + '\n')

    return train_acc, train_pre, train_rec, train_f1, train_mcc, train_auc, test_acc, test_pre, test_rec, test_f1, test_mcc, test_auc

#randomforest_cv(feature, label)

def xgboost_cv(data_x, data_y):
    np.random.seed(100)
    train_x, test_x, train_y, test_y = train_test_split(data_x, data_y, test_size=0.1, random_state=42)
    #cross-validation
    clf = xgb.XGBClassifier(max_depth=8,
                        learning_rate=0.05,
                        n_estimators=20,
                        objective='binary:logistic',
                        eval_metric = "logloss",
                        nthread=4,
                        gamma=0,
                        min_child_weight=5,
                        max_delta_step=0,
                        subsample=0.8,
                        colsample_bytree=0.7,
                        colsample_bylevel=1,
                        reg_alpha=1e-5,
                        reg_lambda=1,
                        scale_pos_weight=1,
                        seed=27,
                        missing=None)

    #clf = xgb.XGBClassifier(n_estimators=30)
    #kfold = KFold(n_splits=10, shuffle=True)
    kfold = StratifiedKFold(n_splits=5)

    #training
    train_acc = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='accuracy')
    train_pre = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='precision')
    train_rec = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='recall')
    train_f1 = cross_val_score(clf, train_x, train_y, cv=kfold, scoring='f1')
    train_pred = cross_val_predict(clf, train_x, train_y, cv=kfold)
    train_mcc = matthews_corrcoef(train_y, train_pred)

    #testing
    clf = xgb.XGBClassifier().fit(train_x, train_y)
    test_acc = accuracy_score(test_y, clf.predict(test_x))
    test_pre = precision_score(test_y, clf.predict(test_x))
    test_rec = recall_score(test_y, clf.predict(test_x))
    test_f1 = f1_score(test_y, clf.predict(test_x))
    test_mcc = matthews_corrcoef(test_y, clf.predict(test_x))

    # print('training set:')
    # print("train_acc: %f" %train_acc.mean() + '\n')
    # print("train_pre: %f" %train_pre.mean() + '\n')
    # print("train_rec: %f" %train_rec.mean() + '\n')
    # print("train_f1 : %f " %train_f1.mean() + '\n')
    # print("train_mcc: %f" %np.mean(train_mcc) + '\n')

    # print('testing set:')
    # print("test_acc: %f" %test_acc.mean() + '\n')
    # print("test_pre: %f" %test_pre.mean() + '\n')
    # print("test_rec: %f" %test_rec.mean() + '\n')
    # print("test_f1 : %f " %test_f1.mean() + '\n')
    # print("test_mcc: %f" %np.mean(test_mcc) + '\n')

    return train_acc, train_pre, train_rec, train_f1, train_mcc, test_acc, test_pre, test_rec, test_f1, test_mcc
#xgboost_cv(feature_imputation, label)

def patient_variable_influence(model, data_x, data_y, patient_info, variable):
    #cross-validation on raw features + patient_level variables (after one-hot coding)
    if variable == 'raw':
        feature_new = data_x
    elif variable == 'age_symptom_onset':
        patient_variable = patient_info.iloc[:, 0:1]
        feature_new = pd.concat([data_x, patient_variable], axis=1)
    elif variable == 'current_age':
        patient_variable = patient_info.iloc[:, 1:2]
        feature_new = pd.concat([data_x, patient_variable], axis=1)
    elif variable == 'ethnicity':
        patient_variable = patient_info.iloc[:, 2:6]
        feature_new = pd.concat([data_x, patient_variable], axis=1)
    elif variable == 'gender':
        patient_variable = patient_info.iloc[:, 6:8]
        feature_new = pd.concat([data_x, patient_variable], axis=1)
    elif variable == 'race':
        patient_variable = patient_info.iloc[:, 8:14]
        feature_new = pd.concat([data_x, patient_variable], axis=1)
    elif variable == 'all':  #all features
        feature_new = pd.concat([data_x, patient_info], axis=1)

    #np.random.seed(100)

    if model == 'rf':
        train_acc, train_pre, train_rec, train_f1, train_mcc, train_auc, test_acc, test_pre, test_rec, test_f1, test_mcc, test_auc = randomforest_cv(feature_new, data_y)

    elif model == 'xgboost':
        train_acc, train_pre, train_rec, train_f1, train_mcc, train_auc, test_acc, test_pre, test_rec, test_f1, test_mcc, test_auc = xgboost_cv(feature_new, data_y)

    print(model + ' + ' + variable + ':')
    #print("test_acc: %f" %test_acc.mean() + '\n')
    #print("test_pre: %f" %test_pre.mean() + '\n')
    #print("test_rec: %f" %test_rec.mean() + '\n')
    #print("test_f1 : %f " %test_f1.mean() + '\n')
    #print("test_mcc: %f" %np.mean(test_mcc) + '\n')
    print("test_auc : %f " %test_auc.mean() + '\n')

def multi_experiment(model, time, train_x, train_y):
    train_acc_all = []
    train_pre_all = []
    train_rec_all = []
    train_f1_all = []
    train_mcc_all = []
    train_auc_all = []

    test_acc_all = []
    test_pre_all = []
    test_rec_all = []
    test_f1_all = []
    test_mcc_all = []
    test_auc_all = []

    for i in range(0, time):
        if model == 'rf':
            if sum(train_x.isnull().sum().values) != 0:
                imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
                x_new = pd.DataFrame(imp.fit_transform(train_x))
                train_x = x_new
            train_acc, train_pre, train_rec, train_f1, train_mcc, train_auc, test_acc, test_pre, test_rec, test_f1, test_mcc, test_auc = randomforest_cv(train_x, train_y)
        elif model == 'lr':
            if sum(train_x.isnull().sum().values) != 0:
                imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
                x_new = pd.DataFrame(imp.fit_transform(train_x))
                train_x = x_new
            train_acc, train_pre, train_rec, train_f1, train_mcc, test_acc, test_pre, test_rec, test_f1, test_mcc = logistic_regression_cv(train_x, train_y)
        elif model == 'xgboost':
            train_acc, train_pre, train_rec, train_f1, train_mcc, test_acc, test_pre, test_rec, test_f1, test_mcc = xgboost_cv(train_x, train_y)

        train_acc_all.append(train_acc)
        train_pre_all.append(train_pre)
        train_rec_all.append(train_rec)
        train_f1_all.append(train_f1)
        train_mcc_all.append(train_mcc)
        train_auc_all.append(train_auc)

        test_acc_all.append(test_acc)
        test_pre_all.append(test_pre)
        test_rec_all.append(test_rec)
        test_f1_all.append(test_f1)
        test_mcc_all.append(test_mcc)
        test_auc_all.append(test_auc)

    train_acc_mean = np.mean(train_acc_all)
    train_acc_sd = np.std(train_acc_all, ddof = 1)
    train_pre_mean = np.mean(train_pre_all)
    train_pre_sd = np.std(train_pre_all, ddof = 1)
    train_rec_mean = np.mean(train_rec_all)
    train_rec_sd = np.std(train_rec_all, ddof = 1)
    train_f1_mean = np.mean(train_f1_all)
    train_f1_sd = np.std(train_f1_all, ddof = 1)
    train_mcc_mean = np.mean(train_mcc_all)
    train_mcc_sd = np.std(train_mcc_all, ddof = 1)
    train_auc_mean = np.mean(train_auc_all)
    train_auc_sd = np.std(train_auc_all, ddof = 1)

    test_acc_mean = np.mean(test_acc_all)
    test_acc_sd = np.std(test_acc_all, ddof = 1)
    test_pre_mean = np.mean(test_pre_all)
    test_pre_sd = np.std(test_pre_all, ddof = 1)
    test_rec_mean = np.mean(test_rec_all)
    test_rec_sd = np.std(test_rec_all, ddof = 1)
    test_f1_mean = np.mean(test_f1_all)
    test_f1_sd = np.std(test_f1_all, ddof = 1)
    test_mcc_mean = np.mean(test_mcc_all)
    test_mcc_sd = np.std(test_mcc_all, ddof = 1)
    test_auc_mean = np.mean(test_auc_all)
    test_auc_sd = np.std(test_auc_all, ddof = 1)

    print(model + '_train')
    print("accuracy:%f" %train_acc_mean + ' ' + "sd:%f"%train_acc_sd + '\n')
    print("precision:%f" %train_pre_mean + ' ' + "sd:%f"%train_pre_sd + '\n')
    print("recall:%f" %train_rec_mean + ' ' + "sd:%f"%train_rec_sd + '\n')
    print("f1:%f" %train_f1_mean + ' ' + "sd:%f"%train_f1_sd + '\n')
    print("auc:%f" %train_auc_mean + ' ' + "sd:%f"%train_auc_sd + '\n')
    #print("mcc:%f" %train_mcc_mean + ' ' + "sd:%f"%train_mcc_sd + '\n')

    print(model + '_test')
    print("accuracy:%f" %test_acc_mean + ' ' + "sd:%f"%test_acc_sd + '\n')
    print("precision:%f" %test_pre_mean + ' ' + "sd:%f"%test_pre_sd + '\n')
    print("recall:%f" %test_rec_mean + ' ' + "sd:%f"%test_rec_sd + '\n')
    print("f1:%f" %test_f1_mean + ' ' + "sd:%f"%test_f1_sd + '\n')
    print("auc:%f" %test_auc_mean + ' ' + "sd:%f"%test_auc_sd + '\n')
    #print("mcc:%f" %test_mcc_mean + ' ' + "sd:%f"%test_mcc_sd + '\n')

#cross-validation train + independent test data to draw roc curve
def draw_roc_train_test_split(data_x, data_y, syn_data, model):
    np.random.seed(100)
    if model == 'lr':
        if sum(data_x.isnull().sum().values) != 0:
            imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
            data_x = pd.DataFrame(imp.fit_transform(data_x))
        clf = LogisticRegression()
    elif model == 'rf':
        if sum(data_x.isnull().sum().values) != 0:
            imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
            data_x = pd.DataFrame(imp.fit_transform(data_x))
        clf = ensemble.RandomForestClassifier(n_estimators=340, max_depth=5, max_features=7, min_samples_split=30,
                                  min_samples_leaf=10, random_state=42)
    elif model == 'xgboost':
        clf = xgb.XGBClassifier(max_depth=8, learning_rate=0.05, n_estimators=20, objective='binary:logistic', nthread=4, gamma=0, min_child_weight=5,
                                max_delta_step=0, subsample=0.8, colsample_bytree=0.7, colsample_bylevel=1, reg_alpha=1e-5, reg_lambda=1,
                                scale_pos_weight=1, seed=27, missing=None)
    else:
        print('error')
    train_x, test_x, train_y, test_y = train_test_split(data_x, data_y, test_size=0.2)

    #plot the roc curve
    probas_ = clf.fit(train_x, train_y).predict_proba(test_x)
    # Compute ROC curve and area the curve
    fpr, tpr, thresholds = roc_curve(test_y, probas_[:, 1])
    roc_auc = auc(fpr, tpr)
    plt.plot([0, 1], [0, 1], linestyle='--', lw=1, color='grey', alpha=.8)
    plt.plot(fpr, tpr, label= 'Real set'  + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=1, alpha=.8)

    #plot the roc curve for synthetic test set
    syn_test_index = syn_data.index
    syn_data_split = split_list(syn_test_index, group_num=10, shuffle=True, retain_left=False)

    for key, value in syn_data_split.items():
        #print(key)
        syn_test = syn_data.iloc[value]
        syn_test_x = syn_test[syn_test.columns[0:-1]]
        syn_test_y = syn_test[syn_test.columns[-1]]

        probas_ = clf.fit(train_x, train_y).predict_proba(syn_test_x)
        fpr, tpr, thresholds = roc_curve(syn_test_y, probas_[:, 1])
        roc_auc = auc(fpr, tpr)
        plt.plot([0, 1], [0, 1], linestyle='--', lw=1, alpha=.8)
        plt.plot(fpr, tpr, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=1, alpha=.8)

    plt.xlim([-0.00, 1.00])
    plt.ylim([-0.00, 1.00])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC curve')
    plt.legend(loc="lower right")
    #plt.show()

#cross validation to draw roc curve
def draw_roc_curve(train_x, train_y, syn_data, model):
    np.random.seed(100)
    if model == 'lr':
        if sum(train_x.isnull().sum().values) != 0:
            imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
            train_x = pd.DataFrame(imp.fit_transform(train_x))
        clf = LogisticRegression()
    elif model == 'rf':
        if sum(train_x.isnull().sum().values) != 0:
            imp = SimpleImputer(missing_values=np.nan, strategy="constant", fill_value=0)
            train_x = pd.DataFrame(imp.fit_transform(train_x))
        clf = ensemble.RandomForestClassifier()
    elif model == 'xgboost':
        clf = xgb.XGBClassifier()
    else:
        print('error')

    #plot the roc curve
    tprs = []
    aucs = []
    mean_fpr = np.linspace(0, 1, 100)

    i = 0
    kfold = StratifiedKFold(n_splits=10, shuffle=True)
    for train, test in kfold.split(train_x, train_y):
        probas_ = clf.fit(train_x.iloc[train], train_y.iloc[train]).predict_proba(train_x.iloc[test])
        # Compute ROC curve and area the curve
        fpr, tpr, thresholds = roc_curve(train_y[test], probas_[:, 1])
        tprs.append(interp(mean_fpr, fpr, tpr))
        tprs[-1][0] = 0.0
        roc_auc = auc(fpr, tpr)
        aucs.append(roc_auc)
        #plt.plot(fpr, tpr, lw=1, alpha=0.3, label='ROC fold %d (AUC = %0.2f)' % (i, roc_auc))
        i = i + 1
    plt.plot([0, 1], [0, 1], linestyle='--', lw=1, color='grey', alpha=.8)

    mean_tpr = np.mean(tprs, axis=0)
    mean_tpr[-1] = 1.0
    mean_auc = auc(mean_fpr, mean_tpr)
    std_auc = np.std(aucs)
    plt.plot(mean_fpr, mean_tpr, label= 'Real set' + ' ' + r'(AUC = %0.3f $\pm$ %0.3f)' % (mean_auc, std_auc), lw=1, alpha=.8)

    #std_tpr = np.std(tprs, axis=0)
    #tprs_upper = np.minimum(mean_tpr + std_tpr, 1)
    #tprs_lower = np.maximum(mean_tpr - std_tpr, 0)
    #plt.fill_between(mean_fpr, tprs_lower, tprs_upper, color='grey', alpha=.2, label=r'$\pm$ 1 std. dev.')

    #plot the roc curve for synthetic test set
    syn_test_index = syn_data.index
    syn_data_split = split_list(syn_test_index, group_num=5, shuffle=True, retain_left=True)

    for key, value in syn_data_split.items():
        #print(key)
        syn_test = syn_data.iloc[value]
        syn_test_x = syn_test[syn_test.columns[0:-1]]
        syn_test_y = syn_test[syn_test.columns[-1]]

        probas_ = clf.fit(train_x, train_y).predict_proba(syn_test_x)
        fpr, tpr, thresholds = roc_curve(syn_test_y, probas_[:, 1])
        roc_auc = auc(fpr, tpr)
        plt.plot([0, 1], [0, 1], linestyle='--', lw=1, alpha=.8)
        plt.plot(fpr, tpr, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=1, alpha=.8)

    plt.xlim([-0.00, 1.00])
    plt.ylim([-0.00, 1.00])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC curve')
    plt.legend(loc="lower right")
    #plt.show()

def subset(alist, idxs):
    '''
        fetch the subset of alist based on idxs
        alist: list
        idxs: list
    '''
    sub_list = []
    for idx in idxs:
        sub_list.append(alist[idx])

    return sub_list

def split_list(alist, group_num=5, shuffle=True, retain_left=False):
    '''
        alist divide into sublist by group that for each subgroup we have len(alist)//group elements
        shuffle: if random, True by default
        retain_left: the remaining elements will be taken as one group
    '''

    index = list(range(len(alist)))

    if shuffle == True:
        random.shuffle(index)

    elem_num = len(alist) // group_num
    sub_lists = {}

    for idx in range(group_num):
        start, end = idx*elem_num, (idx+1)*elem_num
        sub_lists['set'+str(idx+1)] = subset(alist, index[start:end])

    if retain_left and group_num * elem_num != len(index):
        sub_lists['set'+str(idx+2)] = subset(alist, index[end:])

    return sub_lists

def draw_pr_roc(real_x, real_y, syn_data):
    clf = ensemble.RandomForestClassifier(n_estimators=340, max_depth=5, max_features=7, min_samples_split=30,
                                  min_samples_leaf=10, random_state=42)

    train_x, test_x, train_y, test_y = train_test_split(real_x, real_y, test_size=0.1, random_state=0)

    #plot PR curve for synthetic test
    syn_test_index = syn_data.index
    syn_data_split = split_list(syn_test_index, group_num=5, shuffle=True, retain_left=False)

    #plot the roc curve
    probas_real = clf.fit(train_x, train_y).predict_proba(test_x)
    # Compute ROC curve and area the curve
    fpr, tpr, thresholds = roc_curve(test_y, probas_real[:, 1])
    roc_auc = auc(fpr, tpr)
    #print(syn_data_split)
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(22, 11))
    i = -1
    color = ['r','b', 'y', 'g', 'c']
    for key, value in syn_data_split.items():
        #print(key)
        i = i + 1
        syn_test = syn_data.iloc[value]
        syn_test_x = syn_test[syn_test.columns[0:-1]]
        syn_test_y = syn_test[syn_test.columns[-1]]

        #calculate roc
        probas_syn = clf.fit(train_x, train_y).predict_proba(syn_test_x)
        fpr, tpr, thresholds = roc_curve(syn_test_y, probas_syn[:, 1])
        roc_auc = auc(fpr, tpr)
        ax1.plot([0, 1], [0, 1], linestyle='--', lw=1, alpha=.8)
        #roc_display = RocCurveDisplay(fpr=fpr, tpr=tpr).plot(fpr, tpr, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=2, alpha=.8)

        ax1.plot(fpr, tpr, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (roc_auc), color=color[i], lw=2, alpha=.8)

        #calculate the precision and recall
        probas_syn = clf.fit(train_x, train_y).predict_proba(syn_test_x)[:, 1]
        precision, recall, thresholds = precision_recall_curve(syn_test_y, probas_syn)
        #calculate the precision score
        pr_auc = auc(recall, precision)

        #plt.plot([1, 0], [1, 0], linestyle='--', lw=1, alpha=.8)
        #pr_display = PrecisionRecallDisplay(precision=precision, recall=recall).plot(recall, precision, label= 'Syn' + ' ' + key + ' ' + r'(PRAUC = %0.3f)' % (pr_auc), lw=2, alpha=.8)

        ax2.plot(recall, precision, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (pr_auc), color=color[i], lw=2, alpha=.8)

    plt.xlim([-0.00, 1.05])
    plt.ylim([-0.00, 1.05])
    plt.xticks(size = 10)
    plt.yticks(size = 10)
    ax1.title.set_text('ROC Curve')
    ax1.set_xlabel('False Positive Rate')
    ax1.set_ylabel('True Positive Rate')
    ax1.legend(loc="lower right", prop = {'size': 20})
    ax1.title.set_size(25)
    ax1.xaxis.label.set_size(20)
    ax1.yaxis.label.set_size(20)

    plt.xlim([-0.00, 1.05])
    plt.ylim([-0.00, 1.05])
    plt.xticks(size = 10)
    plt.yticks(size = 10)
    ax2.title.set_text('PR Curve')
    ax2.set_xlabel('Recall')
    ax2.set_ylabel('Precision')
    ax2.legend(loc='lower right', prop = {'size': 20})
    ax2.title.set_size(25)
    ax2.xaxis.label.set_size(20)
    ax2.yaxis.label.set_size(20)
    plt.savefig('result/figure/roc_pr.eps', dpi=300)

#cross-validation train + independent test data + external validation data to draw roc curve
def draw_roc_syn_test(real_x, real_y, syn_data, model):
    #np.random.seed(100)
    if model == 'lr':
        clf = LogisticRegression()
    elif model == 'rf':
        clf = ensemble.RandomForestClassifier(n_estimators=340, max_depth=5, max_features=7, min_samples_split=30,
                                  min_samples_leaf=10, random_state=42)
    elif model == 'xgboost':
        clf = xgb.XGBClassifier(max_depth=8, learning_rate=0.05, n_estimators=20, objective='binary:logistic', nthread=4, gamma=0, min_child_weight=5,
                                max_delta_step=0, subsample=0.8, colsample_bytree=0.7, colsample_bylevel=1, reg_alpha=1e-5, reg_lambda=1,
                                scale_pos_weight=1, seed=27, missing=None)
    else:
        print('error')

    train_x, test_x, train_y, test_y = train_test_split(real_x, real_y, test_size=0.1, random_state=0)

    #plot the roc curve
    probas_real = clf.fit(train_x, train_y).predict_proba(test_x)
    # Compute ROC curve and area the curve
    fpr, tpr, thresholds = roc_curve(test_y, probas_real[:, 1])
    roc_auc = auc(fpr, tpr)
    plt.plot([0, 1], [0, 1], linestyle='--', lw=2, color='grey', alpha=.8)
    #plt.plot(fpr, tpr, label= 'Real set' + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=2, alpha=.8)

    #plot the roc curve for synthetic test set
    syn_test_index = syn_data.index
    syn_data_split = split_list(syn_test_index, group_num=5, shuffle=True, retain_left=False)
    #print(syn_data_split)
    for key, value in syn_data_split.items():
        #print(key)
        syn_test = syn_data.iloc[value]
        syn_test_x = syn_test[syn_test.columns[0:-1]]
        syn_test_y = syn_test[syn_test.columns[-1]]

        probas_syn = clf.fit(train_x, train_y).predict_proba(syn_test_x)
        # prediction_syn = clf.fit(train_x, train_y).predict(syn_test_x)
        # print(key)
        # print(accuracy_score(syn_test_y, prediction_syn))
        # print(precision_score(syn_test_y, prediction_syn))
        # print(recall_score(syn_test_y, prediction_syn))
        fpr, tpr, thresholds = roc_curve(syn_test_y, probas_syn[:, 1])
        roc_auc = auc(fpr, tpr)
        plt.plot([0, 1], [0, 1], linestyle='--', lw=1, alpha=.8)
        plt.plot(fpr, tpr, label= 'Syn' + ' ' + key + ' ' + r'(AUC = %0.3f)' % (roc_auc), lw=2, alpha=.8)

    plt.xlim([-0.00, 1.05])
    plt.ylim([-0.00, 1.05])
    plt.xlabel('False Positive Rate')
    plt.ylabel('True Positive Rate')
    plt.title('ROC curve')
    plt.legend(loc="lower right")
    plt.savefig('result/figure/roc.eps', dpi=300)

def draw_pr_curve(real_x, real_y, syn_data):
    clf = ensemble.RandomForestClassifier(n_estimators=340, max_depth=5, max_features=7, min_samples_split=30,
                                  min_samples_leaf=10, random_state=42)

    train_x, test_x, train_y, test_y = train_test_split(real_x, real_y, test_size=0.1, random_state=0)

    #plot PR curve for synthetic test
    syn_test_index = syn_data.index
    syn_data_split = split_list(syn_test_index, group_num=5, shuffle=True, retain_left=False)

    plt.add_shape(type='line', line=dict(dash='dash'), x0=0, x1=1, y0=1, y1=0)
    #print(syn_data_split)
    for key, value in syn_data_split.items():
        #print(key)
        syn_test = syn_data.iloc[value]
        syn_test_x = syn_test[syn_test.columns[0:-1]]
        syn_test_y = syn_test[syn_test.columns[-1]]

        #calculate the precision and recall
        probas_syn = clf.fit(train_x, train_y).predict_proba(syn_test_x)[:, 1]
        precision, recall, thresholds = precision_recall_curve(syn_test_y, probas_syn)
        #calculate the precision score
        pr_auc = auc(recall, precision)

        #plt.plot([1, 0], [1, 0], linestyle='--', lw=1, alpha=.8)
        plt.plot(recall, precision, label= 'Syn' + ' ' + key + ' ' + r'(PRAUC = %0.3f)' % (pr_auc), lw=2, alpha=.8)

    plt.xlim([-0.00, 1.05])
    plt.ylim([-0.00, 1.05])
    plt.title('PR Curve')
    plt.xlabel('Recall')
    plt.ylabel('Precision')
    plt.legend(loc='lower left')
    plt.savefig('result/figure/pr.eps', dpi=300)

def merge_data(real_feature, real_label, syn_feature, syn_label):
    feature_combine = pd.concat([real_feature, syn_feature], axis=0)
    label_combine = pd.concat([real_label, syn_label], axis=0)
    return feature_combine, label_combine

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

max_stars_repo_name: hms-dbmi/VarPPUD

max_issues_repo_name: hms-dbmi/VarPPUD

max_forks_repo_name: hms-dbmi/VarPPUD


text:
import os
import pandas as pd
import torch

from datetime import datetime
from scipy.stats.mstats import gmean
from dateutil.relativedelta import relativedelta

from server.models.portfolio.stats import *
from server.models.portfolio.cost import costs
from server.models.portfolio.tiingo import get_data
from server.models.portfolio.bl import bl, get_mkt_cap
from server.models.portfolio.rs import fama_french, regime_switch, current_regime, business_days, expected_returns, \
    covariance
from server.models.stock.config import SYMBOLS, COLLECTION, START_DATE, END_DATE
from server.models.portfolio.sent_analysis import SentimentAnalysis, predict

def prepare():
    #tickers = list(pd.read_csv(os.getcwd() + r'/data/tickers.csv')['Tickers'])
    tickers = SYMBOLS
    end_date = datetime.now().strftime("%Y-%m-%d")
    start_date = (datetime.strptime(end_date, "%Y-%m-%d") - relativedelta(years=6)).strftime("%Y-%m-%d")

    # target rebalance date ... based on calendar days ... need to adjust for trading days
    rebalance_date = (datetime.strptime(end_date, "%Y-%m-%d") + relativedelta(months=6, days=1)).strftime("%Y-%m-%d")
    rebalance_date = datetime.strftime(pd.bdate_range(end_date, rebalance_date)[-1], "%Y-%m-%d")

    np.set_printoptions(precision=4)

    ## *********************************************************************************************************************
    #  pull the factors and asset prices
    ## *********************************************************************************************************************

    print('\n\n********************************************************************')
    print('\tretrieving asset prices')
    print('********************************************************************')

    prices = get_data(tickers, 'adjClose', start_date, end_date, save=False, fail_safe=True)
    print(prices.tail(10))

    print('\n\n********************************************************************')
    print('\tretrieving factor data')
    print('********************************************************************')

    factors = fama_french(start_date, end_date, save=True)
    print(factors.tail(10))

    print('\n\n********************************************************************')
    print('risk adjusted returns')
    print('********************************************************************')

    returns = (prices / prices.shift(1) - 1).dropna()[:len(factors)]
    R = returns.values
    print(pd.DataFrame(R, index=factors.index, columns=prices.columns).tail(10))

    ## *********************************************************************************************************************
    #  factor model
    ## *********************************************************************************************************************

    factors.drop('RF', axis=1, inplace=True)
    F = np.hstack((np.atleast_2d(np.ones(factors.shape[0])).T, factors))

    print('\n\n********************************************************************')
    print('\tfitting the factor model')
    print('********************************************************************')

    transmat, loadings, covarainces = regime_switch(R, F, tickers)

    print('\n\ntransition matrix')
    print(transmat)

    for i in range(2):
        print('\n\nregime %d' % (i + 1))

        print('\nfactor loadings matrix')
        print(loadings[i])

        print("\ncovariance matrix")
        print(covarainces[i])

    # see what regime best fits the most recent data (# days = baseline)
    baseline = 30
    regime = current_regime(R, F, loadings, baseline)
    print('\nbased on the last %d trading days, the best fitted regime is %d' % (baseline, regime))

    ## *********************************************************************************************************************
    #  expeceted mean and variances from the factor model
    ## *********************************************************************************************************************

    print('\n\n********************************************************************')
    print('\tcalculating period one estimates')
    print('********************************************************************')

    # get the number of days until the next scheduled rebalance
    days = business_days((datetime.strptime(end_date, "%Y-%m-%d") + relativedelta(days=1)).strftime("%Y-%m-%d"),
                         rebalance_date)

    # get the estimate returns and covariances from the factor model
    mu_rsfm = pd.DataFrame(days * expected_returns(F, transmat, loadings, regime), index=tickers, columns=['returns'])
    cov_rsfm = pd.DataFrame(days * covariance(R, F, transmat, loadings, covarainces, regime), index=tickers,
                            columns=tickers)

    # write estimates to a csv file
    #mu_rsfm.to_csv(os.getcwd() + r'/data/mu_rsfm.csv')
    #cov_rsfm.to_csv(os.getcwd() + r'/data/cov_rsfm.csv')

    print('\nexpected returns from the factor model')
    print(mu_rsfm)

    print('\nexpected covariance from the factor model')
    print(cov_rsfm)

    ## *********************************************************************************************************************
    #  black litterman for period one returns
    ## *********************************************************************************************************************

    mktcap = get_mkt_cap(tickers, save=True)

    print("\nmarket cap data")
    print(mktcap)

    # calculate the market coefficient
    l = (gmean(factors.iloc[-days:, :]['MKT'] + 1, axis=0) - 1) / factors.iloc[-days:, :]['MKT'].var()

    mu_bl1, cov_bl1 = bl(tickers=tickers,
                         l=l, tau=1,
                         mktcap=mktcap,
                         Sigma=returns.iloc[-days:, :].cov().values * days,
                         P=np.identity(len(tickers)),
                         Omega=np.diag(np.diag(cov_rsfm)),
                         q=mu_rsfm.values,
                         adjust=False)

    print('\nperiod one returns')
    print(mu_bl1)

    print('\nperiod one covariances')
    print(cov_bl1)

    ## *********************************************************************************************************************
    #  AIDAN ML ... pass along a new returns dataframe, mu_ml, with the same format as mu_rsfm
    ## *********************************************************************************************************************

    print('\n\n********************************************************************')
    print('\tcalculating period two estimates')
    print('********************************************************************')

    # temp mu_ml
    #mu_ml = mu_bl1.mul(
    #    pd.DataFrame(1 + np.random.uniform(-0.05, 0.1, len(tickers)), index=mu_bl1.index, columns=mu_bl1.columns))
    df = pd.DataFrame()
    df['prices'] = prices.apply(lambda x: ','.join(x.astype(str)), axis=1)
    df['prices'] = df.prices.apply(lambda x: [float(y) for y in x.split(',')])
    #df.prices = df.prices.apply(lambda x: ast.literal_eval(x))
    mu_ml = predict(torch.FloatTensor(df['prices'].iloc[-1:].values.tolist()).transpose(0, 1), check_ml=mu_bl1, tickers = tickers)
    ## *********************************************************************************************************************
    #  black litterman for period two returns
    ## *********************************************************************************************************************

    mu_bl2, cov_bl2 = bl(tickers=tickers,
                         l=l, tau=1,
                         mktcap=mktcap,
                         Sigma=returns.iloc[-days:, :].cov().values * days,
                         P=np.identity(len(tickers)),
                         Omega=np.diag(np.diag(cov_rsfm)),
                         q=mu_ml.values,
                         adjust=True)

    print('\nperiod two returns')
    print(mu_bl2)

    print('\nperiod two covariances')
    print(cov_bl2)

    ## *********************************************************************************************************************
    #  calculate transaction cost coefficients
    ## *********************************************************************************************************************

    print('\n\n********************************************************************')
    print('\tcalculating cost coefficients')
    print('********************************************************************')

    cost = costs(tickers=tickers,
                 cov=cov_rsfm,
                 prices=prices.iloc[-2, :] if prices.iloc[-1, :].isnull().values.any() else prices.iloc[-1, :],
                 start_date=(datetime.strptime(end_date, "%Y-%m-%d") - relativedelta(years=1)).strftime("%Y-%m-%d"),
                 end_date=end_date,
                 alpha=5)

    print('\ncost coefficients')
    print(cost)

    return (mu_bl1.values.ravel(), mu_bl2.values.ravel()), (cov_bl1.values, cov_bl2.values), cost, prices.iloc[-2, :].values if prices.iloc[-1, :].isnull().values.any() else prices.iloc[-1, :].values

reason: from scipy