text:
import pandas as pd
import numpy as np
from pathlib import Path
from data_params import Data

import matplotlib.pyplot as plt
import matplotlib
matplotlib.font_manager.findSystemFonts(fontpaths=None, fontext='ttf')
import seaborn as sns

sns.set(style="darkgrid")

from data_utils import (
    preprocess_train_df,
    save_transformed_stats)

def summarize_stats(csv_name, dir_to_data, stat_name_select):

    path_to_input_df = Path(dir_to_data, csv_name)
    path_to_output_df = Path(dir_to_data, "funded_" + csv_name)

    df = pd.read_csv(path_to_input_df)

    df_funded = df.loc[df["STATUS"] == "funded"].drop(columns=["STATUS"])
    df_expired = df.loc[df["STATUS"] == "expired"].drop(columns=["STATUS"])

    df_funded_described = df_funded.describe(include="all")
    df_funded_described
    df_expired_described = df_expired.describe(include="all")
    df_expired_described

    funded_stat = df_funded_described.loc[stat_name_select]
    funded_stat = funded_stat.to_frame()
    funded_stat.reset_index(level=0, inplace=True)
    funded_stat.columns = ["FeatureName", "FeatureValueFunded"]

    expired_stat = df_expired_described.loc[stat_name_select]
    expired_stat = expired_stat.to_frame()
    expired_stat.reset_index(level=0, inplace=True)
    expired_stat.columns = ["FeatureName", "FeatureValueExpired"]

    funded_stat = funded_stat.merge(expired_stat, on="FeatureName")

    funded_stat.to_csv(path_to_output_df, index=False)

    return funded_stat

# %% ###########
data_par = Data()
cols_process = data_par.cols_process
cols_output = data_par.cols_output
valid_status = data_par.valid_status
dir_to_saved_data = data_par.dir_to_saved_data
dir_to_query_data = data_par.dir_to_query_data
path_to_training_data = data_par.path_to_training_data
stat_name_select = data_par.stat_name_select

predict = False

# csv_name_tags = "stats_tags_df.csv"
# csv_name_loanuse = "stats_loanuse_df.csv"
# csv_name_desc = "stats_desc_df.csv"

# df = pd.read_csv(
#     path_to_training_data, usecols=cols_process
#     , skiprows=lambda i: i % 10 != 0)
df = pd.read_csv(
    path_to_training_data, usecols=cols_process)

df = preprocess_train_df(df, valid_status, cols_output, predict)

df = df[["STATUS", "WAS_TRANSLATED", "ORIGINAL_LANGUAGE"]]
df = df[df["ORIGINAL_LANGUAGE"] != "None"]
df = df[df["WAS_TRANSLATED"] != np.nan]

df_status = df[["STATUS"]]
df_status = df_status.groupby(["STATUS"])["STATUS"].count()
df_status = pd.DataFrame(df_status)
df_status.columns = ["count"]
df_status = df_status.reset_index()
df_status["percent"] = 100 * df_status["count"] / df_status["count"].sum()
df_status

df_lang = df[["STATUS", "ORIGINAL_LANGUAGE"]].copy()
df_lang = df_lang.groupby(["STATUS", "ORIGINAL_LANGUAGE"])["STATUS"].count()
df_lang = df_lang.groupby(level=[1]).apply(lambda x: x / x.sum())
df_lang = df_lang.reset_index(level=[1])
df_lang.columns = ["ORIGINAL_LANGUAGE", "ratio"]
df_lang = df_lang.reset_index()
df_lang.columns = ["STATUS", "ORIGINAL_LANGUAGE", "ratio"]
df_lang = df_lang[df_lang["STATUS"] == "funded"]
df_lang = df_lang[["ORIGINAL_LANGUAGE", "ratio"]]
df_lang.columns = ["ORIGINAL_LANGUAGE", "percent_not_funded"]
df_lang["percent_not_funded"] = 100 * (1 - df_lang["percent_not_funded"])
df_lang = df_lang.reset_index().drop(columns=["index"])
df_lang = df_lang[df_lang["percent_not_funded"] != 0.0]
df_lang = df_lang.sort_values("percent_not_funded", ascending=False)
df_lang

df_trans = df[["STATUS", "WAS_TRANSLATED"]].copy()
df_trans = df_trans.groupby(["STATUS", "WAS_TRANSLATED"])["STATUS"].count()
df_trans = df_trans.groupby(level=[1]).apply(lambda x: x / x.sum())
df_trans = df_trans.reset_index(level=[1])
df_trans.columns = ["WAS_TRANSLATED", "ratio"]
df_trans = df_trans.reset_index()
df_trans.columns = ["STATUS", "WAS_TRANSLATED", "ratio"]
df_trans = df_trans[df_trans["STATUS"] == "funded"]
df_trans = df_trans[["WAS_TRANSLATED", "ratio"]]
df_trans.columns = ["WAS_TRANSLATED", "percent_expired"]
df_trans.columns = ["WAS_TRANSLATED", "percent_not_funded"]
df_trans["percent_not_funded"] = 100 * (1 - df_trans["percent_not_funded"])
# df_trans["WAS_TRANSLATED"] = int(df_trans["WAS_TRANSLATED"])
df_trans = df_trans.replace({0: "no", 1: "yes"})
df_trans = df_trans.reset_index().drop(columns=["index"])
df_trans

sns.set_style(
    "white",
    {
        "axes.spines.bottom": False,
        "axes.spines.top": False,
        "axes.spines.right": False,
        "axes.spines.left": False,
        "axes.grid": False,},)
fig_dpi = 400
fig_size = (4.5, 2.5)

plt.figure(figsize=fig_size)
plt.rcParams['font.family'] = 'sans-serif'
plt.rcParams['font.sans-serif'] = 'Arial'
cnt_plot1 = sns.barplot(y="ORIGINAL_LANGUAGE", x="percent_not_funded", data=df_lang, palette=("Greens_r"))
for p in cnt_plot1.patches:
             cnt_plot1.annotate("%.1f" % p.get_width(), (p.get_width(), p.get_y() + p.get_height() / 2.),
                 ha='center', va='center', fontsize=11, color='white', xytext=(-11, 0),
                 textcoords='offset points', weight='semibold')

cnt_plot1.set(xticklabels=[])
cnt_plot1.set(xlabel="")
cnt_plot1.set(ylabel="")
cnt_plot1.set(title="Campaigns Not Funded (%)")
cnt_plot1.figure.tight_layout()
fig = cnt_plot1.get_figure()
fig.savefig("language_not_funded.png", dpi=fig_dpi)
plt.close()

plt.figure(figsize=(1.75, 1.75))
cnt_plot2 = sns.barplot(x="WAS_TRANSLATED", y="percent_not_funded", data=df_trans)
# cnt_plot2.set_xticklabels(cnt_plot2.get_xticklabels(), rotation=30)
cnt_plot2.figure.tight_layout()
cnt_plot2.set(xlabel="Translated?")
cnt_plot2.set(ylabel="Campaigns Not Funded (%)")
fig = cnt_plot2.get_figure()
fig.savefig("translated_not_funded.png", dpi=fig_dpi)
plt.close()

max_stars_repo_name: far-from-normal/Kiva-Insight

max_issues_repo_name: far-from-normal/Kiva-Insight

max_forks_repo_name: far-from-normal/Kiva-Insight


text:
# encoding: utf-8
#----------------------------------------------------------------------------
#  Copyright (C) 2008-2014, IPython Development Team and Enthought, Inc.
#  Distributed under the terms of the BSD License.  See COPYING.rst.
#----------------------------------------------------------------------------

"""
Load the seismic volume, created with create_volume.py, and do stuff with it.

Usage:
    $ python load_volume.py --file <filename> --key <keyname> --hdf5 --dnpy

The --file argument specifies the input filename, default 'seismic.hdf5'.
The --key argument specifies the name of the data block in the HDF5 file, default 'seismic'.
The --hdf5 argument specifies to read from an HDF5 file, default True.
The --dnpy argument specifies to read from .dnpy files, default False.
Only one of --hdf5 or --dnpy can be specified.

The volume is loaded in parallel as a DistArray.
It can be read from either a single HDF5 file or a set of .dnpy files.
Each trace (a z-slice constant in x and y) has statistics calculated for it.
Next we apply a couple of filters to each trace.
Next we extract some slices from the volume and create plots for them.
We also extract slices from one of the filtered results.
Finally, we save the 3-point-average filtered volume, both as an HDF5 file,
and as a set of .dnpy files.
"""

from __future__ import print_function

import argparse
import os.path

import numpy
import h5py
from matplotlib import pyplot

from distarray.globalapi import Context, Distribution
from distarray.globalapi.distarray import DistArray

# Load/save of seismic data from HDF5 for .dnpy files.

def get_hdf5_dataset_shape(pathname, key):
    ''' Get the shape of the array in the HDF5 file. '''
    with h5py.File(pathname, 'r') as f:
        dataset = f[key]
        shape = dataset.shape
    return shape

def load_hdf5_distarray(context, filename, key, dist):
    ''' Create a distarray from the specified section of the HDF5 file. '''
    # Filename for load_hdf5() needs the full path.
    pathname = os.path.abspath(filename)
    # Get array shape.
    print('Getting array shape...')
    array_shape = get_hdf5_dataset_shape(pathname, key)
    # Create distribution.
    print('Creating distribution...')
    distribution = Distribution(context, array_shape, dist=dist)
    # Load HDF5 file into DistArray.
    print('Loading HDF5 file...')
    distarray = context.load_hdf5(filename=pathname,
                                  distribution=distribution,
                                  key=key)
    print('Loaded.')
    return distarray

def save_hdf5_distarray(filename, key, distarray):
    ''' Write a distarray to an HDF5 file. '''
    pathname = os.path.abspath(filename)
    context = distarray.context
    context.save_hdf5(pathname, distarray, key=key, mode='w')

def load_dnpy_distarray(context, filename):
    ''' Create a distarray from the .dnpy files. '''
    # Remove extension from filename.
    filename = os.path.splitext(filename)[0]
    pathname = os.path.abspath(filename)
    print('Loading .dnpy files...')
    distarray = context.load_dnpy(pathname)
    print('Loaded.')
    return distarray

def save_dnpy_distarray(filename, distarray):
    ''' Write a distarray to .dnpy files. '''
    pathname = os.path.abspath(filename)
    context = distarray.context
    context.save_dnpy(pathname, distarray)

# Print some stuff about the array.
# (Mostly only practical for small ones.)

def dump_distarray_info(da):
    ''' Print some stuff about the array. '''
    print('Local Shapes:')
    localshapes = da.localshapes()
    print(localshapes)
    print('Arrays Per Process:')
    ndarrays = da.get_ndarrays()
    for i, ndarray in enumerate(ndarrays):
        print('Process:', i)
        print(ndarray)
    print('Full Array:')
    ndarray = da.toarray()
    print(ndarray)

# Calculate statistics on each trace [z slice].
# We get the same results when calculating via DistArray,
# and via NumPy on the full array.

def calc_distributed_stats(distarray):
    '''Calculate some statistics on each trace [z slices] in the distarray.'''
    trace_min = distarray.min(axis=2)
    trace_max = distarray.max(axis=2)
    trace_mean = distarray.mean(axis=2)
    trace_var = distarray.var(axis=2)
    trace_std = distarray.std(axis=2)
    # Collect into dict.
    distarray_stats = {
        'min': trace_min,
        'max': trace_max,
        'mean': trace_mean,
        'var': trace_var,
        'std': trace_std,
    }
    return distarray_stats

def calc_undistributed_stats(ndarray):
    ''' Calculate some statistics on each trace [z slices] in the ndarray. '''
    trace_min = ndarray.min(axis=2)
    trace_max = ndarray.max(axis=2)
    trace_mean = ndarray.mean(axis=2)
    trace_var = ndarray.var(axis=2)
    trace_std = ndarray.std(axis=2)
    # Collect into dict.
    undistributed_stats = {
        'min': trace_min,
        'max': trace_max,
        'mean': trace_mean,
        'var': trace_var,
        'std': trace_std,
    }
    return undistributed_stats

def compare_stats(distarray_stats, numpy_stats, verbose=False):
    ''' Compare the statistics made on DistArray vs NumPy array. '''

    stat_keys = ['min', 'max', 'mean', 'var', 'std']

    def print_stats(title, stat_dict):
        ''' Print out the statistics. '''
        print(title)
        for stat in stat_keys:
            print(stat)
            print(stat_dict[stat])

    # Convert DistArray to NumPy array to allow comparison.
    for stat in stat_keys:
        distarray_stats[stat] = distarray_stats[stat].toarray()

    # The difference is ideally zero.
    if verbose:
        print('Calculating difference...')
    difference_stats = {}
    for stat in stat_keys:
        diff = distarray_stats[stat] - numpy_stats[stat]
        difference_stats[stat] = diff

    # Print statistics.
    if verbose:
        print_stats('Distributed stats:', distarray_stats)
        print_stats('Undistributed stats:', numpy_stats)
        print_stats('Difference:', difference_stats)

    # Test that results are effectively the same with allclose().
    for stat in stat_keys:
        distributed_stat = distarray_stats[stat]
        undistributed_stat = numpy_stats[stat]
        # For 512x512x1024, single precision, var has errors which exceed
        # the default allclose() bounds, many about 0.004.
        is_close = numpy.allclose(distributed_stat, undistributed_stat)
        print(stat, 'is_close:', is_close)

def analyze_statistics(distarray, compare=True, verbose=False):
    ''' Calculate statistics on each trace, via DistArray and NumPy.

    The results should match within numerical precision.
    '''
    if verbose:
        print('Calculating statistics...')
    # Using DistArray methods.
    if verbose:
        print('DistArray statistics...')
    distributed_stats = calc_distributed_stats(distarray)
    if compare:
        # Convert to NumPy array and use NumPy methods.
        if verbose:
            print('NumPy array statistics...')
        ndarray = distarray.tondarray()
        if verbose:
            print('Converted to ndarray...')
        undistributed_stats = calc_undistributed_stats(ndarray)
        # Compare.
        compare_stats(distributed_stats, undistributed_stats, verbose=verbose)

# 3-point averaging filter. This uses a 2-point average on the ends.

def filter_avg3(a):
    ''' Filter a numpy array via 3-point average over z axis.
    A 2-point average is used at the ends. '''
    from numpy import empty_like
    b = empty_like(a)
    b[:, :, 0] = (a[:, :, 0] + a[:, :, 1]) / 2.0
    b[:, :, 1:-1] = (a[:, :, :-2] + a[:, :, 1:-1] + a[:, :, 2:]) / 3.0
    b[:, :, -1] = (a[:, :, -2] + a[:, :, -1]) / 2.0
    return b

def local_filter_avg3(la):
    ''' Filter a local array via 3-point average over z axis. '''

    def filter_avg3(a):
        ''' Filter a local array via 3-point average over z axis.
        A 2-point average is used at the ends. '''
        from numpy import empty_like
        b = empty_like(a)
        b[:, :, 0] = (a[:, :, 0] + a[:, :, 1]) / 2.0
        b[:, :, 1:-1] = (a[:, :, :-2] + a[:, :, 1:-1] + a[:, :, 2:]) / 3.0
        b[:, :, -1] = (a[:, :, -2] + a[:, :, -1]) / 2.0
        return b

    from distarray.localapi import LocalArray
    a = la.ndarray
    b = filter_avg3(a)
    res = LocalArray(la.distribution, buf=b)
    rtn = proxyize(res)
    return rtn

# 3-point maximum window filter.

def filter_max3(a):
    ''' Filter a numpy array via a 3-element window maximum. '''
    from numpy import empty_like
    b = empty_like(a)
    shape = a.shape
    for k in xrange(shape[2]):
        k0 = max(k - 1, 0)
        k1 = min(k + 1, shape[2] - 1)
        b[:, :, k] = a[:, :, k0:k1 + 1].max(axis=2)
    return b

def local_filter_max3(la):
    ''' Filter a local array via 3-point average over z axis. '''

    def filter_max3(a):
        ''' Filter a numpy array via a 3-element window maximum. '''
        from numpy import empty_like
        b = empty_like(a)
        shape = a.shape
        for k in xrange(shape[2]):
            k0 = max(k - 1, 0)
            k1 = min(k + 1, shape[2] - 1)
            b[:, :, k] = a[:, :, k0:k1 + 1].max(axis=2)
        return b

    from distarray.localapi import LocalArray
    a = la.ndarray
    b = filter_max3(a)
    res = LocalArray(la.distribution, buf=b)
    rtn = proxyize(res)
    return rtn

def distributed_filter(distarray, local_filter):
    ''' Filter a DistArray, returning a new DistArray. '''
    context = distarray.context
    filtered_key = context.apply(local_filter, (distarray.key,))
    filtered_da = DistArray.from_localarrays(filtered_key[0], context=context)
    return filtered_da

def undistributed_filter(ndarray, numpy_filter):
    ''' Filter a NumPy array, returning a new NumPy array. '''
    filtered_nd = numpy_filter(ndarray)
    return filtered_nd

def analyze_filter(da, local_filter, numpy_filter, compare, verbose):
    ''' Apply the filter both via DistArray methods and via NumPy methods. '''
    # Via DistArray.
    result_distarray = distributed_filter(da, local_filter)
    if verbose:
        # Print results of averaging.
        print('Original:')
        print(da.toarray())
        print('Filtered:')
        print(result_distarray.toarray())
    if compare:
        # Filter via NumPy array.
        ndarray = da.toarray()
        result_numpy = undistributed_filter(ndarray, numpy_filter)
        # Difference between DistArray and NumPy results.
        distributed_filtered = result_distarray.toarray()
        undistributed_filtered = result_numpy
        diff = distributed_filtered - undistributed_filtered
        print('Difference of DistArray - NumPy filters:')
        if verbose:
            print(diff)
        is_close = numpy.allclose(distributed_filtered, undistributed_filtered)
        print('is_close:', is_close)
    # Return the filtered distarray.
    return result_distarray

# Slicing examples.

def plot_slice(distarray_slice, filename, title, x_label, y_label):
    ''' Create an array plot of the 2D slice. '''
    print('Visualizing slice', title)
    num_dim = len(distarray_slice.shape)
    if (num_dim != 2):
        raise ValueError('Slice must be 2D for plotting.')
    # Convert to ndarray for plotting.
    slice_nd = distarray_slice.toarray()
    # Transpose to make depth look like depth.
    slice_nd = slice_nd.transpose()
    # Plot.
    cmap = 'RdBu'
    pyplot.matshow(slice_nd, cmap=cmap)
    if False:
        # This looks bad without more effort.
        pyplot.colorbar()
    pyplot.title(title)
    pyplot.xlabel(x_label)
    pyplot.ylabel(y_label)
    pyplot.savefig(filename, dpi=100)

def slice_volume(distarray, base_filename='slice_'):
    ''' Slice the volume three different ways and plot result. '''
    shape = distarray.shape
    # Choose (arbitrary) indices for slicing.
    i_index = shape[0] // 4
    j_index = (2 * shape[1]) // 3
    k_index = shape[2] // 2
    # Slice and visualize result.
    print('Taking I-Slice...')
    i_slice = distarray[i_index, :, :]
    filename = base_filename + 'i_plot.png'
    title = 'Vertical Slice [i = %d]' % (i_index)
    plot_slice(i_slice, filename, title, 'j', 'depth')
    print('Taking J-Slice...')
    filename = base_filename + 'j_plot.png'
    title = 'Vertical Slice [j = %d]' % (j_index)
    j_slice = distarray[:, j_index, :]
    plot_slice(j_slice, filename, title, 'i', 'depth')
    print('Taking K-Slice...')
    filename = base_filename + 'k_plot.png'
    title = 'Horizontal Slice [k = %d]' % (k_index)
    k_slice = distarray[:, :, k_index]
    plot_slice(k_slice, filename, title, 'i', 'j')

# Main processing functions.

def load_seismic_volume(filename, key, dist, use_hdf5):
    ''' Load the seismic volume, from HDF5 or .dnpy files. '''
    # Create context.
    context = Context()
    if use_hdf5:
        print('Loading from .hdf5 file...')
        da = load_hdf5_distarray(context, filename, key, dist)
    else:
        print('Loading from .dnpy files...')
        da = load_dnpy_distarray(context, filename)
    # Print some stuff about the array.
    if False:
        dump_distarray_info(da)
    return da

def process_seismic_volume(da, key, compare=True, verbose=False):
    ''' Do some processsing with the seismic volume as a DistArray.

    Also do the same calculations on the global NumPy array,
    to confirm that the distributed methods give the same results.
    '''
    # Slicing.
    print('Slicing...')
    slice_volume(da, base_filename='slice_')
    # Statistics per-trace
    print('Analyzing statistics...')
    analyze_statistics(da, compare=compare, verbose=verbose)
    # 3-point average filter.
    print('Filtering avg3...')
    filtered_avg3_da = analyze_filter(da,
                                      local_filter_avg3,
                                      filter_avg3,
                                      compare=compare,
                                      verbose=verbose)
    # 3-point maximum filter.
    print('Filtering max3...')
    filtered_max3_da = analyze_filter(da,
                                      local_filter_max3,
                                      filter_max3,
                                      compare=compare,
                                      verbose=verbose)
    # Slice the filtered array.
    print('Slicing filtered array...')
    slice_volume(filtered_avg3_da, base_filename='filtered_')
    # Save filtered output to .dnpy files.
    print('Saving .dnpy files...')
    output_dnpy_filename = 'filtered_avg3'
    save_dnpy_distarray(output_dnpy_filename, filtered_avg3_da)
    # Save filtered distarray to new HDF5 file.
    print('Saving .hdf5 file...')
    output_filename = 'filtered_avg3.hdf5'
    save_hdf5_distarray(output_filename, key, filtered_avg3_da)

def main():
    # Parse arguments:
    #     --file <filename>
    #     --key <keyname>
    #     --hdf5
    #     --dnpy
    parser = argparse.ArgumentParser()
    parser.add_argument('--file',
                        default='seismic.hdf5',
                        help='Name of input file.')
    parser.add_argument('--key',
                        default='seismic',
                        help='Name of HDF5 key for data.')
    parser.add_argument('--hdf5',
                        action='store_true',
                        help='Read input from an HDF5 file.')
    parser.add_argument('--dnpy',
                        action='store_true',
                        help='Read input from .dnpy files.')
    args = parser.parse_args()
    # Extract arguments and complain about invalid ones.
    filename = args.file
    key = args.key
    use_hdf5 = args.hdf5
    use_dnpy = args.dnpy
    # Pick either HDF5 or .dnpy
    if (use_hdf5 == False) and (use_dnpy == False):
        use_hdf5 = True
    if (use_hdf5 == True) and (use_dnpy == True):
        raise ValueError('Can only specify one of --hdf5 or --dnpy.')
    # Desired distribution method.
    dist = ('b', 'b', 'n')
    # Processing options.
    compare = True
    verbose = False
    # Load the seismic volume.
    da = load_seismic_volume(filename, key, dist, use_hdf5)
    # Process the seismic volume.
    process_seismic_volume(da,
                           key=key,
                           compare=compare,
                           verbose=verbose)

if __name__ == '__main__':
    main()
    print('Done.')

ext: py

max_stars_repo_name: sjperkins/distarray

max_issues_repo_name: sjperkins/distarray

max_forks_repo_name: sjperkins/distarray

reason: import numpy,from numpy


text:
#!/usr/bin/env python3

import argparse
import csv
import logging
import math
import sys

import numpy as np
import torch
import torch.hub
import torch.nn.functional as F
import torchvision as tv
import tqdm
from torch.utils.data import DataLoader
from torchvision.io import VideoReader
from torchvision.transforms.functional import F_t
from torch.optim.lr_scheduler import OneCycleLR

from datasets import (AlgaeClassificationDataset, AlgaeUnlabeledDataset)

BACKBONES = [
    'vgg16', 'densenet161', 'shufflenet_v2_x1_0', 'mobilenet_v2',
    'mobilenet_v3_large', 'mobilenet_v3_small', 'resnet18', 'resnet34',
    'resnet50', 'resnet101', 'resnet152', 'efficientnet_b0', 'efficientnet_b1',
    'efficientnet_b2', 'efficientnet_b3', 'efficientnet_b4', 'efficientnet_b5',
    'efficientnet_b6', 'efficientnet_b7', 'fpn_resnet18', 'fpn_resnet34',
    'fpn_resnet50', 'fpn_resnet101', 'fpn_resnet152', 'fpn_efficientnet_b0',
    'fpn_efficientnet_b1', 'fpn_efficientnet_b2', 'fpn_efficientnet_b3',
    'fpn_efficientnet_b4', 'fpn_efficientnet_b5', 'fpn_efficientnet_b6',
    'fpn_efficientnet_b7'
]

AVIRIS_TO_SENTINEL2 = [10, 15, 22, 33, 37, 41, 45, 50, 62, 107, 132, 193]

SENTINEL2_TO_4B = [0, 2, 3, 8]

def cli_parser():
    parser = argparse.ArgumentParser()
    parser.add_argument('--backbone', required=True, type=str, choices=BACKBONES)
    parser.add_argument('--batch-size', required=False, type=int, default=128)
    parser.add_argument('--classification-savezs', required=True, type=str, nargs='+')
    parser.add_argument('--epochs1', required=False, type=int, default=33)
    parser.add_argument('--epochs2', required=False, type=int, default=2003)
    parser.add_argument('--epochs3', required=False, type=int, default=33)
    parser.add_argument('--lr1', required=False, type=float, default=1e-4)
    parser.add_argument('--lr2', required=False, type=float, default=1e-5)
    parser.add_argument('--num-workers', required=False, type=int, default=0)
    parser.add_argument('--prescale', required=False, type=int, default=1)
    parser.add_argument('--pth-cheaplab-donor', required=False, type=str, default=None)
    parser.add_argument('--pth-load', required=False, type=str, default=None)
    parser.add_argument('--pth-save', required=False, type=str, default=None)
    parser.add_argument('--unlabeled-epoch-size', required=False, type=int, default=1e6)
    parser.add_argument('--unlabeled-savezs', required=False, default=[], type=str, nargs='+')
    parser.add_argument('--w0', required=False, type=float, default=1.0)
    parser.add_argument('--w1', required=False, type=float, default=0.0)
    parser.add_argument('--w2', required=False, type=float, default=0.5)
    parser.add_argument('--w3', required=False, type=float, default=0.5)

    parser.add_argument('--ndwi-mask', required=False, dest='ndwi_mask', action='store_true')
    parser.set_defaults(ndwi_mask=False)

    parser.add_argument('--no-schedule', dest='schedule', action='store_false')
    parser.set_defaults(schedule=True)

    parser.add_argument('--no-pretrained', dest='pretrained', action='store_false')
    parser.set_defaults(pretrained=True)

    return parser

def worker_init_fn(x):
    np.random.seed(42 + x)

dataloader1_cfg = {
    'batch_size': None,
    'num_workers': None,
    'shuffle': True,
    'worker_init_fn': worker_init_fn
}

dataloader2_cfg = {
    'batch_size': None,
    'num_workers': None,
    'shuffle': True,
    'worker_init_fn': worker_init_fn
}

# https://discuss.pytorch.org/t/how-to-freeze-bn-layers-while-training-the-rest-of-network-mean-and-var-wont-freeze/89736/9
def freeze_bn(m):
    for (name, child) in m.named_children():
        if isinstance(child, torch.nn.BatchNorm2d):
            for param in child.parameters():
                param.requires_grad = False
            child.eval()
        else:
            freeze_bn(child)

def unfreeze_bn(m):
    for (name, child) in m.named_children():
        if isinstance(child, torch.nn.BatchNorm2d):
            for param in child.parameters():
                param.requires_grad = True
            child.train()
        else:
            unfreeze_bn(child)

def freeze(m):
    for p in m.parameters():
        p.requires_grad = False

def unfreeze(m):
    for p in m.parameters():
        p.requires_grad = True

def greens_function(x, x0, eps=1e-3):
    out = -((x - x0) * (x - x0)) / (eps * eps)
    out = torch.exp(out)
    return out

def entropy_function(x):
    s = torch.sigmoid(x)
    pno_narrow = torch.mean(greens_function(s, 0.25, 1.0 / 16))
    pno_wide = torch.mean(greens_function(s, 0.25, 1.0 / 8))
    pyes_narrow = torch.mean(greens_function(s, 0.75, 1.0 / 16))
    pyes_wide = torch.mean(greens_function(s, 0.75, 1.0 / 8))
    return ((pyes_narrow * torch.log(pyes_wide)) +
            (pno_narrow * torch.log(pno_wide)))

if __name__ == '__main__':

    args = cli_parser().parse_args()
    logging.basicConfig(stream=sys.stderr, level=logging.INFO, format='%(asctime)-15s %(message)s')
    log = logging.getLogger()

    dataloader1_cfg['batch_size'] = args.batch_size
    dataloader1_cfg['num_workers'] = args.num_workers

    dataloader2_cfg['batch_size'] = args.batch_size
    dataloader2_cfg['num_workers'] = args.num_workers

    device = torch.device('cuda')
    model = torch.hub.load(
        'jamesmcclain/algae-classifier:df888fa9c383c976faecada5bef7844afe53cba7',
        'make_algae_model',
        in_channels=[4, 12, 224],
        prescale=args.prescale,
        backbone_str=args.backbone,
        pretrained=args.pretrained)

    if args.pth_cheaplab_donor:
        state = torch.load(args.pth_cheaplab_donor)
        for key in list(state.keys()):
            if 'cheaplab' not in key:
                state.pop(key)
        model.load_state_dict(state, strict=False)

    for cheaplab in model.cheaplab.values():
        cheaplab.to(device)
    model.to(device)

    opt1a = torch.optim.AdamW(model.parameters(), lr=args.lr1)
    opt1b = torch.optim.AdamW(model.parameters(), lr=args.lr2)
    if args.epochs2 > 0:
        sched1 = OneCycleLR(opt1b, max_lr=args.lr2, total_steps=args.epochs2)

    opt2a = torch.optim.SGD(model.parameters(), lr=args.lr1, momentum=0.9)
    opt2b = torch.optim.SGD(model.parameters(), lr=args.lr2, momentum=0.9)
    if args.epochs2 > 0:
        sched2 = OneCycleLR(opt2b, max_lr=args.lr2, total_steps=args.epochs2)

    obj1 = torch.nn.BCEWithLogitsLoss().to(device)
    obj2 = torch.nn.CrossEntropyLoss(ignore_index=2).to(device)

    log.info(f'backbone={args.backbone}')
    log.info(f'batch-size={args.batch_size}')
    log.info(f'classification-savezs={args.classification_savezs}')
    log.info(f'epochs1={args.epochs1}')
    log.info(f'epochs2={args.epochs2}')
    log.info(f'epochs3={args.epochs3}')
    log.info(f'ndwi-mask={args.ndwi_mask}')
    log.info(f'num-workers={args.num_workers}')
    log.info(f'prescale={args.prescale}')
    log.info(f'pretrained={args.pretrained}')
    log.info(f'pth-load={args.pth_load}')
    log.info(f'pth-save={args.pth_save}')
    log.info(f'schedule={args.schedule}')
    log.info(f'unlabeled-epoch-size={args.unlabeled_epoch_size}')
    log.info(f'unlabeled-savezs={args.unlabeled_savezs}')

    log.info(f'parameter lr1: {args.lr1}')
    log.info(f'parameter lr2: {args.lr2}')
    log.info(f'parameter w0:   {args.w0}')
    log.info(f'parameter w1:   {args.w1}')
    log.info(f'parameter w2:   {args.w2}')
    log.info(f'parameter w3:   {args.w3}')

    classification_dls = []
    classification_batches = 0
    for savez in args.classification_savezs:
        log.info(f'loading {savez}')
        ad1 = AlgaeClassificationDataset(savezs=[savez],
                                         ndwi_mask=args.ndwi_mask,
                                         augment=True)
        dl1 = DataLoader(ad1, **dataloader1_cfg)
        classification_batches += len(dl1)
        classification_dls.append(dl1)

    unlabeled_dls = []
    for savez in args.unlabeled_savezs:
        log.info(f'loading {savez}')
        ad2 = AlgaeUnlabeledDataset(savezs=[savez],
                                    ndwi_mask=args.ndwi_mask,
                                    augment=True)
        dl2 = DataLoader(ad2, **dataloader2_cfg)
        unlabeled_dls.append(dl2)

    if args.pth_load is None:
        log.info('Training everything')
        unfreeze(model)
        freeze_bn(model)
        for epoch in range(0, args.epochs1):
            constraints1 = []
            entropies1 = []
            entropies2 = []
            losses1 = []
            losses2 = []

            for dl in classification_dls:
                for (i, batch) in enumerate(dl):
                    imgs = batch[0].float()
                    if i == 0:
                        freeze_bn(model.cheaplab)
                        unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                    while imgs is not None:
                        out = model(imgs.to(device))
                        constraint = obj1(out.get('class').squeeze(), batch[1].float().to(device))
                        if 'seg' in out.keys():
                            constraint += obj2(out.get('seg'), batch[2].to(device))
                        if args.w1 != 0.0:
                            entropy = entropy_function(out.get('class').squeeze())
                            loss = args.w0 * constraint + args.w1 * entropy
                            entropies1.append(entropy.item())
                        else:
                            loss = args.w0 * constraint
                        losses1.append(loss.item())
                        constraints1.append(constraint.item())
                        loss.backward()
                        opt1a.step()
                        opt1a.zero_grad()
                        if imgs.shape[1] == 224:
                            imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                        elif imgs.shape[1] == 12:
                            imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                        else:
                            imgs = None

            for dl in unlabeled_dls:
                for (i, batch) in enumerate(dl):
                    if i > args.unlabeled_epoch_size:
                        break
                    imgs = batch[0].float()
                    if i == 0:
                        freeze_bn(model.cheaplab)
                        unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                    while imgs is not None:
                        out = model(imgs.to(device))
                        entropy = entropy_function(out.get('class').squeeze())
                        entropies2.append(entropy.item())
                        loss = args.w3 * entropy
                        losses2.append(loss.item())
                        loss.backward()
                        opt2a.step()
                        opt2a.zero_grad()
                        if imgs.shape[1] == 224:
                            imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                        elif imgs.shape[1] == 12:
                            imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                        else:
                            imgs = None

            mean_constraint1 = np.mean(constraints1)
            mean_entropy1 = np.mean(entropies1)
            mean_entropy2 = np.mean(entropies2)
            mean_loss1 = np.mean(losses1)
            mean_loss2 = np.mean(losses2)

            log.info(f'epoch={epoch:<3d} loss={mean_loss1:+1.5f} entropy={mean_entropy1:+1.5f} constraint={mean_constraint1:+1.5f}')
            log.info(f'          loss={mean_loss2:+1.5f} entropy={mean_entropy2:+1.5f}')

        log.info('Training input and output layers')
        freeze(model)
        unfreeze(model.cheaplab)
        unfreeze(model.first)
        unfreeze(model.last)
        freeze_bn(model)
        for epoch in range(0, args.epochs1):
            constraints1 = []
            entropies1 = []
            losses1 = []

            for dl in classification_dls:
                for (i, batch) in enumerate(dl):
                    imgs = batch[0].float()
                    if i == 0:
                        freeze_bn(model.cheaplab)
                        unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                    while imgs is not None:
                        out = model(imgs.to(device)).get('class').squeeze()
                        constraint = obj1(out, batch[1].float().to(device))
                        if args.w1 != 0.0:
                            entropy = entropy_function(out)
                            loss = args.w0 * constraint + args.w1 * entropy
                            entropies1.append(entropy.item())
                        else:
                            loss = args.w0 * constraint
                        losses1.append(loss.item())
                        constraints1.append(constraint.item())
                        loss.backward()
                        opt1a.step()
                        opt1a.zero_grad()
                        if imgs.shape[1] == 224:
                            imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                        elif imgs.shape[1] == 12:
                            imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                        else:
                            imgs = None

            mean_constraint1 = np.mean(constraints1)
            mean_entropy1 = np.mean(entropies1)
            mean_loss1 = np.mean(losses1)

            log.info(f'epoch={epoch:<3d} loss={mean_loss1:+1.5f} entropy={mean_entropy1:+1.5f} constraint={mean_constraint1:+1.5f}')
    else:
        log.info(f'Loading model from {args.pth_load}')
        model.load_state_dict(torch.load(args.pth_load))
        for cheaplab in model.cheaplab.values():
            cheaplab.to(device)
        model.to(device)

    log.info('Training everything')
    unfreeze(model.backbone)
    freeze_bn(model)
    for epoch in range(0, args.epochs2):
        constraints1 = []
        constraints2 = []
        entropies1 = []
        entropies2 = []
        losses1 = []
        losses2 = []

        for dl in classification_dls:
            for (i, batch) in enumerate(dl):
                imgs = batch[0].float()
                if i == 0:
                    freeze_bn(model.cheaplab)
                    unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                while imgs is not None:
                    out = model(imgs.to(device))
                    constraint = obj1(out.get('class').squeeze(), batch[1].float().to(device))
                    if 'seg' in out.keys():
                        constraint += obj2(out.get('seg'), batch[2].to(device))
                    if args.w1 != 0.0:
                        entropy = entropy_function(out)
                        loss = args.w0 * constraint + args.w1 * entropy
                        entropies1.append(entropy.item())
                    else:
                        loss = args.w0 * constraint
                    losses1.append(loss.item())
                    constraints1.append(constraint.item())
                    loss.backward()
                    opt1b.step()
                    opt1b.zero_grad()
                    if imgs.shape[1] == 224:
                        imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                    elif imgs.shape[1] == 12:
                        imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                    else:
                        imgs = None

        for dl in unlabeled_dls:
            for (i, batch) in enumerate(dl):
                if i > args.unlabeled_epoch_size:
                    break
                imgs = batch[0].float()
                if i == 0:
                    freeze_bn(model.cheaplab)
                    unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                while imgs is not None:
                    out = model(imgs.to(device))
                    entropy = entropy_function(out.get('class').squeeze())
                    entropies2.append(entropy.item())
                    if 'seg' in out.keys():
                        constraint = obj2(out.get('seg'), batch[1].to(device))
                        loss = (i/args.unlabeled_epoch_size) * args.w2 * constraint + args.w3 * entropy
                        constraints2.append(constraint.item())
                    else:
                        loss = args.w3 * entropy
                    losses2.append(loss.item())
                    loss.backward()
                    opt2b.step()
                    opt2b.zero_grad()
                    if imgs.shape[1] == 224:
                        imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                    elif imgs.shape[1] == 12:
                        imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                    else:
                        imgs = None

        if args.schedule:
            sched1.step()
            sched2.step()
        if epoch % 107 == 0:
            log.info(f'Saving checkpoint to /tmp/checkpoint.pth')
            torch.save(model.state_dict(), '/tmp/checkpoint.pth')

        mean_constraint1 = np.mean(constraints1)
        mean_constraint2 = np.mean(constraints2)
        mean_entropy1 = np.mean(entropies1)
        mean_entropy2 = np.mean(entropies2)
        mean_loss1 = np.mean(losses1)
        mean_loss2 = np.mean(losses2)

        log.info(f'epoch={epoch:<3d} loss={mean_loss1:+1.5f} entropy={mean_entropy1:+1.5f} constraint={mean_constraint1:+1.5f}')
        log.info(f'          loss={mean_loss2:+1.5f} entropy={mean_entropy2:+1.5f} constraint={mean_constraint2:+1.5f}')

    log.info(f'Saving checkpoint to /tmp/checkpoint.pth')
    torch.save(model.state_dict(), '/tmp/checkpoint.pth')

    log.info('Reconciling CheapLabs')
    freeze(model)
    unfreeze(model.cheaplab['224'])
    unfreeze(model.cheaplab['4'])
    freeze_bn(model)
    for epoch in range(0, args.epochs3):
        losses1 = []

        for dl in unlabeled_dls:
            for (i, batch) in enumerate(dl):
                if i > args.unlabeled_epoch_size or batch[0].shape[1] < 224:
                    break
                imgs = batch[0].float()
                if i == 0:
                    freeze_bn(model.cheaplab)
                    unfreeze_bn(model.cheaplab[str(imgs.shape[1])])
                outs = []
                outs.append(model(imgs.to(device)).get('class'))
                if imgs.shape[1] == 224:
                    imgs = imgs[:, AVIRIS_TO_SENTINEL2, :, :]
                    outs.append(model(imgs.to(device)).get('class'))
                if imgs.shape[1] == 12:
                    imgs = imgs[:, SENTINEL2_TO_4B, :, :]
                    outs.append(model(imgs.to(device)).get('class'))

                if len(outs[0].shape) == 2:
                    label = torch.stack([o for o in outs], axis=2)
                    label = torch.round(torch.sigmoid(torch.mean(label, axis=2))).float()
                else:
                    label = torch.stack([o for o in outs], axis=4)
                    label = torch.round(torch.sigmoid(torch.mean(label, axis=4))).float()
                loss = obj1(outs[0], label) + obj1(outs[2], label)
                losses1.append(loss.item())
                loss.backward()
                opt1a.step()
                opt1a.zero_grad()

        mean_loss1 = np.mean(losses1)

        log.info(f'epoch={epoch:<3d} loss={mean_loss1:+1.5f}')

    if args.pth_save is not None:
        log.info(f'Saving model to {args.pth_save}')
        torch.save(model.state_dict(), args.pth_save)

    model.eval()
    tp = 0.0
    tn = 0.0
    fp = 0.0
    fn = 0.0
    with torch.no_grad():
        for (i, batch) in enumerate(dl1):
            pred = torch.sigmoid(model(batch[0].float().to(device)).get('class')).squeeze()
            pred = (pred > 0.5).detach().cpu().numpy().astype(np.uint8)
            gt = batch[1].detach().cpu().numpy().astype(np.uint8)
            tp += np.sum((pred == 1) * (gt == 1))
            tn += np.sum((pred == 0) * (gt == 0))
            fp += np.sum((pred == 1) * (gt == 0))
            fn += np.sum((pred == 0) * (gt == 1))
    total = tp + tn + fp + fn
    log.info(f'tpr={tp/(tp+fn)} tnr={tn/(tn+fp)}')

max_stars_repo_name: jamesmcclain/algae-model

max_issues_repo_name: jamesmcclain/algae-model

max_forks_repo_name: jamesmcclain/algae-model

reason: import numpy


text:
import numpy
from QDTK.Operator import Operator

def get_operator_matrix(operator: Operator) -> numpy.ndarray:
    nprod = operator.countProducts()
    grid_points = operator.pdim[0]

    matrix = numpy.zeros((grid_points, grid_points), dtype=numpy.complex128)
    unit_matrix = numpy.diag(numpy.ones(grid_points))

    for k in range(1, nprod + 1):
        # pylint: disable=protected-access
        term = operator._table.get((k, 1))
        coefficient = operator._table.get(k)

        if coefficient is None:
            coefficient = 1.0
        else:
            coefficient = coefficient.coef

        # unit operator
        if term is None:
            matrix += coefficient * unit_matrix
            continue

        # matrix operator
        if term.type == "matrix":
            if term.dim != grid_points:
                raise ValueError(
                    "invalid dimension {} of matrix term, expected {}".format(
                        term.dim, grid_points
                    )
                )
            matrix += coefficient * term.oterm.astype(numpy.complex128)
            continue

        # vector operator
        if term.type == "diag":
            if term.dim != grid_points:
                raise ValueError(
                    "invalid dimension {} of vector term, expected {}".format(
                        term.dim, grid_points
                    )
                )
            matrix += coefficient * numpy.diag(term.oterm.astype(numpy.complex128))
            continue

        # invalid type
        raise ValueError(f"invalid term type {term.type}")

    return matrix

max_stars_repo_name: f-koehler/mlxtk

max_issues_repo_name: f-koehler/mlxtk

max_forks_repo_name: f-koehler/mlxtk


text:
import numpy as np
import pandas as pd
import random
import pystan

#random.seed(121)
#output = pd.read_csv('../data/output_1000stars.csv', usecols=range(1,4))
output = pd.read_csv('~/y4project/data/output_back_filesremoved.csv', usecols=['kic','numax', 'numax_err'])
output = output.loc[(np.abs(output['numax'] - 75) < 30.0)].reset_index(drop=True)
# Select a random number of stars from the file, recording their IDs
IDs = [random.choice(output['kic']) for i in range(6)]
IDs = [10802837,12003742,2283075,4760954,11358669,11521197]
# Create a list of the numax values corresponding to each ID
Numax = [(output.loc[(output['kic'] == IDs[i])]).iloc[0]['numax'] for i in range(len(IDs))]
Numax_err = [(output.loc[(output['kic'] == IDs[i])]).iloc[0]['numax_err'] for i in range(len(IDs))]

modesID = [pd.read_csv('~/y4project/data/rgbmodes/modes_'+str(IDs[i])+'.csv', usecols=['f0', 'f0_err', 'A0'])
           for i in np.arange(0,len(IDs),1)]
lenmodes = [len(modesID[i]) for i in range(len(IDs))]
maxmodes = max(lenmodes)
arr_n = np.zeros([len(IDs),maxmodes])
arr_freq = np.zeros([len(IDs),maxmodes])
arr_freqerr = np.zeros([len(IDs),maxmodes])

dnu_avgID = []
for i in np.arange(0,len(IDs),1):
    modesID[i] = modesID[i].sort_values(by=['f0'])
    modesID[i] = modesID[i].set_index(np.arange(0,len(modesID[i]),1))
    modesID[i]['dnu'] = (modesID[i]['f0'].diff(2).shift(-1))/2
    dnu_avg = (np.mean(modesID[i]['dnu']))
    dnu_avgID.append(dnu_avg)

    n_min = int(modesID[i]['f0'].min() / dnu_avg)
    n = np.arange(n_min, n_min+len(modesID[i]), 1)
    modesID[i].insert(loc=0, column='n', value=n)

    # Loop to ensure all arrays are the same length
    if lenmodes[i] < maxmodes:
        l = lenmodes[i]
        while l < maxmodes:
            newrow = {'n': int(np.max(modesID[i]['n'])+1),
                      'f0': np.max(modesID[i]['f0'])+dnu_avgID[i],
                      'f0_err': 100000}
            modesID[i] = modesID[i].append(newrow, ignore_index=True)
            l += 1

    arr_n[i,:] = modesID[i]['n']
    arr_freq[i,:] = modesID[i]['f0']
    arr_freqerr[i,:] = modesID[i]['f0_err']

    #display(modesID[i])

dnu_avgarr = np.asarray(dnu_avgID)

#Define start parameters
eps = []
for i in range(len(IDs)):
    epsilon1 = np.median((modesID[i].f0 % dnu_avg) / dnu_avg)
    eps.append(epsilon1)
epsilon = np.asarray(eps)

numax_obs = Numax
nmax = numax_obs/dnu_avgarr - epsilon
alpha = 0.015*dnu_avgarr**(-0.32)
A = 0.06*dnu_avgarr**(-0.88)
G = np.ones([len(IDs)])*3.08
tau = np.ones([len(IDs)])*200
phi = np.ones([len(IDs)])*1.71

import pickle
iterations = int(input("How many iterations would you like to use? " ))
sm = pickle.load(open('no_tau_stan.pkl', 'rb'))

stan_data = {'N': len(IDs),
         'M': maxmodes,
         'n': arr_n,
         'freq': arr_freq,
         'freq_err': arr_freqerr,
         'numax_obs': Numax,
         'numax_err': Numax_err,
         'dnu_guess': dnu_avgarr
        }
start = {'dnu': dnu_avgarr,
         'numax': Numax,
         'eps_std': 0.01,
         'al_std': 0.01,
         'A_err': 0.01,
         #'G_err': 0.65,
         #'epsilon': epsilon,
         #'alpha': alpha,
         #'A': A,
         'G': G,
         'phi': phi,
         'epsA': 0.601,
         'epsB': 0.632,
         'alA': 0.015,
         'alB': 0.32,
         'AA': 0.06,
         'AB': 0.88,
         #'GA': 3.08
         #'tau': tau
    }
nchains = 4

fit = sm.sampling(data=stan_data, iter=iterations, chains=nchains, init=[start for n in range(nchains)])

print(fit)

params = np.zeros([len(IDs), 8])

for i in np.arange(0,len(IDs),1):
    params[i] = [np.mean(fit['dnu'], axis=0)[i],
           np.mean(fit['numax'], axis=0)[i],
           np.mean(fit['epsilon'], axis=0)[i],
           np.mean(fit['alpha'], axis=0)[i],
           np.mean(fit['A'], axis=0)[i],
           np.mean(fit['G'], axis=0)[i],
           np.mean(fit['phi'], axis=0)[i],
                np.nan]

np.savetxt('no_tau_models.csv', params, delimiter=',')


text:
#!/usr/bin/env python

from __future__ import division
import numpy
import sys

def DistMat(cord):
    '''
    Equivalent to Matlab pdsit function. It takes input a NX3 numpy array (cord)
    and returns a NXN matrix where the ijth cell is the distance between the ith
    and jth point in NX3 matrix
    '''
    d0 = numpy.subtract.outer(cord[:,0], cord[:,0])**2
    d1 = numpy.subtract.outer(cord[:,1], cord[:,1])**2
    d2 = numpy.subtract.outer(cord[:,2], cord[:,2])**2
    dist = numpy.sqrt(d0+d1+d2)
    return dist

def Rotate(cord, axis=None, thetha=None):
    '''
    This function applies rotation matrix along the axis specified in 'axis'
    on the input NX3 numpy array (cord) and returns the new set of coordnates as NX3 array.
    'axis' can be either of x, y, or z in which case one single rotation is performed, or it
    can be a list of axis (for ex: [z,x,y,x...]) in which case multiple rotations are performed
    in the order of axis values in 'axis'.
    'thetha' is the angle of rotation in degrees, It can be a single value or a list of values only if axis
    is also a list
    '''
    try:
        assert(axis != None and thetha !=None) #Check if axis is passed
    except AssertionError:
        sys.exit('*** AXIS and/or THETHA CANNOT BE NONE ***. \nValid input to axis is one of: x, y, z or [list of axis in the order to perform multiple rotations].')

    if not type(axis) is list:
        axis = [axis]
    if not type(thetha) is list:
        thetha = [thetha]

    try:
        assert(len(axis) == len(thetha))
    except AssertionError:
        sys.exit('*** Number of input axis does not match number of thetha values ***')

    cos_term = numpy.cos(numpy.radians(numpy.array(thetha)))
    sine_term = numpy.sin(numpy.radians(numpy.array(thetha)))
    newcord = numpy.transpose(cord)

    for rot_axis, cos_thetha, sin_thetha in zip(axis,cos_term,sine_term):
        if rot_axis.lower() == 'x':
            rot_mat = numpy.array([[1, 0, 0], [0, cos_thetha, -sin_thetha], [0, sin_thetha, cos_thetha]])
        elif rot_axis.lower() == 'y':
            rot_mat = numpy.array([[cos_thetha, 0, sin_thetha], [0, 1, 0], [-sin_thetha, 0, cos_thetha]])
        elif rot_axis.lower() == 'z':
            rot_mat = numpy.array([[cos_thetha, -sin_thetha, 0], [sin_thetha, cos_thetha, 0], [0, 0, 1]])
        else:
            sys.exit('*** INVALID AXIS ARGUMENT %s ***' % rot_axis+'. \nValid input to axis is one of: x, y, z or [list of x,y,z ].')

        newcord = numpy.dot(rot_mat, newcord)
    return numpy.transpose(newcord)


text:
'''
This file contains code to characterize detected hand movement. These features are used to characterize bradykinesia.
The features calculated are: 1. hand movement amplitude, 2. smoothness of hand movement (jerk measure)
'''
import pandas as pd
import numpy as np
from signal_preprocessing import preprocess
from features import signal_features as sf

def compute_rms(data_array):
    '''
    Compute RMS of data.

    :param data_array: np.array of data
    :return: rms
    '''
    return np.sqrt(np.mean(np.square(data_array)))

def calculate_amplitude_and_smoothness_features(raw_accelerometer_data_df, fs):
    '''
    Function to calculate hand movement amplitude and smoothness of hand movement (jerk metric) from accelerometer data
    collected from a wrist worn wearable device.

    :param raw_accelerometer_data_df: Pandas DataFrame of raw accelerometer data. Columns = ['ts', 'x', 'y', 'z']
    :param fs: Sampling rate of raw accelerometer data (Float)
    :return: Computed hand movement amplitude (list) and smoothness of hand movement (jerk metric) (list) in 3 second
    windows
    '''

    # Pre-process data
    # Bandpass filter between 0.25-3.5
    filtered_data_df = preprocess.band_pass_filter(raw_accelerometer_data_df, fs, [0.25, 3.5], 4,
                                        channels=['x', 'y', 'z'])
    bp_headers = ['x_bp_filt_[0.25, 3.5]', 'y_bp_filt_[0.25, 3.5]', 'z_bp_filt_[0.25, 3.5]']
    filtered_data_df = filtered_data_df[bp_headers]

    avg_acc_per_window = []
    jerk_per_window = []

    # Segment into 3 second windows
    total_samples = filtered_data_df.shape[0]
    window_samples = fs * 3.0
    total_windows = round(total_samples / float(window_samples))
    for win in range(int(total_windows)):
        current_win_start = int(window_samples * win + 1)
        current_win_end = int(current_win_start + window_samples - 1)
        if current_win_end > filtered_data_df.shape[0]:
            window_data_df = filtered_data_df.loc[current_win_start:, :]
        else:
            window_data_df = filtered_data_df.loc[current_win_start:current_win_end, :]
        window_data_df.reset_index(drop=True, inplace=True)
        window_data_df = window_data_df[bp_headers]

        # Compute Avg RMS -> Amplitude of hand movement
        window_data_df['mag'] = np.sqrt(
            (window_data_df['x_bp_filt_[0.25, 3.5]'] ** 2) + (window_data_df['y_bp_filt_[0.25, 3.5]'] ** 2) + (
            window_data_df['z_bp_filt_[0.25, 3.5]'] ** 2))
        combined_amplitude = compute_rms(window_data_df.mag.tolist())
        avg_acc_per_window.append(combined_amplitude)

        # Compute Jerk -> Smoothness of hand movement
        computed_jerk_df = sf.jerk_metric(window_data_df, fs, ['mag'])
        computed_jerk = computed_jerk_df.mag_jerk_ratio.values[0]
        jerk_per_window.append(computed_jerk)

    return avg_acc_per_window, jerk_per_window

if __name__ == "__main__":
    '''
    Main runner to extract features of hand movement.
    '''

    raw_data_filepath = ''# Insert filepath to raw accelerometer sensor data from wrist location.
    raw_data_df = pd.read_csv(raw_data_filepath) # Load data into Pandas DataFrame. (Channels = 'ts','x','y','z')

    fs = 128. # Enter sampling rate of raw accelerometer data

    # Build Feature Set
    amplitude, jerk = calculate_amplitude_and_smoothness_features(raw_data_df, fs)


text:
# Copyright 2021 Toyota Research Institute.  All rights reserved.

from multiprocessing import Pool
import os
import shutil

import cv2
import numpy as np

IN = "datasets/sceneflow/raw"
OUT = "datasets/sceneflow/flying_things"

def handle_files(args):
    infile, outfile = args

    if ".png" in infile:
        shutil.copy(infile, outfile)
    else:
        disparity = -cv2.imread(infile, cv2.IMREAD_UNCHANGED)
        with open(outfile, "wb") as out_file:
            np.savez_compressed(out_file, disparity)

def main():
    files = []
    for root, dirnames, filenames in os.walk(IN):
        for filename in filenames:

            split_dir = "train" if "/train/" in root else "val"
            camera = "left" if "/left" in root else "right"

            if ".png" in filename:
                out_dir = os.path.join(OUT, split_dir, camera)
                os.makedirs(out_dir, exist_ok=True)

                files.append((os.path.join(root, filename), os.path.join(out_dir, filename)))
            elif ".pfm" in filename:
                out_dir = os.path.join(OUT, split_dir, camera + "_disparity")
                os.makedirs(out_dir, exist_ok=True)

                files.append((os.path.join(root, filename), os.path.join(out_dir, filename.replace(".pfm", ".npz"))))

    pool = Pool(processes=16)
    pool.map(handle_files, files)

if __name__ == "__main__":
    main()

max_stars_repo_name: kevinleestone/mmstereo

max_issues_repo_name: kevinleestone/mmstereo

max_forks_repo_name: kevinleestone/mmstereo


text:
__author__ = 'gkour'

import matplotlib.pyplot as plt
import numpy as np
from creature_actions import Actions
from configsimulator import ConfigSimulator

class Dashboard:

    def __init__(self):

        self._fig = plt.figure(figsize=(9, 5), dpi=120, facecolor='w')
        self._fig.canvas.set_window_title('Origin Dashboard')
        self._axes_pop = self._fig.add_subplot(221)
        self._axes_pop.set_ylabel('Population Size')
        self._line_pop, = self._axes_pop.plot([], [], '-', label=self._axes_pop.yaxis.label.get_text())
        self._axes_age = self._axes_pop.twinx()
        self._axes_age.set_ylabel('AVG Age')
        self._line_age, = self._axes_age.plot([], [], 'y-', label=self._axes_age.yaxis.label.get_text())
        self._axes_pop.legend([self._line_pop, self._line_age],
                              [self._line_pop.get_label(), self._line_age.get_label()], loc=0)

        self._axes_aiq = self._fig.add_subplot(222)
        self._axes_aiq.set_ylabel('Population AIQ')
        self._line_aiq, = self._axes_aiq.plot([], [], '-', label=self._axes_aiq.yaxis.label.get_text())
        self._axes_aiq.set_ylim(bottom=0, top=1)

        self._fig_creatures_loc = self._fig.add_axes([0.05, 0.1, 0.2, 0.3])
        self._fig_creatures_loc.yaxis.set_major_locator(plt.NullLocator())
        self._fig_creatures_loc.xaxis.set_major_locator(plt.NullLocator())
        self._fig_food_loc = self._fig.add_axes([0.26, 0.1, 0.2, 0.3])
        self._fig_food_loc.yaxis.set_major_locator(plt.NullLocator())
        self._fig_food_loc.xaxis.set_major_locator(plt.NullLocator())
        self._fig_action = self._fig.add_subplot(2, 3, 6)
        self._fig_death = self._fig.add_subplot(4, 6, 16)
        self._fig_races = self._fig.add_subplot(4, 6, 22)

    def update_epoch_dash(self, epoch_stats_df):
        if epoch_stats_df is None or epoch_stats_df.empty:
            return

    @staticmethod
    def make_autopct(values):
        def my_autopct(pct):
            if pct < 0.1:
                return ''
            total = sum(values)
            val = int(round(pct * total / 100.0))
            return '{v:d}'.format(v=val)
        return my_autopct

    @staticmethod
    def make_numpct(values):
        def my_autopct(pct):
            if pct == 0:
                return ''
            total = sum(values)
            val = int(round(pct * total / 100.0))
            return '({v:d}, {p:1.1f})'.format(v=val, p=pct / 100)

        return my_autopct

    def update_step_dash(self, step_stats_df):
        if step_stats_df is None or step_stats_df.empty:
            return
        self._line_pop.set_xdata(step_stats_df['Time'])
        self._line_pop.set_ydata(step_stats_df['Population'])

        self._line_age.set_xdata(step_stats_df['Time'])
        self._line_age.set_ydata(step_stats_df['MeanAge'])

        self._line_aiq.set_xdata(step_stats_df['Time'])
        self._line_aiq.set_ydata(step_stats_df['AIQ'])

        self._fig.canvas.draw()
        self._fig.canvas.flush_events()
        self._axes_pop.relim()
        self._axes_pop.autoscale_view()
        self._axes_age.relim()
        self._axes_age.autoscale_view()
        self._axes_aiq.relim()
        self._axes_aiq.autoscale_view()

        ## Creatures Dist
        creatures_dist = np.asarray(step_stats_df['CreaturesDist'].iloc[-1])
        self._fig_creatures_loc.clear()
        self._fig_creatures_loc.imshow(creatures_dist, cmap="Purples", aspect="auto", vmin=0, vmax=10)
        self._fig_creatures_loc.set_title('Creatures Location')
        self._fig_creatures_loc.yaxis.set_major_locator(plt.NullLocator())
        self._fig_creatures_loc.xaxis.set_major_locator(plt.NullLocator())

        ## Food Supply
        food_supply = np.asarray(step_stats_df['FoodDist'].iloc[-1])
        self._fig_food_loc.clear()
        self._fig_food_loc.imshow(food_supply, cmap="Greens", aspect="auto", vmin=0, vmax=100)
        self._fig_food_loc.set_title('Food Dist')
        self._fig_food_loc.yaxis.set_major_locator(plt.NullLocator())
        self._fig_food_loc.xaxis.set_major_locator(plt.NullLocator())

        ## Action Dist Pie
        actions_dist = np.mean(step_stats_df['ActionDist'].tail(ConfigSimulator.LOGGING_BATCH_SIZE).values, axis=0)
        self._fig_action.clear()
        patches, texts, autotexts = self._fig_action.pie(actions_dist,
                                                         startangle=90, autopct=self.make_autopct(actions_dist))
        self._fig_action.legend(patches, labels=Actions.get_actions_str(), loc=(1, 0))

        ## Death Dist Pie
        death_cause = np.mean(step_stats_df['DeathCause'].tail(ConfigSimulator.LOGGING_BATCH_SIZE).values, axis=0)
        self._fig_death.clear()
        self._fig_death.pie(death_cause, labels=['Fatigue', 'Fight', 'Elderly', 'Fall'],
                            startangle=90, autopct=self.make_autopct(death_cause))

        ## races Dist
        races = np.mean(step_stats_df['RacesDist'].tail(ConfigSimulator.LOGGING_BATCH_SIZE).values, axis=0)
        self._fig_races.clear()
        self._fig_races.pie(races, labels=[race.race_name() for race in ConfigSimulator.RACES],
                            startangle=90, autopct=self.make_autopct(races))

    def get_figure(self):
        return self._fig

max_stars_repo_name: kourgeorge/life

max_issues_repo_name: kourgeorge/life

max_forks_repo_name: kourgeorge/life


text:
"""
Bootstrap resampling tools.

Compute estimator bias, variance, confidence intervals with bootstrap resampling.

Several forms of bootstrapping on N-dimensional data are supported: ordinary, balanced,
extended, parametric, and stratified sampling, see :func:`resample` for details.
Parametric bootstrapping fits a user-specified distribution to the data and samples
from the parametric distribution. The distributions are taken from scipy.stats.

Confidence intervals can be computed with the ordinary percentile method and with the
more efficient BCa method, see :func:`confidence_interval` for details.
"""
import typing as _tp

import numpy as np
from scipy import stats

from . import _util
from .empirical import quantile_function_gen
from .jackknife import jackknife

_Kwargs = _tp.Any
_Args = _tp.Any
_ArrayLike = _tp.Collection

def resample(
    sample: _ArrayLike,
    *args: _Args,
    size: int = 100,
    method: str = "balanced",
    strata: _tp.Optional[_ArrayLike] = None,
    random_state: _tp.Optional[_tp.Union[np.random.Generator, int]] = None,
) -> _tp.Generator[np.ndarray, None, None]:
    """
    Return generator of bootstrap samples.

    Parameters
    ----------
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    size : int, optional
        Number of bootstrap samples to generate. Default is 100.
    method : str or None, optional
        How to generate bootstrap samples. Supported are 'ordinary', 'balanced',
        'extended', or a distribution name for a parametric bootstrap.
        Default is 'balanced'. Supported distribution names: 'normal' (also:
        'gaussian', 'norm'), 'student' (also: 't'), 'laplace', 'logistic',
        'F' (also: 'f'), 'beta', 'gamma', 'log-normal' (also: 'lognorm',
        'log-gaussian'), 'inverse-gaussian' (also: 'invgauss'), 'pareto', 'poisson'.
    strata : array-like, optional
        Stratification labels. Must have the same shape as `sample`. Default is None.
    random_state : numpy.random.Generator or int, optional
        Random number generator instance. If an integer is passed, seed the numpy
        default generator with it. Default is to use `numpy.random.default_rng()`.

    Yields
    ------
    ndarray
        Bootstrap sample.

    Examples
    --------
    Compute uncertainty of arithmetic mean.

    >>> from resample.bootstrap import resample
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> fx = np.mean(x)
    >>> fb = []
    >>> for b in resample(x, size=10000, random_state=1):
    ...     fb.append(np.mean(b))
    >>> print(f"f(x) = {fx:.1f} +/- {np.std(fb):.1f}")
    f(x) = 4.5 +/- 0.9

    Compute uncertainty of function applied to multivariate data.

    >>> from resample.bootstrap import resample
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> y = np.arange(10, 20)
    >>> fx = np.mean((x, y))
    >>> fb = []
    >>> for bx, by in resample(x, y, size=10000, random_state=1):
    ...     fb.append(np.mean((bx, by)))
    >>> print(f"f(x, y) = {fx:.1f} +/- {np.std(fb):.1f}")
    f(x, y) = 9.5 +/- 0.9

    Notes
    -----
    Balanced vs. ordinary bootstrap:

    The balanced bootstrap produces more accurate results for the same number of
    bootstrap samples than the ordinary bootstrap, but needs to allocate memory for
    B integers, where B is the number of bootstrap samples. Since values of B larger
    than 10000 are rarely needed, this is usually not an issue.

    Non-parametric vs. parametric bootstrap:

    If you know that the data follow a particular parametric distribution, it is
    better to sample from this parametric distribution, but in most cases it is
    sufficient and more convenient to do a non-parametric bootstrap (using "balanced",
    "ordinary", "extended"). The parametric bootstrap is essential for estimators
    sensitive to the tails of a distribution (for example, a quantile close to 0 or 1).
    In this case, only a parametric bootstrap will give reasonable answers, since the
    non-parametric bootstrap cannot include rare events in the tail if the original
    sample did not have them.

    Extended bootstrap:

    In particle physics and perhaps also in other fields, estimators are used which are
    that are a function of both the size and shape of a sample (for example, fit of a
    peak over smooth background to the mass distribution of decay candidates). In this
    case, the normal bootstrap (parametric or non-parametric) is not correct, since the
    sample size is kept constant. For such cases, one needs the "extended" bootstrap.
    The name alludes to the so-called extended maximum-likelihood (EML) method in
    particle physics. Estimates obtained with the EML need to be bootstrapped with the
    "extended" bootstrap.

    Stratification:

    If the sample consists of several distinct classes, stratification
    ensures that the relative proportions of each class are maintained in each
    replicated sample. This is a stricter constraint than that offered by the
    balanced bootstrap, which only guarantees that classes have the original
    proportions over all replicates, but not within each one replicate.
    """
    sample_np = np.atleast_1d(sample)
    n_sample = len(sample_np)
    args_np: _tp.List[np.ndarray] = []

    if args:
        if not isinstance(args[0], _tp.Collection):
            import warnings

            from numpy import VisibleDeprecationWarning

            warnings.warn(
                "Calling resample with positional instead of keyword parameters is "
                "deprecated",
                VisibleDeprecationWarning,
            )
            kwargs: _tp.Dict[str, _tp.Any] = {
                "size": size,
                "method": method,
                "strata": strata,
                "random_state": random_state,
            }
            if len(args) > len(kwargs):
                raise ValueError("too many arguments")
            for key, val in zip(kwargs, args):
                kwargs[key] = val
            size = kwargs["size"]
            method = kwargs["method"]
            strata = kwargs["strata"]
            random_state = kwargs["random_state"]
            del args
        else:
            args_np.append(sample_np)
            for arg in args:
                arg = np.atleast_1d(arg)
                n_arg = len(arg)
                if n_arg != n_sample:
                    raise ValueError(
                        f"extra argument has wrong length {n_arg} != {n_sample}"
                    )
                args_np.append(arg)

    rng = _util.normalize_rng(random_state)

    if strata is not None:
        strata_np = np.atleast_1d(strata)
        if args_np:
            raise ValueError("Stratified resampling only works with one sample array")
        if len(strata_np) != n_sample:
            raise ValueError("a and strata must have the same shape")
        return _resample_stratified(sample_np, size, method, strata_np, rng)

    if method == "balanced":
        if args_np:
            return _resample_balanced_n(args_np, size, rng)
        else:
            return _resample_balanced_1(sample_np, size, rng)
    if method == "ordinary":
        if args_np:
            return _resample_ordinary_n(args_np, size, rng)
        else:
            return _resample_ordinary_1(sample_np, size, rng)
    if method == "extended":
        if args_np:
            return _resample_extended_n(args_np, size, rng)
        else:
            return _resample_extended_1(sample_np, size, rng)

    if args_np:
        raise ValueError("Parametric resampling only works with one sample array")

    dist = {
        # put aliases here
        "gaussian": stats.norm,
        "normal": stats.norm,
        "log-normal": stats.lognorm,
        "log-gaussian": stats.lognorm,
        "inverse-gaussian": stats.invgauss,
        "student": stats.t,
    }.get(method, None)

    # fallback to scipy.stats name
    if dist is None:
        try:
            dist = getattr(stats, method.lower())
        except AttributeError:
            raise ValueError(f"Invalid family: '{method}'")

    if sample_np.ndim > 1:
        if dist != stats.norm:
            raise ValueError(f"family '{method}' only supports 1D samples")
        dist = stats.multivariate_normal
        if sample_np.ndim > 2:
            raise ValueError("multivariate normal only works with 2D samples")

    return _resample_parametric(sample_np, size, dist, rng)

def bootstrap(
    fn: _tp.Callable, sample: _ArrayLike, *args: _Args, **kwargs: _Kwargs
) -> np.ndarray:
    """
    Calculate function values from bootstrap samples.

    This is equivalent to ``numpy.array([fn(b) for b in resample(sample)])`` and
    implemented for convenience.

    Parameters
    ----------
    fn : Callable
        Bootstrap samples are passed to this function.
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    **kwargs
        Keywords are forwarded to :func:`resample`.

    Returns
    -------
    np.array
        Results of `fn` applied to each bootstrap sample.

    Examples
    --------
    >>> from resample.bootstrap import bootstrap
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> fx = np.mean(x)
    >>> fb = bootstrap(np.mean, x, size=10000, random_state=1)
    >>> print(f"f(x) = {fx:.1f} +/- {np.std(fb):.1f}")
    f(x) = 4.5 +/- 0.9
    """
    gen = resample(sample, *args, **kwargs)
    if args:
        return np.array([fn(*b) for b in gen])
    return np.array([fn(x) for x in gen])

def bias(
    fn: _tp.Callable, sample: _ArrayLike, *args: _Args, **kwargs: _Kwargs
) -> np.ndarray:
    """
    Calculate bias of the function estimate with the bootstrap.

    Parameters
    ----------
    fn : callable
        Function to be bootstrapped.
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    **kwargs
        Keyword arguments forwarded to :func:`resample`.

    Returns
    -------
    ndarray
        Bootstrap estimate of bias (= expectation of estimator - true value).

    Examples
    --------
    Compute bias of numpy.var with and without bias-correction.

    >>> from resample.bootstrap import bias
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> round(bias(np.var, x, size=10000, random_state=1), 1)
    -0.8
    >>> round(bias(lambda x: np.var(x, ddof=1), x, size=10000, random_state=1), 1)
    0.0

    Notes
    -----
    This function has special space requirements, it needs to hold `size` replicates of
    the original sample in memory at once. The balanced bootstrap is recommended over
    the ordinary bootstrap for bias estimation, it tends to converge faster.
    """
    thetas = []
    if args:
        replicates: _tp.List[_tp.List] = [[] for _ in range(len(args) + 1)]
        for b in resample(sample, *args, **kwargs):
            for ri, bi in zip(replicates, b):
                ri.append(bi)
            thetas.append(fn(*b))
        population_theta = fn(*(np.concatenate(r) for r in replicates))
    else:
        replicates = []
        for b in resample(sample, *args, **kwargs):
            replicates.append(b)
            thetas.append(fn(b))
        population_theta = fn(np.concatenate(replicates))
    return np.mean(thetas, axis=0) - population_theta

def bias_corrected(
    fn: _tp.Callable, sample: _ArrayLike, *args: _Args, **kwargs: _Kwargs
) -> np.ndarray:
    """
    Calculate bias-corrected estimate of the function with the bootstrap.

    Parameters
    ----------
    fn : callable
        Estimator. Can be any mapping ℝⁿ → ℝᵏ, where n is the sample size
        and k is the length of the output array.
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    **kwargs
        Keyword arguments forwarded to :func:`resample`.

    Returns
    -------
    ndarray
        Estimate with some bias removed.

    Examples
    --------
    Compute bias-corrected estimate of numpy.var.

    >>> from resample.bootstrap import bias_corrected
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> round(np.var(x), 1)
    8.2
    >>> round(bias_corrected(np.var, x, size=10000, random_state=1), 1)
    9.1
    """
    return fn(sample, *args) - bias(fn, sample, *args, **kwargs)

def variance(
    fn: _tp.Callable, sample: _ArrayLike, *args: _Args, **kwargs: _Kwargs
) -> np.ndarray:
    """
    Calculate bootstrap estimate of variance.

    Parameters
    ----------
    fn : callable
        Estimator. Can be any mapping ℝⁿ → ℝᵏ, where n is the sample size
        and k is the length of the output array.
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    **kwargs
        Keyword arguments forwarded to :func:`resample`.

    Returns
    -------
    ndarray
        Bootstrap estimate of variance.

    Examples
    --------
    Compute variance of arithmetic mean.

    >>> from resample.bootstrap import variance
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> round(variance(np.mean, x, size=10000, random_state=1), 1)
    0.8
    """
    thetas = bootstrap(fn, sample, *args, **kwargs)
    return np.var(thetas, ddof=1, axis=0)

def confidence_interval(
    fn: _tp.Callable,
    sample: _ArrayLike,
    *args: _Args,
    cl: float = 0.95,
    ci_method: str = "bca",
    **kwargs: _Kwargs,
) -> _tp.Tuple[float, float]:
    """
    Calculate bootstrap confidence intervals.

    Parameters
    ----------
    fn : callable
        Function to be bootstrapped.
    sample : array-like
        Original sample.
    *args : array-like
        Optional additional arrays of the same length to resample.
    cl : float, default : 0.95
        Confidence level. Asymptotically, this is the probability that the interval
        contains the true value.
    ci_method : str, {'bca', 'percentile'}, optional
        Confidence interval method. Default is 'bca'. See notes for details.
    **kwargs
        Keyword arguments forwarded to :func:`resample`.

    Returns
    -------
    (float, float)
        Upper and lower confidence limits.

    Examples
    --------
    Compute confidence interval for arithmetic mean.

    >>> from resample.bootstrap import confidence_interval
    >>> import numpy as np
    >>> x = np.arange(10)
    >>> a, b = confidence_interval(np.mean, x, size=10000, random_state=1)
    >>> round(a, 1), round(b, 1)
    (2.6, 6.2)

    Notes
    -----
    Both the 'percentile' and 'bca' methods produce intervals that are invariant to
    monotonic transformations of the data values, a desirable and consistent property.

    The 'percentile' method is straightforward and useful as a fallback. The 'bca'
    method is 2nd order accurate (to O(1/n) where n is the sample size) and generally
    preferred. It computes a jackknife estimate in addition to the bootstrap, which
    increases the number of function evaluations in a direct comparison to
    'percentile'. However the increase in accuracy should compensate for this, with the
    result that less bootstrap replicas are needed overall to achieve the same accuracy.
    """
    if args:
        if not isinstance(args[0], _tp.Collection):
            import warnings

            from numpy import VisibleDeprecationWarning

            warnings.warn(
                "Calling confidence_interval with positional instead of keyword "
                "arguments is deprecated",
                VisibleDeprecationWarning,
            )

            if len(args) == 1:
                (cl,) = args
            elif len(args) == 2:
                cl, ci_method = args
            else:
                raise ValueError("too many arguments")
            args = ()

    if not 0 < cl < 1:
        raise ValueError("cl must be between zero and one")

    thetas = bootstrap(fn, sample, *args, **kwargs)
    alpha = 1 - cl

    if ci_method == "percentile":
        return _confidence_interval_percentile(thetas, alpha / 2)

    if ci_method == "bca":
        theta = fn(sample, *args)
        j_thetas = jackknife(fn, sample, *args)
        return _confidence_interval_bca(theta, thetas, j_thetas, alpha / 2)

    raise ValueError(
        f"ci_method must be 'percentile' or 'bca', but '{ci_method}' was supplied"
    )

def _resample_stratified(
    sample: np.ndarray,
    size: int,
    method: str,
    strata: np.ndarray,
    rng: np.random.Generator,
) -> _tp.Generator[np.ndarray, None, None]:
    # call resample on sub-samples and merge the replicates
    sub_samples = [sample[strata == x] for x in np.unique(strata)]
    for sub_replicates in zip(
        *[resample(s, size=size, method=method, random_state=rng) for s in sub_samples]
    ):
        yield np.concatenate(sub_replicates, axis=0)

def _resample_ordinary_1(
    sample: np.ndarray, size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    # i.i.d. sampling from empirical cumulative distribution of sample
    n = len(sample)
    for _ in range(size):
        yield rng.choice(sample, size=n, replace=True)

def _resample_ordinary_n(
    samples: _tp.List[np.ndarray], size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    n = len(samples[0])
    indices = np.arange(n)
    for _ in range(size):
        m = rng.choice(indices, size=n, replace=True)
        yield tuple(s[m] for s in samples)

def _resample_balanced_1(
    sample: np.ndarray, size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    # effectively computes a random permutation of `size` concatenated
    # copies of `sample` and returns `size` equal chunks of that
    n = len(sample)
    indices = rng.permutation(n * size)
    for i in range(size):
        m = indices[i * n : (i + 1) * n] % n
        yield sample[m]

def _resample_balanced_n(
    samples: _tp.List[np.ndarray], size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    n = len(samples[0])
    indices = rng.permutation(n * size)
    for i in range(size):
        m = indices[i * n : (i + 1) * n] % n
        yield tuple(s[m] for s in samples)

def _resample_extended_1(
    sample: np.ndarray, size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    # randomly generates the sample size from a Poisson distribution
    n = len(sample)
    for i in range(size):
        k = rng.poisson(1, size=n)
        yield np.repeat(sample, k)

def _resample_extended_n(
    samples: _tp.List[np.ndarray], size: int, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    n = len(samples[0])
    for i in range(size):
        k = rng.poisson(1, size=n)
        yield tuple(np.repeat(s, k) for s in samples)

def _fit_parametric_family(
    dist: stats.rv_continuous, sample: np.ndarray
) -> _tp.Tuple[float, ...]:
    if dist == stats.multivariate_normal:
        # has no fit method...
        return np.mean(sample, axis=0), np.cov(sample.T, ddof=1)

    if dist in {stats.f, stats.beta}:
        fit_kwargs = {"floc": 0, "fscale": 1}
    elif dist in {stats.gamma, stats.lognorm, stats.invgauss, stats.pareto}:
        fit_kwargs = {"floc": 0}
    else:
        fit_kwargs = {}

    return dist.fit(sample, **fit_kwargs)  # type: ignore

def _resample_parametric(
    sample: np.ndarray, size: int, dist: stats.rv_continuous, rng: np.random.Generator
) -> _tp.Generator[np.ndarray, None, None]:
    n = len(sample)

    # fit parameters by maximum likelihood and sample from that
    if dist == stats.poisson:
        # - poisson has no fit method and there is no scale parameter
        # - random number generation for poisson distribution in scipy seems to be buggy
        mu = np.mean(sample)
        for _ in range(size):
            yield rng.poisson(mu, size=n)
    else:
        args = _fit_parametric_family(dist, sample)
        dist = dist(*args)
        for _ in range(size):
            yield dist.rvs(size=n, random_state=rng)

def _confidence_interval_percentile(
    thetas: np.ndarray, alpha_half: float
) -> _tp.Tuple[float, float]:
    quant = quantile_function_gen(thetas)
    return quant(alpha_half), quant(1 - alpha_half)

def _confidence_interval_bca(
    theta: float, thetas: np.ndarray, j_thetas: np.ndarray, alpha_half: float
) -> _tp.Tuple[float, float]:
    norm = stats.norm

    # bias correction; implementation notes:
    # - if prop_less is zero, z_naught would become -inf;
    #   we set z_naught to zero then (no bias)
    prop_less = np.mean(thetas < theta)  # proportion of replicates less than obs
    z_naught = norm.ppf(prop_less) if prop_less > 0 else 0.0

    # acceleration; implementation notes:
    # - np.mean returns float even if j_thetas are int,
    #   must convert type explicity to make -= operator work
    # - it is possible that all j_thetas are zero, it then follows
    #   that den and num are zero; we set acc to zero then (no acceleration)
    j_mean = np.mean(j_thetas)
    j_thetas = j_thetas.astype(j_mean.dtype, copy=False)
    j_thetas -= j_mean
    num = np.sum((-j_thetas) ** 3)
    den = np.sum(j_thetas**2)
    acc = num / (6 * den**1.5) if den > 0 else 0.0

    z_low = z_naught + norm.ppf(alpha_half)
    z_high = z_naught + norm.ppf(1 - alpha_half)

    p_low = norm.cdf(z_naught + z_low / (1 - acc * z_low))
    p_high = norm.cdf(z_naught + z_high / (1 - acc * z_high))

    quant = quantile_function_gen(thetas)
    return quant(p_low), quant(p_high)

max_stars_repo_name: scikit-hep/resample

max_issues_repo_name: scikit-hep/resample

max_forks_repo_name: scikit-hep/resample

reason: import numpy,from numpy,from scipy


text:
from rest_framework.response import Response
from rest_framework import status
from .models import FileModel
from .serializers import  FileSerializer
from rest_framework.decorators import api_view
from stl import mesh
import os
import numpy as np
from .utils import MeshSolve
from rest_framework.exceptions import NotFound
from django.core.exceptions import ValidationError
from rest_framework.exceptions import NotFound,APIException,ParseError

def validate_extension_function(name:str):
    """ Function, which check  file extentions
    Parameters
    ----------
    name : str
        the name of file
    Returns
    -------
    Boolean
        if file extention is allowed return True
    """
    accepted_extensions=['.stl']
    file_extension = os.path.splitext(name)[1].lower()
    if file_extension not in accepted_extensions:
        raise ParseError('File extention is not supported', code=400)

    return True

def get_mesh(fileUrl:str):
    """ Get file content.
    Parameters
    ----------
    fileUrl : str
        The relative path of the uploaded file.
    Returns
    -------
    Mesh
        Mesh is a network that is formed of cells and points.
    """
    BASE_DIR = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
    full_path = os.path.join(BASE_DIR,fileUrl.strip("/"))
    response=None

    try:
        return mesh.Mesh.from_file(full_path)
    except Exception as my:
        raise NotFound('The file is not existing or too large, maximum  number of triangles of 100000000 ')

def delete_file(id:int):
    """ After data extraction delete the object
    Parameters
    ----------
    id : int
        The id of uploaded file
    Returns
    -------
    Boolean
        if the file is found and deleted, return true; otherweise false.

    """
    try:
        file=FileModel.objects.get(id=id)
        file.delete()

        return True
    except:
        raise NotFound('The file  not found')

@api_view(['POST'])
def file_exract(request):
    if request.method == 'POST':
        #Check the validity of  file
        validate_extension_function(request.FILES['file'].name)

        serializer = FileSerializer(data=request.data)
        if serializer.is_valid():
            serializer.save()
            mesh=get_mesh(serializer.data['file'])
            if not mesh:
                message={'detail':'Uploaded file is defect'}
                return  Response(message, status=status.HTTP_400_BAD_REQUEST)

            meshObject=MeshSolve(mesh)
            res=meshObject.solve()

            delete_response=delete_file(serializer.data['id'])
            if not delete_response:
                message = {'detail':str(delete_response['status'])}
                return Response(message, status=status.HTTP_404_NOT_FOUND)

            return Response(res, status=status.HTTP_200_OK)
    return Response(serializer.errors, status=status.HTTP_400_BAD_REQUEST)


text:
    # -*- coding: utf-8 -*-
"""
Main function that trains the model from input data, using a c3d pretrained net.

This function gathers the data and feeds tensors into the network, training it.
"""

from utils.model import c3d_base
from utils.dataio import make_crops, BrainSequence
import pandas as pd
import configparser
import time
import argparse
import os
import numpy as np
from sklearn.cross_validation import train_test_split
from sklearn.metrics import log_loss
from keras.optimizers import SGD, Adam
from keras.utils import to_categorical
from keras.callbacks import TensorBoard
import bids.layout

# Parser
def get_parser():
    """Parse data for main function."""
    parser = argparse.ArgumentParser(description='Training 3D conv network.')
    parser.add_argument("--config_file",
                        type=str, nargs=1, required=True, help='config file')
    parser.add_argument("--output_directory_name", type=str, nargs=1,
                        required=True, help='directory where the output will be')
    parser.add_argument("--crop", action="store_true",
                        help='wether to make new cropping or not.')

    return parser

def train(config_file, out_dir_name, crop):
    """
    Execute Main function for training.

    Trains the model with a given dataset.
    """
    t0 = time.time()
    # Load configuration
    # Load the configuration of a given experiment.
    config = configparser.ConfigParser()
    config.read(config_file)

    # Create output directory to store results
    out_dir = (config["folders"]["EXPERIMENTS"] +
               out_dir_name + os.sep)
    # Create out directory
    if not os.path.exists(out_dir):
        os.makedirs(out_dir)

    # Load tranining parameters
    batch_size = int(config["model"]["batch_size"])
    new_img_size = (int(x) for x in config["general"]["new_size"].split(','))
    epochs = int(config["model"]["epochs"])
    weights_file = config["model"]["weights"]
    mean_file = config["model"]["mean"]
    metadata_file = config["data"]["metadata"]
    bids_folder = config["data"]["bids_folder"]

    crops_dir = out_dir + '/' + 'CROPS' + '/'
    if not os.path.exists(crops_dir):
        os.makedirs(crops_dir)
    # The crops folder is a parent directory where all the
    # images are stored. The metadata file associates each
    # image to its label
    crop_metadata_file = out_dir + config["data"]["metadata_crops"]
    # Create crops (optional), if already done there is no need.
    if crop:
        # Load dataset and labels using BIDS
        # The metadata file contains information about
        # all the samples we will use. All the subjects mentioned in the
        # metadata must have a corresponding MRI image present
        # in the BIDS folder.
        bids_folder = config["data"]["bids_folder"]
        metadata_file = config["data"]["metadata"]
        new_size = config["general"]["new_size"].split(',')
        make_crops(bids_folder, metadata_file,
                   crops_dir, crop_metadata_file, new_size)

    # Divide between test, train and validation
    df_metadata = pd.read_csv(crop_metadata_file)

    # Get list of unique subjects
    subj = df_metadata.PTID.values
    dx = df_metadata.DX.values
    s = list(set(zip(subj, dx)))
    x, y = zip(*s)

    # Get train/test/val
    # THOSE ARE SUBJECT NAMES
    rd_seed = 1714
    S_train, S_test, DX_train, DX_test = train_test_split(x, y, test_size=.2, random_state=rd_seed)
    S_train, S_val, DX_train, DX_val = train_test_split(S_train, DX_train, test_size=.2, random_state=rd_seed)

    # Preprocess labels
    label_dict = dict(zip(["NL", "MCI", "AD"], range(0, 3)))

    # GET CORRESPONDING DX AND PATHS OF SAID SUBJECTS
    X_train = df_metadata[df_metadata["PTID"].isin(S_train)].path.values
    Y_train = df_metadata[df_metadata["PTID"].isin(S_train)].DX.map(label_dict, na_action='ignore').values

    X_valid = df_metadata[df_metadata["PTID"].isin(S_val)].path.values
    Y_valid = df_metadata[df_metadata["PTID"].isin(S_val)].DX.map(label_dict, na_action='ignore').values

    # Load neural network parameters
    batch_size = int(config["model"]["batch_size"])
    epochs = int(config["model"]["epochs"])
    weights_file = config["model"]["weights"]
    mean_file = config["model"]["mean"]
    mean = np.load(mean_file)
    # Crop mean
    mean = mean[:, :, 8:120, 30:142]
    print("mean shape")
    print(mean.shape)
    # Convert mean to black and white
    mean = np.mean(mean, axis=0)

    # Create generator File
    BrainSeq = BrainSequence(X_train, to_categorical(Y_train), batch_size, norm='mean',
                             norm_param=mean, color=True, train=True, crop=False, img_aug=False)
    BrainSeq_val = BrainSequence(X_valid, to_categorical(Y_valid), batch_size, norm='mean',
                             norm_param=mean, color=True, train=False, crop=False, img_aug=False)

    # Initialize model
    print(new_img_size)
    x, y, z = new_img_size
    model = c3d_base(input_shape=(x, y, z, 3), weights=True, weights_file=weights_file)
    # sgd = SGD(lr=1e-3, decay=1e-6, momentum=0.9, nesterov=True)
    adam = Adam(lr=0.001, amsgrad=False)

    model.compile(optimizer=adam, loss='categorical_crossentropy', metrics=['accuracy'])

    # Train
    callb = TensorBoard(log_dir=out_dir + 'logs/', histogram_freq=0, batch_size=batch_size,
                        write_graph=True, write_grads=True, write_images=False,
                        embeddings_freq=0, embeddings_layer_names=None,
                        embeddings_metadata=None, embeddings_data=None)

    # use early stopping callback

    model.fit_generator(BrainSeq,
                        steps_per_epoch=None,
                        epochs=epochs,
                        shuffle=True,
                        callbacks=[callb],
                        verbose=1,
                        validation_data=BrainSeq_val)

    # Validate the model with a custom predictive function

    X_test = df_metadata[df_metadata["subj"].isin(S_test)].path.values
    Y_test = df_metadata[df_metadata["subj"].isin(S_test)].DX.map(label_dict, na_action='ignore').values

    BrainSeq_test = BrainSequence(X_test, to_categorical(Y_test), batch_size, norm='mean',
                             norm_param=mean, color=True, train=False, crop=False)

    # evaluate
    score = model.evaluate_generator(BrainSeq_test)
    print(score)

    print('Proces finished.')
    t1 = time.time()
    print('Time to compute the script: ', t1 - t0)

if __name__ == '__main__':
    parser = get_parser()
    args = parser.parse_args()
    train(args.config_file[0], args.output_directory_name[0], args.crop)


text:
"""A module containing the core message-passing functions for belief propagation"""

import functools

import jax
import jax.numpy as jnp

import pgmax.bp.bp_utils as bp_utils

NEG_INF = -100000.0  # A large negative value to use as -inf to avoid NaN's

@jax.jit
def pass_var_to_fac_messages(
    ftov_msgs: jnp.array,
    evidence: jnp.array,
    var_states_for_edges: jnp.array,
) -> jnp.array:
    """Passes messages from Variables to Factors.

    The update works by first summing the evidence and neighboring factor to variable messages for
    each variable. Next, it subtracts messages from the correct elements of this sum to yield the
    correct updated messages.

    Args:
        ftov_msgs: Array of shape (num_edge_state,). This holds all the flattened factor to variable
            messages.
        evidence: Array of shape (num_var_states,) representing the flattened evidence for each variable
        var_states_for_edges: Array of shape (num_edge_states,)
            Global variable state indices for each edge state
    Returns:
        Array of shape (num_edge_state,). This holds all the flattened variable to factor messages.
    """
    var_sums_arr = evidence.at[var_states_for_edges].add(ftov_msgs)
    vtof_msgs = var_sums_arr[var_states_for_edges] - ftov_msgs
    return vtof_msgs

@functools.partial(jax.jit, static_argnames=("num_val_configs", "temperature"))
def pass_fac_to_var_messages(
    vtof_msgs: jnp.ndarray,
    factor_configs_edge_states: jnp.ndarray,
    log_potentials: jnp.ndarray,
    num_val_configs: int,
    temperature: float,
) -> jnp.ndarray:

    """Passes messages from Factors to Variables.

    The update is performed in two steps. First, a "summary" array is generated that has an entry for every valid
    configuration for every factor. The elements of this array are simply the sums of messages across each valid
    config. Then, the info from factor_configs_edge_states is used to apply the scattering operation and
    generate a flat set of output messages.

    Args:
        vtof_msgs: Array of shape (num_edge_state,). This holds all the flattened variable to factor messages.
        factor_configs_edge_states: Array of shape (num_factor_configs, 2)
            factor_configs_edge_states[ii] contains a pair of global factor_config and edge_state indices
            factor_configs_edge_states[ii, 0] contains the global factor config index
            factor_configs_edge_states[ii, 1] contains the corresponding global edge_state index
        log_potentials: Array of shape (num_val_configs, ). An entry at index i is the log potential
            function value for the configuration with global factor config index i.
        num_val_configs: the total number of valid configurations for factors in the factor graph.
        temperature: Temperature for loopy belief propagation.
            1.0 corresponds to sum-product, 0.0 corresponds to max-product.

    Returns:
        Array of shape (num_edge_state,). This holds all the flattened factor to variable messages.
    """
    fac_config_summary_sum = (
        jnp.zeros(shape=(num_val_configs,))
        .at[factor_configs_edge_states[..., 0]]
        .add(vtof_msgs[factor_configs_edge_states[..., 1]])
    ) + log_potentials
    max_factor_config_summary_for_edge_states = (
        jnp.full(shape=(vtof_msgs.shape[0],), fill_value=NEG_INF)
        .at[factor_configs_edge_states[..., 1]]
        .max(fac_config_summary_sum[factor_configs_edge_states[..., 0]])
    )
    ftov_msgs = max_factor_config_summary_for_edge_states - vtof_msgs
    if temperature != 0.0:
        ftov_msgs = ftov_msgs + (
            temperature
            * jnp.log(
                jnp.full(shape=(vtof_msgs.shape[0],), fill_value=jnp.exp(NEG_INF))
                .at[factor_configs_edge_states[..., 1]]
                .add(
                    jnp.exp(
                        (
                            fac_config_summary_sum[factor_configs_edge_states[..., 0]]
                            - max_factor_config_summary_for_edge_states[
                                factor_configs_edge_states[..., 1]
                            ]
                        )
                        / temperature
                    )
                )
            )
        )

    return ftov_msgs

@functools.partial(jax.jit, static_argnames=("max_msg_size"))
def normalize_and_clip_msgs(
    msgs: jnp.ndarray,
    edges_num_states: jnp.ndarray,
    max_msg_size: int,
) -> jnp.ndarray:
    """Performs normalization and clipping of flattened messages

    Normalization is done by subtracting the maximum value of every message from every element of every message,
    clipping is done to keep every message value in the range [-1000, 0].

    Args:
        msgs: Array of shape (num_edge_state,). This holds all the flattened factor to variable messages.
        edges_num_states: Array of shape (num_edges,). Number of states for the variables connected to each edge
        max_msg_size: the max of edges_num_states

    Returns:
        Array of shape (num_edge_state,). This holds all the flattened factor to variable messages
            after normalization and clipping
    """
    msgs = msgs - jnp.repeat(
        bp_utils.segment_max_opt(msgs, edges_num_states, max_msg_size),
        edges_num_states,
        total_repeat_length=msgs.shape[0],
    )
    # Clip message values to be always greater than -1000
    msgs = jnp.clip(msgs, -1000, None)
    return msgs

max_stars_repo_name: antoine-dedieu/PGMax

max_issues_repo_name: antoine-dedieu/PGMax

max_forks_repo_name: antoine-dedieu/PGMax

reason: import jax